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 #include <math.h>
00023 #include <float.h>
00024 #include "../defs.h"
00025 #include "../athena.h"
00026 #include "../globals.h"
00027 #include "prototypes.h"
00028 #include "../prototypes.h"
00029
00030 #ifdef VISCOSITY
00031
00032
00033
00034 typedef struct ViscFlux_t{
00035 Real Mx;
00036 Real My;
00037 Real Mz;
00038 #ifndef BAROTROPIC
00039 Real E;
00040 #endif
00041 }ViscFluxS;
00042
00043 static ViscFluxS ***x1Flux=NULL, ***x2Flux=NULL, ***x3Flux=NULL;
00044 static Real3Vect ***Vel=NULL;
00045 static Real ***divv=NULL;
00046
00047
00048
00049
00050
00051
00052
00053 void ViscStress_iso(DomainS *pD);
00054 void ViscStress_aniso(DomainS *pD);
00055
00056
00057
00058
00059
00060
00061
00062 void viscosity(DomainS *pD)
00063 {
00064 GridS *pG = (pD->Grid);
00065 int i, is = pG->is, ie = pG->ie;
00066 int j, jl, ju, js = pG->js, je = pG->je;
00067 int k, kl, ku, ks = pG->ks, ke = pG->ke;
00068 Real x1,x2,x3,dtodx1=pG->dt/pG->dx1, dtodx2=0.0, dtodx3=0.0;
00069
00070 if (pG->Nx[1] > 1){
00071 jl = js - 2;
00072 ju = je + 2;
00073 dtodx2 = pG->dt/pG->dx2;
00074 } else {
00075 jl = js;
00076 ju = je;
00077 }
00078 if (pG->Nx[2] > 1){
00079 kl = ks - 2;
00080 ku = ke + 2;
00081 dtodx3 = pG->dt/pG->dx3;
00082 } else {
00083 kl = ks;
00084 ku = ke;
00085 }
00086
00087
00088
00089 for (k=kl; k<=ku; k++) {
00090 for (j=jl; j<=ju; j++) {
00091 for (i=is-2; i<=ie+2; i++) {
00092
00093 x1Flux[k][j][i].Mx = 0.0;
00094 x1Flux[k][j][i].My = 0.0;
00095 x1Flux[k][j][i].Mz = 0.0;
00096 #ifndef BAROTROPIC
00097 x1Flux[k][j][i].E = 0.0;
00098 #endif
00099
00100 x2Flux[k][j][i].Mx = 0.0;
00101 x2Flux[k][j][i].My = 0.0;
00102 x2Flux[k][j][i].Mz = 0.0;
00103 #ifndef BAROTROPIC
00104 x2Flux[k][j][i].E = 0.0;
00105 #endif
00106
00107 x3Flux[k][j][i].Mx = 0.0;
00108 x3Flux[k][j][i].My = 0.0;
00109 x3Flux[k][j][i].Mz = 0.0;
00110 #ifndef BAROTROPIC
00111 x3Flux[k][j][i].E = 0.0;
00112 #endif
00113
00114 Vel[k][j][i].x = pG->U[k][j][i].M1/pG->U[k][j][i].d;
00115 Vel[k][j][i].y = pG->U[k][j][i].M2/pG->U[k][j][i].d;
00116 #ifdef FARGO
00117 cc_pos(pG,i,j,k,&x1,&x2,&x3);
00118 Vel[k][j][i].y -= qshear*Omega_0*x1;
00119 #endif
00120 Vel[k][j][i].z = pG->U[k][j][i].M3/pG->U[k][j][i].d;
00121 }
00122 }}
00123
00124 for (k=kl; k<=ku; k++) {
00125 for (j=jl; j<=ju; j++) {
00126 for (i=is-1; i<=ie+1; i++) {
00127 divv[k][j][i] = (Vel[k][j][i+1].x - Vel[k][j][i-1].x)/(2.0*pG->dx1);
00128 }
00129 }}
00130
00131 if (pG->Nx[1] > 1) {
00132 for (k=kl; k<=ku; k++) {
00133 for (j=js-1; j<=je+1; j++) {
00134 for (i=is-1; i<=ie+1; i++) {
00135 divv[k][j][i] += (Vel[k][j+1][i].y - Vel[k][j-1][i].y)/(2.0*pG->dx2);
00136 }
00137 }}
00138 }
00139
00140 if (pG->Nx[2] > 1) {
00141 for (k=ks-1; k<=ke+1; k++) {
00142 for (j=js-1; j<=je+1; j++) {
00143 for (i=is-1; i<=ie+1; i++) {
00144 divv[k][j][i] += (Vel[k+1][j][i].z - Vel[k-1][j][i].z)/(2.0*pG->dx3);
00145 }
00146 }}
00147 }
00148
00149
00150
00151
00152 if (nu_iso > 0.0) ViscStress_iso(pD);
00153 if (nu_aniso > 0.0) ViscStress_aniso(pD);
00154
00155
00156
00157 for (k=ks; k<=ke; k++) {
00158 for (j=js; j<=je; j++) {
00159 for (i=is; i<=ie; i++) {
00160 pG->U[k][j][i].M1 += dtodx1*(x1Flux[k][j][i+1].Mx-x1Flux[k][j][i].Mx);
00161 pG->U[k][j][i].M2 += dtodx1*(x1Flux[k][j][i+1].My-x1Flux[k][j][i].My);
00162 pG->U[k][j][i].M3 += dtodx1*(x1Flux[k][j][i+1].Mz-x1Flux[k][j][i].Mz);
00163 #ifndef BAROTROPIC
00164 pG->U[k][j][i].E += dtodx1*(x1Flux[k][j][i+1].E -x1Flux[k][j][i].E );
00165 #endif
00166 }
00167 }}
00168
00169
00170
00171 if (pG->Nx[1] > 1){
00172 for (k=ks; k<=ke; k++) {
00173 for (j=js; j<=je; j++) {
00174 for (i=is; i<=ie; i++) {
00175 pG->U[k][j][i].M1 += dtodx2*(x2Flux[k][j+1][i].Mx-x2Flux[k][j][i].Mx);
00176 pG->U[k][j][i].M2 += dtodx2*(x2Flux[k][j+1][i].My-x2Flux[k][j][i].My);
00177 pG->U[k][j][i].M3 += dtodx2*(x2Flux[k][j+1][i].Mz-x2Flux[k][j][i].Mz);
00178 #ifndef BAROTROPIC
00179 pG->U[k][j][i].E += dtodx2*(x2Flux[k][j+1][i].E -x2Flux[k][j][i].E );
00180 #endif
00181 }
00182 }}
00183 }
00184
00185
00186
00187 if (pG->Nx[2] > 1){
00188 for (k=ks; k<=ke; k++) {
00189 for (j=js; j<=je; j++) {
00190 for (i=is; i<=ie; i++) {
00191 pG->U[k][j][i].M1 += dtodx3*(x3Flux[k+1][j][i].Mx-x3Flux[k][j][i].Mx);
00192 pG->U[k][j][i].M2 += dtodx3*(x3Flux[k+1][j][i].My-x3Flux[k][j][i].My);
00193 pG->U[k][j][i].M3 += dtodx3*(x3Flux[k+1][j][i].Mz-x3Flux[k][j][i].Mz);
00194 #ifndef BAROTROPIC
00195 pG->U[k][j][i].E += dtodx3*(x3Flux[k+1][j][i].E -x3Flux[k][j][i].E );
00196 #endif
00197 }
00198 }}
00199 }
00200
00201 return;
00202 }
00203
00204
00205
00206
00207
00208
00209 void ViscStress_iso(DomainS *pD)
00210 {
00211 GridS *pG = (pD->Grid);
00212 ViscFluxS VStress;
00213 int i, is = pG->is, ie = pG->ie;
00214 int j, js = pG->js, je = pG->je;
00215 int k, ks = pG->ks, ke = pG->ke;
00216 Real nud;
00217
00218
00219
00220 for (k=ks; k<=ke; k++) {
00221 for (j=js; j<=je; j++) {
00222 for (i=is; i<=ie+1; i++) {
00223 VStress.Mx = 2.0*(Vel[k][j][i].x - Vel[k][j][i-1].x)/pG->dx1
00224 - ONE_3RD*(divv[k][j][i] + divv[k][j][i-1]);
00225
00226 VStress.My = (Vel[k][j][i].y - Vel[k][j][i-1].y)/pG->dx1;
00227 if (pG->Nx[1] > 1) {
00228 VStress.My += (0.25/pG->dx2)*((Vel[k][j+1][i].x + Vel[k][j+1][i-1].x) -
00229 (Vel[k][j-1][i].x + Vel[k][j-1][i-1].x));
00230 }
00231
00232 VStress.Mz = (Vel[k][j][i].z - Vel[k][j][i-1].z)/pG->dx1;
00233 if (pG->Nx[2] > 1) {
00234 VStress.Mz += (0.25/pG->dx3)*((Vel[k+1][j][i].x + Vel[k+1][j][i-1].x) -
00235 (Vel[k-1][j][i].x + Vel[k-1][j][i-1].x));
00236 }
00237
00238 nud = nu_iso*0.5*(pG->U[k][j][i].d + pG->U[k][j][i-1].d);
00239 x1Flux[k][j][i].Mx += nud*VStress.Mx;
00240 x1Flux[k][j][i].My += nud*VStress.My;
00241 x1Flux[k][j][i].Mz += nud*VStress.Mz;
00242
00243 #ifndef BAROTROPIC
00244 x1Flux[k][j][i].E +=
00245 0.5*(Vel[k][j][i-1].x + Vel[k][j][i].x)*VStress.Mx +
00246 0.5*(Vel[k][j][i-1].y + Vel[k][j][i].y)*VStress.My +
00247 0.5*(Vel[k][j][i-1].z + Vel[k][j][i].z)*VStress.Mz;
00248 #endif
00249 }
00250 }}
00251
00252
00253
00254 if (pG->Nx[1] > 1) {
00255 for (k=ks; k<=ke; k++) {
00256 for (j=js; j<=je+1; j++) {
00257 for (i=is; i<=ie; i++) {
00258 VStress.Mx = (Vel[k][j][i].x - Vel[k][j-1][i].x)/pG->dx2
00259 + ((Vel[k][j][i+1].y + Vel[k][j-1][i+1].y) -
00260 (Vel[k][j][i-1].y + Vel[k][j-1][i-1].y))/(4.0*pG->dx1);
00261
00262 VStress.My = 2.0*(Vel[k][j][i].y - Vel[k][j-1][i].y)/pG->dx2
00263 - ONE_3RD*(divv[k][j][i] + divv[k][j-1][i]);
00264
00265 VStress.Mz = (Vel[k][j][i].z - Vel[k][j-1][i].z)/pG->dx2;
00266 if (pG->Nx[2] > 1) {
00267 VStress.Mz += (0.25/pG->dx3)*((Vel[k+1][j][i].y + Vel[k+1][j-1][i].y) -
00268 (Vel[k-1][j][i].y + Vel[k-1][j-1][i].y));
00269 }
00270
00271 nud = nu_iso*0.5*(pG->U[k][j][i].d + pG->U[k][j-1][i].d);
00272 x2Flux[k][j][i].Mx += nud*VStress.Mx;
00273 x2Flux[k][j][i].My += nud*VStress.My;
00274 x2Flux[k][j][i].Mz += nud*VStress.Mz;
00275
00276 #ifndef BAROTROPIC
00277 x2Flux[k][j][i].E +=
00278 0.5*(Vel[k][j-1][i].x + Vel[k][j][i].x)*VStress.Mx +
00279 0.5*(Vel[k][j-1][i].y + Vel[k][j][i].y)*VStress.My +
00280 0.5*(Vel[k][j-1][i].z + Vel[k][j][i].z)*VStress.Mz;
00281 #endif
00282 }
00283 }}
00284 }
00285
00286
00287
00288 if (pG->Nx[2] > 1) {
00289 for (k=ks; k<=ke+1; k++) {
00290 for (j=js; j<=je; j++) {
00291 for (i=is; i<=ie; i++) {
00292 VStress.Mx = (Vel[k][j][i].x - Vel[k-1][j][i].x)/pG->dx3
00293 + ((Vel[k][j][i+1].z + Vel[k-1][j][i+1].z) -
00294 (Vel[k][j][i-1].z + Vel[k-1][j][i-1].z))/(4.0*pG->dx1);
00295
00296 VStress.My = (Vel[k][j][i].y - Vel[k-1][j][i].y)/pG->dx3
00297 + ((Vel[k][j+1][i].z + Vel[k-1][j+1][i].z) -
00298 (Vel[k][j-1][i].z + Vel[k-1][j-1][i].z))/(4.0*pG->dx2);
00299
00300 VStress.Mz = 2.0*(Vel[k][j][i].z - Vel[k-1][j][i].z)/pG->dx3
00301 - ONE_3RD*(divv[k][j][i] + divv[k-1][j][i]);
00302
00303 nud = nu_iso*0.5*(pG->U[k][j][i].d + pG->U[k-1][j][i].d);
00304 x3Flux[k][j][i].Mx += nud*VStress.Mx;
00305 x3Flux[k][j][i].My += nud*VStress.My;
00306 x3Flux[k][j][i].Mz += nud*VStress.Mz;
00307
00308 #ifndef BAROTROPIC
00309 x3Flux[k][j][i].E +=
00310 0.5*(Vel[k-1][j][i].x + Vel[k][j][i].x)*VStress.Mx +
00311 0.5*(Vel[k-1][j][i].y + Vel[k][j][i].y)*VStress.My +
00312 0.5*(Vel[k-1][j][i].z + Vel[k][j][i].z)*VStress.Mz;
00313 #endif
00314 }
00315 }}
00316 }
00317
00318 return;
00319 }
00320
00321
00322
00323
00324
00325
00326 void ViscStress_aniso(DomainS *pD)
00327 {
00328 GridS *pG = (pD->Grid);
00329 ViscFluxS VStress;
00330 int i, is = pG->is, ie = pG->ie;
00331 int j, js = pG->js, je = pG->je;
00332 int k, ks = pG->ks, ke = pG->ke;
00333 Real Bx,By,Bz,B02,dVc,dVl,dVr,lim_slope,BBdV,divV,qa,nud;
00334 Real dVxdx,dVydx,dVzdx,dVxdy,dVydy,dVzdy,dVxdz,dVydz,dVzdz;
00335
00336 if (pD->Nx[1] == 1) return;
00337
00338
00339
00340 for (k=ks; k<=ke; k++) {
00341 for (j=js; j<=je; j++) {
00342 for (i=is; i<=ie+1; i++) {
00343
00344
00345 dVr = 0.5*((Vel[k][j+1][i-1].x + Vel[k][j+1][i].x) -
00346 (Vel[k][j ][i-1].x + Vel[k][j ][i].x));
00347 dVl = 0.5*((Vel[k][j ][i-1].x + Vel[k][j ][i].x) -
00348 (Vel[k][j-1][i-1].x + Vel[k][j-1][i].x));
00349 dVc = dVr + dVl;
00350
00351 dVxdy = 0.0;
00352 if (dVl*dVr > 0.0) {
00353 lim_slope = MIN(fabs(dVl),fabs(dVr));
00354 dVxdy = SIGN(dVc)*MIN(0.5*fabs(dVc),2.0*lim_slope)/pG->dx2;
00355 }
00356
00357
00358 dVr = 0.5*((Vel[k][j+1][i-1].y + Vel[k][j+1][i].y) -
00359 (Vel[k][j ][i-1].y + Vel[k][j ][i].y));
00360 dVl = 0.5*((Vel[k][j ][i-1].y + Vel[k][j ][i].y) -
00361 (Vel[k][j-1][i-1].y + Vel[k][j-1][i].y));
00362 dVc = dVr + dVl;
00363
00364 dVydy = 0.0;
00365 if (dVl*dVr > 0.0) {
00366 lim_slope = MIN(fabs(dVl),fabs(dVr));
00367 dVydy = SIGN(dVc)*MIN(0.5*fabs(dVc),2.0*lim_slope)/pG->dx2;
00368 }
00369
00370
00371 dVr = 0.5*((Vel[k][j+1][i-1].z + Vel[k][j+1][i].z) -
00372 (Vel[k][j ][i-1].z + Vel[k][j ][i].z));
00373 dVl = 0.5*((Vel[k][j ][i-1].z + Vel[k][j ][i].z) -
00374 (Vel[k][j-1][i-1].z + Vel[k][j-1][i].z));
00375 dVc = dVr + dVl;
00376
00377 dVzdy = 0.0;
00378 if (dVl*dVr > 0.0) {
00379 lim_slope = MIN(fabs(dVl),fabs(dVr));
00380 dVzdy = SIGN(dVc)*MIN(0.5*fabs(dVc),2.0*lim_slope)/pG->dx2;
00381 }
00382
00383
00384 if (pD->Nx[2] > 1) {
00385 dVr = 0.5*((Vel[k+1][j][i-1].x + Vel[k+1][j][i].x) -
00386 (Vel[k ][j][i-1].x + Vel[k ][j][i].x));
00387 dVl = 0.5*((Vel[k ][j][i-1].x + Vel[k ][j][i].x) -
00388 (Vel[k-1][j][i-1].x + Vel[k-1][j][i].x));
00389 dVc = dVr + dVl;
00390
00391 dVxdz = 0.0;
00392 if (dVl*dVr > 0.0) {
00393 lim_slope = MIN(fabs(dVl),fabs(dVr));
00394 dVxdz = SIGN(dVc)*MIN(0.5*fabs(dVc),2.0*lim_slope)/pG->dx3;
00395 }
00396 }
00397
00398
00399 if (pD->Nx[2] > 1) {
00400 dVr = 0.5*((Vel[k+1][j][i-1].y + Vel[k+1][j][i].y) -
00401 (Vel[k ][j][i-1].y + Vel[k ][j][i].y));
00402 dVl = 0.5*((Vel[k ][j][i-1].y + Vel[k ][j][i].y) -
00403 (Vel[k-1][j][i-1].y + Vel[k-1][j][i].y));
00404 dVc = dVr + dVl;
00405
00406 dVydz = 0.0;
00407 if (dVl*dVr > 0.0) {
00408 lim_slope = MIN(fabs(dVl),fabs(dVr));
00409 dVydz = SIGN(dVc)*MIN(0.5*fabs(dVc),2.0*lim_slope)/pG->dx3;
00410 }
00411 }
00412
00413
00414 if (pD->Nx[2] > 1) {
00415 dVr = 0.5*((Vel[k+1][j][i-1].z + Vel[k+1][j][i].z) -
00416 (Vel[k ][j][i-1].z + Vel[k ][j][i].z));
00417 dVl = 0.5*((Vel[k ][j][i-1].z + Vel[k ][j][i].z) -
00418 (Vel[k-1][j][i-1].z + Vel[k-1][j][i].z));
00419 dVc = dVr + dVl;
00420
00421 dVzdz = 0.0;
00422 if (dVl*dVr > 0.0) {
00423 lim_slope = MIN(fabs(dVl),fabs(dVr));
00424 dVzdz = SIGN(dVc)*MIN(0.5*fabs(dVc),2.0*lim_slope)/pG->dx3;
00425 }
00426 }
00427
00428
00429
00430 Bx = pG->B1i[k][j][i];
00431 By = 0.5*(pG->U[k][j][i].B2c + pG->U[k][j][i-1].B2c);
00432 Bz = 0.5*(pG->U[k][j][i].B3c + pG->U[k][j][i-1].B3c);
00433 B02 = Bx*Bx + By*By + Bz*Bz;
00434 B02 = MAX(B02,TINY_NUMBER);
00435
00436
00437
00438 if (pD->Nx[2] == 1) {
00439 BBdV =
00440 Bx*(Bx*(Vel[k][j][i].x-Vel[k][j][i-1].x)/pG->dx1 + By*dVxdy) +
00441 By*(Bx*(Vel[k][j][i].y-Vel[k][j][i-1].y)/pG->dx1 + By*dVydy) +
00442 Bz*(Bx*(Vel[k][j][i].z-Vel[k][j][i-1].z)/pG->dx1 + By*dVzdy);
00443 BBdV /= B02;
00444
00445 divV = (Vel[k][j][i].x-Vel[k][j][i-1].x)/pG->dx1 + dVydy;
00446
00447 } else {
00448 BBdV =
00449 Bx*(Bx*(Vel[k][j][i].x-Vel[k][j][i-1].x)/pG->dx1 +By*dVxdy +Bz*dVxdz)+
00450 By*(Bx*(Vel[k][j][i].y-Vel[k][j][i-1].y)/pG->dx1 +By*dVydy +Bz*dVydz)+
00451 Bz*(Bx*(Vel[k][j][i].z-Vel[k][j][i-1].z)/pG->dx1 +By*dVzdy +Bz*dVzdz);
00452 BBdV /= B02;
00453
00454 divV = (Vel[k][j][i].x-Vel[k][j][i-1].x)/pG->dx1 + dVydy + dVzdz;
00455 }
00456
00457
00458
00459 nud = nu_aniso*0.5*(pG->U[k][j][i].d + pG->U[k][j][i-1].d);
00460 qa = nud*(BBdV - ONE_3RD*divV);
00461
00462 VStress.Mx = qa*(3.0*Bx*Bx/B02 - 1.0);
00463 VStress.My = qa*(3.0*By*Bx/B02);
00464 VStress.Mz = qa*(3.0*Bz*Bx/B02);
00465
00466 x1Flux[k][j][i].Mx += VStress.Mx;
00467 x1Flux[k][j][i].My += VStress.My;
00468 x1Flux[k][j][i].Mz += VStress.Mz;
00469
00470 #ifndef BAROTROPIC
00471 x1Flux[k][j][i].E =
00472 0.5*(Vel[k][j][i-1].x + Vel[k][j][i].x)*VStress.Mx +
00473 0.5*(Vel[k][j][i-1].y + Vel[k][j][i].y)*VStress.My +
00474 0.5*(Vel[k][j][i-1].z + Vel[k][j][i].z)*VStress.Mz;
00475 #endif
00476 }
00477 }}
00478
00479
00480
00481 for (k=ks; k<=ke; k++) {
00482 for (j=js; j<=je+1; j++) {
00483 for (i=is; i<=ie; i++) {
00484
00485
00486 dVr = 0.5*((Vel[k][j-1][i+1].x + Vel[k][j][i+1].x) -
00487 (Vel[k][j-1][i ].x + Vel[k][j][i ].x));
00488 dVl = 0.5*((Vel[k][j-1][i ].x + Vel[k][j][i ].x) -
00489 (Vel[k][j-1][i-1].x + Vel[k][j][i-1].x));
00490 dVc = dVr + dVl;
00491
00492 dVxdx = 0.0;
00493 if (dVl*dVr > 0.0) {
00494 lim_slope = MIN(fabs(dVl),fabs(dVr));
00495 dVxdx = SIGN(dVc)*MIN(0.5*fabs(dVc),2.0*lim_slope)/pG->dx1;
00496 }
00497
00498
00499 dVr = 0.5*((Vel[k][j-1][i+1].y + Vel[k][j][i+1].y) -
00500 (Vel[k][j-1][i ].y + Vel[k][j][i ].y));
00501 dVl = 0.5*((Vel[k][j-1][i ].y + Vel[k][j][i ].y) -
00502 (Vel[k][j-1][i-1].y + Vel[k][j][i-1].y));
00503 dVc = dVr + dVl;
00504
00505 dVydx = 0.0;
00506 if (dVl*dVr > 0.0) {
00507 lim_slope = MIN(fabs(dVl),fabs(dVr));
00508 dVydx = SIGN(dVc)*MIN(0.5*fabs(dVc),2.0*lim_slope)/pG->dx1;
00509 }
00510
00511
00512 dVr = 0.5*((Vel[k][j-1][i+1].z + Vel[k][j][i+1].z) -
00513 (Vel[k][j-1][i ].z + Vel[k][j][i ].z));
00514 dVl = 0.5*((Vel[k][j-1][i ].z + Vel[k][j][i ].z) -
00515 (Vel[k][j-1][i-1].z + Vel[k][j][i-1].z));
00516 dVc = dVr + dVl;
00517
00518 dVzdx = 0.0;
00519 if (dVl*dVr > 0.0) {
00520 lim_slope = MIN(fabs(dVl),fabs(dVr));
00521 dVzdx = SIGN(dVc)*MIN(0.5*fabs(dVc),2.0*lim_slope)/pG->dx1;
00522 }
00523
00524
00525 if (pD->Nx[2] > 1) {
00526 dVr = 0.5*((Vel[k+1][j-1][i].x + Vel[k+1][j][i].x) -
00527 (Vel[k ][j-1][i].x + Vel[k ][j][i].x));
00528 dVl = 0.5*((Vel[k ][j-1][i].x + Vel[k ][j][i].x) -
00529 (Vel[k-1][j-1][i].x + Vel[k-1][j][i].x));
00530 dVc = dVr + dVl;
00531
00532 dVxdz = 0.0;
00533 if (dVl*dVr > 0.0) {
00534 lim_slope = MIN(fabs(dVl),fabs(dVr));
00535 dVxdz = SIGN(dVc)*MIN(0.5*fabs(dVc),2.0*lim_slope)/pG->dx3;
00536 }
00537 }
00538
00539
00540 if (pD->Nx[2] > 1) {
00541 dVr = 0.5*((Vel[k+1][j-1][i].y + Vel[k+1][j][i].y) -
00542 (Vel[k ][j-1][i].y + Vel[k ][j][i].y));
00543 dVl = 0.5*((Vel[k ][j-1][i].y + Vel[k ][j][i].y) -
00544 (Vel[k-1][j-1][i].y + Vel[k-1][j][i].y));
00545 dVc = dVr + dVl;
00546
00547 dVydz = 0.0;
00548 if (dVl*dVr > 0.0) {
00549 lim_slope = MIN(fabs(dVl),fabs(dVr));
00550 dVydz = SIGN(dVc)*MIN(0.5*fabs(dVc),2.0*lim_slope)/pG->dx3;
00551 }
00552 }
00553
00554
00555 if (pD->Nx[2] > 1) {
00556 dVr = 0.5*((Vel[k+1][j-1][i].z + Vel[k+1][j][i].z) -
00557 (Vel[k ][j-1][i].z + Vel[k ][j][i].z));
00558 dVl = 0.5*((Vel[k ][j-1][i].z + Vel[k ][j][i].z) -
00559 (Vel[k-1][j-1][i].z + Vel[k-1][j][i].z));
00560 dVc = dVr + dVl;
00561
00562 dVzdz = 0.0;
00563 if (dVl*dVr > 0.0) {
00564 lim_slope = MIN(fabs(dVl),fabs(dVr));
00565 dVzdz = SIGN(dVc)*MIN(0.5*fabs(dVc),2.0*lim_slope)/pG->dx3;
00566 }
00567 }
00568
00569
00570
00571 Bx = 0.5*(pG->U[k][j][i].B1c + pG->U[k][j-1][i].B1c);
00572 By = pG->B2i[k][j][i];
00573 Bz = 0.5*(pG->U[k][j][i].B3c + pG->U[k][j-1][i].B3c);
00574 B02 = Bx*Bx + By*By + Bz*Bz;
00575 B02 = MAX(B02,TINY_NUMBER);
00576
00577
00578
00579 if (pD->Nx[2] == 1) {
00580 BBdV =
00581 Bx*(Bx*dVxdx + By*(Vel[ks][j][i].x-Vel[ks][j-1][i].x)/pG->dx2) +
00582 By*(Bx*dVydx + By*(Vel[ks][j][i].y-Vel[ks][j-1][i].y)/pG->dx2) +
00583 Bz*(Bx*dVzdx + By*(Vel[ks][j][i].z-Vel[ks][j-1][i].z)/pG->dx2);
00584 BBdV /= B02;
00585
00586 divV = dVxdx + (Vel[ks][j][i].y-Vel[ks][j-1][i].y)/pG->dx2;
00587
00588 } else {
00589 BBdV =
00590 Bx*(Bx*dVxdx +By*(Vel[k][j][i].x-Vel[k][j-1][i].x)/pG->dx2 +Bz*dVxdz)+
00591 By*(Bx*dVydx +By*(Vel[k][j][i].y-Vel[k][j-1][i].y)/pG->dx2 +Bz*dVydz)+
00592 Bz*(Bx*dVzdx +By*(Vel[k][j][i].z-Vel[k][j-1][i].z)/pG->dx2 +Bz*dVzdz);
00593 BBdV /= B02;
00594
00595 divV = dVxdx + (Vel[k][j][i].y-Vel[k][j-1][i].y)/pG->dx2 + dVzdz;
00596 }
00597
00598
00599
00600 nud = nu_aniso*0.5*(pG->U[k][j][i].d + pG->U[k][j-1][i].d);
00601 qa = nud*(BBdV - ONE_3RD*divV);
00602
00603 VStress.Mx = qa*(3.0*Bx*By/B02);
00604 VStress.My = qa*(3.0*By*By/B02 - 1.0);
00605 VStress.Mz = qa*(3.0*Bz*By/B02);
00606
00607 x2Flux[k][j][i].Mx += VStress.Mx;
00608 x2Flux[k][j][i].My += VStress.My;
00609 x2Flux[k][j][i].Mz += VStress.Mz;
00610
00611 #ifndef BAROTROPIC
00612 x2Flux[k][j][i].E +=
00613 0.5*(Vel[k][j-1][i].x + Vel[k][j][i].x)*VStress.Mx +
00614 0.5*(Vel[k][j-1][i].y + Vel[k][j][i].y)*VStress.My +
00615 0.5*(Vel[k][j-1][i].z + Vel[k][j][i].z)*VStress.Mz;
00616 #endif
00617 }
00618 }}
00619
00620
00621
00622 if (pD->Nx[2] > 1) {
00623 for (k=ks; k<=ke+1; k++) {
00624 for (j=js; j<=je; j++) {
00625 for (i=is; i<=ie; i++) {
00626
00627
00628 dVr = 0.5*((Vel[k-1][j][i+1].x + Vel[k][j][i+1].x) -
00629 (Vel[k-1][j][i ].x + Vel[k][j][i ].x));
00630 dVl = 0.5*((Vel[k-1][j][i ].x + Vel[k][j][i ].x) -
00631 (Vel[k-1][j][i-1].x + Vel[k][j][i-1].x));
00632 dVc = dVr + dVl;
00633
00634 dVxdx = 0.0;
00635 if (dVl*dVr > 0.0) {
00636 lim_slope = MIN(fabs(dVl),fabs(dVr));
00637 dVxdx = SIGN(dVc)*MIN(0.5*fabs(dVc),2.0*lim_slope)/pG->dx1;
00638 }
00639
00640
00641 dVr = 0.5*((Vel[k-1][j][i+1].y + Vel[k][j][i+1].y) -
00642 (Vel[k-1][j][i ].y + Vel[k][j][i ].y));
00643 dVl = 0.5*((Vel[k-1][j][i ].y + Vel[k][j][i ].y) -
00644 (Vel[k-1][j][i-1].y + Vel[k][j][i-1].y));
00645 dVc = dVr + dVl;
00646
00647 dVydx = 0.0;
00648 if (dVl*dVr > 0.0) {
00649 lim_slope = MIN(fabs(dVl),fabs(dVr));
00650 dVydx = SIGN(dVc)*MIN(0.5*fabs(dVc),2.0*lim_slope)/pG->dx1;
00651 }
00652
00653
00654 dVr = 0.5*((Vel[k-1][j][i+1].z + Vel[k][j][i+1].z) -
00655 (Vel[k-1][j][i ].z + Vel[k][j][i ].z));
00656 dVl = 0.5*((Vel[k-1][j][i ].z + Vel[k][j][i ].z) -
00657 (Vel[k-1][j][i-1].z + Vel[k][j][i-1].z));
00658 dVc = dVr + dVl;
00659
00660 dVzdx = 0.0;
00661 if (dVl*dVr > 0.0) {
00662 lim_slope = MIN(fabs(dVl),fabs(dVr));
00663 dVzdx = SIGN(dVc)*MIN(0.5*fabs(dVc),2.0*lim_slope)/pG->dx1;
00664 }
00665
00666
00667 dVr = 0.5*((Vel[k-1][j+1][i].x + Vel[k][j+1][i].x) -
00668 (Vel[k-1][j ][i].x + Vel[k][j ][i].x));
00669 dVl = 0.5*((Vel[k-1][j ][i].x + Vel[k][j ][i].x) -
00670 (Vel[k-1][j-1][i].x + Vel[k][j-1][i].x));
00671 dVc = dVr + dVl;
00672
00673 dVxdy = 0.0;
00674 if (dVl*dVr > 0.0) {
00675 lim_slope = MIN(fabs(dVl),fabs(dVr));
00676 dVxdy = SIGN(dVc)*MIN(0.5*fabs(dVc),2.0*lim_slope)/pG->dx2;
00677 }
00678
00679
00680 dVr = 0.5*((Vel[k-1][j+1][i].y + Vel[k][j+1][i].y) -
00681 (Vel[k-1][j ][i].y + Vel[k][j ][i].y));
00682 dVl = 0.5*((Vel[k-1][j ][i].y + Vel[k][j ][i].y) -
00683 (Vel[k-1][j-1][i].y + Vel[k][j-1][i].y));
00684 dVc = dVr + dVl;
00685
00686 dVydy = 0.0;
00687 if (dVl*dVr > 0.0) {
00688 lim_slope = MIN(fabs(dVl),fabs(dVr));
00689 dVydy = SIGN(dVc)*MIN(0.5*fabs(dVc),2.0*lim_slope)/pG->dx2;
00690 }
00691
00692
00693 dVr = 0.5*((Vel[k-1][j+1][i].z + Vel[k][j+1][i].z) -
00694 (Vel[k-1][j ][i].z + Vel[k][j ][i].z));
00695 dVl = 0.5*((Vel[k-1][j ][i].z + Vel[k][j ][i].z) -
00696 (Vel[k-1][j-1][i].z + Vel[k][j-1][i].z));
00697 dVc = dVr + dVl;
00698
00699 dVzdy = 0.0;
00700 if (dVl*dVr > 0.0) {
00701 lim_slope = MIN(fabs(dVl),fabs(dVr));
00702 dVzdy = SIGN(dVc)*MIN(0.5*fabs(dVc),2.0*lim_slope)/pG->dx2;
00703 }
00704
00705
00706
00707 Bx = 0.5*(pG->U[k][j][i].B1c + pG->U[k-1][j][i].B1c);
00708 By = 0.5*(pG->U[k][j][i].B2c + pG->U[k-1][j][i].B2c);
00709 Bz = pG->B3i[k][j][i];
00710 B02 = Bx*Bx + By*By + Bz*Bz;
00711 B02 = MAX(B02,TINY_NUMBER);
00712
00713
00714
00715 BBdV =
00716 Bx*(Bx*dVxdx +By*dVxdy +Bz*(Vel[k][j][i].x-Vel[k-1][j][i].x)/pG->dx3)+
00717 By*(Bx*dVydx +By*dVydy +Bz*(Vel[k][j][i].y-Vel[k-1][j][i].y)/pG->dx3)+
00718 Bz*(Bx*dVzdx +By*dVzdy +Bz*(Vel[k][j][i].z-Vel[k-1][j][i].z)/pG->dx3);
00719 BBdV /= B02;
00720
00721 divV = dVxdx + dVydy + (Vel[k][j][i].z-Vel[k-1][j][i].z)/pG->dx3;
00722
00723
00724
00725 nud = nu_aniso*0.5*(pG->U[k][j][i].d + pG->U[k-1][j][i].d);
00726 qa = nud*(BBdV - ONE_3RD*divV);
00727
00728 VStress.Mx = qa*(3.0*Bx*Bz/B02);
00729 VStress.My = qa*(3.0*By*Bz/B02);
00730 VStress.Mz = qa*(3.0*Bz*Bz/B02 - 1.0);
00731
00732 x3Flux[k][j][i].Mx += VStress.Mx;
00733 x3Flux[k][j][i].My += VStress.My;
00734 x3Flux[k][j][i].Mz += VStress.Mz;
00735
00736 #ifndef BAROTROPIC
00737 x3Flux[k][j][i].E +=
00738 0.5*(Vel[k-1][j][i].x + Vel[k][j][i].x)*VStress.Mx +
00739 0.5*(Vel[k-1][j][i].y + Vel[k][j][i].y)*VStress.My +
00740 0.5*(Vel[k-1][j][i].z + Vel[k][j][i].z)*VStress.Mz;
00741 #endif
00742 }
00743 }}
00744 }
00745
00746 return;
00747 }
00748
00749
00750
00751
00752
00753
00754 void viscosity_init(MeshS *pM)
00755 {
00756 int nl,nd,size1=0,size2=0,size3=0,Nx1,Nx2,Nx3;
00757
00758
00759 for (nl=0; nl<(pM->NLevels); nl++){
00760 for (nd=0; nd<(pM->DomainsPerLevel[nl]); nd++){
00761 if (pM->Domain[nl][nd].Grid != NULL) {
00762 if (pM->Domain[nl][nd].Grid->Nx[0] > size1){
00763 size1 = pM->Domain[nl][nd].Grid->Nx[0];
00764 }
00765 if (pM->Domain[nl][nd].Grid->Nx[1] > size2){
00766 size2 = pM->Domain[nl][nd].Grid->Nx[1];
00767 }
00768 if (pM->Domain[nl][nd].Grid->Nx[2] > size3){
00769 size3 = pM->Domain[nl][nd].Grid->Nx[2];
00770 }
00771 }
00772 }
00773 }
00774
00775 Nx1 = size1 + 2*nghost;
00776
00777 if (pM->Nx[1] > 1){
00778 Nx2 = size2 + 2*nghost;
00779 } else {
00780 Nx2 = size2;
00781 }
00782
00783 if (pM->Nx[2] > 1){
00784 Nx3 = size3 + 2*nghost;
00785 } else {
00786 Nx3 = size3;
00787 }
00788
00789 if ((x1Flux = (ViscFluxS***)calloc_3d_array(Nx3,Nx2,Nx1, sizeof(ViscFluxS)))
00790 == NULL) goto on_error;
00791 if ((x2Flux = (ViscFluxS***)calloc_3d_array(Nx3,Nx2,Nx1, sizeof(ViscFluxS)))
00792 == NULL) goto on_error;
00793 if ((x3Flux = (ViscFluxS***)calloc_3d_array(Nx3,Nx2,Nx1, sizeof(ViscFluxS)))
00794 == NULL) goto on_error;
00795 if ((Vel = (Real3Vect***)calloc_3d_array(Nx3,Nx2,Nx1, sizeof(Real3Vect)))
00796 == NULL) goto on_error;
00797 if ((divv = (Real***)calloc_3d_array(Nx3,Nx2,Nx1, sizeof(Real))) == NULL)
00798 goto on_error;
00799 return;
00800
00801 on_error:
00802 viscosity_destruct();
00803 ath_error("[viscosity_init]: malloc returned a NULL pointer\n");
00804 return;
00805 }
00806
00807
00808
00809
00810
00811
00812 void viscosity_destruct(void)
00813 {
00814 if (x1Flux != NULL) free_3d_array(x1Flux);
00815 if (x2Flux != NULL) free_3d_array(x2Flux);
00816 if (x3Flux != NULL) free_3d_array(x3Flux);
00817 if (Vel != NULL) free_3d_array(Vel);
00818 if (divv != NULL) free_3d_array(divv);
00819 return;
00820 }
00821 #endif