00001 #include "copyright.h"
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021 #include <stdio.h>
00022 #include <math.h>
00023 #include "defs.h"
00024 #include "athena.h"
00025 #include "globals.h"
00026 #include "prototypes.h"
00027
00028
00029
00030
00031
00032 void new_dt(MeshS *pM)
00033 {
00034 GridS *pGrid;
00035 #ifndef SPECIAL_RELATIVITY
00036 int i,j,k;
00037 Real di,v1,v2,v3,qsq,asq,cf1sq,cf2sq,cf3sq;
00038 #ifdef ADIABATIC
00039 Real p;
00040 #endif
00041 #ifdef MHD
00042 Real b1,b2,b3,bsq,tsum,tdif;
00043 #endif
00044 #ifdef PARTICLES
00045 long q;
00046 #endif
00047 #endif
00048 #ifdef MPI_PARALLEL
00049 double dt, my_dt;
00050 int ierr;
00051 #endif
00052 int nl,nd;
00053 Real tlim,max_v1=0.0,max_v2=0.0,max_v3=0.0,max_dti = 0.0;
00054 Real x1,x2,x3;
00055
00056
00057
00058 for (nl=0; nl<(pM->NLevels); nl++){
00059 for (nd=0; nd<(pM->DomainsPerLevel[nl]); nd++){
00060
00061 if (pM->Domain[nl][nd].Grid != NULL) {
00062 pGrid=(pM->Domain[nl][nd].Grid);
00063
00064
00065 #ifdef SPECIAL_RELATIVITY
00066 max_v1 = max_v2 = max_v3 = 1.0;
00067 #else
00068
00069 for (k=pGrid->ks; k<=pGrid->ke; k++) {
00070 for (j=pGrid->js; j<=pGrid->je; j++) {
00071 for (i=pGrid->is; i<=pGrid->ie; i++) {
00072 di = 1.0/(pGrid->U[k][j][i].d);
00073 v1 = pGrid->U[k][j][i].M1*di;
00074 v2 = pGrid->U[k][j][i].M2*di;
00075 v3 = pGrid->U[k][j][i].M3*di;
00076 qsq = v1*v1 + v2*v2 + v3*v3;
00077
00078 #ifdef MHD
00079
00080
00081 b1 = pGrid->U[k][j][i].B1c
00082 + fabs((double)(pGrid->B1i[k][j][i] - pGrid->U[k][j][i].B1c));
00083 b2 = pGrid->U[k][j][i].B2c
00084 + fabs((double)(pGrid->B2i[k][j][i] - pGrid->U[k][j][i].B2c));
00085 b3 = pGrid->U[k][j][i].B3c
00086 + fabs((double)(pGrid->B3i[k][j][i] - pGrid->U[k][j][i].B3c));
00087 bsq = b1*b1 + b2*b2 + b3*b3;
00088
00089 #ifdef ADIABATIC
00090 p = MAX(Gamma_1*(pGrid->U[k][j][i].E - 0.5*pGrid->U[k][j][i].d*qsq
00091 - 0.5*bsq), TINY_NUMBER);
00092 asq = Gamma*p*di;
00093 #elif defined ISOTHERMAL
00094 asq = Iso_csound2;
00095 #endif
00096
00097 tsum = bsq*di + asq;
00098 tdif = bsq*di - asq;
00099 cf1sq = 0.5*(tsum + sqrt(tdif*tdif + 4.0*asq*(b2*b2+b3*b3)*di));
00100 cf2sq = 0.5*(tsum + sqrt(tdif*tdif + 4.0*asq*(b1*b1+b3*b3)*di));
00101 cf3sq = 0.5*(tsum + sqrt(tdif*tdif + 4.0*asq*(b1*b1+b2*b2)*di));
00102
00103 #else
00104
00105
00106 #ifdef ADIABATIC
00107 p = MAX(Gamma_1*(pGrid->U[k][j][i].E - 0.5*pGrid->U[k][j][i].d*qsq),
00108 TINY_NUMBER);
00109 asq = Gamma*p*di;
00110 #elif defined ISOTHERMAL
00111 asq = Iso_csound2;
00112 #endif
00113
00114 cf1sq = asq;
00115 cf2sq = asq;
00116 cf3sq = asq;
00117
00118 #endif
00119
00120
00121 if (pGrid->Nx[0] > 1)
00122 max_v1 = MAX(max_v1,fabs(v1)+sqrt((double)cf1sq));
00123 if (pGrid->Nx[1] > 1)
00124 #ifdef CYLINDRICAL
00125 cc_pos(pGrid,i,j,k,&x1,&x2,&x3);
00126 max_v2 = MAX(max_v2,(fabs(v2)+sqrt((double)cf2sq))/x1);
00127 #else
00128 max_v2 = MAX(max_v2,fabs(v2)+sqrt((double)cf2sq));
00129 #endif
00130 if (pGrid->Nx[2] > 1)
00131 max_v3 = MAX(max_v3,fabs(v3)+sqrt((double)cf3sq));
00132
00133 }
00134 }}
00135
00136 #endif
00137
00138
00139 #ifdef PARTICLES
00140 for (q=0; q<pGrid->nparticle; q++) {
00141 if (pGrid->Nx[0] > 1)
00142 max_v1 = MAX(max_v1, pGrid->particle[q].v1);
00143 if (pGrid->Nx[1] > 1)
00144 max_v2 = MAX(max_v2, pGrid->particle[q].v2);
00145 if (pGrid->Nx[2] > 1)
00146 max_v3 = MAX(max_v3, pGrid->particle[q].v3);
00147 }
00148 #endif
00149
00150
00151 if (pGrid->Nx[0] > 1)
00152 max_dti = MAX(max_dti, max_v1/pGrid->dx1);
00153 if (pGrid->Nx[1] > 1)
00154 max_dti = MAX(max_dti, max_v2/pGrid->dx2);
00155 if (pGrid->Nx[2] > 1)
00156 max_dti = MAX(max_dti, max_v3/pGrid->dx3);
00157
00158 }}}
00159
00160
00161
00162 if (pM->nstep == 0) {
00163 pM->dt = CourNo/max_dti;
00164 } else {
00165 pM->dt = MIN(2.0*pM->dt, CourNo/max_dti);
00166 }
00167
00168
00169
00170 #ifdef MPI_PARALLEL
00171 my_dt = pM->dt;
00172 ierr = MPI_Allreduce(&my_dt, &dt, 1, MPI_DOUBLE, MPI_MIN, MPI_COMM_WORLD);
00173 pM->dt = dt;
00174 #endif
00175
00176
00177
00178 tlim = par_getd("time","tlim");
00179 if ((pM->time < tlim) && ((tlim - pM->time) < pM->dt))
00180 pM->dt = tlim - pM->time;
00181
00182
00183
00184 for (nl=0; nl<=(pM->NLevels)-1; nl++){
00185 for (nd=0; nd<=(pM->DomainsPerLevel[nl])-1; nd++){
00186 if (pM->Domain[nl][nd].Grid != NULL) {
00187 pM->Domain[nl][nd].Grid->dt = pM->dt;
00188 }
00189 }
00190 }
00191
00192 return;
00193 }