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 #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
00039 static Real ***Temp=NULL;
00040 static Real3Vect ***Q=NULL;
00041
00042
00043
00044
00045
00046
00047
00048 void HeatFlux_iso(DomainS *pD);
00049 void HeatFlux_aniso(DomainS *pD);
00050
00051
00052
00053
00054
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
00082
00083
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
00105
00106
00107 if (kappa_iso > 0.0) HeatFlux_iso(pD);
00108 if (kappa_aniso > 0.0) HeatFlux_aniso(pD);
00109
00110
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
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
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
00146
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
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
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
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
00192
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;
00205
00206
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
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
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
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);
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
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);
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
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
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
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
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);
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
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);
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
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
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
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
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);
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
00369
00370 return;
00371 }
00372
00373
00374
00375
00376
00377
00378
00379 void conduction_init(MeshS *pM)
00380 {
00381 int nl,nd,size1=0,size2=0,size3=0,Nx1,Nx2,Nx3;
00382
00383
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
00426
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