JSOC_Documentation: proj/mag/nlfff/apps/relax.c Source File

00001 /*
00002  * Relaxation step for NLFFF model code
00003  * Adapted from Wiegelmann's relax1.c
00004  * Takes in initial guess and model parameters
00005  * Computes B
00006  * Adapted by X. Sun
00007  * Version 1.0 Feb 18 2010
00008  * Version 1.1 Mar 17 2010
00009  *  Took out green(), removed B_PF
00010  * Version 2.0 Mar 26 2010
00011  *  Move boundary substitution out to main
00012  *  Move macros to main
00013  *  Work with multigrid now
00014  */
00015 
00016 #ifndef PI
00017 #include <math.h>
00018 #define PI (M_PI)
00019 #endif 
00020 
00021 /* Macros */
00022 #ifndef Macro
00023 #define Macro
00024 #define MCENTERGRAD(f,id) ((f[i+id]-f[i-id])/(2*h))
00025 #define MLEFTGRAD(f,id)   ((-3*f[i]+4*f[i+id]-f[i+2*id])/(2*h))
00026 #define MRIGHTGRAD(f,id)  ((+3*f[i]-4*f[i-id]+f[i-2*id])/(2*h))
00027 #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))))
00028 #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))))
00029 #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))))
00030 #endif
00031 
00032 /*===========================================*/
00033 /*=============== Prototypes ================*/
00034 /*===========================================*/
00035 
00036 
00037 /* Originally in relax.c */
00038 double zeit();
00039 
00040 /* Originally in optimization.c */
00041 double calculateL(double *Bx, double *By, double *Bz, 
00042                   float *DivB, 
00043                   float *oxa, float *oya, float *oza, 
00044                   float *oxb, float *oyb, float *ozb, 
00045                   float *oxbx, float *oxby, float *oxbz, 
00046                   float *oxjx, float *oxjy, float *oxjz, 
00047                   float *odotb, 
00048                   double *wa0, double *wb0, 
00049                   float *wax, float *way, float *waz, 
00050                   float *wbx, float *wby, float *wbz, 
00051                   float *Fx, float *Fy, float *Fz, 
00052                   int nx, int ny, int nz, int nynz, 
00053                   double dx, double dy, double dz);
00054                   
00055 double calculateLi(double *Bx, double *By, double *Bz, 
00056                    float *oxa, float *oya, float *oza, 
00057                    float *oxb, float *oyb, float *ozb, 
00058                    int nx, int ny, int nz, int nynz, int nd, 
00059                    double dx, double dy, double dz);
00060 
00061 
00062 /*=======================================================*/
00063 /*===================== NLFFF Main ======================*/
00064 /*=======================================================*/
00065 
00066 void relax(double *Bx, double *By, double *Bz, 
00067            int nx, int ny, int nz, 
00068            int nd, int maxit, int verb)
00069 {
00070     /* Global Variables */
00071   //FILE *streamw, *initfile;
00072   //double Cmn[102][102], lmn[102][102], rmn[102][102];
00073     double Lx, Ly, Lz;  //, alpha, pi, q1;
00074   //double *t, *xt, *yt, *zt, *absB, *t2, *xyzB, *refB;
00075   //double *x, *y, *z;  // unused
00076   //double *Bx, *By, *Bz;  // ---- PASSED IN
00077     /* double *B; */
00078     double *Bx1, *By1, *Bz1, *Bx2, *By2, *Bz2;
00079   //float *VxBx, *VxBy, *VxBz, *Rotvx, *Rotvy, *Rotvz;  // ---- COMMENTED OUT IN RELAX1
00080   //float *Jx, *Jy, *Jz, *Vx, *Vy, *Vz; // ---- COMMENTED OUT IN RELAX1
00081   //float *Jx1, *Jy1, *Jz1; // ---- COMMENTED OUT IN RELAX1
00082     float *ox, *oy, *oz, *oxbx, *oxby, *oxbz, *odotb, *DivB;
00083     float *oxjx, *oxjy, *oxjz, *Fx, *Fy, *Fz;
00084     /* double *Bxana, *Byana, *Bzana; */
00085     float *oxa, *oya, *oza, *oxb, *oyb, *ozb; /*, *ux, *uy, *uz, *u; */
00086     double *wa0, *wb0;
00087     float *wax, *way, *waz, *wbx, *wby, *wbz;
00088   //double dt, length;
00089   //double C;
00090     double dx, dy, dz, h;
00091     double mue;     //, mu0, rmu0, gausT, rmu0gausT, Ls;
00092   //double Lxyave;  //, Bave;
00093   //int maxm, maxn;
00094   //int maxsteps, max2;
00095   //int nx, ny, nz; // ---- PASSED IN
00096     int nynz, nxnynz, nxmax, nymax, nzmax;  //, nd ---- PASSED IN
00097   //int calcb, maxit;   //, firsttime, m0, n0;
00098     int diagstep;//, nxnynz1, nxnynz2;  // not used
00099     
00100     int i, it, ix, iy, iz;
00101   //int zaehlv; // not used
00102 //    double dummy, dummy1, dummy2;
00103     double L, Li;   //, nu;
00104     double time1;//, timetot;
00105     double nave, Bave;  //, mue2;
00106     double oldL;    //, w2, w, Llimit;
00107     double prevL, newL, gradL;
00108     float *prof, *xprof, *yprof, *zprof;
00109     int statcount;
00110 //    char leer[25];
00111     
00112     L = oldL = 0.0; // = w= w2 = Llimit = 0.0;
00113     statcount = 0;
00114 
00115 /*    
00116     timetot = 0.0;
00117     q1 = 1.0;
00118     mu0 = 1.2566e-6;
00119     rmu0 = 1. / mu0;
00120     Ls = 1.0;   // Length in Mm
00121     gausT = 0.0001;  // Gauss to Tesla
00122     rmu0gausT = rmu0 * gausT;
00123  */
00124 
00125     nynz = ny * nz;
00126     nxnynz = nx * ny * nz;
00127 //  nxnynz1 = (nx - 0) * (ny - 0) * (nz - 0);   // ???
00128 //  nxnynz2 = (nx - 0) * (ny - 0) * (nz - 0);
00129 //  zaehlv = 0;     // not used
00130 
00131     nxmax = nx; nymax = ny; nzmax = nz;
00132     nave = sqrt(1.0 * (nx - 2 * nd - 1) * (ny - 2 * nd - 1));
00133 
00134     Lx = 1.0 * (nx - 1) / nave;
00135     Ly = 1.0 * (ny - 1) / nave;
00136     Lz = 1.0 * (nz - 1) / nave;
00137     dx = dy = dz = Lx / (nx - 1);
00138     
00139 /* === Unused ===
00140     x = (double *) calloc(nx, sizeof(double));
00141     y = (double *) calloc(ny, sizeof(double));
00142     z = (double *) calloc(nz, sizeof(double));
00143  */
00144 /* ======== Bxyz PASSED IN AS PARAMETER ========
00145     Bx = (double *) calloc(nxnynz, sizeof(double));
00146     By = (double *) calloc(nxnynz, sizeof(double));
00147     Bz = (double *) calloc(nxnynz, sizeof(double));
00148  */
00149 
00150     Bx1 = (double *) calloc(nxnynz, sizeof(double));
00151     By1 = (double *) calloc(nxnynz, sizeof(double));
00152     Bz1 = (double *) calloc(nxnynz, sizeof(double));
00153     Bx2 = (double *) calloc(nxnynz, sizeof(double));
00154     By2 = (double *) calloc(nxnynz, sizeof(double));
00155     Bz2 = (double *) calloc(nxnynz, sizeof(double));
00156     ox = (float *) calloc(nxnynz, sizeof(float));
00157     oy = (float *) calloc(nxnynz, sizeof(float));
00158     oz = (float *) calloc(nxnynz, sizeof(float));
00159     oxbx = (float *) calloc(nxnynz, sizeof(float));
00160     oxby = (float *) calloc(nxnynz, sizeof(float));
00161     oxbz = (float *) calloc(nxnynz, sizeof(float));
00162     odotb = (float *) calloc(nxnynz, sizeof(float));
00163     DivB = (float *) calloc(nxnynz, sizeof(float));
00164     oxjx = (float *) calloc(nxnynz, sizeof(float));
00165     oxjy = (float *) calloc(nxnynz, sizeof(float));
00166     oxjz = (float *) calloc(nxnynz, sizeof(float));
00167     Fx = (float *) calloc(nxnynz, sizeof(float));
00168     Fy = (float *) calloc(nxnynz, sizeof(float));
00169     Fz = (float *) calloc(nxnynz, sizeof(float));
00170     oxa = (float *) calloc(nxnynz, sizeof(float));
00171     oya = (float *) calloc(nxnynz, sizeof(float));
00172     oza = (float *) calloc(nxnynz, sizeof(float));
00173     oxb = (float *) calloc(nxnynz, sizeof(float));
00174     oyb = (float *) calloc(nxnynz, sizeof(float));
00175     ozb = (float *) calloc(nxnynz, sizeof(float));
00176     wax = (float *) calloc(nxnynz, sizeof(float));
00177     way = (float *) calloc(nxnynz, sizeof(float));
00178     waz = (float *) calloc(nxnynz, sizeof(float));
00179     wa0 = (double *) calloc(nxnynz, sizeof(double));
00180     wbx = (float *) calloc(nxnynz, sizeof(float));
00181     wby = (float *) calloc(nxnynz, sizeof(float));
00182     wbz = (float *) calloc(nxnynz, sizeof(float));
00183     wb0 = (double *) calloc(nxnynz, sizeof(double));
00184 
00185 #ifdef _OPENMP
00186 #pragma omp parallel for private (i)
00187 #endif
00188     for (i = 0; i < nxnynz; i++) {
00189         ox[i] = oy[i] = oz[i] = oxbx[i] = oxby[i] = oxbz[i] = odotb[i] = DivB[i] = 0.0;
00190         oxjx[i] = oxjy[i] = oxjz[i] = Fx[i] = Fy[i] = Fz[i] = 0.0;
00191         // Copy Bxyz_PF to Bxyz, originally "load B0.bin"
00192 //        Bx[i] = Bx_PF[i]; By[i] = By_PF[i]; Bz[i] = Bz_PF[i]; // Now Done in wrapper 
00193     }
00194     // alpha = 0.0; // only for linear force free field
00195 
00196     // Substitute bottom with vector magnetograms, "loadvectormag()"
00197 /*
00198     for (iy = 0; iy < ny; iy++)
00199     for (ix = 0; ix < nx; ix++) {
00200         i = nynz * ix + nz * iy;    // iz = 0
00201         Bx[i] = Bx0[ix * ny + iy];
00202         By[i] = By0[ix * ny + iy];
00203         Bz[i] = Bz0[ix * ny + iy];
00204     }
00205     if (verb) printf("PF bottom layer replaced with vector magnetogram\n");
00206 */
00207 /*
00208 printf("\n********************************************\n");
00209 printf("Bz[0]=%f, Bz[72]=%f\n", Bz[0], Bz[72]);
00210 printf("********************************************\n");
00211 */
00212 
00213     iz = 0; Bave = 0.0;
00214     for (ix = nd; ix < nx - nd; ix++)
00215     for (iy = nd; iy < ny - nd; iy++) {
00216         i = nynz * ix + nz * iy + iz;
00217         Bave = Bave + sqrt(Bx[i] * Bx[i] + By[i] * By[i] + Bz[i]*Bz[i]);
00218     }
00219     Bave = Bave / (nx - 2 * nd) / (ny - 2 * nd);
00220 
00221 #ifdef _OPENMP
00222 #pragma omp parallel for private (i)
00223 #endif
00224     for (i = 0; i < nxnynz; i++) {
00225         Bx[i] = Bx[i] / Bave;
00226         By[i] = By[i] / Bave;
00227         Bz[i] = Bz[i] / Bave;
00228     }
00229 
00230 #ifdef _OPENMP
00231 #pragma omp parallel for private (i)
00232 #endif
00233     for (i = 0; i < nxnynz; i++) {
00234         Bx1[i] = Bx2[i] = Bx[i];
00235         By1[i] = By2[i] = By[i];
00236         Bz1[i] = Bz2[i] = Bz[i];
00237 
00238     }
00239 
00240     /* ================= Weighting Function ================= */
00241 
00242     h = dx;  // FOR GRAD!
00243 
00244     if (verb) printf("Calculate wa and wb\n");
00245 
00246 #ifdef _OPENMP
00247 #pragma omp parallel for private (i)
00248 #endif
00249     for (i = 0; i < nxnynz; i++) {
00250         wa0[i] = wb0[i] = 1.0;
00251     }
00252 
00253     if (nd > 1) {
00254         prof = (float*) calloc(nd, sizeof(float));
00255         xprof = (float *) calloc(nx, sizeof(float));
00256         yprof = (float *) calloc(ny, sizeof(float));
00257         zprof = (float *) calloc(nz, sizeof(float));
00258         for (i = 0; i < nd; i++) prof[i] = 0.5 + 0.5 * cos(1.0 * i / (nd - 1) * PI);
00259         for (ix = 0; ix < nx; ix++) xprof[ix] = 1.0;
00260         for (iy = 0; iy < ny; iy++) yprof[iy] = 1.0;
00261         for (iz = 0; iz < nz; iz++) zprof[iz] = 1.0;
00262         for (i = 0; i < nd; i++) {
00263             xprof[nx - nd + i] = prof[i];
00264             xprof[nd - i - 1] = prof[i]; 
00265             yprof[ny - nd + i] = prof[i];
00266             yprof[nd - i - 1] = prof[i];
00267             zprof[nz - nd + i] = prof[i];
00268         }
00269 
00270 #ifdef _OPENMP
00271 #pragma omp parallel for private (i,ix,iy,iz)
00272 #endif
00273         for (ix = 0; ix < nx; ix++)
00274         for (iy = 0; iy < ny; iy++)
00275         for (iz = 0; iz < nz; iz++) {
00276             i = ix * nynz + iy * nz + iz;
00277             wa0[i] = xprof[ix] * yprof[iy] * zprof[iz];
00278             wb0[i] = xprof[ix] * yprof[iy] * zprof[iz];
00279         }
00280     }
00281 
00282 #ifdef _OPENMP
00283 #pragma omp parallel for private (i)
00284 #endif
00285     for (i = 0; i < nxnynz; i++) {
00286         wax[i] = GRADX(wa0,i);
00287         way[i] = GRADY(wa0,i);
00288         waz[i] = GRADZ(wa0,i);
00289         wbx[i] = GRADX(wb0,i);
00290         wby[i] = GRADY(wb0,i);
00291         wbz[i] = GRADZ(wb0,i);
00292     }
00293 
00294 //  nxnynz1 = (nx - 2 * nd) * (ny - 2 * nd) * (nz - nd);
00295   //printf("\n Boundary layer nd= %i ",nd);
00296   //printf("\n Inner region = %.2f",(double)(nxnynz1)/(nxnynz));
00297 
00298     /* ================= Relaxation Loop ================= */
00299 
00300     it = -1;
00301     mue = 0.1 * dx * dx;
00302 
00303     // Added by xudong mar 18 2010
00304     // When L > oldL we need to back up one step, B, B1, B2 all get restored
00305     // and L gets recomputed in the next step, using the restored B
00306     // Ideally we'll get L = oldL, but rounding errors may cause this to fail
00307     // Added a restore flag
00308     int restore = 0;
00309 
00310     while (it < maxit && statcount < 10 && mue > (1.0e-7 * dx * dx)) 
00311     {
00312         it++;
00313         L = calculateL(Bx, By, Bz, 
00314                        DivB, 
00315                        oxa, oya, oza, 
00316                        oxb, oyb, ozb, 
00317                        oxbx, oxby, oxbz, 
00318                        oxjx, oxjy, oxjz, 
00319                        odotb, 
00320                        wa0, wb0, 
00321                        wax, way, waz, 
00322                        wbx, wby, wbz, 
00323                        Fx, Fy, Fz, 
00324                        nx, ny, nz, nynz, 
00325                        dx, dy, dz);
00326 
00327         if (it == 0)
00328             oldL = L;
00329 
00330         if (restore) L = oldL;
00331 
00332         if (it > 0 && L > oldL) {
00333             restore = 1;
00334             mue = mue / 2.0;
00335             if (verb) {
00336                 printf("mue reduced, mue= %lf \t mue/dx^2 = %lf\n", mue, mue / (dx * dx));
00337                 printf("oldL= %lf \t L=%lf\n", oldL, L);
00338             }
00339             it--;
00340 #ifdef _OPENMP
00341 #pragma omp parallel for private (i)
00342 #endif
00343             for (i = 0; i < nxnynz; i++) {  // step too large, refused
00344                 Bx[i] = Bx1[i] = Bx2[i];
00345                 By[i] = By1[i] = By2[i];
00346                 Bz[i] = Bz1[i] = Bz2[i];
00347             }
00348         } else {
00349             restore = 0;
00350             mue = mue * 1.01;
00351             oldL = L;
00352         }
00353         
00354 //        if (oldL >= L) {
00355         if (restore == 0) {
00356         /* L minimize */
00357           //nu = 0.1e12;
00358 #ifdef _OPENMP
00359 #pragma omp parallel for private (i,ix,iy,iz)
00360 #endif
00361             for (ix = 1; ix < nx - 1; ix++)
00362             for (iy = 1; iy < ny - 1; iy++)
00363             for (iz = 1; iz < nz - 1; iz++) {
00364                 i = ix * nynz + iy * nz + iz;
00365                 Bx1[i] = Bx[i] + mue * Fx[i];
00366                 By1[i] = By[i] + mue * Fy[i];
00367                 Bz1[i] = Bz[i] + mue * Fz[i];
00368             }
00369         }
00370         
00371         /* Some tests with (Anti) Symmetric Boundary-conditions */
00372         diagstep = 10;
00373         if (it % diagstep == 0) /* diagnostic */
00374         {
00375             time1 = zeit();
00376           //if (time1 > 0.0) timetot += time1;
00377             Li = calculateLi(Bx, By, Bz, 
00378                              oxa, oya, oza, 
00379                              oxb, oyb, ozb, 
00380                              nx, ny, nz, nynz, nd, 
00381                              dx, dy, dz);
00382             if (verb) printf("%i L= %.4f, L_i= %.4f,  cpu %.1f s", it, L, Li, time1);
00383          // dummy = 0.0;
00384          // dummy1 = dummy2 = 0.0;
00385             /* NEW: calc gradient (dL/dt)L for stopping rule */
00386             if (it==0) {
00387                 prevL = 2.0 * L;
00388                 newL = L;
00389             } else {
00390                 prevL = newL;
00391                 newL = L;
00392             }
00393             gradL = fabs((newL - prevL) / newL);
00394             if (verb) printf(",  gradL/L= %lf\n", gradL);
00395             if (gradL < 0.0001) {
00396                 statcount++;
00397                 if (verb) printf("*** STATIONARY STATE count: %i *** grad L/L= %lf \n", statcount, gradL);
00398             }
00399             if (gradL > 0.0001)
00400                 statcount = 0;
00401          // dummy1 = 0.0;
00402          // dummy2 = 1.0;
00403          // dummy = sqrt(dummy1 / dummy2);
00404         }   /* End Diagnostic */
00405 
00406 #ifdef _OPENMP
00407 #pragma omp parallel for private (i)
00408 #endif
00409         for (i = 0; i < nxnynz; i++) {
00410             Bx2[i] = Bx[i];
00411             By2[i] = By[i];
00412             Bz2[i] = Bz[i];
00413         }
00414 #ifdef _OPENMP
00415 #pragma omp parallel for private (i)
00416 #endif
00417         for (i = 0; i < nxnynz; i++) {
00418             Bx[i] = Bx1[i];
00419             By[i] = By1[i];
00420             Bz[i] = Bz1[i];
00421         }
00422     }   // while
00423 
00424     /* Normalize B */
00425 #ifdef _OPENMP
00426 #pragma omp parallel for private (i)
00427 #endif
00428     for (i = 0; i < nxnynz; i++) {
00429         Bx[i] = Bx[i] * Bave;
00430         By[i] = By[i] * Bave;
00431         Bz[i] = Bz[i] * Bave;
00432     }
00433 
00434     /* Clean up */
00435     free(Bx1); free(By1); free(Bz1);
00436     free(Bx2); free(By2); free(Bz2);
00437     free(ox); free(oy); free(oz);
00438     free(oxbx); free(oxby); free(oxbz);
00439     free(odotb);
00440     free(oxjx); free(oxjy); free(oxjz);
00441     free(Fx); free(Fy); free(Fz);
00442     free(oxa); free(oya); free(oza);
00443     free(oxb); free(oyb); free(ozb);
00444     free(wa0); free(wb0);
00445     free(wax); free(way); free(waz);
00446     free(wbx); free(wby); free(wbz);
00447     if (nd > 1) {
00448         free(prof); free(xprof); free(yprof); free(zprof);
00449     }
00450 }
00451 
00452 
00453 /* ========================================== */
00454 /* ====== Routine fuer CPU Zeitmessung ====== */
00455 /* ========================================== */
00456 
00457 double zeit()
00458 {
00459     static double tima = 0.;
00460     double tim, t;
00461     tim = (double)clock() / CLOCKS_PER_SEC;
00462     t = tim - tima;
00463     tima = tim;
00464     return t;
00465 }
00466 
00467 
00468 
00469 /* ========================================== */
00470 /* ============== Optimization ============== */
00471 /* ========================================== */
00472 
00473 double calculateL(double *Bx, double *By, double *Bz, 
00474                   float *DivB, 
00475                   float *oxa, float *oya, float *oza, 
00476                   float *oxb, float *oyb, float *ozb, 
00477                   float *oxbx, float *oxby, float *oxbz, 
00478                   float *oxjx, float *oxjy, float *oxjz, 
00479                   float *odotb, 
00480                   double *wa0, double *wb0, 
00481                   float *wax, float *way, float *waz, 
00482                   float *wbx, float *wby, float *wbz, 
00483                   float *Fx, float *Fy, float *Fz, 
00484                   int nx, int ny, int nz, int nynz,
00485                   double dx, double dy, double dz)
00486 {
00487     double h;
00488     double bx, by, bz, cbx, cby, cbz, fx, fy, fz;
00489     double divB, b2, helpL;//, helpf;
00490     double o2a, o2b, term1x, term2x, term3x, term4x, term5ax, term5bx;
00491     double term1y, term2y, term3y, term4y, term5ay, term5by;
00492     double term1z, term2z, term3z, term4z, term5az, term5bz;
00493     double term6x, term6y, term6z, term7x, term7y, term7z;
00494 //    double dummy1, dummy2, dummy3, dummy4, dummy5, dummy6;
00495     int ix, iy, iz, i;
00496 //    int zaehl;
00497 
00498     h = dx;
00499     helpL = 0.0;
00500 //    helpf = 0.0;
00501 //    zaehl = 0;
00502 
00503 #ifdef _OPENMP
00504 #pragma omp parallel for private (i,ix,iy,iz,bx,by,bz,cbx,cby,cbz,fx,fy,fz,divB,b2,o2a,o2b) reduction(+:helpL)
00505 #endif
00506     for (ix = 0; ix < nx; ix++)
00507     for (iy = 0; iy < ny; iy++)
00508     for (iz = 0; iz < nz; iz++)
00509     {
00510         i = ix * nynz + iy * nz + iz;
00511         bx = Bx[i];
00512         by = By[i];
00513         bz = Bz[i];
00514         b2 = (bx * bx + by * by + bz * bz);
00515 
00516         cbx = GRADY(Bz,i) - GRADZ(By,i);
00517         cby = GRADZ(Bx,i) - GRADX(Bz,i);
00518         cbz = GRADX(By,i) - GRADY(Bx,i);
00519 
00520         fx = cby * bz - cbz * by;
00521         fy = cbz * bx - cbx * bz;
00522         fz = cbx * by - cby * bx;
00523         divB = GRADX(Bx,i) + GRADY(By,i) + GRADZ(Bz,i);
00524         DivB[i] = divB;
00525 
00526         if (b2 > 0.0) {
00527             oxa[i] = (1.0 / b2) * fx;
00528             oya[i] = (1.0 / b2) * fy;
00529             oza[i] = (1.0 / b2) * fz;
00530             oxb[i] = (1.0 / b2) * (divB * bx);
00531             oyb[i] = (1.0 / b2) * (divB * by);
00532             ozb[i] = (1.0 / b2) * (divB * bz);
00533             o2a = oxa[i] * oxa[i] + oya[i] * oya[i] + oza[i] * oza[i];
00534             o2b = oxb[i] * oxb[i] + oyb[i] * oyb[i] + ozb[i] * ozb[i];
00535             helpL = helpL + (wa0[i] * b2 * o2a + wb0[i] * b2 * o2b);
00536         }
00537 
00538         oxbx[i] = oya[i] * bz - oza[i] * by;
00539         oxby[i] = oza[i] * bx - oxa[i] * bz;
00540         oxbz[i] = oxa[i] * by - oya[i] * bx;
00541         odotb[i] = oxb[i] * bx + oyb[i] * by + ozb[i] * bz;
00542         oxjx[i] = oya[i] * cbz - oza[i] * cby;
00543         oxjy[i] = oza[i] * cbx - oxa[i] * cbz;
00544         oxjz[i] = oxa[i] * cby - oya[i] * cbx;
00545     }   // ix, iy, iz
00546 
00547     helpL = helpL * (dx * dy * dz);
00548     
00549     /*10.07.02 diagnostik */
00550 //    dummy1 = dummy2 = dummy3 = dummy4 = dummy5 = dummy6 = 0.0;
00551 #ifdef _OPENMP
00552 #pragma omp parallel for private (i, ix, iy, iz, o2a, o2b, term1x, term2x, term3x, term4x, term5ax, term5bx, term6x, term7x, term1y, term2y, term3y, term4y, term5ay, term5by, term6y, term7y, term1z, term2z, term3z, term4z, term5az, term5bz, term6z, term7z)
00553 #endif
00554     for (ix = 0; ix < nx; ix++)
00555     for (iy = 0; iy < ny; iy++)
00556     for (iz = 0; iz < nz; iz++)
00557     {
00558         i = ix * nynz + iy * nz + iz;
00559         term1x = GRADY(oxbz,i) - GRADZ(oxby,i);
00560         term1y = GRADZ(oxbx,i) - GRADX(oxbz,i);
00561         term1z = GRADX(oxby,i) - GRADY(oxbx,i);
00562 
00563         term2x = oxjx[i];
00564         term2y = oxjy[i];
00565         term2z = oxjz[i];
00566 
00567         term3x = GRADX(odotb,i);
00568         term3y = GRADY(odotb,i);
00569         term3z = GRADZ(odotb,i);
00570 
00571         term4x = oxb[i] * DivB[i];
00572         term4y = oyb[i] * DivB[i];
00573         term4z = ozb[i] * DivB[i];
00574 
00575         o2a = oxa[i] * oxa[i] + oya[i] * oya[i] + oza[i] * oza[i];
00576         o2b = oxb[i] * oxb[i] + oyb[i] * oyb[i] + ozb[i] * ozb[i];
00577         term5ax = Bx[i] * o2a;
00578         term5ay = By[i] * o2a;
00579         term5az = Bz[i] * o2a;
00580         term5bx = Bx[i] * o2b;
00581         term5by = By[i] * o2b;
00582         term5bz = Bz[i] * o2b;
00583 
00584         /* Terms regarding weighting function */
00585         term6x = oxby[i] * waz[i] - oxbz[i] * way[i];
00586         term6y = oxbz[i] * wax[i] - oxbx[i] * waz[i];
00587         term6z = oxbx[i] * way[i] - oxby[i] * wax[i];
00588         term7x = odotb[i] * wbx[i];
00589         term7y = odotb[i] * wby[i];
00590         term7z = odotb[i] * wbz[i];
00591 
00592         Fx[i] = wa0[i] * (term1x - term2x + term5ax) + wb0[i] * (term3x - term4x + term5bx) + term6x + term7x;
00593         Fy[i] = wa0[i] * (term1y - term2y + term5ay) + wb0[i] * (term3y - term4y + term5by) + term6y + term7y;
00594         Fz[i] = wa0[i] * (term1z - term2z + term5az) + wb0[i] * (term3z - term4z + term5bz) + term6z + term7z;
00595     }
00596 
00597     return helpL; 
00598 }
00599 
00600 
00601 /* ========================================== */
00602 /* ======== Optimization (inner box) ======== */
00603 /* ========================================== */
00604 
00605 double calculateLi(double *Bx, double *By, double *Bz, 
00606                    float *oxa, float *oya, float *oza, 
00607                    float *oxb, float *oyb, float *ozb, 
00608                    int nx, int ny, int nz, int nynz, int nd, 
00609                    double dx, double dy, double dz)
00610 {
00611     double h;
00612     double bx, by, bz;
00613     double b2, helpL;
00614     double o2a, o2b, o2;
00615     int ix, iy, iz, i;
00616 
00617     h = dx;
00618     helpL = 0.0;
00619 
00620 #ifdef _OPENMP
00621 #pragma omp parallel for private (i,ix,iy,iz,bx,by,bz,b2,o2a,o2b,o2) reduction(+:helpL)
00622 #endif
00623     for (ix = nd; ix < nx - nd; ix++)
00624     for (iy = nd; iy < ny - nd; iy++)
00625     for (iz = 0; iz < nz - nd; iz++)
00626     {
00627         i = ix * nynz + iy * nz + iz;
00628         bx = Bx[i];
00629         by = By[i];
00630         bz = Bz[i];
00631         b2 = (bx * bx + by * by + bz * bz);
00632         o2a = oxa[i] * oxa[i] + oya[i] * oya[i] + oza[i] * oza[i];
00633         o2b = oxb[i] * oxb[i] + oyb[i] * oyb[i] + ozb[i] * ozb[i];
00634         o2 = o2a + o2b;
00635 
00636         helpL = helpL + b2*o2;
00637     }
00638 
00639     helpL = helpL * (dx * dy * dz);
00640     return helpL;
00641 }