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