Athena: microphysics/conduction.c Source File

Go to the documentation of this file.
00001 #include "../copyright.h"
00002 /*============================================================================*/
00003 /*! \file conduction.c
00004  *  \brief Adds explicit thermal conduction term to the energy equation,
00005  *      dE/dt = Div(Q)
00006  *
00007  *   where 
00008  *    - Q = kappa_iso Grad(T) + kappa_aniso([b Dot Grad(T)]b) = heat flux
00009  *    -    T = (P/d)*(mbar/k_B) = temperature
00010  *    -    b = magnetic field unit vector
00011  *
00012  *   Here 
00013  *    - kappa_iso   is the coeffcient for   isotropic conduction
00014  *    -   kappa_aniso is the coeffcient for anisotropic conduction
00015  *
00016  *   The heat flux Q is calculated by calls to HeatFlux_* functions.
00017  *
00018  * CONTAINS PUBLIC FUNCTIONS:
00019  * - conduction() - updates energy equation with thermal conduction
00020  * - conduction_init() - allocates memory needed
00021  * - conduction_destruct() - frees memory used */
00022 /*============================================================================*/
00023 
00024 #include <math.h>
00025 #include <float.h>
00026 #include "../defs.h"
00027 #include "../athena.h"
00028 #include "../globals.h"
00029 #include "prototypes.h"
00030 #include "../prototypes.h"
00031 
00032 #ifdef THERMAL_CONDUCTION
00033 
00034 #ifdef BAROTROPIC
00035 #error : Thermal conduction requires an adiabatic EOS
00036 #endif
00037 
00038 /* Arrays for the temperature and heat fluxes */
00039 static Real ***Temp=NULL;
00040 static Real3Vect ***Q=NULL;
00041 
00042 /*==============================================================================
00043  * PRIVATE FUNCTION PROTOTYPES:
00044  *   HeatFlux_iso   - computes   isotropic heat flux
00045  *   HeatFlux_aniso - computes anisotropic heat flux
00046  *============================================================================*/
00047 
00048 void HeatFlux_iso(DomainS *pD);
00049 void HeatFlux_aniso(DomainS *pD);
00050 
00051 /*=========================== PUBLIC FUNCTIONS ===============================*/
00052 /*----------------------------------------------------------------------------*/
00053 /*! \fn void conduction(DomainS *pD)
00054  *  \brief Explicit thermal conduction
00055  */
00056 void conduction(DomainS *pD)
00057 {
00058   GridS *pG = (pD->Grid);
00059   int i, is = pG->is, ie = pG->ie;
00060   int j, jl, ju, js = pG->js, je = pG->je;
00061   int k, kl, ku, ks = pG->ks, ke = pG->ke;
00062   Real dtodx1=pG->dt/pG->dx1, dtodx2=0.0, dtodx3=0.0;
00063 
00064   if (pG->Nx[1] > 1){
00065     jl = js - 1;
00066     ju = je + 1;
00067     dtodx2 = pG->dt/pG->dx2;
00068   } else {
00069     jl = js;
00070     ju = je;
00071   }
00072   if (pG->Nx[2] > 1){
00073     kl = ks - 1;
00074     ku = ke + 1;
00075     dtodx3 = pG->dt/pG->dx3;
00076   } else {
00077     kl = ks;
00078     ku = ke;
00079   }
00080 
00081 /* Zero heat heat flux array; compute temperature at cell centers.  Temperature
00082  * is dimensionless (in units of [mbar/k_B]).  For cgs units, the value of
00083  * kappa_iso or kappa_aniso must contain a factor (be multiplied by) mbar/k_B.
00084  */
00085 
00086   for (k=kl; k<=ku; k++) {
00087   for (j=jl; j<=ju; j++) {
00088   for (i=is-1; i<=ie+1; i++) {
00089 
00090     Q[k][j][i].x = 0.0;
00091     Q[k][j][i].y = 0.0;
00092     Q[k][j][i].z = 0.0;
00093 
00094     Temp[k][j][i] = pG->U[k][j][i].E - (0.5/pG->U[k][j][i].d)*
00095       (SQR(pG->U[k][j][i].M1) +SQR(pG->U[k][j][i].M2) +SQR(pG->U[k][j][i].M3));
00096 #ifdef MHD
00097     Temp[k][j][i] -= (0.5)*(SQR(pG->U[k][j][i].B1c) +
00098       SQR(pG->U[k][j][i].B2c) + SQR(pG->U[k][j][i].B3c));
00099 #endif
00100     Temp[k][j][i] *= (Gamma_1/pG->U[k][j][i].d);
00101 
00102   }}}
00103 
00104 /* Compute isotropic and anisotropic heat fluxes.  Heat fluxes and temperature
00105  * are global variables in this file. */
00106 
00107   if (kappa_iso > 0.0)   HeatFlux_iso(pD);
00108   if (kappa_aniso > 0.0) HeatFlux_aniso(pD);
00109 
00110 /* Update energy using x1-fluxes */
00111 
00112   for (k=ks; k<=ke; k++) {
00113   for (j=js; j<=je; j++) {
00114     for (i=is; i<=ie; i++) {
00115       pG->U[k][j][i].E += dtodx1*(Q[k][j][i+1].x - Q[k][j][i].x);
00116     }
00117   }}
00118 
00119 /* Update energy using x2-fluxes */
00120 
00121   if (pG->Nx[1] > 1){
00122     for (k=ks; k<=ke; k++) {
00123     for (j=js; j<=je; j++) {
00124       for (i=is; i<=ie; i++) {
00125         pG->U[k][j][i].E += dtodx2*(Q[k][j+1][i].y - Q[k][j][i].y);
00126       }
00127     }}
00128   }
00129 
00130 /* Update energy using x3-fluxes */
00131 
00132   if (pG->Nx[2] > 1){
00133     for (k=ks; k<=ke; k++) {
00134     for (j=js; j<=je; j++) {
00135       for (i=is; i<=ie; i++) {
00136         pG->U[k][j][i].E += dtodx3*(Q[k+1][j][i].z - Q[k][j][i].z);
00137       }
00138     }}
00139   }
00140 
00141   return;
00142 }
00143 
00144 /*----------------------------------------------------------------------------*/
00145 /*! \fn void HeatFlux_iso(DomainS *pD)
00146  *  \brief Calculate heat fluxes with isotropic conduction
00147  */
00148 
00149 void HeatFlux_iso(DomainS *pD)
00150 { 
00151   GridS *pG = (pD->Grid);
00152   int i, is = pG->is, ie = pG->ie;
00153   int j, js = pG->js, je = pG->je;
00154   int k, ks = pG->ks, ke = pG->ke;
00155 
00156 /* Add heat fluxes in 1-direction */ 
00157 
00158   for (k=ks; k<=ke; k++) {
00159   for (j=js; j<=je; j++) {
00160     for (i=is; i<=ie+1; i++) {
00161       Q[k][j][i].x += kappa_iso*(Temp[k][j][i] - Temp[k][j][i-1])/pG->dx1;
00162     }
00163   }}
00164 
00165 /* Add heat fluxes in 2-direction */ 
00166 
00167   if (pG->Nx[1] > 1) {
00168     for (k=ks; k<=ke; k++) {
00169     for (j=js; j<=je+1; j++) {
00170       for (i=is; i<=ie; i++) {
00171         Q[k][j][i].y += kappa_iso*(Temp[k][j][i] - Temp[k][j-1][i])/pG->dx2;
00172       }
00173     }}
00174   }
00175 
00176 /* Add heat fluxes in 3-direction */
00177 
00178   if (pG->Nx[2] > 1) {
00179     for (k=ks; k<=ke+1; k++) {
00180     for (j=js; j<=je; j++) {
00181       for (i=is; i<=ie; i++) {
00182         Q[k][j][i].z += kappa_iso*(Temp[k][j][i] - Temp[k-1][j][i])/pG->dx3;
00183       }
00184     }}
00185   }
00186 
00187   return;
00188 }
00189 
00190 /*----------------------------------------------------------------------------*/
00191 /*! \fn void HeatFlux_aniso(DomainS *pD)
00192  *  \brief Calculate heat fluxes with anisotropic conduction
00193  */
00194 
00195 void HeatFlux_aniso(DomainS *pD)
00196 {
00197   GridS *pG = (pD->Grid);
00198   int i, is = pG->is, ie = pG->ie;
00199   int j, js = pG->js, je = pG->je;
00200   int k, ks = pG->ks, ke = pG->ke;
00201   Real Bx,By,Bz,B02,dTc,dTl,dTr,lim_slope,dTdx,dTdy,dTdz,bDotGradT;
00202 
00203 #ifdef MHD
00204   if (pD->Nx[1] == 1) return;  /* problem must be at least 2D */
00205 
00206 /* Compute heat fluxes in 1-direction  --------------------------------------*/
00207 
00208   for (k=ks; k<=ke; k++) {
00209   for (j=js; j<=je; j++) {
00210     for (i=is; i<=ie+1; i++) {
00211 
00212 /* Monotonized temperature difference dT/dy */
00213       dTr = 0.5*((Temp[k][j+1][i-1] + Temp[k][j+1][i]) -
00214                  (Temp[k][j  ][i-1] + Temp[k][j  ][i]));
00215       dTl = 0.5*((Temp[k][j  ][i-1] + Temp[k][j  ][i]) -
00216                  (Temp[k][j-1][i-1] + Temp[k][j-1][i]));
00217       dTc = dTr + dTl;
00218 
00219       dTdy = 0.0;
00220       if (dTl*dTr > 0.0) {
00221         lim_slope = MIN(fabs(dTl),fabs(dTr));
00222         dTdy = SIGN(dTc)*MIN(0.5*fabs(dTc),2.0*lim_slope)/pG->dx2;
00223       }
00224 
00225 /* Monotonized temperature difference dT/dz, 3D problem ONLY */
00226       if (pD->Nx[2] > 1) {
00227         dTr = 0.5*((Temp[k+1][j][i-1] + Temp[k+1][j][i]) -
00228                    (Temp[k  ][j][i-1] + Temp[k  ][j][i]));
00229         dTl = 0.5*((Temp[k  ][j][i-1] + Temp[k  ][j][i]) -
00230                    (Temp[k-1][j][i-1] + Temp[k-1][j][i]));
00231         dTc = dTr + dTl; 
00232 
00233         dTdz = 0.0;
00234         if (dTl*dTr > 0.0) {
00235           lim_slope = MIN(fabs(dTl),fabs(dTr));
00236           dTdz = SIGN(dTc)*MIN(0.5*fabs(dTc),2.0*lim_slope)/pG->dx3;
00237         }
00238       }
00239 
00240 /* Add flux at x1-interface, 2D PROBLEM */
00241 
00242       if (pD->Nx[2] == 1) {
00243         By = 0.5*(pG->U[k][j][i-1].B2c + pG->U[k][j][i].B2c);
00244         B02 = SQR(pG->B1i[k][j][i]) + SQR(By);
00245         B02 = MAX(B02,TINY_NUMBER); /* limit in case B=0 */
00246         bDotGradT = pG->B1i[k][j][i]*(Temp[k][j][i]-Temp[k][j][i-1])/pG->dx1
00247            + By*dTdy;
00248         Q[k][j][i].x += kappa_aniso*(pG->B1i[k][j][i]*bDotGradT)/B02;
00249 
00250 /* Add flux at x1-interface, 3D PROBLEM */
00251 
00252       } else {
00253         By = 0.5*(pG->U[k][j][i-1].B2c + pG->U[k][j][i].B2c);
00254         Bz = 0.5*(pG->U[k][j][i-1].B3c + pG->U[k][j][i].B3c);
00255         B02 = SQR(pG->B1i[k][j][i]) + SQR(By) + SQR(Bz);
00256         B02 = MAX(B02,TINY_NUMBER); /* limit in case B=0 */
00257         bDotGradT = pG->B1i[k][j][i]*(Temp[k][j][i]-Temp[k][j][i-1])/pG->dx1
00258            + By*dTdy + Bz*dTdz;
00259         Q[k][j][i].x += kappa_aniso*(pG->B1i[k][j][i]*bDotGradT)/B02;
00260       }
00261     }
00262   }}
00263 
00264 /* Compute heat fluxes in 2-direction  --------------------------------------*/
00265 
00266   for (k=ks; k<=ke; k++) {
00267   for (j=js; j<=je+1; j++) {
00268     for (i=is; i<=ie; i++) {
00269 
00270 /* Monotonized temperature difference dT/dx */
00271       dTr = 0.5*((Temp[k][j-1][i+1] + Temp[k][j][i+1]) -
00272                  (Temp[k][j-1][i  ] + Temp[k][j][i  ]));
00273       dTl = 0.5*((Temp[k][j-1][i  ] + Temp[k][j][i  ]) -
00274                  (Temp[k][j-1][i-1] + Temp[k][j][i-1]));
00275       dTc = dTr + dTl;
00276 
00277       dTdx = 0.0;
00278       if (dTl*dTr > 0.0) {
00279         lim_slope = MIN(fabs(dTl),fabs(dTr));
00280         dTdx = SIGN(dTc)*MIN(0.5*fabs(dTc),2.0*lim_slope)/pG->dx1;
00281       }
00282 
00283 /* Monotonized temperature difference dT/dz, 3D problem ONLY */
00284       if (pD->Nx[2] > 1) {
00285         dTr = 0.5*((Temp[k+1][j-1][i] + Temp[k+1][j][i]) -
00286                    (Temp[k  ][j-1][i] + Temp[k  ][j][i]));
00287         dTl = 0.5*((Temp[k  ][j-1][i] + Temp[k  ][j][i]) -
00288                    (Temp[k-1][j-1][i] + Temp[k-1][j][i]));
00289         dTc = dTr + dTl;
00290 
00291         dTdz = 0.0;
00292         if (dTl*dTr > 0.0) {
00293           lim_slope = MIN(fabs(dTl),fabs(dTr));
00294           dTdz = SIGN(dTc)*MIN(0.5*fabs(dTc),2.0*lim_slope)/pG->dx3;
00295         }
00296       }
00297 
00298 /* Add flux at x2-interface, 2D PROBLEM */
00299 
00300       if (pD->Nx[2] == 1) {
00301         Bx = 0.5*(pG->U[k][j-1][i].B1c + pG->U[k][j][i].B1c);
00302         B02 = SQR(Bx) + SQR(pG->B2i[k][j][i]);
00303         B02 = MAX(B02,TINY_NUMBER); /* limit in case B=0 */
00304 
00305         bDotGradT = pG->B2i[k][j][i]*(Temp[k][j][i]-Temp[k][j-1][i])/pG->dx2
00306            + Bx*dTdx;
00307         Q[k][j][i].y += kappa_aniso*(pG->B2i[k][j][i]*bDotGradT)/B02;
00308 
00309 /* Add flux at x2-interface, 3D PROBLEM */
00310 
00311       } else {
00312         Bx = 0.5*(pG->U[k][j-1][i].B1c + pG->U[k][j][i].B1c);
00313         Bz = 0.5*(pG->U[k][j-1][i].B3c + pG->U[k][j][i].B3c);
00314         B02 = SQR(Bx) + SQR(pG->B2i[k][j][i]) + SQR(Bz);
00315         B02 = MAX(B02,TINY_NUMBER); /* limit in case B=0 */
00316         bDotGradT = pG->B2i[k][j][i]*(Temp[k][j][i]-Temp[k][j-1][i])/pG->dx2
00317            + Bx*dTdx + Bz*dTdz;
00318         Q[k][j][i].y += kappa_aniso*(pG->B2i[k][j][i]*bDotGradT)/B02;
00319       }
00320     }
00321   }}
00322 
00323 /* Compute heat fluxes in 3-direction, 3D problem ONLY  ---------------------*/
00324 
00325   if (pD->Nx[2] > 1) {
00326     for (k=ks; k<=ke+1; k++) {
00327     for (j=js; j<=je; j++) {
00328       for (i=is; i<=ie; i++) {
00329 
00330 /* Monotonized temperature difference dT/dx */
00331         dTr = 0.5*((Temp[k-1][j][i+1] + Temp[k][j][i+1]) -
00332                    (Temp[k-1][j][i  ] + Temp[k][j][i  ]));
00333         dTl = 0.5*((Temp[k-1][j][i  ] + Temp[k][j][i  ]) -
00334                    (Temp[k-1][j][i-1] + Temp[k][j][i-1]));
00335         dTc = dTr + dTl;
00336 
00337         dTdx = 0.0; 
00338         if (dTl*dTr > 0.0) {
00339           lim_slope = MIN(fabs(dTl),fabs(dTr));
00340           dTdx = SIGN(dTc)*MIN(0.5*fabs(dTc),2.0*lim_slope)/pG->dx1;
00341         }
00342 
00343 /* Monotonized temperature difference dT/dy */
00344         dTr = 0.5*((Temp[k-1][j+1][i] + Temp[k][j+1][i]) -
00345                    (Temp[k-1][j  ][i] + Temp[k][j  ][i]));
00346         dTl = 0.5*((Temp[k-1][j  ][i] + Temp[k][j  ][i]) -
00347                    (Temp[k-1][j-1][i] + Temp[k][j-1][i]));
00348         dTc = dTr + dTl;
00349 
00350         dTdy = 0.0; 
00351         if (dTl*dTr > 0.0) {
00352           lim_slope = MIN(fabs(dTl),fabs(dTr));
00353           dTdy = SIGN(dTc)*MIN(0.5*fabs(dTc),2.0*lim_slope)/pG->dx2;
00354         }
00355 
00356 /* Add flux at x3-interface, 3D PROBLEM */
00357 
00358         Bx = 0.5*(pG->U[k-1][j][i].B1c + pG->U[k][j][i].B1c);
00359         By = 0.5*(pG->U[k-1][j][i].B2c + pG->U[k][j][i].B2c);
00360         B02 = SQR(Bx) + SQR(By) + SQR(pG->B3i[k][j][i]);
00361         B02 = MAX(B02,TINY_NUMBER); /* limit in case B=0 */
00362         bDotGradT = pG->B3i[k][j][i]*(Temp[k][j][i]-Temp[k-1][j][i])/pG->dx3
00363            + Bx*dTdx + By*dTdy;
00364         Q[k][j][i].z += kappa_aniso*(pG->B3i[k][j][i]*bDotGradT)/B02;
00365       }
00366     }}
00367   }
00368 #endif /* MHD */
00369 
00370   return;
00371 }
00372 
00373 
00374 /*----------------------------------------------------------------------------*/
00375 /*! \fn void conduction_init(MeshS *pM) 
00376  *  \brief Allocate temporary arrays
00377  */
00378 
00379 void conduction_init(MeshS *pM)
00380 {
00381   int nl,nd,size1=0,size2=0,size3=0,Nx1,Nx2,Nx3;
00382 
00383 /* Cycle over all Grids on this processor to find maximum Nx1, Nx2, Nx3 */
00384   for (nl=0; nl<(pM->NLevels); nl++){
00385     for (nd=0; nd<(pM->DomainsPerLevel[nl]); nd++){
00386       if (pM->Domain[nl][nd].Grid != NULL) {
00387         if (pM->Domain[nl][nd].Grid->Nx[0] > size1){
00388           size1 = pM->Domain[nl][nd].Grid->Nx[0];
00389         }
00390         if (pM->Domain[nl][nd].Grid->Nx[1] > size2){
00391           size2 = pM->Domain[nl][nd].Grid->Nx[1];
00392         }
00393         if (pM->Domain[nl][nd].Grid->Nx[2] > size3){
00394           size3 = pM->Domain[nl][nd].Grid->Nx[2];
00395         }
00396       }
00397     }
00398   }
00399 
00400   Nx1 = size1 + 2*nghost;
00401 
00402   if (pM->Nx[1] > 1){
00403     Nx2 = size2 + 2*nghost;
00404   } else {
00405     Nx2 = size2;
00406   }
00407 
00408   if (pM->Nx[2] > 1){
00409     Nx3 = size3 + 2*nghost;
00410   } else {
00411     Nx3 = size3;
00412   }
00413   if ((Temp = (Real***)calloc_3d_array(Nx3,Nx2,Nx1, sizeof(Real))) == NULL)
00414     goto on_error;
00415   if ((Q = (Real3Vect***)calloc_3d_array(Nx3,Nx2,Nx1,sizeof(Real3Vect)))==NULL)
00416     goto on_error;
00417   return;
00418 
00419   on_error:
00420   conduction_destruct();
00421   ath_error("[conduct_init]: malloc returned a NULL pointer\n");
00422 }
00423 
00424 /*----------------------------------------------------------------------------*/
00425 /*! \fn void conduction_destruct(void)
00426  *  \brief Free temporary arrays
00427  */
00428 
00429 void conduction_destruct(void)
00430 {
00431   if (Temp != NULL) free_3d_array(Temp);
00432   if (Q != NULL) free_3d_array(Q);
00433   return;
00434 }
00435 #endif /* THERMAL_CONDUCTION */