Go to the documentation of this file.00001 #include "../copyright.h"
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018 #include <stdio.h>
00019 #include <math.h>
00020 #include "../defs.h"
00021 #include "../athena.h"
00022 #include "../globals.h"
00023 #include "prototypes.h"
00024 #include "../prototypes.h"
00025
00026
00027
00028
00029 Real diff_dt(MeshS *pM)
00030 {
00031 int irefine, ir;
00032 Real dtmin_diffusion=(HUGE_NUMBER);
00033 Real dxmin,qa,fac;
00034 #ifdef RESISTIVITY
00035 int i,j,k,nl,nd;
00036 int il,iu,jl,ju,kl,ku;
00037 Real qb;
00038 GridS *pG;
00039 #endif
00040 #ifdef MPI_PARALLEL
00041 double my_dt, dt;
00042 int ierr;
00043 #endif
00044
00045
00046
00047 irefine = 1;
00048 for (ir=1; ir<(pM->NLevels); ir++) irefine *= 2;
00049
00050 dxmin = pM->dx[0]/(Real)(irefine);
00051 if (pM->Nx[1] > 1) dxmin = MIN( dxmin, (pM->dx[1]/(Real)(irefine)) );
00052 if (pM->Nx[2] > 1) dxmin = MIN( dxmin, (pM->dx[2]/(Real)(irefine)) );
00053
00054 qa = CourNo*(dxmin*dxmin)/2.0;
00055
00056 fac = 1.0;
00057 if (pM->Nx[1] > 1) fac = 2.0;
00058 if (pM->Nx[2] > 1) fac = 3.0;
00059
00060 qa = qa / fac;
00061
00062 #ifdef THERMAL_CONDUCTION
00063 dtmin_diffusion = MIN(dtmin_diffusion,(qa/(kappa_iso + kappa_aniso)));
00064 #endif
00065 #ifdef VISCOSITY
00066 dtmin_diffusion = MIN(dtmin_diffusion,(qa/(nu_iso + nu_aniso)));
00067 #endif
00068
00069 #ifdef RESISTIVITY
00070 qb = 0.25*qa;
00071 for (nl=pM->NLevels-1; nl>=0; nl--){
00072 qb *= 4.0;
00073 for (nd=0; nd<(pM->DomainsPerLevel[nl]); nd++){
00074 if (pM->Domain[nl][nd].Grid != NULL){
00075
00076 pG = pM->Domain[nl][nd].Grid;
00077
00078 il = pG->is - 2; iu = pG->ie + 2;
00079 if (pG->Nx[1] > 1){
00080 jl = pG->js - 2; ju = pG->je + 2;
00081 } else {
00082 jl = pG->js; ju = pG->je;
00083 }
00084 if (pG->Nx[2] > 1){
00085 kl = pG->ks - 2; ku = pG->ke + 2;
00086 } else {
00087 kl = pG->ks; ku = pG->ke;
00088 }
00089
00090 for (k=kl; k<=ku; k++) {
00091 for (j=jl; j<=ju; j++) {
00092 for (i=il; i<=iu; i++) {
00093
00094 if (eta_Ohm > 0.0)
00095 dtmin_diffusion = MIN(dtmin_diffusion,(qb/pG->eta_Ohm[k][j][i]));
00096
00097 if (Q_Hall > 0.0)
00098 dtmin_diffusion = MIN(dtmin_diffusion,
00099 (0.5*fac*qb/pG->eta_Hall[k][j][i]));
00100
00101 if (Q_AD > 0.0)
00102 dtmin_diffusion = MIN(dtmin_diffusion, (qb/pG->eta_AD[k][j][i]));
00103 }}}
00104 }
00105 }
00106 }
00107 #endif
00108
00109
00110 #ifdef MPI_PARALLEL
00111 my_dt = dtmin_diffusion;
00112 ierr = MPI_Allreduce(&my_dt, &dt, 1, MPI_DOUBLE, MPI_MIN, MPI_COMM_WORLD);
00113 dtmin_diffusion = dt;
00114 #endif
00115
00116 return dtmin_diffusion;
00117 }