00001 #include "../copyright.h"
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
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
00038
00039
00040 static Real3Vect ***J=NULL, ***emf=NULL, ***EnerFlux=NULL;
00041
00042
00043 static Real3Vect ***Jcor=NULL, ***Bcor=NULL, ***emfp=NULL;
00044
00045
00046
00047
00048
00049
00050
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
00062
00063
00064
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
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
00118
00119
00120
00121
00122
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
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
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
00160
00161
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
00167
00168
00169
00170
00171
00172
00173
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
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
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
00234
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
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
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
00265
00266
00267
00268
00269
00270
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
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
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
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
00302 pG->B3i[ks][j][i] = pG->U[ks][j][i].B3c;
00303 }
00304 }
00305 }
00306
00307
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
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
00350
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
00366
00367
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
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
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
00429
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
00459 hyper_diffusion6(pD);
00460
00461
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
00477 EField_Hall_sub(pD, Bcor, J, emfp, 3);
00478
00479
00480
00481
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
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
00518
00519
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
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
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
00556 EField_Hall_sub(pD, Bcor, Jcor, emf, 1);
00557
00558 return;
00559
00560 }
00561
00562
00563
00564
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
00579
00580
00581
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
00593
00594
00595
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
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
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
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
00632
00633
00634
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
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
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
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
00685
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
00701
00702
00703
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
00722
00723
00724
00725
00726 if (ndim == 2){
00727 for (j=js; j<=je+1; j++) {
00728 for (i=is; i<=ie+1; i++) {
00729
00730
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
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
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
00786
00787
00788
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
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
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
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
00863
00864
00865
00866
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
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
00905 if (ndim == 2) {
00906 for (j=js; j<=je+1; j++) {
00907 for (i=is; i<=ie+1; i++) {
00908
00909
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
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
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
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
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
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
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
00988
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
00997 mycase = par_geti_def("problem","CASE",1);
00998
00999 switch (mycase)
01000 {
01001
01002 case 1: get_myeta = eta_single_const; break;
01003
01004
01005 case 2: get_myeta = eta_single_user; break;
01006
01007
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
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
01071
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