JSOC_Documentation: proj/mag/d4vm/apps/test

00001 /*
00002  * Module name:     test_d4vm.c
00003  *
00004  * Description:
00005  *          DAVE4VM C version
00006  *
00007  * Calling:
00008  *                  
00009  *
00010  * Original by:     Peter Schuck, Jacob Hageman
00011  * Adapted by:      Xudong Sun (xudongs@stanford.edu)
00012  *
00013  * Version:
00014  *          v1.0        Mar 01 2010
00015  *          v1.1        Mar 15 2010
00016  *
00017 
00018  * Issues:
00019  *          v1.0
00020  *          1. Designed to deal with a continuous time series, computing result
00021  *          with each pair of neighbouring records.
00022  *          2. Added "cadence" keyword to stride over a minimum fixed period between
00023  *          two records, so the result can be computed at a lower cadence.
00024  *          3. Commented out #infdef USE_INTEL_LIB and DEBUG, seems not properly set up
00025  *          4. Will have to make sure the stack size is big enough (at least ~70MB)
00026  *          run "limit s u" to set unlimited
00027  *          v1.1
00028  *          Updated according to Jacob's new code
00029  *
00030  *
00031  * Example:
00032  *          test_d4vm_old 'in=su_xudong.d4vmin[2004.08.31_23:55:00_TAI/10m]'
00033  */
00034 
00035 
00036 #include <jsoc_main.h>
00037 #include <stdio.h>
00038 #include <stdlib.h>
00039 #include <math.h>
00040 #include <sys/time.h>
00041 #include <sys/resource.h>
00042 
00043 #define PI  (M_PI)
00044 #define DTOR    (PI / 180.)
00045 
00046 #define ARRLENGTH(ARR) (sizeof(ARR) / sizeof(ARR[0]))
00047 #define DIE(msg) {fflush(stdout); fprintf(stderr, "%s, status=%d\n", msg, status); return(status);}
00048 #define SHOW(msg) {printf("%s", msg); fflush(stdout);}
00049 
00050 #include "timing.c"
00051 
00052 
00053 
00054 /* ################## Prototypes ################## */
00055 
00056 extern void d4vm_preproc_(float *bx0, float *bx1, 
00057                           float *by0, float *by1, 
00058                           float *bz0, float *bz1, 
00059                           double *dx, double *dy, double *dt, 
00060                           float *bx, float *bxx, float *bxy, 
00061                           float *by, float *byx, float *byy, 
00062                           float *bz, float *bzx, float *bzy, float *bzt,
00063                           long *xsize, long *ysize);
00064 
00065 void d4vm_kernel_(double *dx, double *dy, 
00066                   long *ksize, float *k_th,
00067                   float *k_x, float *k_y, 
00068                   float *k_xx, float *k_yy, float *k_xy);
00069 
00070 void d4vm_matrix_(float *a, 
00071                   float *bx, float *bxx, float *bxy, 
00072                   float *by, float *byx, float *byy, 
00073                   float *bz, float *bzx, float *bzy, float *bzt, 
00074                   float *k_th, float *k_x, float *k_y,
00075                   float *k_xx, float *k_yy, float *k_xy,
00076                   long *xsize, long *ysize, long *ksize);
00077 
00078 void d4vm_solver_(float *a, 
00079                   float *u0, float *v0, float *w0, 
00080                   float *ux, float *vx, float *wx, 
00081                   float *uy, float *vy, float *wy,
00082                   long *xsize, long *ysize, double *threshold);
00083 
00084 
00085 char *module_name = "test_d4vm";    /* Module name */
00086 
00087 ModuleArgs_t module_args[] =
00088 {
00089     {ARG_STRING, "in", "hmi.sharp_cea_720s", "Input data series."},
00090     {ARG_STRING, "out", "su_xudong.d4vm",  "Output data series."},
00091     {ARG_INT, "ksz0", "19", "Width of window for convolution, must be odd."},
00092     {ARG_INT, "ksz1", "19", "Height of window for convolution, must be odd."},
00093     {ARG_DOUBLE, "dx", "364.43", "X scaling."},
00094     {ARG_DOUBLE, "dy", "364.43", "Y scaling."},
00095     {ARG_DOUBLE, "threshold", "0.", "Threshold for resolving aperture problem."},
00096     {ARG_DOUBLE, "cadence", "240.0", "Minimum stride between two records."},
00097     {ARG_INT, "VERB", "2", "Level of verbosity: 0=errors/warnings; 1=status; 2=all"},
00098     {ARG_END}
00099 };
00100 
00101 
00102 /* ################## Main Module ################## */
00103 
00104 int DoIt(void)
00105 {
00106     int status = DRMS_SUCCESS;
00107     char *inQuery, *outQuery;
00108     DRMS_RecordSet_t *inRS;
00109     int irec, nrecs;
00110     DRMS_Record_t *inRec0, *inRec1, *outRec;
00111     DRMS_Segment_t *inSeg0, *inSeg1;
00112     DRMS_Segment_t *outSegV, *outSegDVDx, *outSegDVDy;
00113     DRMS_Array_t *inArrayX0, *inArrayY0, *inArrayZ0, *inArrayX1, *inArrayY1, *inArrayZ1;
00114     DRMS_Array_t *outArrayV, *outArrayDVDx, *outArrayDVDy;
00115     float *outDataV, *outDataDVDx, *outDataDVDy;
00116     int harpnum;
00117     
00118     float *Bx0, *Bx1, *By0, *By1, *Bz0, *Bz1;               // ancillary
00119     float *Bx, *Bxx, *Bxy, *By, *Byx, *Byy, *Bz, *Bzx, *Bzy, *Bzt;  // processed
00120     float *k_th, *k_x, *k_y, *k_xx, *k_yy, *k_xy;           // convolution kernels
00121     float *a;                               // matrix
00122     float *Vx, *Vy, *Vz;                        // outputs
00123     float *DVxDx, *DVyDx, *DVzDx;
00124     float *DVxDy, *DVyDy, *DVzDy;
00125 
00126     char *inchecklist[] = {"T_REC","T_OBS","HARPNUM"};      // For now
00127     char *outchecklist[] = {"T_REC", "T_REC1", "K0", "K1", "THRESH","HARPNUM"};
00128     DRMS_Record_t *inRec_tmp, *outRec_tmp;
00129     DRMS_Keyword_t *inkeytest, *outkeytest;
00130 
00131     int outDims[3] = {0, 0, 3};
00132     long nx, ny, nxny, ksize[2];    // Seems compatible with long in Fortran?
00133     int verbflag;
00134     double threshold, cadence;
00135     TIME t0, t1;    // Need to check format
00136     double dx, dy, dt;
00137 
00138     int i, l, m, itest;
00139     int nsuccess = 0;
00140     
00141     // Time measuring
00142     double wt0, wt1, wt;
00143     double ut0, ut1, ut;
00144     double st0, st1, st;
00145     double ct0, ct1, ct;
00146     wt0 = getwalltime();
00147     ct0 = getcputime(&ut0, &st0);
00148 
00149     /* Get parameters */
00150     inQuery = (char *)params_get_str(&cmdparams, "in");
00151     outQuery = (char *)params_get_str(&cmdparams, "out");
00152     ksize[0] = params_get_int(&cmdparams, "ksz0");
00153     ksize[1] = params_get_int(&cmdparams, "ksz1");
00154     dx = params_get_double(&cmdparams, "dx");
00155     dy = params_get_double(&cmdparams, "dy");
00156     threshold = params_get_double(&cmdparams, "threshold");
00157     cadence = params_get_double(&cmdparams, "cadence");
00158     verbflag = params_get_int(&cmdparams, "VERB");
00159     for (i = 0; i < 2; i++) {
00160         if (ksize[i] % 2 != 1) {
00161             DIE("Kernel size must be odd");
00162         }
00163     }
00164 
00165     /* Open input */
00166     inRS = drms_open_records(drms_env, inQuery, &status);
00167     if (status || inRS->n == 0) DIE("No input data found");
00168     if (inRS->n < 2) DIE("At least 2 records are needed");
00169     nrecs = inRS->n;
00170 //printf("nrecs=%d\n", nrecs);
00171 
00172     /* Check for essential input keywords */
00173     inRec_tmp = inRS->records[0];
00174     for (itest = 0; itest < ARRLENGTH(inchecklist); itest++) {
00175         inkeytest = hcon_lookup_lower(&inRec_tmp->keywords, inchecklist[itest]);
00176         if (!inkeytest) {
00177             fprintf(stderr, "ERROR: required input keyword %s is missing\n", inchecklist[itest]);
00178             return 0;
00179         }
00180     }
00181 
00182     /* Check for essential output keywords */
00183     outRec_tmp = drms_create_record(drms_env, outQuery, DRMS_TRANSIENT, &status);
00184     if (status) DIE("Unable to open record in output series");
00185     for (itest = 0; itest < ARRLENGTH(outchecklist); itest++) {
00186         outkeytest = hcon_lookup_lower(&outRec_tmp->keywords, outchecklist[itest]);
00187         if (!outkeytest || outkeytest->info->islink || outkeytest->info->recscope == 1) {
00188             fprintf(stderr, "Output keyword %s is missing/constant/a link.\n", outchecklist[itest]);
00189             return 0;
00190         }
00191     }
00192     drms_close_record(outRec_tmp, DRMS_FREE_RECORD);
00193     
00194     /* Read in the first record and save it */
00195     inRec0 = inRS->records[0];
00196     t0 = drms_getkey_time(inRec0, "T_OBS", NULL);   // NEED TO CHECK!!!!!
00197     harpnum = drms_getkey_int(inRec0, "HARPNUM", NULL);
00198     
00199     inSeg0 = drms_segment_lookup(inRec0, "Bp");
00200     inArrayX0 = drms_segment_read(inSeg0, DRMS_TYPE_FLOAT, &status);
00201     Bx0 = (float *)inArrayX0->data;
00202     nx = inArrayX0->axis[0]; ny = inArrayX0->axis[1]; nxny = nx * ny;
00203     outDims[0] = nx; outDims[1] = ny;
00204     inSeg0 = drms_segment_lookup(inRec0, "Bt");
00205     inArrayY0 = drms_segment_read(inSeg0, DRMS_TYPE_FLOAT, &status);
00206     By0 = (float *)inArrayY0->data;
00207     for (int ii = 0; ii < nxny; ii++) By0[ii] *= -1.;
00208     inSeg0 = drms_segment_lookup(inRec0, "Br");
00209     inArrayZ0 = drms_segment_read(inSeg0, DRMS_TYPE_FLOAT, &status);
00210     Bz0 = (float *)inArrayZ0->data;
00211 
00212     /* Do this for record #1 to #nrecs-1 */
00213     for (irec = 1; irec < nrecs; irec++)
00214     {
00215         /* Input record and data */
00216         inRec1 = inRS->records[irec];
00217         t1 = drms_getkey_time(inRec1, "T_OBS", NULL);   // NEED TO CHECK!!!!!
00218         dt = t1 - t0;
00219         if (dt < cadence) {
00220             printf("Cadence higher then limit, skipping current record. \n");
00221             drms_free_array(inArrayX1);
00222             continue;
00223         }
00224         if (harpnum != drms_getkey_time(inRec1, "HARPNUM", NULL)) {
00225             printf("HARP number incorrect, skipping current record. \n");
00226             drms_free_array(inArrayX1);
00227             continue;
00228         }
00229 
00230         inSeg1 = drms_segment_lookup(inRec1, "Bp"); // Single segment
00231         inArrayX1 = drms_segment_read(inSeg1, DRMS_TYPE_FLOAT, &status);
00232 //printf("t0=%f, t1=%f\n", t0, t1);
00233 //printf("dt=%f\n",t1-t0);
00234         Bx1 = (float *)inArrayX1->data;
00235         if (inArrayX1->axis[0] != nx || inArrayX1->axis[1] != ny)
00236             DIE("Wrong data dimension for current record, stop. \n"); 
00237         inSeg1 = drms_segment_lookup(inRec1, "Bt");
00238         inArrayY1 = drms_segment_read(inSeg1, DRMS_TYPE_FLOAT, &status);
00239         By1 = (float *)inArrayY1->data;
00240         for (int ii = 0; ii < nxny; ii++) By1[ii] *= -1.;
00241         inSeg1 = drms_segment_lookup(inRec1, "Br");
00242         inArrayZ1 = drms_segment_read(inSeg1, DRMS_TYPE_FLOAT, &status);
00243         Bz1 = (float *)inArrayZ1->data;
00244         
00245         /* Time measure */
00246         if (verbflag) {
00247             wt1 = getwalltime();
00248             ct1 = getcputime(&ut1, &st1);
00249             printf("processing record %d...\n", nsuccess);
00250         }
00251 
00252         /* Output data */
00253         outArrayV = drms_array_create(DRMS_TYPE_FLOAT, 3, outDims, NULL, &status);
00254         outArrayDVDx = drms_array_create(DRMS_TYPE_FLOAT, 3, outDims, NULL, &status);
00255         outArrayDVDy = drms_array_create(DRMS_TYPE_FLOAT, 3, outDims, NULL, &status);
00256         outDataV = (float *)outArrayV->data;
00257         outDataDVDx = (float *)outArrayDVDx->data;
00258         outDataDVDy = (float *)outArrayDVDy->data;
00259 
00260         /* ======================== */
00261         /* This is the working part.*/
00262         /* ======================== */
00263 
00264         // Preproc
00265         Bx = (float *)calloc(sizeof(float), nxny);
00266         Bxx = (float *)calloc(sizeof(float), nxny);
00267         Bxy = (float *)calloc(sizeof(float), nxny);
00268         By = (float *)calloc(sizeof(float), nxny);
00269         Byx = (float *)calloc(sizeof(float), nxny);
00270         Byy = (float *)calloc(sizeof(float), nxny);
00271         Bz = (float *)calloc(sizeof(float), nxny);
00272         Bzx = (float *)calloc(sizeof(float), nxny);
00273         Bzy = (float *)calloc(sizeof(float), nxny);
00274         Bzt = (float *)calloc(sizeof(float), nxny);
00275         // Kernel
00276         k_th = (float *)calloc(sizeof(float), ksize[0] * ksize[1]);
00277         k_x = (float *)calloc(sizeof(float), ksize[0] * ksize[1]);
00278         k_y = (float *)calloc(sizeof(float), ksize[0] * ksize[1]);
00279         k_xx = (float *)calloc(sizeof(float), ksize[0] * ksize[1]);
00280         k_yy = (float *)calloc(sizeof(float), ksize[0] * ksize[1]);
00281         k_xy = (float *)calloc(sizeof(float), ksize[0] * ksize[1]);
00282         // Matrix
00283         a = (float *)calloc(sizeof(float), 10 * 10 * nxny);
00284         // Output
00285         Vx = outDataV; Vy = outDataV + nxny; Vz = outDataV + nxny * 2;
00286         DVxDx = outDataDVDx; DVyDx = outDataDVDx + nxny; DVzDx = outDataDVDx + nxny * 2;
00287         DVxDy = outDataDVDy; DVyDy = outDataDVDy + nxny; DVzDy = outDataDVDy + nxny * 2;
00288 //SHOW("111\n");
00289         // Preprocess data
00290         d4vm_preproc_(Bx0, Bx1, By0, By1, Bz0, Bz1, 
00291                       &dx, &dy, &dt, 
00292                       Bx, Bxx, Bxy, By, Byx, Byy, Bz, Bzx, Bzy, Bzt,
00293                       &nx, &ny);
00294 //SHOW("112\n");
00295         // Generate kernel
00296         d4vm_kernel_(&dx, &dy, ksize,
00297                      k_th, k_x, k_y, k_xx, k_yy, k_xy);
00298 //SHOW("113\n");
00299 
00300         // Process data
00301         d4vm_matrix_(a, 
00302                      Bx, Bxx, Bxy, By, Byx, Byy, Bz, Bzx, Bzy, Bzt,
00303                      k_th, k_x, k_y, k_xx, k_yy, k_xy,
00304                      &nx, &ny, ksize);
00305 //SHOW("114\n");
00306         // Solve matrix
00307         d4vm_solver_(a, 
00308                      Vx, Vy, Vz, 
00309                      DVxDx, DVyDx, DVzDx, 
00310                      DVxDy, DVyDy, DVzDy, 
00311                      &nx, &ny, &threshold);
00312 //SHOW("115\n");
00313 
00314         /* Output record, create one at a time */
00315         outRec = drms_create_record(drms_env, outQuery, DRMS_PERMANENT, &status);
00316         if (status) DIE("Output record not created");
00317         outSegV = drms_segment_lookupnum(outRec, 0);
00318         outSegDVDx = drms_segment_lookupnum(outRec, 1);
00319         outSegDVDy = drms_segment_lookupnum(outRec, 2);
00320         for (i = 0; i < 3; i++) {
00321             outSegV->axis[i] = outArrayV->axis[i];  // For variable dimensions
00322             outSegDVDx->axis[i] = outArrayDVDx->axis[i];
00323             outSegDVDy->axis[i] = outArrayDVDy->axis[i];
00324         }
00325         outArrayV->israw = 0;
00326         outArrayV->bzero = 0.0;
00327         outArrayV->bscale = 0.0001;
00328         outArrayDVDx->israw = 0;
00329         outArrayDVDx->bzero = 0.0;
00330         outArrayDVDx->bscale = 0.000001;
00331         outArrayDVDy->israw = 0;
00332         outArrayDVDy->bzero = 0.0;
00333         outArrayDVDy->bscale = 0.000001;
00334 
00335         /* Set keywords */
00336         drms_copykey(outRec, inRec0, "HARPNUM");
00337         drms_copykey(outRec, inRec0, "T_REC");
00338         drms_copykey(outRec, inRec0, "T_OBS");
00339         drms_setkey_time(outRec, "T_REC1", drms_getkey_time(inRec1, "T_REC", NULL));
00340         drms_setkey_time(outRec, "T_OBS1", drms_getkey_time(inRec1, "T_OBS", NULL));
00341         drms_setkey_int(outRec, "K0", ksize[0]);
00342         drms_setkey_int(outRec, "K1", ksize[1]);
00343         drms_setkey_double(outRec, "THRESH", threshold);
00344 
00345         /* Result writing */
00346         status = drms_segment_write(outSegV, outArrayV, 0);
00347         if (status) DIE("Problem writing file V");
00348         status = drms_segment_write(outSegDVDx, outArrayDVDx, 0);
00349         if (status) DIE("Problem writing file DVDx"); 
00350         status = drms_segment_write(outSegDVDy, outArrayDVDy, 0);
00351         if (status) DIE("Problem writing file DVDy"); 
00352         drms_close_record(outRec, DRMS_INSERT_RECORD);
00353         drms_free_array(outArrayV);
00354         drms_free_array(outArrayDVDx);
00355         drms_free_array(outArrayDVDy);
00356 
00357         /* Clean up, move to next record */
00358         drms_free_array(inArrayX0);
00359         drms_free_array(inArrayY0);
00360         drms_free_array(inArrayZ0);
00361         free(Bx); free(Bxx); free(Bxy);
00362         free(By); free(Byx); free(Byy);
00363         free(Bz); free(Bzx); free(Bzy); free(Bzt);
00364         free(a);
00365         free(k_th); free(k_x); free(k_y);
00366         free(k_xx); free(k_yy); free(k_xy);
00367         
00368         inRec0 = inRec1;
00369         inSeg0 = inSeg1;
00370         inArrayX0 = inArrayX1;
00371         inArrayY0 = inArrayY1;
00372         inArrayZ0 = inArrayZ1;
00373         Bx0 = Bx1;
00374         By0 = By1;
00375         Bz0 = Bz1;
00376         t0 = t1;
00377         
00378         /* Time measure */
00379         if (verbflag) {
00380             wt = getwalltime();
00381             ct = getcputime(&ut, &st);
00382             printf("record %d done, %.2f ms wall time, %.2f ms cpu time\n", 
00383                      nsuccess, wt - wt1, ct - ct1);
00384         }
00385 
00386         nsuccess++;
00387     }
00388 
00389     drms_free_array(inArrayX1);
00390     drms_free_array(inArrayY1);
00391     drms_free_array(inArrayZ1);
00392     drms_close_records(inRS, DRMS_FREE_RECORD);
00393     return(DRMS_SUCCESS);
00394 }
proj/mag/d4vm/apps/test_d4vm.c
