Athena: microphysics/viscosity.c Source File

Go to the documentation of this file.
00001 #include "../copyright.h"
00002 /*============================================================================*/
00003 /*! \file viscosity.c
00004  *  \brief Adds explicit viscosity terms to the momentum and energy equations
00005  *
00006  * PURPOSE: Adds explicit viscosity terms to the momentum and energy equations,
00007  *  -   dM/dt = Div(T)    
00008  *  -   dE/dt = Div(v.T)
00009  *   where 
00010  *  - T = nu_iso Grad(V) + T_Brag = TOTAl viscous stress tensor
00011  *
00012  *   Note T contains contributions from both isotropic (Navier-Stokes) and
00013  *   anisotropic (Braginskii) viscosity.  These contributions are computed in
00014  *   calls to ViscStress_* functions.
00015  *
00016  * CONTAINS PUBLIC FUNCTIONS:
00017  *- viscosity() - updates momentum and energy equations with viscous terms
00018  *- viscosity_init() - allocates memory needed
00019  *- viscosity_destruct() - frees memory used */
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 /*! \struct ViscFluxS
00033  *  \brief Structure to contain 4-components of the viscous fluxes */
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  * PRIVATE FUNCTION PROTOTYPES:
00049  *   ViscStress_iso   - computes   isotropic viscous fluxes
00050  *   ViscStress_aniso - computes anisotropic viscous fluxes
00051  *============================================================================*/
00052 
00053 void ViscStress_iso(DomainS *pD);
00054 void ViscStress_aniso(DomainS *pD);
00055 
00056 /*=========================== PUBLIC FUNCTIONS ===============================*/
00057 /*----------------------------------------------------------------------------*/
00058 /*! \fn void viscosity(DomainS *pD)
00059  *  \brief Adds explicit viscosity terms to the momentum and energy equations
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 /* Zero viscous fluxes; compute vel and div(v) at cell centers. */
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 /* Compute isotropic and anisotropic viscous fluxes.  Fluxes, V and div(V)
00150  * are global variables in this file. */
00151 
00152   if (nu_iso > 0.0)   ViscStress_iso(pD);
00153   if (nu_aniso > 0.0) ViscStress_aniso(pD);
00154 
00155 /* Update momentum and energy using x1-fluxes (dM/dt = Div(T)) */
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 /* BAROTROPIC */
00166     }
00167   }}
00168 
00169 /* Update momentum and energy using x2-fluxes (dM/dt = Div(T)) */
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 /* BAROTROPIC */
00181       }
00182     }}
00183   }
00184 
00185 /* Update momentum and energy using x3-fluxes (dM/dt = Div(T)) */
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 /* BAROTROPIC */
00197       }
00198     }}
00199   }
00200 
00201   return;
00202 }
00203 
00204 /*----------------------------------------------------------------------------*/
00205 /*! \fn void ViscStress_iso(DomainS *pD)
00206  *  \brief Calculate viscous stresses with isotropic (NS) viscosity
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 /* Add viscous fluxes in 1-direction */
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 /* BAROTROPIC */
00249     }
00250   }}
00251 
00252 /* Add viscous fluxes in 2-direction */
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 /* BAROTROPIC */
00282       }
00283     }}
00284   }
00285 
00286 /* Add viscous fluxes in 3-direction */
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 /* BAROTROPIC */
00314       }
00315     }}
00316   }
00317 
00318   return;
00319 }
00320 
00321 /*----------------------------------------------------------------------------*/
00322 /*! \fn void ViscStress_aniso(DomainS *pD) 
00323  *  \brief Calculate viscous stresses with anisotropic (Braginskii)
00324  *  viscosity */
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;  /* problem must be at least 2D */
00337 
00338 /* Compute viscous fluxes in 1-direction, centered at x1-Faces --------------- */
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 /* Monotonized Velocity gradient dVx/dy */
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 /* Monotonized Velocity gradient dVy/dy */
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 /* Monotonized Velocity gradient dVz/dy */
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 /* Monotonized Velocity gradient dVx/dz, 3D problem ONLY */
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 /* Monotonized Velocity gradient dVy/dz */
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 /* Monotonized Velocity gradient dVz/dz */
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 /* Compute field components at x1-interface */
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);  /* limit in case B=0 */
00435 
00436 /* compute BBdV and div(V) */
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 /* Add fluxes */
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 /* BAROTROPIC */
00476     }
00477   }}
00478 
00479 /* Compute viscous fluxes in 2-direction, centered at X2-Faces ---------------*/
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 /* Monotonized Velocity gradient dVx/dx */
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 /* Monotonized Velocity gradient dVy/dx */
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 /* Monotonized Velocity gradient dVz/dx */
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 /* Monotonized Velocity gradient dVx/dz */
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 /* Monotonized Velocity gradient dVy/dz */
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 /* Monotonized Velocity gradient dVz/dz */
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 /* Compute field components at x2-interface */
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); /* limit in case B=0 */
00576 
00577 /* compute BBdV and div(V) */
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 /* Add fluxes */
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 /* BAROTROPIC */
00617     }
00618   }}
00619 
00620 /* Compute viscous fluxes in 3-direction, centered at x3-Faces ---------------*/
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 /* Monotonized Velocity gradient dVx/dx */
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 /* Monotonized Velocity gradient dVy/dx */
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 /* Monotonized Velocity gradient dVz/dx */
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 /* Monotonized Velocity gradient dVx/dy */
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 /* Monotonized Velocity gradient dVy/dy */
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 /* Monotonized Velocity gradient dVz/dy */
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 /* Compute field components at x3-interface */
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); /* limit in case B=0 */
00712 
00713 /* compute BBdV and div(V) */
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 /* Add fluxes */
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 /* BAROTROPIC */
00742       }
00743     }}
00744   }
00745 
00746   return;
00747 }
00748 
00749 /*----------------------------------------------------------------------------*/
00750 /*! \fn void viscosity_init(MeshS *pM)
00751  *  \brief Allocate temporary arrays
00752  */
00753 
00754 void viscosity_init(MeshS *pM)
00755 {
00756   int nl,nd,size1=0,size2=0,size3=0,Nx1,Nx2,Nx3;
00757 
00758 /* Cycle over all Grids on this processor to find maximum Nx1, Nx2, Nx3 */
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 /*! \fn void viscosity_destruct(void) 
00809  *  \brief Free temporary arrays
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 /* VISCOSITY */