JSOC_Documentation: proj/mag/nlfff/apps/test

00001 /*
00002  * Module name:     test_localpf.c
00003  *
00004  * Description:     Local PF using Wiegelmann's Green function code
00005  *
00006  * Calling:     
00007  *
00008  * Original source: C NLFFF model by Thomas Wiegelmann (weigelmann@linmpi.mg.de)
00009  *
00010  * Written by:      Xudong Sun (xudongs@stanford.edu)
00011  *
00012  * Version:
00013  *          v1.0        Oct 13 2010
00014  *
00015  * Issues:
00016  *          v1.0
00017  *
00018  * Example:
00019  *          test_preproc "in=su_xudong.nlfff_test_in[2010]" "out=su_xudong.nlfff_preproc" "test=1"
00020  *          test_preproc "in=su_xudong.mercator_vec_S5[2010.07.01_11:48:00_TAI]" "out=su_xudong.nlfff_preproc" "test=0"
00021  */
00022 
00023 
00024 #include <jsoc_main.h>
00025 #include <stdio.h>
00026 #include <stdlib.h>
00027 #include <math.h>
00028 #include <string.h>
00029 #include <sys/time.h>
00030 #include <sys/resource.h>
00031 
00032 #define PI  (M_PI)
00033 #define DTOR    (PI / 180.)
00034 
00035 #define ARRLENGTH(ARR) (sizeof(ARR) / sizeof(ARR[0]))
00036 #define DIE(msg) {fflush(stdout); fprintf(stderr, "%s, status=%d\n", msg, status); return(status);}
00037 #define SHOW(msg) {printf("%s", msg); fflush(stdout);}
00038 
00039 #define Macro
00040 #define MCENTERGRAD(f,id) ((f[i+id]-f[i-id])/(2*h))
00041 #define MLEFTGRAD(f,id)   ((-3*f[i]+4*f[i+id]-f[i+2*id])/(2*h))
00042 #define MRIGHTGRAD(f,id)  ((+3*f[i]-4*f[i-id]+f[i-2*id])/(2*h))
00043 #define GRADX(f,i) ((ix>0 && ix<nx-1) ? (MCENTERGRAD(f,nynz)) : ((ix==0) ? (MLEFTGRAD(f,nynz)) : ((ix==nx-1) ? (MRIGHTGRAD(f,nynz)) : (0.0))))
00044 #define GRADY(f,i) ((iy>0 && iy<ny-1) ? (MCENTERGRAD(f,nz)) : ((iy==0) ? (MLEFTGRAD(f,nz)) : ((iy==ny-1) ? (MRIGHTGRAD(f,nz)) : (0.0))))
00045 #define GRADZ(f,i) ((iz>0 && iz<nz-1) ? (MCENTERGRAD(f,1)) : ((iz==0) ? (MLEFTGRAD(f,1)) : ((iz==nz-1) ? (MRIGHTGRAD(f,1)) : (0.0))))
00046 
00047 // Working part
00048 // ======== OpenMP =========
00049 #ifdef _OPENMP 
00050 #include <omp.h>
00051 #endif 
00052 
00053 #include "rebin.c"
00054 #include "preproc.c"
00055 #include "green.c"
00056 #include "relax.c"
00057 
00058 
00059 /* Timing by T. P. Larson */
00060 
00061 double getwalltime(void)
00062 {
00063     struct timeval tv;
00064     gettimeofday(&tv, NULL);
00065     return tv.tv_sec * 1000.0 + tv.tv_usec/1000.0;
00066 }
00067 
00068 double getcputime(double *utime, double *stime)
00069 {
00070     
00071     struct rusage ru;
00072     getrusage(RUSAGE_SELF, &ru);
00073     *utime = ru.ru_utime.tv_sec * 1000.0 + ru.ru_utime.tv_usec/1000.0;
00074     *stime = ru.ru_stime.tv_sec * 1000.0 + ru.ru_stime.tv_usec/1000.0;
00075     return *utime + *stime;
00076 }
00077 
00078 
00079 
00080 
00081 char *module_name = "test_preproc"; /* Module name */
00082 
00083 ModuleArgs_t module_args[] =
00084 {
00085     {ARG_STRING, "in", NULL, "Input data series."},
00086     {ARG_STRING, "out", NULL,  "Output data series."},
00087     {ARG_DOUBLE, "mu1", "1.0", "For vector magnetogram preprocessing."},
00088     {ARG_DOUBLE, "mu2", "1.0", "For vector magnetogram preprocessing."},
00089     {ARG_DOUBLE, "mu3", "0.001", "For vector magnetogram preprocessing."},
00090     {ARG_DOUBLE, "mu4", "0.01", "For vector magnetogram preprocessing."},
00091     {ARG_INT, "maxit", "10000", "Maximum itertation number."},
00092     {ARG_INT, "test", "0", "1 for test data with 3 layer input cube rather than 3 images."},
00093     {ARG_INT, "VERB", "2", "Level of verbosity: 0=errors/warnings; 1=status; 2=all"},
00094     {ARG_STRING, "COMMENT", " ", "Comments."}, 
00095     {ARG_END}
00096 };
00097 
00098 
00099 /* ################## Main Module ################## */
00100 
00101 int DoIt(void)
00102 {
00103     int status = DRMS_SUCCESS;
00104     
00105 #ifdef _OPENMP
00106     printf("Compiled by OpenMP\n");
00107 #endif
00108     
00109     char *inQuery, *outQuery;
00110     DRMS_RecordSet_t *inRS, *outRS;
00111     int irec, nrecs;
00112     DRMS_Record_t *inRec, *outRec;
00113     DRMS_Segment_t *inSegBx, *inSegBy, *inSegBz;
00114     DRMS_Segment_t *outSegBx, *outSegBy, *outSegBz;
00115     DRMS_Array_t *inArrayBx, *inArrayBy, *inArrayBz;
00116     DRMS_Array_t *outArrayBx, *outArrayBy, *outArrayBz;
00117     float *inDataBx, *inDataBy, *inDataBz;
00118     float *outDataBx, *outDataBy, *outDataBz;
00119     
00120     double *bx0, *by0, *bz0;
00121     double *bx0_now, *by0_now, *bz0_now;
00122     double *Bx, *By, *Bz, *Bx_tmp, *By_tmp, *Bz_tmp;
00123     
00124     int verbflag, prepflag;
00125     int outDims[3];
00126     
00127     int i, j, k, itest, dpt, dpt0;
00128     int nx, ny, nz, nxny, nynz, nxnynz;
00129     int nd, nd_now, maxit;
00130     double mu1, mu2, mu3, mu4;
00131     int multi, multi_now;
00132     int nx_now, ny_now, nz_now, nxny_now, nynz_now, nxnynz_now;
00133     int test;
00134     char *comment;
00135     
00136     // Time measuring
00137     double wt0, wt1, wt;
00138     double ut0, ut1, ut;
00139     double st0, st1, st;
00140     double ct0, ct1, ct;
00141     wt0 = getwalltime();
00142     ct0 = getcputime(&ut0, &st0);
00143     
00144     /* Get parameters */
00145     inQuery = (char *)params_get_str(&cmdparams, "in");
00146     outQuery = (char *)params_get_str(&cmdparams, "out");
00147     comment = (char *)params_get_str(&cmdparams, "COMMENT");
00148     verbflag = params_get_int(&cmdparams, "VERB");
00149     maxit = params_get_int(&cmdparams, "maxit");
00150     mu1 = params_get_double(&cmdparams, "mu1");
00151     mu2 = params_get_double(&cmdparams, "mu2");
00152     mu3 = params_get_double(&cmdparams, "mu3");
00153     mu4 = params_get_double(&cmdparams, "mu4");
00154     test = params_get_int(&cmdparams, "test");
00155     multi = 1;
00156     
00157     /* Open input */
00158     inRS = drms_open_records(drms_env, inQuery, &status);
00159     if (status || inRS->n == 0) DIE("No input data found");
00160     nrecs = inRS->n;
00161     
00162     /* Create output */
00163     outRS = drms_create_records(drms_env, nrecs, outQuery, DRMS_PERMANENT, &status);
00164     if (status) DIE("Output recordset not created");
00165     
00166     /* Do this for each record */
00167     for (irec = 0; irec < nrecs; irec++)
00168     {
00169         
00170         if (verbflag) {
00171             wt1 = getwalltime();
00172             ct1 = getcputime(&ut1, &st1);
00173             printf("processing record %d...\n", irec);
00174         }
00175         
00176         /* Input record and data */
00177         inRec = inRS->records[irec];
00178         if (!test) {
00179             inSegBx = drms_segment_lookup(inRec, "Bx");
00180             inSegBy = drms_segment_lookup(inRec, "By");
00181             inSegBz = drms_segment_lookup(inRec, "Bz");
00182             inArrayBx = drms_segment_read(inSegBx, DRMS_TYPE_FLOAT, &status);
00183             if (status) DIE("No Bx data file found. \n");
00184             inArrayBy = drms_segment_read(inSegBy, DRMS_TYPE_FLOAT, &status);
00185             if (status) DIE("No By data file found. \n");
00186             inArrayBz = drms_segment_read(inSegBz, DRMS_TYPE_FLOAT, &status);
00187             if (status) DIE("No Bz data file found. \n");
00188             inDataBx = (float *)inArrayBx->data;
00189             inDataBy = (float *)inArrayBy->data;
00190             inDataBz = (float *)inArrayBz->data;
00191             nx = inArrayBx->axis[0]; ny = inArrayBx->axis[1];
00192             if (inArrayBy->axis[0] != nx || inArrayBz->axis[0] != nx || 
00193                 inArrayBy->axis[1] != ny || inArrayBz->axis[1] != ny)
00194                 DIE("Dimension error. \n");
00195             nxny = nx * ny; // nynz = ny * nz; nxnynz = nxny * nz;
00196         } else {
00197             inSegBz = drms_segment_lookupnum(inRec, 0);
00198             inArrayBz = drms_segment_read(inSegBz, DRMS_TYPE_FLOAT, &status);
00199             if (status) DIE("No data file found. \n");
00200             if (inArrayBz->naxis != 3) DIE("Wrong data dimension. \n");
00201             nx = inArrayBz->axis[0]; ny = inArrayBz->axis[1];
00202             nxny = nx * ny; // nynz = ny * nz; nxnynz = nxny * nz;
00203             inDataBx = (float *)inArrayBz->data;
00204             inDataBy = inDataBx + nxny;
00205             inDataBz = inDataBx + 2 * nxny;
00206         }
00207         
00208         /* Output data */
00209         outDims[0] = nx; outDims[1] = ny;
00210         
00211         outArrayBx = drms_array_create(DRMS_TYPE_FLOAT, 2, outDims, NULL, &status);
00212         outArrayBy = drms_array_create(DRMS_TYPE_FLOAT, 2, outDims, NULL, &status);
00213         outArrayBz = drms_array_create(DRMS_TYPE_FLOAT, 2, outDims, NULL, &status);
00214         outDataBx = (float *)outArrayBx->data;
00215         outDataBy = (float *)outArrayBy->data;
00216         outDataBz = (float *)outArrayBz->data;
00217         
00218         
00219         /* ======================== */
00220         /* This is the working part.*/
00221         /* ======================== */
00222         
00223         // Copy in magnetogram--------------- CHECKED AGAINST WIEGELMANN'S CODE, COLUMN FIRST
00224         Bx = (double *)(calloc(nxny, sizeof(double)));
00225         By = (double *)(calloc(nxny, sizeof(double)));
00226         Bz = (double *)(calloc(nxny, sizeof(double)));
00227         dpt = 0;
00228         for (i = 0; i < nx; i++)    // See green()
00229             for (j = 0; j < ny; j++) {
00230                 dpt0 = j * nx + i;
00231                 Bx[dpt] = inDataBx[dpt0];
00232                 By[dpt] = inDataBy[dpt0];
00233                 Bz[dpt] = inDataBz[dpt0];
00234                 dpt++;
00235             }
00236         
00237         // Main
00238         preproc(Bx, By, Bz, nx, ny, 
00239                 mu1, mu2, mu3, mu4, 0.0, maxit, verbflag);
00240         
00241         // Copy out result -------------- CHECKED AGAINST WIEGELMANN'S CODE
00242         for (j = 0; j < ny; j++) {
00243         for (i = 0; i < nx; i++) {
00244             outDataBx[j * nx + i] = Bx[i * ny + j];
00245             outDataBy[j * nx + i] = By[i * ny + j];
00246             outDataBz[j * nx + i] = Bz[i * ny + j];
00247         }}
00248         
00249         /* Output record */
00250         outRec = outRS->records[irec];
00251         outSegBx = drms_segment_lookupnum(outRec, 0);
00252         outSegBy = drms_segment_lookupnum(outRec, 1);
00253         outSegBz = drms_segment_lookupnum(outRec, 2);
00254         
00255         for (i = 0; i < 2; i++) {   // For variable dimensions
00256             outSegBx->axis[i] = outArrayBx->axis[i];
00257             outSegBy->axis[i] = outArrayBy->axis[i];
00258             outSegBz->axis[i] = outArrayBz->axis[i];
00259         }
00260         
00261         outArrayBx->parent_segment = outSegBx;
00262         outArrayBy->parent_segment = outSegBy;
00263         outArrayBz->parent_segment = outSegBz;
00264         
00265         /* Result writing */
00266         status = drms_segment_write(outSegBx, outArrayBx, 0);
00267         if (status) DIE("Problem writing file Bx");
00268         drms_free_array(outArrayBx);
00269         status = drms_segment_write(outSegBy, outArrayBy, 0);
00270         if (status) DIE("Problem writing file By");
00271         drms_free_array(outArrayBy);
00272         status = drms_segment_write(outSegBz, outArrayBz, 0);
00273         if (status) DIE("Problem writing file Bz");
00274         drms_free_array(outArrayBz);
00275         
00276         /* Set keywords */
00277         drms_setkey_string(outRec, "BLD_VERS", jsoc_version);
00278         drms_setkey_time(outRec, "DATE", CURRENT_SYSTEM_TIME);
00279         if (test) {
00280             drms_copykey(outRec, inRec, "NUM");
00281         } else {
00282             drms_copykey(outRec, inRec, "T_REC");
00283             drms_copykey(outRec, inRec, "T_OBS");
00284             drms_copykey(outRec, inRec, "PNUM");
00285             drms_setkey_string(outRec, "COMMENT", comment);
00286         }
00287         
00288         drms_setkey_int(outRec, "NUM", 100);
00289         
00290         /* Clean up */
00291         if (!test) {
00292             drms_free_array(inArrayBx);
00293             drms_free_array(inArrayBy);
00294             drms_free_array(inArrayBz);
00295         } else {
00296             drms_free_array(inArrayBz);
00297         }
00298         free(Bx); free(By); free(Bz);
00299         
00300         /* Time measure */
00301         if (verbflag) {
00302             wt = getwalltime();
00303             ct = getcputime(&ut, &st);
00304             printf("record %d done, %.2f ms wall time, %.2f ms cpu time\n", 
00305                    irec, wt - wt1, ct - ct1);
00306         }
00307         
00308     }
00309     
00310     drms_close_records(inRS, DRMS_FREE_RECORD);
00311     drms_close_records(outRS, DRMS_INSERT_RECORD);
00312     
00313     if (verbflag) {
00314         wt = getwalltime();
00315         ct = getcputime(&ut, &st);
00316         printf("total time spent: %.2f ms wall time, %.2f ms cpu time\n", 
00317                wt - wt0, ct - ct0);
00318     }
00319     
00320     return(DRMS_SUCCESS);
00321 }
00322
proj/mag/nlfff/apps/test_preproc.c
