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  *          v3.0        Oct 19 2017
00015  *
00016  * Issues:
00017  *          v1.0
00018  *          v3.0
00019  *          exclude nd layer from energy calculation
00020  *          add -n flag to suppress output
00021  *          add RUNNUM keyword for multiple runs
00022  *          save ND keyword
00023  *
00024  * Example:
00025  *          test_localpf "in=su_xudong.nlfff_test_in[2010]" "out=su_xudong.localpf_test_out" "nz=256" "test=1" "COMMENT=04080510_Preprocessed"
00026  */
00027 
00028 
00029 #include <jsoc_main.h>
00030 #include <stdio.h>
00031 #include <stdlib.h>
00032 #include <math.h>
00033 #include <string.h>
00034 #include <sys/time.h>
00035 #include <sys/resource.h>
00036 
00037 #define PI  (M_PI)
00038 #define DTOR    (PI / 180.)
00039 
00040 #define ARRLENGTH(ARR) (sizeof(ARR) / sizeof(ARR[0]))
00041 #define DIE(msg) {fflush(stdout); fprintf(stderr, "%s, status=%d\n", msg, status); return(status);}
00042 #define SHOW(msg) {printf("%s", msg); fflush(stdout);}
00043 
00044 #define Macro
00045 #define MCENTERGRAD(f,id) ((f[i+id]-f[i-id])/(2*h))
00046 #define MLEFTGRAD(f,id)   ((-3*f[i]+4*f[i+id]-f[i+2*id])/(2*h))
00047 #define MRIGHTGRAD(f,id)  ((+3*f[i]-4*f[i-id]+f[i-2*id])/(2*h))
00048 #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))))
00049 #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))))
00050 #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))))
00051 
00052 // Working part
00053 // ======== OpenMP =========
00054 #ifdef _OPENMP 
00055 #include <omp.h>
00056 #endif 
00057 
00058 #include "rebin.c"
00059 #include "preproc.c"
00060 #include "green.c"
00061 #include "relax.c"
00062 
00063 
00064 /* Timing by T. P. Larson */
00065 
00066 double getwalltime(void)
00067 {
00068   struct timeval tv;
00069   gettimeofday(&tv, NULL);
00070   return tv.tv_sec * 1000.0 + tv.tv_usec/1000.0;
00071 }
00072 
00073 double getcputime(double *utime, double *stime)
00074 {
00075 
00076   struct rusage ru;
00077   getrusage(RUSAGE_SELF, &ru);
00078   *utime = ru.ru_utime.tv_sec * 1000.0 + ru.ru_utime.tv_usec/1000.0;
00079   *stime = ru.ru_stime.tv_sec * 1000.0 + ru.ru_stime.tv_usec/1000.0;
00080   return *utime + *stime;
00081 }
00082 
00083 
00084 
00085 
00086 char *module_name = "test_localpf"; /* Module name */
00087 
00088 ModuleArgs_t module_args[] =
00089 {
00090     {ARG_STRING, "in", NULL, "Input data series."},
00091     {ARG_STRING, "out", NULL,  "Output data series."},
00092     {ARG_INT, "nz", "256", "Number of grids in z axis."},
00093     {ARG_INT, "nd", "16", "Number of grids at boundary."},
00094     {ARG_INT, "PREP", "0", "Preprocess vector magnetogram?"},
00095     {ARG_INT, "RUNNUM", "0", "ID of run"},
00096     {ARG_DOUBLE, "mu1", "1.0", "For vector magnetogram preprocessing."},
00097     {ARG_DOUBLE, "mu2", "1.0", "For vector magnetogram preprocessing."},
00098     {ARG_DOUBLE, "mu3", "0.001", "For vector magnetogram preprocessing."},
00099     {ARG_DOUBLE, "mu4", "0.01", "For vector magnetogram preprocessing."},
00100     {ARG_INT, "maxit", "10000", "Maximum itertation number."},
00101     {ARG_INT, "test", "0", "1 for test data with 3 layer input cube rather than 3 images."},
00102     {ARG_INT, "VERB", "2", "Level of verbosity: 0=errors/warnings; 1=status; 2=all"},
00103     {ARG_FLAG, "n", "", "Suppress output"},
00104     {ARG_STRING, "COMMENT", " ", "Comments."}, 
00105     {ARG_END}
00106 };
00107 
00108 
00109 /* ################## Main Module ################## */
00110 
00111 int DoIt(void)
00112 {
00113     int status = DRMS_SUCCESS;
00114 
00115     #ifdef _OPENMP
00116     printf("Compiled by OpenMP\n");
00117     #endif
00118 
00119     char *inQuery, *outQuery;
00120     DRMS_RecordSet_t *inRS, *outRS;
00121     int irec, nrecs;
00122     DRMS_Record_t *inRec, *outRec;
00123     DRMS_Segment_t *inSegBx, *inSegBy, *inSegBz;
00124     DRMS_Segment_t *outSegBx, *outSegBy, *outSegBz;
00125     DRMS_Array_t *inArrayBx, *inArrayBy, *inArrayBz;
00126     DRMS_Array_t *outArrayBx, *outArrayBy, *outArrayBz;
00127     float *inDataBx, *inDataBy, *inDataBz;
00128     float *outDataBx, *outDataBy, *outDataBz;
00129     
00130     double *bx0, *by0, *bz0;
00131     double *bx0_now, *by0_now, *bz0_now;
00132     double *Bx, *By, *Bz, *Bx_tmp, *By_tmp, *Bz_tmp;
00133     double energy;
00134     
00135     int verbflag, prepflag, noOutput;
00136     int outDims[3];
00137 
00138     int i, j, k, itest, dpt, dpt0;
00139     int nx, ny, nz, nxny, nynz, nxnynz;
00140     int nd, nd_now, maxit;
00141     double mu1, mu2, mu3, mu4;
00142     int multi, multi_now;
00143     int nx_now, ny_now, nz_now, nxny_now, nynz_now, nxnynz_now;
00144     int test, runnum;
00145     char *comment;
00146 
00147     // Time measuring
00148     double wt0, wt1, wt;
00149     double ut0, ut1, ut;
00150     double st0, st1, st;
00151     double ct0, ct1, ct;
00152     wt0 = getwalltime();
00153     ct0 = getcputime(&ut0, &st0);
00154 
00155     /* Get parameters */
00156     inQuery = (char *)params_get_str(&cmdparams, "in");
00157     outQuery = (char *)params_get_str(&cmdparams, "out");
00158     comment = (char *)params_get_str(&cmdparams, "COMMENT");
00159     prepflag = params_get_int(&cmdparams, "PREP");
00160     verbflag = params_get_int(&cmdparams, "VERB");
00161     nz = outDims[2] = params_get_int(&cmdparams, "nz");     // z nodes fixed for dataset
00162     nd = params_get_int(&cmdparams, "nd");
00163     maxit = params_get_int(&cmdparams, "maxit");
00164     mu1 = params_get_double(&cmdparams, "mu1");
00165     mu2 = params_get_double(&cmdparams, "mu2");
00166     mu3 = params_get_double(&cmdparams, "mu3");
00167     mu4 = params_get_double(&cmdparams, "mu4");
00168     test = params_get_int(&cmdparams, "test");
00169     multi = 1;
00170     runnum = params_get_int(&cmdparams, "RUNNUM");
00171     
00172     noOutput = params_isflagset(&cmdparams, "n");
00173     printf("runnum=%d, noOutput=%d\n", runnum, noOutput);
00174 
00175     /* Open input */
00176     inRS = drms_open_records(drms_env, inQuery, &status);
00177     if (status || inRS->n == 0) DIE("No input data found");
00178     nrecs = inRS->n;
00179 
00180     /* Create output */
00181     outRS = drms_create_records(drms_env, nrecs, outQuery, DRMS_PERMANENT, &status);
00182     if (status) DIE("Output recordset not created");
00183 
00184     /* Do this for each record */
00185     for (irec = 0; irec < nrecs; irec++)
00186     {
00187 
00188         if (verbflag) {
00189             wt1 = getwalltime();
00190             ct1 = getcputime(&ut1, &st1);
00191             printf("processing record %d...\n", irec);
00192         }
00193         
00194         /* Input record and data */
00195         inRec = inRS->records[irec];
00196         if (!test) {
00197             inSegBx = drms_segment_lookup(inRec, "Bx");
00198             inSegBy = drms_segment_lookup(inRec, "By");
00199             inSegBz = drms_segment_lookup(inRec, "Bz");
00200             inArrayBx = drms_segment_read(inSegBx, DRMS_TYPE_FLOAT, &status);
00201             if (status) DIE("No Bx data file found. \n");
00202             inArrayBy = drms_segment_read(inSegBy, DRMS_TYPE_FLOAT, &status);
00203             if (status) DIE("No By data file found. \n");
00204             inArrayBz = drms_segment_read(inSegBz, DRMS_TYPE_FLOAT, &status);
00205             if (status) DIE("No Bz data file found. \n");
00206             inDataBx = (float *)inArrayBx->data;
00207             inDataBy = (float *)inArrayBy->data;
00208             inDataBz = (float *)inArrayBz->data;
00209             nx = inArrayBx->axis[0]; ny = inArrayBx->axis[1];
00210             if (inArrayBy->axis[0] != nx || inArrayBz->axis[0] != nx || 
00211                 inArrayBy->axis[1] != ny || inArrayBz->axis[1] != ny)
00212                 DIE("Dimension error. \n");
00213             nxny = nx * ny; nynz = ny * nz; nxnynz = nxny * nz;
00214         } else {
00215            inSegBz = drms_segment_lookupnum(inRec, 0);
00216            inArrayBz = drms_segment_read(inSegBz, DRMS_TYPE_FLOAT, &status);
00217            if (status) DIE("No data file found. \n");
00218            if (inArrayBz->naxis != 3) DIE("Wrong data dimension. \n");
00219            nx = inArrayBz->axis[0]; ny = inArrayBz->axis[1];
00220            nxny = nx * ny; nynz = ny * nz; nxnynz = nxny * nz;
00221            inDataBx = (float *)inArrayBz->data;
00222            inDataBy = inDataBx + nxny;
00223            inDataBz = inDataBx + 2 * nxny;
00224         }
00225 
00226         /* Test if multigrid can be performed */
00227         nx_now = nx; ny_now = ny; nz_now = nz; nd_now = nd;
00228         
00229         multi = 1;      // constant, so no multi-grid
00230 
00231         /* Output data */
00232         outDims[0] = nx; outDims[1] = ny;
00233 
00234         outArrayBx = drms_array_create(DRMS_TYPE_FLOAT, 3, outDims, NULL, &status);
00235         outArrayBy = drms_array_create(DRMS_TYPE_FLOAT, 3, outDims, NULL, &status);
00236         outArrayBz = drms_array_create(DRMS_TYPE_FLOAT, 3, outDims, NULL, &status);
00237         outDataBx = (float *)outArrayBx->data;
00238         outDataBy = (float *)outArrayBy->data;
00239         outDataBz = (float *)outArrayBz->data;
00240 
00241         /* ======================== */
00242         /* This is the working part.*/
00243         /* ======================== */
00244 
00245         // Copy in magnetogram--------------- CHECKED AGAINST WIEGELMANN'S CODE, COLUMN FIRST
00246         bx0 = (double *)(calloc(nxny, sizeof(double)));
00247         by0 = (double *)(calloc(nxny, sizeof(double)));
00248         bz0 = (double *)(calloc(nxny, sizeof(double)));
00249         dpt = 0;
00250         for (i = 0; i < nx; i++)    // See green()
00251         for (j = 0; j < ny; j++) {
00252             dpt0 = j * nx + i;
00253             bx0[dpt] = inDataBx[dpt0];
00254             by0[dpt] = inDataBy[dpt0];
00255             bz0[dpt] = inDataBz[dpt0];
00256             dpt++;
00257         }
00258 
00259         /* Multigrid processing, essetially disabled for this */
00260 
00261         for (multi_now = multi; multi_now > 0; multi_now--) {
00262 
00263             if (verbflag) {
00264                 printf("multi=%d, nx_now=%d, ny_now=%d, nz_now=%d, nd_now=%d\n", 
00265                         multi_now, nx_now, ny_now, nz_now, nd_now);
00266                 fflush(stdout);
00267             }
00268 
00269             /* Step 0: Grid */
00270             nxny_now = nx_now * ny_now; nynz_now = ny_now * nz_now;
00271             nxnynz_now = nxny_now * nz_now;
00272 
00273             /* Step 1: Get magnetogram at this level */
00274             bx0_now = (double *)(calloc(nxny_now, sizeof(double)));
00275             by0_now = (double *)(calloc(nxny_now, sizeof(double)));
00276             bz0_now = (double *)(calloc(nxny_now, sizeof(double)));
00277             SHOW("here");
00278 printf("%d, %d, %d\n", nx_now, ny_now, multi_now); fflush(stdout);
00279             rebin2d(bx0, by0, bz0, bx0_now, by0_now, bz0_now, nx, ny, multi_now);
00280 
00281             /* Step 2: Preprocessing */
00282             if (prepflag) {
00283                 preproc(bx0_now, by0_now, bz0_now, nx_now, ny_now, 
00284                         mu1, mu2, mu3, mu4, 0.0, maxit, verbflag);
00285             }
00286 printf("preproc"); fflush(stdout);
00287             /* Step 3: Get intial guess */
00288             if (multi_now == multi) {
00289                 Bx = (double *)(calloc(nxnynz_now, sizeof(double)));
00290                 By = (double *)(calloc(nxnynz_now, sizeof(double)));
00291                 Bz = (double *)(calloc(nxnynz_now, sizeof(double)));
00292                 // Potential solver
00293                 printf("ahyaya%d", nxnynz_now); fflush(stdout);
00294                 green(bz0_now, Bx, By, Bz, nx_now, ny_now, nz_now, verbflag);
00295                 printf("%d", multi_now); fflush(stdout);
00296             }
00297                 printf("why"); fflush(stdout);
00298 
00299             free(bx0_now); free(by0_now); free(bz0_now);
00300                 printf("heyhey"); fflush(stdout);
00301 
00302             /* Step 0: Grid */
00303             nx_now *= 2; ny_now *= 2; nz_now *= 2; nd_now *= 2;
00304 
00305         }   // end of multi
00306      
00307         // Copy out result -------------- CHECKED AGAINST WIEGELMANN'S CODE
00308         dpt = 0;
00309         for (k = 0; k < nz; k++) {
00310         for (j = 0; j < ny; j++) {
00311         for (i = 0; i < nx; i++) {
00312             dpt0 = i * nynz + j * nz + k;
00313             outDataBx[dpt] = Bx[dpt0];
00314             outDataBy[dpt] = By[dpt0];
00315             outDataBz[dpt] = Bz[dpt0];
00316             dpt++;
00317         }}}
00318         
00319         // no checking of nd    // Oct 19 2017
00320         energy = 0.;        // Feb 13 2011
00321         for (k = 0; k < nz - nd; k++) {
00322         for (j = nd; j < ny - nd; j++) {
00323         for (i = nd; i < nx - nd; i++) {
00324             dpt0 = i * nynz + j * nz + k;
00325             energy += (Bx[dpt0] * Bx[dpt0] + By[dpt0] * By[dpt0] + Bz[dpt0] * Bz[dpt0]);
00326         }}}
00327 
00328         /* Output record */
00329         outRec = outRS->records[irec];
00330         
00331         if (!noOutput) {
00332             
00333             outSegBx = drms_segment_lookupnum(outRec, 0);
00334             outSegBy = drms_segment_lookupnum(outRec, 1);
00335             outSegBz = drms_segment_lookupnum(outRec, 2);
00336             
00337             for (i = 0; i < 3; i++) {    // For variable dimensions
00338                 outSegBx->axis[i] = outArrayBx->axis[i];
00339                 outSegBy->axis[i] = outArrayBy->axis[i];
00340                 outSegBz->axis[i] = outArrayBz->axis[i];
00341             }
00342             
00343             outArrayBx->parent_segment = outSegBx;
00344             outArrayBy->parent_segment = outSegBy;
00345             outArrayBz->parent_segment = outSegBz;
00346             
00347             /* Result writing */
00348             status = drms_segment_write(outSegBx, outArrayBx, 0);
00349             if (status) DIE("Problem writing file Bx");
00350             drms_free_array(outArrayBx);
00351             status = drms_segment_write(outSegBy, outArrayBy, 0);
00352             if (status) DIE("Problem writing file By");
00353             drms_free_array(outArrayBy);
00354             status = drms_segment_write(outSegBz, outArrayBz, 0);
00355             if (status) DIE("Problem writing file Bz");
00356             drms_free_array(outArrayBz);
00357             
00358         }
00359 
00360         /* Set keywords */
00361         drms_setkey_string(outRec, "COMMENT", comment);
00362         drms_setkey_int(outRec, "MULTI", multi);
00363         drms_setkey_int(outRec, "PREP", prepflag);
00364         drms_setkey_string(outRec, "BLD_VERS", jsoc_version);
00365         drms_setkey_time(outRec, "DATE", CURRENT_SYSTEM_TIME);
00366         drms_setkey_double(outRec, "ENERGY", energy);
00367         drms_copykey(outRec, inRec, "BIN");
00368         if (test) {
00369             drms_copykey(outRec, inRec, "NUM");
00370         } else {
00371             drms_copykey(outRec, inRec, "T_REC");
00372             drms_copykey(outRec, inRec, "T_OBS");
00373             drms_copykey(outRec, inRec, "HARPNUM");
00374             drms_setkey_string(outRec, "COMMENT", comment);
00375         }
00376         drms_setkey_int(outRec, "RUNNUM", runnum);
00377         drms_setkey_int(outRec, "ND", nd);
00378 
00379         /* Clean up */
00380         if (!test) {
00381             drms_free_array(inArrayBx);
00382             drms_free_array(inArrayBy);
00383             drms_free_array(inArrayBz);
00384         } else {
00385             drms_free_array(inArrayBz);
00386         }
00387         free(bx0); free(by0); free(bz0);
00388         free(Bx); free(By); free(Bz);
00389 
00390         /* Time measure */
00391         if (verbflag) {
00392             wt = getwalltime();
00393             ct = getcputime(&ut, &st);
00394             printf("record %d done, %.2f ms wall time, %.2f ms cpu time\n", 
00395                      irec, wt - wt1, ct - ct1);
00396         }
00397 
00398     }
00399 
00400     drms_close_records(inRS, DRMS_FREE_RECORD);
00401     drms_close_records(outRS, DRMS_INSERT_RECORD);
00402 
00403     if (verbflag) {
00404         wt = getwalltime();
00405         ct = getcputime(&ut, &st);
00406         printf("total time spent: %.2f ms wall time, %.2f ms cpu time\n", 
00407                 wt - wt0, ct - ct0);
00408     }
00409 
00410     return(DRMS_SUCCESS);
00411 }
00412
proj/mag/nlfff/apps/test_localpf.c
