Athena: microphysics/resistivity.c Source File

Go to the documentation of this file.
00001 #include "../copyright.h"
00002 /*============================================================================*/
00003 /*! \file resistivity.c
00004  *  \brief Adds explicit resistivity terms to the induction and energy eqns
00005  *
00006  * PURPOSE: Adds explicit resistivity terms to the induction and energy eqns,
00007  *  -   dB/dt = -Curl(E)
00008  *  -   dE/dt = Div(B X E)
00009  *   where 
00010  *  -   E = eta_Ohm J + eta_Hall(J X B)/B + eta_AD J_perp = (emf) 
00011  *  -      J = Curl(B) = current density
00012  *  -      eta_Ohm = Ohmic resistivity
00013  *  -      eta_Hall = Hall diffusion coefficient
00014  *  -      eta_AD = ambipolar diffusion coefficient
00015  *   The induction equation is updated using CT to keep div(B)=0.  The total
00016  *   electric field (resistive EMF) is computed from calls to the EField_*
00017  *   functions.
00018  *
00019  * CONTAINS PUBLIC FUNCTIONS:
00020  *- resistivity() - updates induction and energy eqns with resistive term.
00021  *- resistivity_init() - allocates memory needed
00022  *- resistivity_destruct() - frees memory used                                */
00023 /*============================================================================*/
00024 
00025 #include <math.h>
00026 #include <float.h>
00027 #include "../defs.h"
00028 #include "../athena.h"
00029 #include "../globals.h"
00030 #include "prototypes.h"
00031 #include "../prototypes.h"
00032 
00033 #ifdef RESISTIVITY
00034 
00035 #ifdef HYDRO
00036 #error : resistivity only works for MHD.
00037 #endif /* HYDRO */
00038 
00039 /* current, emf, and energy flux, contained in 3D vector structure */
00040 static Real3Vect ***J=NULL, ***emf=NULL, ***EnerFlux=NULL;
00041 
00042 /* current, Bi in the corrector step & emf in the predict step for Hall MHD */
00043 static Real3Vect ***Jcor=NULL, ***Bcor=NULL, ***emfp=NULL;
00044 
00045 /*==============================================================================
00046  * PRIVATE FUNCTION PROTOTYPES:
00047  *   EField_Ohm  - computes electric field due to Ohmic dissipation
00048  *   EField_Hall - computes electric field due to Hall effect
00049  *   EField_AD   - computes electric field due to ambipolar diffusion
00050  *   hyper_diffusion? -- compute higher-order derivatives of current
00051  *============================================================================*/
00052 
00053 void EField_Ohm(DomainS *pD);
00054 void EField_Hall(DomainS *pD);
00055 void EField_AD(DomainS *pD);
00056 
00057 void EField_Hall_sub(DomainS *pD, Real3Vect ***Bs, Real3Vect ***Js,
00058                                   Real3Vect ***emfs, int noff);
00059 void hyper_diffusion6(DomainS *pD);
00060 
00061 /*=========================== PUBLIC FUNCTIONS ===============================*/
00062 /*----------------------------------------------------------------------------*/
00063 /*! \fn void resistivity(DomainS *pD)
00064  *  \brief Explicit resistivity
00065  */
00066 
00067 void resistivity(DomainS *pD)
00068 {
00069   GridS *pG = (pD->Grid);
00070   int i, is = pG->is, ie = pG->ie;
00071   int j, jl, ju, js = pG->js, je = pG->je;
00072   int k, kl, ku, ks = pG->ks, ke = pG->ke;
00073   int ndim=1;
00074   Real dtodx1 = pG->dt/pG->dx1, dtodx2 = 0.0, dtodx3 = 0.0;
00075 
00076   if (pG->Nx[1] > 1){
00077     jl = js - 3;
00078     ju = je + 4;
00079     dtodx2 = pG->dt/pG->dx2;
00080     ndim++;
00081   } else {
00082     jl = js;
00083     ju = je;
00084   }
00085   if (pG->Nx[2] > 1){
00086     kl = ks - 3;
00087     ku = ke + 4;
00088     dtodx3 = pG->dt/pG->dx3;
00089     ndim++;
00090   } else {
00091     kl = ks;
00092     ku = ke;
00093   }
00094 
00095 /* zero fluxes (electric fields) */
00096 
00097   for (k=kl; k<=ku; k++) {
00098   for (j=jl; j<=ju; j++) {
00099     for (i=is-1; i<=ie+1; i++) {
00100       emf[k][j][i].x = 0.0;
00101       emf[k][j][i].y = 0.0;
00102       emf[k][j][i].z = 0.0;
00103     }
00104   }}
00105 
00106   if (Q_Hall > 0.0) {
00107     for (k=kl; k<=ku; k++) {
00108     for (j=jl; j<=ju; j++) {
00109       for (i=is-3; i<=ie+4; i++) {
00110         emfp[k][j][i].x = 0.0;
00111         emfp[k][j][i].y = 0.0;
00112         emfp[k][j][i].z = 0.0;
00113       }
00114     }}
00115   }
00116 
00117 /*--- Step 1. Compute currents.------------------------------------------------
00118  * Note:  J1 = (dB3/dx2 - dB2/dx3)
00119  *        J2 = (dB1/dx3 - dB3/dx1)
00120  *        J3 = (dB2/dx1 - dB1/dx2) */
00121 
00122 /* 1D PROBLEM */
00123   if (ndim == 1){
00124     for (i=is-3; i<=ie+4; i++) {
00125       J[ks][js][i].x = 0.0;
00126       J[ks][js][i].y = -(pG->U[ks][js][i].B3c - pG->U[ks][js][i-1].B3c)/pG->dx1;
00127       J[ks][js][i].z =  (pG->U[ks][js][i].B2c - pG->U[ks][js][i-1].B2c)/pG->dx1;
00128     }
00129   }
00130 
00131 /* 2D PROBLEM */
00132   if (ndim == 2){
00133     for (j=js-3; j<=je+4; j++) {
00134       for (i=is-3; i<=ie+4; i++) {
00135         J[ks][j][i].x=  (pG->U[ks][j][i].B3c - pG->U[ks][j-1][i  ].B3c)/pG->dx2;
00136         J[ks][j][i].y= -(pG->U[ks][j][i].B3c - pG->U[ks][j  ][i-1].B3c)/pG->dx1;
00137         J[ks][j][i].z=  (pG->B2i[ks][j][i] - pG->B2i[ks][j  ][i-1])/pG->dx1 -
00138                         (pG->B1i[ks][j][i] - pG->B1i[ks][j-1][i  ])/pG->dx2;
00139       }
00140     }
00141   }
00142 
00143 /* 3D PROBLEM */
00144   if (ndim == 3){
00145     for (k=ks-3; k<=ke+4; k++) {
00146      for (j=js-3; j<=je+4; j++) {
00147       for (i=is-3; i<=ie+4; i++) {
00148         J[k][j][i].x = (pG->B3i[k][j][i] - pG->B3i[k  ][j-1][i  ])/pG->dx2 -
00149                        (pG->B2i[k][j][i] - pG->B2i[k-1][j  ][i  ])/pG->dx3;
00150         J[k][j][i].y = (pG->B1i[k][j][i] - pG->B1i[k-1][j  ][i  ])/pG->dx3 -
00151                        (pG->B3i[k][j][i] - pG->B3i[k  ][j  ][i-1])/pG->dx1;
00152         J[k][j][i].z = (pG->B2i[k][j][i] - pG->B2i[k  ][j  ][i-1])/pG->dx1 -
00153                        (pG->B1i[k][j][i] - pG->B1i[k  ][j-1][i  ])/pG->dx2;
00154       }
00155      }
00156     }
00157   }
00158 
00159 /*--- Step 2.  Call functions to compute resistive EMFs ------------------------
00160  * including Ohmic dissipation, the Hall effect, and ambipolar diffusion.
00161  * Current density (J) and emfs are global variables in this file. */
00162   if (eta_Ohm > 0.0)  EField_Ohm(pD);
00163   if (Q_Hall > 0.0) EField_Hall(pD);
00164   if (Q_AD > 0.0)   EField_AD(pD);
00165 #ifndef BAROTROPIC
00166 /*--- Step 3.  Compute energy fluxes -------------------------------------------
00167  * flux of total energy due to resistive diffusion = B X emf
00168  *  EnerFlux.x =  By*emf.z - Bz*emf.y
00169  *  EnerFlux.y =  Bz*emf.x - Bx*emf.z
00170  *  EnerFlux.z =  Bx*emf.y - By*emf.x
00171  */
00172 
00173 /* 1D PROBLEM */
00174   if (ndim == 1){
00175     for (i=is; i<=ie+1; i++) {
00176       EnerFlux[ks][js][i].x =
00177          0.5*(pG->U[ks][js][i].B2c + pG->U[ks][js][i-1].B2c)*emf[ks][js][i].z
00178        - 0.5*(pG->U[ks][js][i].B3c + pG->U[ks][js][i-1].B3c)*emf[ks][js][i].y;
00179     }
00180   } 
00181       
00182 /* 2D PROBLEM */
00183   if (ndim == 2){
00184     for (j=js; j<=je; j++) {
00185     for (i=is; i<=ie+1; i++) {
00186       EnerFlux[ks][j][i].x = 0.25*(pG->U[ks][j][i].B2c + pG->U[ks][j][i-1].B2c)*
00187                             (emf[ks][j][i].z + emf[ks][j+1][i].z)
00188          - 0.5*(pG->U[ks][j][i].B3c + pG->U[ks][j][i-1].B3c)*emf[ks][j][i].y;
00189     }}
00190     
00191     for (j=js; j<=je+1; j++) {
00192     for (i=is; i<=ie; i++) {
00193       EnerFlux[ks][j][i].y =
00194          0.5*(pG->U[ks][j][i].B3c + pG->U[ks][j-1][i].B3c)*emf[ks][j][i].x
00195        - 0.25*(pG->U[ks][j][i].B1c + pG->U[ks][j-1][i].B1c)*
00196                 (emf[ks][j][i].z + emf[ks][j][i+1].z);
00197     }}
00198   }   
00199 
00200 /* 3D PROBLEM */
00201   if (ndim == 3){
00202     for (k=ks; k<=ke; k++) {
00203     for (j=js; j<=je; j++) {
00204       for (i=is; i<=ie+1; i++) {
00205         EnerFlux[k][j][i].x = 0.25*(pG->U[k][j][i].B2c + pG->U[k][j][i-1].B2c)*
00206                              (emf[k][j][i].z + emf[k][j+1][i].z)
00207                             - 0.25*(pG->U[k][j][i].B3c + pG->U[k][j][i-1].B3c)*
00208                              (emf[k][j][i].y + emf[k+1][j][i].y);
00209       }
00210     }}
00211 
00212     for (k=ks; k<=ke; k++) {
00213     for (j=js; j<=je+1; j++) {
00214       for (i=is; i<=ie; i++) {
00215         EnerFlux[k][j][i].y = 0.25*(pG->U[k][j][i].B3c + pG->U[k][j-1][i].B3c)*
00216                              (emf[k][j][i].x + emf[k+1][j][i].x)
00217                             - 0.25*(pG->U[k][j][i].B1c + pG->U[k][j-1][i].B1c)*
00218                              (emf[k][j][i].z + emf[k][j][i+1].z);
00219       }
00220     }}
00221 
00222     for (k=ks; k<=ke+1; k++) {
00223     for (j=js; j<=je; j++) {
00224       for (i=is; i<=ie; i++) {
00225         EnerFlux[k][j][i].z = 0.25*(pG->U[k][j][i].B1c + pG->U[k-1][j][i].B1c)*
00226                              (emf[k][j][i].y + emf[k][j][i+1].y)
00227                             - 0.25*(pG->U[k][j][i].B2c + pG->U[k-1][j][i].B2c)*
00228                              (emf[k][j][i].x + emf[k][j+1][i].x);
00229       }
00230     }}
00231   }
00232 
00233 /*--- Step 4.  Update total energy ---------------------------------------------
00234  * Update energy using x1-fluxes */
00235 
00236   for (k=ks; k<=ke; k++) {
00237   for (j=js; j<=je; j++) {
00238     for (i=is; i<=ie; i++) {
00239       pG->U[k][j][i].E += dtodx1*(EnerFlux[k][j][i+1].x - EnerFlux[k][j][i].x);
00240     }
00241   }}
00242 
00243 /* Update energy using x2-fluxes */
00244 
00245   if (pG->Nx[1] > 1){
00246     for (k=ks; k<=ke; k++) {
00247     for (j=js; j<=je; j++) {
00248       for (i=is; i<=ie; i++) {
00249         pG->U[k][j][i].E += dtodx2*(EnerFlux[k][j+1][i].y -EnerFlux[k][j][i].y);
00250       }
00251     }}
00252   }
00253 
00254 /* Update energy using x3-fluxes */
00255 
00256   if (pG->Nx[2] > 1){
00257     for (k=ks; k<=ke; k++) {
00258     for (j=js; j<=je; j++) {
00259       for (i=is; i<=ie; i++) {
00260         pG->U[k][j][i].E += dtodx3*(EnerFlux[k+1][j][i].z -EnerFlux[k][j][i].z);
00261       }
00262     }}
00263   }
00264 #endif /* BAROTROPIC */
00265 
00266 /*--- Step 5. CT update of magnetic field -------------------------------------
00267  * using total resistive EMFs.  This is identical
00268  * to the CT update in the integrators: dB/dt = -Curl(E) */
00269 
00270 /* 1D PROBLEM: centered differences for B2c and B3c */
00271   if (ndim == 1){
00272     for (i=is; i<=ie; i++) {
00273       pG->U[ks][js][i].B2c += dtodx1*(emf[ks][js][i+1].z - emf[ks][js][i].z);
00274       pG->U[ks][js][i].B3c -= dtodx1*(emf[ks][js][i+1].y - emf[ks][js][i].y);
00275 /* For consistency, set B2i and B3i to cell-centered values. */
00276       pG->B2i[ks][js][i] = pG->U[ks][js][i].B2c;
00277       pG->B3i[ks][js][i] = pG->U[ks][js][i].B3c;
00278     }
00279   }
00280 
00281 /* 2D PROBLEM: CT +  centered differences for B3c */
00282   if (ndim == 2){
00283     for (j=js; j<=je; j++) {
00284       for (i=is; i<=ie; i++) {
00285         pG->B1i[ks][j][i] -= dtodx2*(emf[ks][j+1][i  ].z - emf[ks][j][i].z);
00286         pG->B2i[ks][j][i] += dtodx1*(emf[ks][j  ][i+1].z - emf[ks][j][i].z);
00287 
00288         pG->U[ks][j][i].B3c += dtodx2*(emf[ks][j+1][i  ].x - emf[ks][j][i].x) -
00289                                dtodx1*(emf[ks][j  ][i+1].y - emf[ks][j][i].y);
00290       }
00291       pG->B1i[ks][j][ie+1] -= dtodx2*(emf[ks][j+1][ie+1].z -emf[ks][j][ie+1].z);
00292     }
00293     for (i=is; i<=ie; i++) {
00294       pG->B2i[ks][je+1][i] += dtodx1*(emf[ks][je+1][i+1].z -emf[ks][je+1][i].z);
00295     } 
00296 /* Set cell centered magnetic fields to average of face centered */
00297     for (j=js; j<=je; j++) {
00298       for (i=is; i<=ie; i++) {
00299         pG->U[ks][j][i].B1c = 0.5*(pG->B1i[ks][j][i] + pG->B1i[ks][j][i+1]);
00300         pG->U[ks][j][i].B2c = 0.5*(pG->B2i[ks][j][i] + pG->B2i[ks][j+1][i]);
00301 /* Set the 3-interface magnetic field equal to the cell center field. */
00302         pG->B3i[ks][j][i] = pG->U[ks][j][i].B3c;
00303       }
00304     }
00305   }
00306 
00307 /* 3D PROBLEM: CT */
00308   if (ndim == 3){
00309     for (k=ks; k<=ke; k++) {
00310       for (j=js; j<=je; j++) {
00311         for (i=is; i<=ie; i++) {
00312           pG->B1i[k][j][i] += dtodx3*(emf[k+1][j  ][i  ].y - emf[k][j][i].y) -
00313                               dtodx2*(emf[k  ][j+1][i  ].z - emf[k][j][i].z);
00314           pG->B2i[k][j][i] += dtodx1*(emf[k  ][j  ][i+1].z - emf[k][j][i].z) -
00315                               dtodx3*(emf[k+1][j  ][i  ].x - emf[k][j][i].x);
00316           pG->B3i[k][j][i] += dtodx2*(emf[k  ][j+1][i  ].x - emf[k][j][i].x) -
00317                               dtodx1*(emf[k  ][j  ][i+1].y - emf[k][j][i].y);
00318         }
00319         pG->B1i[k][j][ie+1] +=
00320           dtodx3*(emf[k+1][j  ][ie+1].y - emf[k][j][ie+1].y) -
00321           dtodx2*(emf[k  ][j+1][ie+1].z - emf[k][j][ie+1].z);
00322       }
00323       for (i=is; i<=ie; i++) {
00324         pG->B2i[k][je+1][i] +=
00325           dtodx1*(emf[k  ][je+1][i+1].z - emf[k][je+1][i].z) -
00326           dtodx3*(emf[k+1][je+1][i  ].x - emf[k][je+1][i].x);
00327       }
00328     }
00329     for (j=js; j<=je; j++) {
00330     for (i=is; i<=ie; i++) {
00331       pG->B3i[ke+1][j][i] +=
00332         dtodx2*(emf[ke+1][j+1][i  ].x - emf[ke+1][j][i].x) -
00333         dtodx1*(emf[ke+1][j  ][i+1].y - emf[ke+1][j][i].y);
00334     }}
00335 /* Set cell centered magnetic fields to average of face centered */
00336     for (k=ks; k<=ke; k++) {
00337     for (j=js; j<=je; j++) {
00338     for (i=is; i<=ie; i++) {
00339       pG->U[k][j][i].B1c = 0.5*(pG->B1i[k][j][i] + pG->B1i[k][j][i+1]);
00340       pG->U[k][j][i].B2c = 0.5*(pG->B2i[k][j][i] + pG->B2i[k][j+1][i]);
00341       pG->U[k][j][i].B3c = 0.5*(pG->B3i[k][j][i] + pG->B3i[k+1][j][i]);
00342     }}}
00343   }
00344 
00345   return;
00346 }
00347 
00348 /*----------------------------------------------------------------------------*/
00349 /*! \fn void EField_Ohm(DomainS *pD)
00350  *  \brief Resistive EMF from Ohmic dissipation.   E = \eta_Ohm J
00351  */
00352 
00353 void EField_Ohm(DomainS *pD)
00354 {
00355   GridS *pG = (pD->Grid);
00356   int i, is = pG->is, ie = pG->ie;
00357   int j, js = pG->js, je = pG->je;
00358   int k, ks = pG->ks, ke = pG->ke;
00359   int ndim=1;
00360   Real eta_O;
00361 
00362   if (pG->Nx[1] > 1)    ndim++;
00363   if (pG->Nx[2] > 1)    ndim++;
00364 
00365 /* For Ohmic resistivity, E = \eta_Ohm J  */
00366 
00367 /* 1D PROBLEM: */
00368   if (ndim == 1){
00369     for (i=is; i<=ie+1; i++) {
00370 
00371       eta_O = 0.5*(pG->eta_Ohm[ks][js][i] + pG->eta_Ohm[ks][js][i-1]);
00372 
00373       emf[ks][js][i].y += eta_O * J[ks][js][i].y;
00374       emf[ks][js][i].z += eta_O * J[ks][js][i].z;
00375     }
00376   }
00377 
00378 /* 2D PROBLEM: */
00379   if (ndim == 2){
00380     for (j=js; j<=je+1; j++) {
00381     for (i=is; i<=ie+1; i++) {
00382 
00383       eta_O = 0.5*(pG->eta_Ohm[ks][j][i] + pG->eta_Ohm[ks][j-1][i]);
00384 
00385       emf[ks][j][i].x += eta_O * J[ks][j][i].x;
00386 
00387       eta_O = 0.5*(pG->eta_Ohm[ks][j][i] + pG->eta_Ohm[ks][j][i-1]);
00388 
00389       emf[ks][j][i].y += eta_O * J[ks][j][i].y; 
00390 
00391       eta_O = 0.25*(pG->eta_Ohm[ks][j][i  ] + pG->eta_Ohm[ks][j-1][i  ] +
00392                     pG->eta_Ohm[ks][j][i-1] + pG->eta_Ohm[ks][j-1][i-1]);
00393 
00394       emf[ks][j][i].z += eta_O * J[ks][j][i].z;
00395     }}
00396   }
00397 
00398 
00399 /* 3D PROBLEM: */
00400   
00401   if (ndim == 3){
00402     for (k=ks; k<=ke+1; k++) {
00403     for (j=js; j<=je+1; j++) {
00404       for (i=is; i<=ie+1; i++) {
00405 
00406         eta_O = 0.25*(pG->eta_Ohm[k][j  ][i] + pG->eta_Ohm[k-1][j  ][i] +
00407                       pG->eta_Ohm[k][j-1][i] + pG->eta_Ohm[k-1][j-1][i]);
00408 
00409         emf[k][j][i].x += eta_O * J[k][j][i].x;
00410 
00411         eta_O = 0.25*(pG->eta_Ohm[k][j][i  ] + pG->eta_Ohm[k-1][j][i  ] +
00412                       pG->eta_Ohm[k][j][i-1] + pG->eta_Ohm[k-1][j][i-1]);
00413 
00414         emf[k][j][i].y += eta_O * J[k][j][i].y;
00415 
00416         eta_O = 0.25*(pG->eta_Ohm[k][j][i  ] + pG->eta_Ohm[k][j-1][i  ] +
00417                       pG->eta_Ohm[k][j][i-1] + pG->eta_Ohm[k][j-1][i-1]);
00418 
00419         emf[k][j][i].z += eta_O * J[k][j][i].z;
00420       }
00421     }}
00422   }
00423 
00424   return;
00425 }
00426 
00427 /*----------------------------------------------------------------------------*/
00428 /*! \fn void EField_Hall(DomainS *pD)
00429  *  \brief Resistive EMF from Hall effect.  E = Q_H (J X B)
00430  */
00431 void EField_Hall(DomainS *pD)
00432 {
00433   GridS *pG = (pD->Grid);
00434   int i, il,iu, is = pG->is, ie = pG->ie;
00435   int j, jl,ju, js = pG->js, je = pG->je;
00436   int k, kl,ku, ks = pG->ks, ke = pG->ke;
00437   int ndim=1;
00438   Real eta_H, Bmag;
00439   Real dtodx1 = pG->dt/pG->dx1, dtodx2 = 0.0, dtodx3 = 0.0;
00440 
00441   il = is - 3;  iu = ie + 4;
00442 
00443   if (pG->Nx[1] > 1){
00444     jl = js - 3;    ju = je + 4;
00445     dtodx2 = pG->dt/pG->dx2;
00446     ndim++;
00447   } else {
00448     jl = js;        ju = je;
00449   }
00450   if (pG->Nx[2] > 1){
00451     kl = ks - 3;    ku = ke + 4;
00452     dtodx3 = pG->dt/pG->dx3;
00453     ndim++;
00454   } else {
00455     kl = ks;        ku = ke;
00456   }
00457 
00458 /* Preliminary: add hyper-diffusion */
00459   hyper_diffusion6(pD);
00460 
00461 /* Preliminary: divide eta_Hall by B for convenience */
00462   for (k=kl; k<=ku; k++) {
00463   for (j=jl; j<=ju; j++) {
00464   for (i=il; i<=iu; i++) {
00465 
00466       Bmag = sqrt(SQR(pG->U[k][j][i].B1c)
00467                 + SQR(pG->U[k][j][i].B2c) + SQR(pG->U[k][j][i].B3c));
00468 
00469       pG->eta_Hall[k][j][i] /= Bmag+TINY_NUMBER;
00470 
00471       Bcor[k][j][i].x = pG->B1i[k][j][i];
00472       Bcor[k][j][i].y = pG->B2i[k][j][i];
00473       Bcor[k][j][i].z = pG->B3i[k][j][i];
00474   }}}
00475 
00476 /* Step 1: Predictor step for emf calculation at t=t(n)  */
00477   EField_Hall_sub(pD, Bcor, J, emfp, 3);
00478 
00479 /* Step 2: Calculate Bcor and Jcor at t=t(n+1/2) */
00480 
00481   /* 1D PROBLEM: centered differences for B2c and B3c */
00482   if (ndim == 1){
00483     for (i=is-2; i<=ie+2; i++) {
00484       Bcor[ks][js][i].y += 0.5*dtodx1*(emfp[ks][js][i+1].z - emfp[ks][js][i].z);
00485       Bcor[ks][js][i].z -= 0.5*dtodx1*(emfp[ks][js][i+1].y - emfp[ks][js][i].y);
00486     }
00487   }
00488 
00489   /* 2D PROBLEM: CT +  centered differences for B3c */
00490   if (ndim == 2){
00491     for (j=js-2; j<=je+2; j++) {
00492       for (i=is-2; i<=ie+2; i++) {
00493         Bcor[ks][j][i].x -= 0.5*dtodx2*(emfp[ks][j+1][i  ].z - emfp[ks][j][i].z);
00494         Bcor[ks][j][i].y += 0.5*dtodx1*(emfp[ks][j  ][i+1].z - emfp[ks][j][i].z);
00495 
00496         Bcor[ks][j][i].z += 0.5*dtodx2*(emfp[ks][j+1][i  ].x - emfp[ks][j][i].x) -
00497                             0.5*dtodx1*(emfp[ks][j  ][i+1].y - emfp[ks][j][i].y);
00498       }
00499     }
00500   }
00501   
00502   if (ndim == 3){
00503     for (k=ks-2; k<=ke+2; k++) {
00504       for (j=js-2; j<=je+2; j++) {
00505         for (i=is-2; i<=ie+2; i++) {
00506           Bcor[k][j][i].x += 0.5*dtodx3*(emfp[k+1][j  ][i  ].y - emfp[k][j][i].y) -
00507                              0.5*dtodx2*(emfp[k  ][j+1][i  ].z - emfp[k][j][i].z);
00508           Bcor[k][j][i].y += 0.5*dtodx1*(emfp[k  ][j  ][i+1].z - emfp[k][j][i].z) -
00509                              0.5*dtodx3*(emfp[k+1][j  ][i  ].x - emfp[k][j][i].x);
00510           Bcor[k][j][i].z += 0.5*dtodx2*(emfp[k  ][j+1][i  ].x - emfp[k][j][i].x) -
00511                              0.5*dtodx1*(emfp[k  ][j  ][i+1].y - emfp[k][j][i].y);
00512         }
00513       }
00514     }
00515   }
00516 
00517 /* Step 3: compute current for the corrector step */
00518 
00519 /* 1D PROBLEM */
00520   if (ndim == 1){
00521     for (i=is-1; i<=ie+2; i++) {
00522       Jcor[ks][js][i].x = 0.0;
00523       Jcor[ks][js][i].y = -(Bcor[ks][js][i].z - Bcor[ks][js][i-1].z)/pG->dx1;
00524       Jcor[ks][js][i].z =  (Bcor[ks][js][i].y - Bcor[ks][js][i-1].y)/pG->dx1;
00525     }
00526   }
00527 
00528 /* 2D PROBLEM */
00529   if (ndim == 2){
00530     for (j=js-1; j<=je+2; j++) {
00531       for (i=is-1; i<=ie+2; i++) {
00532         Jcor[ks][j][i].x=  (Bcor[ks][j][i].z - Bcor[ks][j-1][i  ].z)/pG->dx2;
00533         Jcor[ks][j][i].y= -(Bcor[ks][j][i].z - Bcor[ks][j  ][i-1].z)/pG->dx1;
00534         Jcor[ks][j][i].z=  (Bcor[ks][j][i].y - Bcor[ks][j  ][i-1].y)/pG->dx1 -
00535                            (Bcor[ks][j][i].x - Bcor[ks][j-1][i  ].x)/pG->dx2;
00536       }
00537     }
00538   }
00539 
00540 /* 3D PROBLEM */
00541   if (ndim == 3){
00542     for (k=ks-1; k<=ke+2; k++) {
00543     for (j=js-1; j<=je+2; j++) {
00544       for (i=is-1; i<=ie+2; i++) {
00545         Jcor[k][j][i].x = (Bcor[k][j][i].z - Bcor[k  ][j-1][i  ].z)/pG->dx2 -
00546                           (Bcor[k][j][i].y - Bcor[k-1][j  ][i  ].y)/pG->dx3;
00547         Jcor[k][j][i].y = (Bcor[k][j][i].x - Bcor[k-1][j  ][i  ].x)/pG->dx3 -
00548                           (Bcor[k][j][i].z - Bcor[k  ][j  ][i-1].z)/pG->dx1;
00549         Jcor[k][j][i].z = (Bcor[k][j][i].y - Bcor[k  ][j  ][i-1].y)/pG->dx1 -
00550                           (Bcor[k][j][i].x - Bcor[k  ][j-1][i  ].x)/pG->dx2;
00551       }
00552     }}
00553   }
00554 
00555 /* Step 4: Corrector step fir emf  calculation at t=t(n+1/2) */
00556   EField_Hall_sub(pD, Bcor, Jcor, emf, 1);
00557 
00558   return;
00559 
00560 }
00561 
00562 /*! \fn void EField_Hall_sub(DomainS *pD, Real3Vect ***Bs, Real3Vect ***Js,
00563  *                                Real3Vect ***emfs, int noff)
00564  *  \brief Update the Hall emfs (subroutine for the Hall calculation) */
00565 void EField_Hall_sub(DomainS *pD, Real3Vect ***Bs, Real3Vect ***Js,
00566                                   Real3Vect ***emfs, int noff)
00567 {
00568   GridS *pG = (pD->Grid);
00569   int i, is = pG->is, ie = pG->ie;
00570   int j, js = pG->js, je = pG->je;
00571   int k, ks = pG->ks, ke = pG->ke;
00572   int ndim=1;
00573   Real eta_H;
00574 
00575   if (pG->Nx[1] > 1)  ndim++;
00576   if (pG->Nx[2] > 1)  ndim++;
00577 
00578 /* 1D PROBLEM:
00579  *   emf.x =  0.0
00580  *   emf.y =  Jz*Bx
00581  *   emf.z = -Jy*Bx   */
00582 
00583   if (ndim == 1){
00584     for (i=is-noff+1; i<=ie+noff; i++) {
00585       eta_H = 0.5*(pG->eta_Hall[ks][js][i] + pG->eta_Hall[ks][js][i-1]);
00586 
00587       emfs[ks][js][i].y += eta_H*Js[ks][js][i].z * Bs[ks][js][i].x;
00588       emfs[ks][js][i].z -= eta_H*Js[ks][js][i].y * Bs[ks][js][i].x;
00589     }
00590   }
00591 
00592 /* 2D PROBLEM:
00593  *  emf.x =  Jy*Bz - Jz*By
00594  *  emf.y =  Jz*Bx - Jx*Bz
00595  *  emf.z =  Jx*By - Jy*Bx  */
00596 
00597   if (ndim == 2){
00598     for (j=js-noff+1; j<=je+noff; j++) {
00599     for (i=is-noff+1; i<=ie+noff; i++) {
00600 
00601       /* x1 */
00602       eta_H = 0.5*(pG->eta_Hall[ks][j][i] + pG->eta_Hall[ks][j-1][i]);
00603 
00604       emfs[ks][j][i].x += eta_H*(
00605         0.125*(Js[ks][j  ][i].y + Js[ks][j  ][i+1].y
00606              + Js[ks][j-1][i].y + Js[ks][j-1][i+1].y)
00607              *(Bs[ks][j  ][i].z + Bs[ks][j-1][i  ].z) -
00608          0.5*((Js[ks][j  ][i].z + Js[ks][j  ][i+1].z)*Bs[ks][j  ][i].y));
00609 
00610       /* x2 */
00611       eta_H = 0.5*(pG->eta_Hall[ks][j][i] + pG->eta_Hall[ks][j][i-1]);
00612 
00613       emfs[ks][j][i].y += eta_H*(
00614          0.5*((Js[ks][j][i  ].z + Js[ks][j+1][i  ].z)*Bs[ks][j][i  ].x) -
00615         0.125*(Js[ks][j][i  ].x + Js[ks][j+1][i  ].x
00616              + Js[ks][j][i-1].x + Js[ks][j+1][i-1].x)
00617              *(Bs[ks][j][i  ].z + Bs[ks][j  ][i-1].z));
00618 
00619       /* x3 */
00620       eta_H = 0.25*(pG->eta_Hall[ks][j][i  ] + pG->eta_Hall[ks][j-1][i  ] +
00621                     pG->eta_Hall[ks][j][i-1] + pG->eta_Hall[ks][j-1][i-1]);
00622 
00623       emfs[ks][j][i].z += eta_H*(
00624         0.25*(Js[ks][j][i].x + Js[ks][j][i-1].x)
00625             *(Bs[ks][j][i].y + Bs[ks][j][i-1].y) -
00626         0.25*(Js[ks][j][i].y + Js[ks][j-1][i].y)
00627             *(Bs[ks][j][i].x + Bs[ks][j-1][i].x));
00628     }}
00629   }
00630 
00631 /* 3D PROBLEM:
00632  *  emf.x =  Jy*Bz - Jz*By
00633  *  emf.y =  Jz*Bx - Jx*Bz
00634  *  emf.z =  Jx*By - Jy*Bx  */
00635 
00636   if (ndim == 3){
00637     for (k=ks-noff+1; k<=ke+noff; k++) {
00638     for (j=js-noff+1; j<=je+noff; j++) {
00639       for (i=is-noff+1; i<=ie+noff; i++) {
00640 
00641         /* x1 */
00642         eta_H = 0.25*(pG->eta_Hall[k][j  ][i] + pG->eta_Hall[k-1][j  ][i] +
00643                       pG->eta_Hall[k][j-1][i] + pG->eta_Hall[k-1][j-1][i]);
00644 
00645         emfs[k][j][i].x += 0.125*eta_H*(
00646                 (Js[k  ][j  ][i].y + Js[k  ][j  ][i+1].y
00647                + Js[k  ][j-1][i].y + Js[k  ][j-1][i+1].y)
00648                *(Bs[k  ][j  ][i].z + Bs[k  ][j-1][i  ].z)-
00649                 (Js[k  ][j  ][i].z + Js[k  ][j  ][i+1].z
00650                + Js[k-1][j  ][i].z + Js[k-1][j  ][i+1].z)
00651                *(Bs[k  ][j  ][i].y + Bs[k-1][j  ][i  ].y));
00652 
00653         /* x2 */
00654         eta_H = 0.25*(pG->eta_Hall[k][j][i  ] + pG->eta_Hall[k-1][j][i  ] +
00655                       pG->eta_Hall[k][j][i-1] + pG->eta_Hall[k-1][j][i-1]);
00656 
00657         emfs[k][j][i].y += 0.125*eta_H*(
00658                 (Js[k  ][j][i  ].z + Js[k  ][j+1][i  ].z
00659                + Js[k-1][j][i  ].z + Js[k-1][j+1][i  ].z)
00660                *(Bs[k  ][j][i  ].x + Bs[k-1][j  ][i  ].x)-
00661                 (Js[k  ][j][i  ].x + Js[k  ][j+1][i  ].x
00662                + Js[k  ][j][i-1].x + Js[k  ][j+1][i-1].x)
00663                *(Bs[k  ][j][i  ].z + Bs[k  ][j  ][i-1].z));
00664 
00665         /* x3 */
00666         eta_H = 0.25*(pG->eta_Hall[k][j][i  ] + pG->eta_Hall[k][j-1][i  ] +
00667                       pG->eta_Hall[k][j][i-1] + pG->eta_Hall[k][j-1][i-1]);
00668 
00669         emfs[k][j][i].z += 0.125*eta_H*(
00670                 (Js[k][j  ][i  ].x + Js[k+1][j  ][i  ].x
00671                + Js[k][j  ][i-1].x + Js[k+1][j  ][i-1].x)
00672                *(Bs[k][j  ][i  ].y + Bs[k  ][j  ][i-1].y)-
00673                 (Js[k][j  ][i  ].y + Js[k+1][j  ][i  ].y
00674                + Js[k][j-1][i  ].y + Js[k+1][j-1][i  ].y)
00675                *(Bs[k][j  ][i  ].x + Bs[k  ][j-1][i  ].x));
00676       }
00677     }}
00678   }
00679 
00680   return;
00681 }
00682 
00683 /*----------------------------------------------------------------------------*/
00684 /*! \fn void EField_AD(DomainS *pD)
00685  *  \brief Resistive EMF from ambipolar diffusion.  E = Q_AD (J X B) X B
00686  */
00687 void EField_AD(DomainS *pD)
00688 {
00689   GridS *pG = (pD->Grid);
00690   int i, is = pG->is, ie = pG->ie;
00691   int j, js = pG->js, je = pG->je;
00692   int k, ks = pG->ks, ke = pG->ke;
00693   int ndim=1;
00694   Real eta_A;
00695   Real intBx,intBy,intBz,intJx,intJy,intJz,Bsq,JdotB;
00696 
00697   if (pG->Nx[1] > 1) ndim++;
00698   if (pG->Nx[2] > 1) ndim++;
00699 
00700 /* 1D PROBLEM:
00701  *   emf.x = 0
00702  *   emf.y = (J_perp)_y
00703  *   emf.z = (J_perp)_z  */
00704 
00705   if (ndim == 1){
00706     for (i=is; i<=ie+1; i++) {
00707       eta_A = 0.5*(pG->eta_AD[ks][js][i] + pG->eta_AD[ks][js][i-1]);
00708 
00709       intBx = pG->B1i[ks][js][i];
00710       intBy = 0.5*(pG->U[ks][js][i].B2c + pG->U[ks][js][i-1].B2c);
00711       intBz = 0.5*(pG->U[ks][js][i].B3c + pG->U[ks][js][i-1].B3c);
00712 
00713       Bsq = SQR(intBx) + SQR(intBy) + SQR(intBz);
00714       JdotB = J[ks][js][i].y*intBy + J[ks][js][i].z*intBz;
00715 
00716       emf[ks][js][i].y += eta_A*(J[ks][js][i].y - JdotB*intBy/Bsq);
00717       emf[ks][js][i].z += eta_A*(J[ks][js][i].z - JdotB*intBz/Bsq);
00718     }
00719   }
00720 
00721 /* 2D PROBLEM:
00722  *   emf.x = (J_perp)_x
00723  *   emf.y = (J_perp)_y
00724  *   emf.z = (J_perp)_z  */
00725 
00726   if (ndim == 2){
00727     for (j=js; j<=je+1; j++) {
00728     for (i=is; i<=ie+1; i++) {
00729 
00730       /* emf.x */
00731       eta_A = 0.5*(pG->eta_AD[ks][j][i] + pG->eta_AD[ks][j-1][i]);
00732 
00733       intJx = J[ks][j][i].x;
00734       intJy = 0.25*(J[ks][j  ][i].y + J[ks][j  ][i+1].y
00735                   + J[ks][j-1][i].y + J[ks][j-1][i+1].y);
00736       intJz = 0.5 *(J[ks][j][i].z   + J[ks][j][i+1].z);
00737 
00738       intBx = 0.5*(pG->U[ks][j][i].B1c + pG->U[ks][j-1][i].B1c);
00739       intBy = pG->B2i[ks][j][i];
00740       intBz = 0.5*(pG->U[ks][j][i].B3c + pG->U[ks][j-1][i].B3c);
00741 
00742       Bsq = SQR(intBx) + SQR(intBy) + SQR(intBz);
00743       JdotB = intJx*intBx + intJy*intBy + intJz*intBz;
00744 
00745       emf[ks][j][i].x += eta_A*(J[ks][j][i].x - JdotB*intBx/Bsq);
00746 
00747       /* emf.y */
00748       eta_A = 0.5*(pG->eta_AD[ks][j][i] + pG->eta_AD[ks][j][i-1]);
00749 
00750       intJx = 0.25*(J[ks][j][i  ].x + J[ks][j+1][i  ].x
00751                   + J[ks][j][i-1].x + J[ks][j+1][i-1].x);
00752       intJy = J[ks][j][i].y;
00753       intJz = 0.5 *(J[ks][j][i].z   + J[ks][j+1][i].z);
00754 
00755       intBx = pG->B1i[ks][j][i];
00756       intBy = 0.5*(pG->U[ks][j][i].B2c + pG->U[ks][j][i-1].B2c);
00757       intBz = 0.5*(pG->U[ks][j][i].B3c + pG->U[ks][j][i-1].B3c);
00758 
00759       Bsq = SQR(intBx) + SQR(intBy) + SQR(intBz);
00760       JdotB = intJx*intBx + intJy*intBy + intJz*intBz;
00761 
00762       emf[ks][j][i].y += eta_A*(J[ks][j][i].y - JdotB*intBy/Bsq);
00763 
00764      /* emf.z */
00765       eta_A = 0.25*(pG->eta_AD[ks][j  ][i] + pG->eta_AD[ks][j  ][i-1]
00766                   + pG->eta_AD[ks][j-1][i] + pG->eta_AD[ks][j-1][i-1]);
00767 
00768       intJx = 0.5*(J[ks][j][i].x + J[ks][j][i-1].x);
00769       intJy = 0.5*(J[ks][j][i].y + J[ks][j-1][i].y);
00770       intJz = J[ks][j][i].z;
00771 
00772       intBx = 0.5*(pG->B1i[ks][j][i] + pG->B1i[ks][j-1][i]);
00773       intBy = 0.5*(pG->B2i[ks][j][i] + pG->B2i[ks][j][i-1]);
00774       intBz = 0.25*(pG->U[ks][j  ][i].B3c + pG->U[ks][j  ][i-1].B3c
00775                   + pG->U[ks][j-1][i].B3c + pG->U[ks][j-1][i-1].B3c);
00776 
00777       Bsq = SQR(intBx) + SQR(intBy) + SQR(intBz);
00778       JdotB = intJx*intBx + intJy*intBy + intJz*intBz;
00779 
00780       emf[ks][j][i].z += eta_A*(J[ks][j][i].z - JdotB*intBz/Bsq);
00781 
00782     }}
00783   }
00784 
00785 /* 3D PROBLEM:
00786  *   emf.x = (J_perp)_x
00787  *   emf.y = (J_perp)_y
00788  *   emf.z = (J_perp)_z  */
00789 
00790   if (ndim == 3){
00791     for (k=ks; k<=ke+1; k++) {
00792     for (j=js; j<=je+1; j++) {
00793       for (i=is; i<=ie+1; i++) {
00794 
00795         /* emf.x */
00796         eta_A = 0.25*(pG->eta_AD[k][j  ][i] + pG->eta_AD[k-1][j  ][i] +
00797                       pG->eta_AD[k][j-1][i] + pG->eta_AD[k-1][j-1][i]);
00798 
00799         intJx = J[k][j][i].x;
00800         intJy = 0.25*(J[k][j  ][i].y + J[k][j  ][i+1].y
00801                     + J[k][j-1][i].y + J[k][j-1][i+1].y);
00802         intJz = 0.25*(J[k  ][j][i].z + J[k  ][j][i+1].z
00803                     + J[k-1][j][i].z + J[k-1][j][i+1].z);
00804 
00805         intBx = 0.25*(pG->U[k][j  ][i].B1c + pG->U[k-1][j  ][i].B1c +
00806                       pG->U[k][j-1][i].B1c + pG->U[k-1][j-1][i].B1c);
00807         intBy = 0.5*(pG->B2i[k][j][i] + pG->B2i[k-1][j][i]);
00808         intBz = 0.5*(pG->B3i[k][j][i] + pG->B3i[k][j-1][i]);
00809 
00810         Bsq = SQR(intBx) + SQR(intBy) + SQR(intBz);
00811         JdotB = intJx*intBx + intJy*intBy + intJz*intBz;
00812 
00813         emf[k][j][i].x += eta_A*(J[k][j][i].x - JdotB*intBx/Bsq);
00814 
00815         /* emf.y */
00816         eta_A = 0.25*(pG->eta_AD[k][j][i  ] + pG->eta_AD[k-1][j][i  ] +
00817                       pG->eta_AD[k][j][i-1] + pG->eta_AD[k-1][j][i-1]);
00818 
00819         intJx = 0.25*(J[k][j][i  ].x + J[k][j+1][i  ].x
00820                     + J[k][j][i-1].x + J[k][j+1][i-1].x);;
00821         intJy = J[k][j][i].y;
00822         intJz = 0.25*(J[k  ][j][i].z + J[k  ][j+1][i].z
00823                     + J[k-1][j][i].z + J[k-1][j+1][i].z);
00824 
00825         intBx = 0.5*(pG->B1i[k][j][i] + pG->B1i[k-1][j][i]);
00826         intBy = 0.25*(pG->U[k][j][i  ].B2c + pG->U[k-1][j][i  ].B2c +
00827                       pG->U[k][j][i-1].B2c + pG->U[k-1][j][i-1].B2c);
00828         intBz = 0.5*(pG->B3i[k][j][i] + pG->B3i[k][j][i-1]);
00829 
00830         Bsq = SQR(intBx) + SQR(intBy) + SQR(intBz);
00831         JdotB = intJx*intBx + intJy*intBy + intJz*intBz;
00832 
00833         emf[k][j][i].y += eta_A*(J[k][j][i].y - JdotB*intBy/Bsq);
00834 
00835         /* emf.z */
00836         eta_A = 0.25*(pG->eta_AD[k][j][i  ] + pG->eta_AD[k][j-1][i  ] +
00837                       pG->eta_AD[k][j][i-1] + pG->eta_AD[k][j-1][i-1]);
00838 
00839         intJx = 0.25*(J[k][j][i  ].x + J[k+1][j][i  ].x
00840                     + J[k][j][i-1].x + J[k+1][j][i-1].x);;
00841         intJy = 0.25*(J[k][j  ][i].y + J[k+1][j  ][i].y
00842                     + J[k][j-1][i].y + J[k+1][j-1][i].y);
00843         intJz = J[k][j][i].z;
00844 
00845         intBx = 0.5*(pG->B1i[k][j][i] + pG->B1i[k][j-1][i]);
00846         intBy = 0.5*(pG->B2i[k][j][i] + pG->B2i[k][j][i-1]);
00847         intBz = 0.25*(pG->U[k][j  ][i].B3c + pG->U[k][j  ][i-1].B3c +
00848                       pG->U[k][j-1][i].B3c + pG->U[k][j-1][i-1].B3c);
00849 
00850         Bsq = SQR(intBx) + SQR(intBy) + SQR(intBz);
00851         JdotB = intJx*intBx + intJy*intBy + intJz*intBz;
00852 
00853         emf[k][j][i].z += eta_A*(J[k][j][i].z - JdotB*intBz/Bsq);
00854       }
00855     }}
00856   }
00857 
00858   return;
00859 }
00860 
00861 /*----------------------------------------------------------------------------*/
00862 /*! \fn void hyper_diffusion6(DomainS *pD)
00863  *  \brief Calculate the higher-order derivatives of J -- 6th order diffusion
00864  */
00865 
00866 /* 6th order diffusion */
00867 void hyper_diffusion6(DomainS *pD)
00868 {
00869   GridS *pG = (pD->Grid);
00870   int i, is = pG->is, ie = pG->ie;
00871   int j, js = pG->js, je = pG->je;
00872   int k, ks = pG->ks, ke = pG->ke;
00873   int ndim=1;
00874   Real eta_H,eta_6,dx41,dy41=0.0,dz41=0.0;
00875   Real fac,fac2,fac3;
00876 
00877   dx41 = 1.0/SQR(SQR(pG->dx1));
00878   if (pG->Nx[1]>1) {
00879     ndim++;
00880     dy41 = 1.0/SQR(SQR(pG->dx2));
00881     fac2 = SQR(pG->dx1/pG->dx2);
00882   }
00883   if (pG->Nx[2]>1) {
00884     ndim++;
00885     dz41 = 1.0/SQR(SQR(pG->dx3));
00886     fac3 = SQR(pG->dx1/pG->dx3);
00887   }
00888   fac = 2.0*SQR(pG->dt/pG->dx1)*pG->dt;
00889   
00890   /* 1D */ 
00891   if (ndim == 1) {
00892     for (i=is; i<=ie+1; i++) {
00893 
00894       eta_H = 0.5*(pG->eta_Hall[ks][js][i] + pG->eta_Hall[ks][js][i-1]);
00895       eta_6 = SQR(SQR(eta_H)) * fac;
00896 
00897       emf[ks][js][i].y += eta_6 *  (J[ks][js][i-2].y - 4.0*J[ks][js][i-1].y
00898          + 6.0*J[ks][js][i].y - 4.0*J[ks][js][i+1].y + J[ks][js][i+2].y) * dx41;
00899       emf[ks][js][i].z += eta_6 *  (J[ks][js][i-2].z - 4.0*J[ks][js][i-1].z
00900          + 6.0*J[ks][js][i].z - 4.0*J[ks][js][i+1].z + J[ks][js][i+2].z) * dx41;
00901     }
00902   }
00903   
00904   /* 2D */ 
00905   if (ndim == 2) {
00906     for (j=js; j<=je+1; j++) {
00907     for (i=is; i<=ie+1; i++) {
00908 
00909       /* x1 */
00910       eta_H = 0.5*(pG->eta_Hall[ks][j][i] + pG->eta_Hall[ks][j-1][i]);
00911       eta_6 = SQR(SQR(eta_H)) * fac;
00912 
00913       emf[ks][j][i].x += eta_6 * ((J[ks][j][i-2].x - 4.0*J[ks][j][i-1].x
00914          + 6.0*J[ks][j][i].x - 4.0*J[ks][j][i+1].x + J[ks][j][i+2].x)*dx41
00915                          + fac2 * (J[ks][j-2][i].x - 4.0*J[ks][j-1][i].x
00916          + 6.0*J[ks][j][i].x - 4.0*J[ks][j+1][i].x + J[ks][j+2][i].x) * dy41);
00917 
00918       /* x2 */
00919       eta_H = 0.5*(pG->eta_Hall[ks][j][i] + pG->eta_Hall[ks][j][i-1]);
00920       eta_6 = SQR(SQR(eta_H)) * fac;
00921 
00922       emf[ks][j][i].y += eta_6 * ((J[ks][j][i-2].y - 4.0*J[ks][j][i-1].y
00923          + 6.0*J[ks][j][i].y - 4.0*J[ks][j][i+1].y + J[ks][j][i+2].y) * dx41
00924                          + fac2 * (J[ks][j-2][i].y - 4.0*J[ks][j-1][i].y
00925          + 6.0*J[ks][j][i].y - 4.0*J[ks][j+1][i].y + J[ks][j+2][i].y) * dy41);
00926 
00927       /* x3 */
00928       eta_H = 0.25*(pG->eta_Hall[ks][j][i  ] + pG->eta_Hall[ks][j-1][i  ] +
00929                     pG->eta_Hall[ks][j][i-1] + pG->eta_Hall[ks][j-1][i-1]);
00930       eta_6 = SQR(SQR(eta_H)) * fac;
00931 
00932       emf[ks][j][i].z += eta_6 * ((J[ks][j][i-2].z - 4.0*J[ks][j][i-1].z
00933          + 6.0*J[ks][j][i].z - 4.0*J[ks][j][i+1].z + J[ks][j][i+2].z) * dx41
00934                          + fac2 * (J[ks][j-2][i].z - 4.0*J[ks][j-1][i].z
00935          + 6.0*J[ks][j][i].z - 4.0*J[ks][j+1][i].z + J[ks][j+2][i].z) * dy41);
00936     }}
00937   }
00938 
00939   /* 3D */
00940   if (ndim == 3) {
00941     for (k=ks; k<=ke+1; k++) {
00942     for (j=js; j<=je+1; j++) {
00943     for (i=is; i<=ie+1; i++) {
00944 
00945       /* x1 */
00946       eta_H = 0.25*(pG->eta_Hall[k][j  ][i] + pG->eta_Hall[k-1][j  ][i] +
00947                     pG->eta_Hall[k][j-1][i] + pG->eta_Hall[k-1][j-1][i]);
00948       eta_6 = SQR(SQR(eta_H)) * fac;
00949 
00950       emf[k][j][i].x += eta_6 * ((J[k][j][i-2].x - 4.0*J[k][j][i-1].x
00951          + 6.0*J[k][j][i].x - 4.0*J[k][j][i+1].x + J[k][j][i+2].x) * dx41
00952                         + fac2 * (J[k][j-2][i].x - 4.0*J[k][j-1][i].x
00953          + 6.0*J[k][j][i].x - 4.0*J[k][j+1][i].x + J[k][j+2][i].x) * dy41
00954                         + fac3 * (J[k-2][j][i].x - 4.0*J[k-1][j][i].x
00955          + 6.0*J[k][j][i].x - 4.0*J[k+1][j][i].x + J[k+2][j][i].x) * dz41);
00956 
00957       /* x2 */
00958       eta_H = 0.25*(pG->eta_Hall[k][j][i  ] + pG->eta_Hall[k-1][j][i  ] +
00959                     pG->eta_Hall[k][j][i-1] + pG->eta_Hall[k-1][j][i-1]);
00960       eta_6 = SQR(SQR(eta_H)) * fac;
00961 
00962       emf[k][j][i].y += eta_6 * ((J[k][j][i-2].y - 4.0*J[k][j][i-1].y
00963          + 6.0*J[k][j][i].y - 4.0*J[k][j][i+1].y + J[k][j][i+2].y) * dx41
00964                         + fac2 * (J[k][j-2][i].y - 4.0*J[k][j-1][i].y
00965          + 6.0*J[k][j][i].y - 4.0*J[k][j+1][i].y + J[k][j+2][i].y) * dy41
00966                         + fac3 * (J[k-2][j][i].y - 4.0*J[k-1][j][i].y
00967          + 6.0*J[k][j][i].y - 4.0*J[k+1][j][i].y + J[k+2][j][i].y) * dz41);
00968 
00969       /* x3 */
00970       eta_H = 0.25*(pG->eta_Hall[k][j][i  ] + pG->eta_Hall[k][j-1][i  ] +
00971                     pG->eta_Hall[k][j][i-1] + pG->eta_Hall[k][j-1][i-1]);
00972       eta_6 = SQR(SQR(eta_H)) * fac;
00973 
00974       emf[k][j][i].z += eta_6 * ((J[k][j][i-2].z - 4.0*J[k][j][i-1].z
00975          + 6.0*J[k][j][i].z - 4.0*J[k][j][i+1].z + J[k][j][i+2].z) * dx41
00976                         + fac2 * (J[k][j-2][i].z - 4.0*J[k][j-1][i].z
00977          + 6.0*J[k][j][i].z - 4.0*J[k][j+1][i].z + J[k][j+2][i].z) * dy41
00978                         + fac3 * (J[k-2][j][i].z - 4.0*J[k-1][j][i].z
00979          + 6.0*J[k][j][i].z - 4.0*J[k+1][j][i].z + J[k+2][j][i].z) * dz41);
00980     }}}
00981   }
00982 
00983   return;
00984 }
00985 
00986 /*----------------------------------------------------------------------------*/
00987 /*! \fn void resistivity_init(MeshS *pM) 
00988  *  \brief Allocate temporary arrays
00989  */
00990 
00991 void resistivity_init(MeshS *pM)
00992 {
00993   int nl,nd,size1=0,size2=0,size3=0,Nx1,Nx2,Nx3;
00994   int mycase;
00995 
00996 /* Assign the function pointer for diffusivity calculation */
00997   mycase = par_geti_def("problem","CASE",1);
00998 
00999   switch (mycase)
01000   {
01001     /* single-ion prescription with constant coefficients */
01002     case 1:  get_myeta = eta_single_const; break;
01003 
01004     /* single-ion prescription with user defined diffusiviteis */
01005     case 2:  get_myeta = eta_single_user;  break;
01006 
01007     /* general prescription with user defined diffusivities */
01008     case 3:  get_myeta = eta_general_user; break;
01009 
01010     default: ath_error("[resistivity_init]: CASE must equal to 1, 2 or 3!\n");
01011              return;
01012   }
01013 
01014 /* Cycle over all Grids on this processor to find maximum Nx1, Nx2, Nx3 */
01015   for (nl=0; nl<(pM->NLevels); nl++){
01016     for (nd=0; nd<(pM->DomainsPerLevel[nl]); nd++){
01017       if (pM->Domain[nl][nd].Grid != NULL) {
01018         if (pM->Domain[nl][nd].Grid->Nx[0] > size1){
01019           size1 = pM->Domain[nl][nd].Grid->Nx[0];
01020         }
01021         if (pM->Domain[nl][nd].Grid->Nx[1] > size2){
01022           size2 = pM->Domain[nl][nd].Grid->Nx[1];
01023         }
01024         if (pM->Domain[nl][nd].Grid->Nx[2] > size3){
01025           size3 = pM->Domain[nl][nd].Grid->Nx[2];
01026         }
01027       }
01028     }
01029   }
01030 
01031   Nx1 = size1 + 2*nghost;
01032 
01033   if (pM->Nx[1] > 1){
01034     Nx2 = size2 + 2*nghost;
01035   } else {
01036     Nx2 = size2;
01037   }
01038 
01039   if (pM->Nx[2] > 1){
01040     Nx3 = size3 + 2*nghost;
01041   } else {
01042     Nx3 = size3;
01043   }
01044 
01045   if ((J = (Real3Vect***)calloc_3d_array(Nx3,Nx2,Nx1,sizeof(Real3Vect)))==NULL)
01046     goto on_error;
01047   if ((emf=(Real3Vect***)calloc_3d_array(Nx3,Nx2,Nx1,sizeof(Real3Vect)))==NULL)
01048     goto on_error;
01049 #ifndef BAROTROPIC
01050   if ((EnerFlux = (Real3Vect***)calloc_3d_array(Nx3,Nx2,Nx1, sizeof(Real3Vect)))
01051     == NULL) goto on_error;
01052 #endif
01053   if (Q_Hall > 0.0) {
01054     if ((Jcor = (Real3Vect***)calloc_3d_array(Nx3,Nx2,Nx1,sizeof(Real3Vect)))==NULL)
01055       goto on_error;
01056     if ((Bcor = (Real3Vect***)calloc_3d_array(Nx3,Nx2,Nx1,sizeof(Real3Vect)))==NULL)
01057       goto on_error;
01058     if ((emfp = (Real3Vect***)calloc_3d_array(Nx3,Nx2,Nx1,sizeof(Real3Vect)))==NULL)
01059       goto on_error;
01060   }
01061   return;
01062 
01063   on_error:
01064   resistivity_destruct();
01065   ath_error("[resistivity_init]: malloc returned a NULL pointer\n");
01066   return;
01067 }
01068 
01069 /*----------------------------------------------------------------------------*/
01070 /*! \fn void resistivity_destruct()
01071  *  \brief Free temporary arrays
01072  */
01073 void resistivity_destruct()
01074 {
01075   get_myeta = NULL;
01076 
01077   if (J != NULL) free_3d_array(J);
01078   if (emf != NULL) free_3d_array(emf);
01079 #ifndef BAROTROPIC
01080   if (EnerFlux != NULL) free_3d_array(EnerFlux);
01081 #endif
01082 
01083   if (Jcor != NULL) free_3d_array(Jcor);
01084   if (Bcor != NULL) free_3d_array(Bcor);
01085   if (emfp != NULL) free_3d_array(emfp);
01086 
01087   return;
01088 }
01089 
01090 #endif /* RESISTIVITY */