Main Page | Alphabetical List | Data Structures | File List | Data Fields | Globals | Related Pages

timestep.c

Go to the documentation of this file.
00001 #include <stdio.h>
00002 #include <stdlib.h>
00003 #include <string.h>
00004 #include <math.h>
00005 #include <mpi.h>
00006 #include "allvars.h"
00007 #include "proto.h"
00008 
00013 static double fac1, fac2, fac3, hubble_a, atime, a3inv;
00014 static double dt_displacement = 0;
00015 
00016 
00024 void advance_and_find_timesteps(void)
00025 {
00026   int i, j, no, ti_step, ti_min, tend, tstart;
00027   double dt_entr, dt_entr2, dt_gravkick, dt_hydrokick, dt_gravkick2, dt_hydrokick2, t0, t1;
00028   double minentropy, aphys;
00029   FLOAT dv[3];
00030 
00031 #if defined(PSEUDOSYMMETRIC) && !defined(FLEXSTEPS)
00032   double apred, prob;
00033   int ti_step2;
00034 #endif
00035 #ifdef PMGRID
00036   double dt_gravkickA, dt_gravkickB;
00037 #endif
00038 #ifdef MAKEGLASS
00039   double disp, dispmax, globmax, dmean, fac, disp2sum, globdisp2sum;
00040 #endif
00041 
00042   t0 = second();
00043 
00044   if(All.ComovingIntegrationOn)
00045     {
00046       fac1 = 1 / (All.Time * All.Time);
00047       fac2 = 1 / pow(All.Time, 3 * GAMMA - 2);
00048       fac3 = pow(All.Time, 3 * (1 - GAMMA) / 2.0);
00049       hubble_a = All.Omega0 / (All.Time * All.Time * All.Time)
00050         + (1 - All.Omega0 - All.OmegaLambda) / (All.Time * All.Time) + All.OmegaLambda;
00051 
00052       hubble_a = All.Hubble * sqrt(hubble_a);
00053       a3inv = 1 / (All.Time * All.Time * All.Time);
00054       atime = All.Time;
00055     }
00056   else
00057     fac1 = fac2 = fac3 = hubble_a = a3inv = atime = 1;
00058 
00059   if(Flag_FullStep || dt_displacement == 0)
00060     find_dt_displacement_constraint(hubble_a * atime * atime);
00061 
00062 #ifdef PMGRID
00063   if(All.ComovingIntegrationOn)
00064     dt_gravkickB = get_gravkick_factor(All.PM_Ti_begstep, All.Ti_Current) -
00065       get_gravkick_factor(All.PM_Ti_begstep, (All.PM_Ti_begstep + All.PM_Ti_endstep) / 2);
00066   else
00067     dt_gravkickB = (All.Ti_Current - (All.PM_Ti_begstep + All.PM_Ti_endstep) / 2) * All.Timebase_interval;
00068 
00069   if(All.PM_Ti_endstep == All.Ti_Current)       /* need to do long-range kick */
00070     {
00071       /* make sure that we reconstruct the domain/tree next time because we don't kick the tree nodes in this case */
00072       All.NumForcesSinceLastDomainDecomp = 1 + All.TotNumPart * All.TreeDomainUpdateFrequency;
00073     }
00074 #endif
00075 
00076 
00077 #ifdef MAKEGLASS
00078   for(i = 0, dispmax = 0, disp2sum = 0; i < NumPart; i++)
00079     {
00080       for(j = 0; j < 3; j++)
00081         {
00082           P[i].GravPM[j] *= -1;
00083           P[i].GravAccel[j] *= -1;
00084           P[i].GravAccel[j] += P[i].GravPM[j];
00085           P[i].GravPM[j] = 0;
00086         }
00087 
00088       disp = sqrt(P[i].GravAccel[0] * P[i].GravAccel[0] +
00089                   P[i].GravAccel[1] * P[i].GravAccel[1] + P[i].GravAccel[2] * P[i].GravAccel[2]);
00090 
00091       disp *= 2.0 / (3 * All.Hubble * All.Hubble);
00092 
00093       disp2sum += disp * disp;
00094 
00095       if(disp > dispmax)
00096         dispmax = disp;
00097     }
00098 
00099   MPI_Allreduce(&dispmax, &globmax, 1, MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD);
00100   MPI_Allreduce(&disp2sum, &globdisp2sum, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
00101 
00102   dmean = pow(P[0].Mass / (All.Omega0 * 3 * All.Hubble * All.Hubble / (8 * M_PI * All.G)), 1.0 / 3);
00103 
00104   if(globmax > dmean)
00105     fac = dmean / globmax;
00106   else
00107     fac = 1.0;
00108 
00109   if(ThisTask == 0)
00110     {
00111       printf("\nglass-making:  dmean= %g  global disp-maximum= %g  rms= %g\n\n",
00112              dmean, globmax, sqrt(globdisp2sum / All.TotNumPart));
00113       fflush(stdout);
00114     }
00115 
00116   for(i = 0, dispmax = 0; i < NumPart; i++)
00117     {
00118       for(j = 0; j < 3; j++)
00119         {
00120           P[i].Vel[j] = 0;
00121           P[i].Pos[j] += fac * P[i].GravAccel[j] * 2.0 / (3 * All.Hubble * All.Hubble);
00122           P[i].GravAccel[j] = 0;
00123         }
00124     }
00125 #endif
00126 
00127 
00128 
00129 
00130   /* Now assign new timesteps and kick */
00131 
00132 #ifdef FLEXSTEPS
00133   if(All.PresentMinStep < TIMEBASE)
00134     All.PresentMinStep *= 2;
00135 
00136   for(i = 0; i < NumPart; i++)
00137     {
00138       if(P[i].Ti_endstep == All.Ti_Current)
00139         {
00140           ti_step = get_timestep(i, &aphys, 0);
00141 
00142           /* make it a power 2 subdivision */
00143           ti_min = TIMEBASE;
00144           while(ti_min > ti_step)
00145             ti_min >>= 1;
00146           ti_step = ti_min;
00147 
00148           if(ti_step < All.PresentMinStep)
00149             All.PresentMinStep = ti_step;
00150         }
00151     }
00152 
00153   ti_step = All.PresentMinStep;
00154   MPI_Allreduce(&ti_step, &All.PresentMinStep, 1, MPI_INT, MPI_MIN, MPI_COMM_WORLD);
00155 #endif
00156 
00157 
00158   for(i = 0; i < NumPart; i++)
00159     {
00160       if(P[i].Ti_endstep == All.Ti_Current)
00161         {
00162           ti_step = get_timestep(i, &aphys, 0);
00163 
00164 #ifdef FLEXSTEPS
00165           ti_step = (ti_step / All.PresentMinStep) * All.PresentMinStep;
00166 
00167           if((((P[i].Ti_endstep + ti_step) / All.PresentMinStep) * All.PresentMinStep - P[i].Ti_endstep)
00168              >= All.PresentMinStep)
00169             ti_step =
00170               ((P[i].Ti_endstep + ti_step) / All.PresentMinStep) * All.PresentMinStep - P[i].Ti_endstep;
00171 #else
00172           /* make it a power 2 subdivision */
00173           ti_min = TIMEBASE;
00174           while(ti_min > ti_step)
00175             ti_min >>= 1;
00176           ti_step = ti_min;
00177 
00178 #ifdef PSEUDOSYMMETRIC
00179           if(P[i].Type != 0)
00180             {
00181               if(P[i].Ti_endstep > P[i].Ti_begstep)
00182                 {
00183                   apred = aphys + ((aphys - P[i].AphysOld) / (P[i].Ti_endstep - P[i].Ti_begstep)) * ti_step;
00184                   if(fabs(apred - aphys) < 0.5 * aphys)
00185                     {
00186                       ti_step2 = get_timestep(i, &apred, -1);
00187                       ti_min = TIMEBASE;
00188                       while(ti_min > ti_step2)
00189                         ti_min >>= 1;
00190                       ti_step2 = ti_min;
00191 
00192                       if(ti_step2 < ti_step)
00193                         {
00194                           get_timestep(i, &apred, ti_step);
00195                           prob =
00196                             ((apred - aphys) / (aphys - P[i].AphysOld) * (P[i].Ti_endstep -
00197                                                                           P[i].Ti_begstep)) / ti_step;
00198                           if(prob < get_random_number(P[i].ID))
00199                             ti_step /= 2;
00200                         }
00201                       else if(ti_step2 > ti_step)
00202                         {
00203                           get_timestep(i, &apred, 2 * ti_step);
00204                           prob =
00205                             ((apred - aphys) / (aphys - P[i].AphysOld) * (P[i].Ti_endstep -
00206                                                                           P[i].Ti_begstep)) / ti_step;
00207                           if(prob < get_random_number(P[i].ID + 1))
00208                             ti_step *= 2;
00209                         }
00210                     }
00211                 }
00212               P[i].AphysOld = aphys;
00213             }
00214 #endif
00215 
00216 #ifdef SYNCHRONIZATION
00217           if(ti_step > (P[i].Ti_endstep - P[i].Ti_begstep))     /* timestep wants to increase */
00218             {
00219               if(((TIMEBASE - P[i].Ti_endstep) % ti_step) > 0)
00220                 ti_step = P[i].Ti_endstep - P[i].Ti_begstep;    /* leave at old step */
00221             }
00222 #endif
00223 #endif /* end of FLEXSTEPS */
00224 
00225           if(All.Ti_Current == TIMEBASE)        /* we here finish the last timestep. */
00226             ti_step = 0;
00227 
00228           if((TIMEBASE - All.Ti_Current) < ti_step)     /* check that we don't run beyond the end */
00229             ti_step = TIMEBASE - All.Ti_Current;
00230 
00231           tstart = (P[i].Ti_begstep + P[i].Ti_endstep) / 2;     /* midpoint of old step */
00232           tend = P[i].Ti_endstep + ti_step / 2; /* midpoint of new step */
00233 
00234           if(All.ComovingIntegrationOn)
00235             {
00236               dt_entr = (tend - tstart) * All.Timebase_interval;
00237               dt_entr2 = (tend - P[i].Ti_endstep) * All.Timebase_interval;
00238               dt_gravkick = get_gravkick_factor(tstart, tend);
00239               dt_hydrokick = get_hydrokick_factor(tstart, tend);
00240               dt_gravkick2 = get_gravkick_factor(P[i].Ti_endstep, tend);
00241               dt_hydrokick2 = get_hydrokick_factor(P[i].Ti_endstep, tend);
00242             }
00243           else
00244             {
00245               dt_entr = dt_gravkick = dt_hydrokick = (tend - tstart) * All.Timebase_interval;
00246               dt_gravkick2 = dt_hydrokick2 = dt_entr2 = (tend - P[i].Ti_endstep) * All.Timebase_interval;
00247             }
00248 
00249           P[i].Ti_begstep = P[i].Ti_endstep;
00250           P[i].Ti_endstep = P[i].Ti_begstep + ti_step;
00251 
00252 
00253           /* do the kick */
00254 
00255           for(j = 0; j < 3; j++)
00256             {
00257               dv[j] = P[i].GravAccel[j] * dt_gravkick;
00258               P[i].Vel[j] += dv[j];
00259             }
00260 
00261           if(P[i].Type == 0)    /* SPH stuff */
00262             {
00263               for(j = 0; j < 3; j++)
00264                 {
00265                   dv[j] += SphP[i].HydroAccel[j] * dt_hydrokick;
00266                   P[i].Vel[j] += SphP[i].HydroAccel[j] * dt_hydrokick;
00267 
00268                   SphP[i].VelPred[j] =
00269                     P[i].Vel[j] - dt_gravkick2 * P[i].GravAccel[j] - dt_hydrokick2 * SphP[i].HydroAccel[j];
00270 #ifdef PMGRID
00271                   SphP[i].VelPred[j] += P[i].GravPM[j] * dt_gravkickB;
00272 #endif
00273                 }
00274 
00275               /* In case of cooling, we prevent that the entropy (and
00276                  hence temperature decreases by more than a factor 0.5 */
00277 
00278               if(SphP[i].DtEntropy * dt_entr > -0.5 * SphP[i].Entropy)
00279                 SphP[i].Entropy += SphP[i].DtEntropy * dt_entr;
00280               else
00281                 SphP[i].Entropy *= 0.5;
00282 
00283               if(All.MinEgySpec)
00284                 {
00285                   minentropy = All.MinEgySpec * GAMMA_MINUS1 / pow(SphP[i].Density * a3inv, GAMMA_MINUS1);
00286                   if(SphP[i].Entropy < minentropy)
00287                     {
00288                       SphP[i].Entropy = minentropy;
00289                       SphP[i].DtEntropy = 0;
00290                     }
00291                 }
00292 
00293               /* In case the timestep increases in the new step, we
00294                  make sure that we do not 'overcool' when deriving
00295                  predicted temperatures. The maximum timespan over
00296                  which prediction can occur is ti_step/2, i.e. from
00297                  the middle to the end of the current step */
00298 
00299               dt_entr = ti_step / 2 * All.Timebase_interval;
00300               if(SphP[i].Entropy + SphP[i].DtEntropy * dt_entr < 0.5 * SphP[i].Entropy)
00301                 SphP[i].DtEntropy = -0.5 * SphP[i].Entropy / dt_entr;
00302             }
00303 
00304 
00305           /* if tree is not going to be reconstructed, kick parent nodes dynamically.
00306            */
00307           if(All.NumForcesSinceLastDomainDecomp < All.TotNumPart * All.TreeDomainUpdateFrequency)
00308             {
00309               no = Father[i];
00310               while(no >= 0)
00311                 {
00312                   for(j = 0; j < 3; j++)
00313                     Extnodes[no].vs[j] += dv[j] * P[i].Mass / Nodes[no].u.d.mass;
00314 
00315                   no = Nodes[no].u.d.father;
00316                 }
00317             }
00318         }
00319     }
00320 
00321 
00322 
00323 #ifdef PMGRID
00324   if(All.PM_Ti_endstep == All.Ti_Current)       /* need to do long-range kick */
00325     {
00326       ti_step = TIMEBASE;
00327       while(ti_step > (dt_displacement / All.Timebase_interval))
00328         ti_step >>= 1;
00329 
00330       if(ti_step > (All.PM_Ti_endstep - All.PM_Ti_begstep))     /* PM-timestep wants to increase */
00331         {
00332           /* we only increase if an integer number of steps will bring us to the end */
00333           if(((TIMEBASE - All.PM_Ti_endstep) % ti_step) > 0)
00334             ti_step = All.PM_Ti_endstep - All.PM_Ti_begstep;    /* leave at old step */
00335         }
00336 
00337       if(All.Ti_Current == TIMEBASE)    /* we here finish the last timestep. */
00338         ti_step = 0;
00339 
00340       tstart = (All.PM_Ti_begstep + All.PM_Ti_endstep) / 2;
00341       tend = All.PM_Ti_endstep + ti_step / 2;
00342 
00343       if(All.ComovingIntegrationOn)
00344         dt_gravkick = get_gravkick_factor(tstart, tend);
00345       else
00346         dt_gravkick = (tend - tstart) * All.Timebase_interval;
00347 
00348       All.PM_Ti_begstep = All.PM_Ti_endstep;
00349       All.PM_Ti_endstep = All.PM_Ti_begstep + ti_step;
00350 
00351       if(All.ComovingIntegrationOn)
00352         dt_gravkickB = -get_gravkick_factor(All.PM_Ti_begstep, (All.PM_Ti_begstep + All.PM_Ti_endstep) / 2);
00353       else
00354         dt_gravkickB =
00355           -((All.PM_Ti_begstep + All.PM_Ti_endstep) / 2 - All.PM_Ti_begstep) * All.Timebase_interval;
00356 
00357       for(i = 0; i < NumPart; i++)
00358         {
00359           for(j = 0; j < 3; j++)        /* do the kick */
00360             P[i].Vel[j] += P[i].GravPM[j] * dt_gravkick;
00361 
00362           if(P[i].Type == 0)
00363             {
00364               if(All.ComovingIntegrationOn)
00365                 {
00366                   dt_gravkickA = get_gravkick_factor(P[i].Ti_begstep, All.Ti_Current) -
00367                     get_gravkick_factor(P[i].Ti_begstep, (P[i].Ti_begstep + P[i].Ti_endstep) / 2);
00368                   dt_hydrokick = get_hydrokick_factor(P[i].Ti_begstep, All.Ti_Current) -
00369                     get_hydrokick_factor(P[i].Ti_begstep, (P[i].Ti_begstep + P[i].Ti_endstep) / 2);
00370                 }
00371               else
00372                 dt_gravkickA = dt_hydrokick =
00373                   (All.Ti_Current - (P[i].Ti_begstep + P[i].Ti_endstep) / 2) * All.Timebase_interval;
00374 
00375               for(j = 0; j < 3; j++)
00376                 SphP[i].VelPred[j] = P[i].Vel[j]
00377                   + P[i].GravAccel[j] * dt_gravkickA
00378                   + SphP[i].HydroAccel[j] * dt_hydrokick + P[i].GravPM[j] * dt_gravkickB;
00379             }
00380         }
00381     }
00382 #endif
00383 
00384   t1 = second();
00385   All.CPU_TimeLine += timediff(t0, t1);
00386 }
00387 
00388 
00389 
00390 
00400 int get_timestep(int p,         
00401                  double *aphys, 
00402                  int flag        )
00404 {
00405   double ax, ay, az, ac, csnd;
00406   double dt = 0, dt_courant = 0, dt_accel;
00407   int ti_step;
00408 
00409 #ifdef CONDUCTION
00410   double dt_cond;
00411 #endif
00412 
00413   if(flag == 0)
00414     {
00415       ax = fac1 * P[p].GravAccel[0];
00416       ay = fac1 * P[p].GravAccel[1];
00417       az = fac1 * P[p].GravAccel[2];
00418 
00419 #ifdef PMGRID
00420       ax += fac1 * P[p].GravPM[0];
00421       ay += fac1 * P[p].GravPM[1];
00422       az += fac1 * P[p].GravPM[2];
00423 #endif
00424 
00425       if(P[p].Type == 0)
00426         {
00427           ax += fac2 * SphP[p].HydroAccel[0];
00428           ay += fac2 * SphP[p].HydroAccel[1];
00429           az += fac2 * SphP[p].HydroAccel[2];
00430         }
00431 
00432       ac = sqrt(ax * ax + ay * ay + az * az);   /* this is now the physical acceleration */
00433       *aphys = ac;
00434     }
00435   else
00436     ac = *aphys;
00437 
00438   if(ac == 0)
00439     ac = 1.0e-30;
00440 
00441   switch (All.TypeOfTimestepCriterion)
00442     {
00443     case 0:
00444       if(flag > 0)
00445         {
00446           dt = flag * All.Timebase_interval;
00447           dt /= hubble_a;       /* convert dloga to physical timestep  */
00448           ac = 2 * All.ErrTolIntAccuracy * atime * All.SofteningTable[P[p].Type] / (dt * dt);
00449           *aphys = ac;
00450           return flag;
00451         }
00452       dt = dt_accel = sqrt(2 * All.ErrTolIntAccuracy * atime * All.SofteningTable[P[p].Type] / ac);
00453 #ifdef ADAPTIVE_GRAVSOFT_FORGAS
00454       if(P[p].Type == 0)
00455         dt = dt_accel = sqrt(2 * All.ErrTolIntAccuracy * atime * SphP[p].Hsml / 2.8 / ac);
00456 #endif
00457       break;
00458     default:
00459       endrun(888);
00460       break;
00461     }
00462 
00463   if(P[p].Type == 0)
00464     {
00465       csnd = sqrt(GAMMA * SphP[p].Pressure / SphP[p].Density);
00466 
00467       if(All.ComovingIntegrationOn)
00468         dt_courant = 2 * All.CourantFac * All.Time * SphP[p].Hsml / (fac3 * SphP[p].MaxSignalVel);
00469       else
00470         dt_courant = 2 * All.CourantFac * SphP[p].Hsml / SphP[p].MaxSignalVel;
00471 
00472       if(dt_courant < dt)
00473         dt = dt_courant;
00474     }
00475 
00476   /* convert the physical timestep to dloga if needed. Note: If comoving integration has not been selected,
00477      hubble_a=1.
00478    */
00479   dt *= hubble_a;
00480 
00481   if(dt >= All.MaxSizeTimestep)
00482     dt = All.MaxSizeTimestep;
00483 
00484   if(dt >= dt_displacement)
00485     dt = dt_displacement;
00486 
00487   if(dt < All.MinSizeTimestep)
00488     {
00489 #ifndef NOSTOP_WHEN_BELOW_MINTIMESTEP
00490       printf("warning: Timestep wants to be below the limit `MinSizeTimestep'\n");
00491 
00492       if(P[p].Type == 0)
00493         {
00494           printf
00495             ("Part-ID=%d  dt=%g dtc=%g ac=%g xyz=(%g|%g|%g)  hsml=%g  maxsignalvel=%g dt0=%g eps=%g\n",
00496              (int) P[p].ID, dt, dt_courant * hubble_a, ac, P[p].Pos[0], P[p].Pos[1], P[p].Pos[2],
00497              SphP[p].Hsml, SphP[p].MaxSignalVel,
00498              sqrt(2 * All.ErrTolIntAccuracy * atime * All.SofteningTable[P[p].Type] / ac) * hubble_a,
00499              All.SofteningTable[P[p].Type]);
00500         }
00501       else
00502         {
00503           printf("Part-ID=%d  dt=%g ac=%g xyz=(%g|%g|%g)\n", (int) P[p].ID, dt, ac, P[p].Pos[0], P[p].Pos[1],
00504                  P[p].Pos[2]);
00505         }
00506       fflush(stdout);
00507       endrun(888);
00508 #endif
00509       dt = All.MinSizeTimestep;
00510     }
00511 
00512   ti_step = dt / All.Timebase_interval;
00513 
00514   if(!(ti_step > 0 && ti_step < TIMEBASE))
00515     {
00516       printf("\nError: A timestep of size zero was assigned on the integer timeline!\n"
00517              "We better stop.\n"
00518              "Task=%d Part-ID=%d dt=%g tibase=%g ti_step=%d ac=%g xyz=(%g|%g|%g) tree=(%g|%g%g)\n\n",
00519              ThisTask, (int) P[p].ID, dt, All.Timebase_interval, ti_step, ac,
00520              P[p].Pos[0], P[p].Pos[1], P[p].Pos[2], P[p].GravAccel[0], P[p].GravAccel[1], P[p].GravAccel[2]);
00521 #ifdef PMGRID
00522       printf("pm_force=(%g|%g|%g)\n", P[p].GravPM[0], P[p].GravPM[1], P[p].GravPM[2]);
00523 #endif
00524       if(P[p].Type == 0)
00525         printf("hydro-frc=(%g|%g|%g)\n", SphP[p].HydroAccel[0], SphP[p].HydroAccel[1], SphP[p].HydroAccel[2]);
00526 
00527       fflush(stdout);
00528       endrun(818);
00529     }
00530 
00531   return ti_step;
00532 }
00533 
00534 
00543 void find_dt_displacement_constraint(double hfac  )
00544 {
00545   int i, j, type, *temp;
00546   int count[6];
00547   long long count_sum[6];
00548   double v[6], v_sum[6], mim[6], min_mass[6];
00549   double dt, dmean, asmth = 0;
00550 
00551   dt_displacement = All.MaxSizeTimestep;
00552 
00553   if(All.ComovingIntegrationOn)
00554     {
00555       for(type = 0; type < 6; type++)
00556         {
00557           count[type] = 0;
00558           v[type] = 0;
00559           mim[type] = 1.0e30;
00560         }
00561 
00562       for(i = 0; i < NumPart; i++)
00563         {
00564           v[P[i].Type] += P[i].Vel[0] * P[i].Vel[0] + P[i].Vel[1] * P[i].Vel[1] + P[i].Vel[2] * P[i].Vel[2];
00565           if(mim[P[i].Type] > P[i].Mass)
00566             mim[P[i].Type] = P[i].Mass;
00567           count[P[i].Type]++;
00568         }
00569 
00570       MPI_Allreduce(v, v_sum, 6, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
00571       MPI_Allreduce(mim, min_mass, 6, MPI_DOUBLE, MPI_MIN, MPI_COMM_WORLD);
00572 
00573       temp = malloc(NTask * 6 * sizeof(int));
00574       MPI_Allgather(count, 6, MPI_INT, temp, 6, MPI_INT, MPI_COMM_WORLD);
00575       for(i = 0; i < 6; i++)
00576         {
00577           count_sum[i] = 0;
00578           for(j = 0; j < NTask; j++)
00579             count_sum[i] += temp[j * 6 + i];
00580         }
00581       free(temp);
00582 
00583       for(type = 0; type < 6; type++)
00584         {
00585           if(count_sum[type] > 0)
00586             {
00587               if(type == 0)
00588                 dmean =
00589                   pow(min_mass[type] / (All.OmegaBaryon * 3 * All.Hubble * All.Hubble / (8 * M_PI * All.G)),
00590                       1.0 / 3);
00591               else
00592                 dmean =
00593                   pow(min_mass[type] /
00594                       ((All.Omega0 - All.OmegaBaryon) * 3 * All.Hubble * All.Hubble / (8 * M_PI * All.G)),
00595                       1.0 / 3);
00596 
00597               dt = All.MaxRMSDisplacementFac * hfac * dmean / sqrt(v_sum[type] / count_sum[type]);
00598 
00599 #ifdef PMGRID
00600               asmth = All.Asmth[0];
00601 #ifdef PLACEHIGHRESREGION
00602               if(((1 << type) & (PLACEHIGHRESREGION)))
00603                 asmth = All.Asmth[1];
00604 #endif
00605               if(asmth < dmean)
00606                 dt = All.MaxRMSDisplacementFac * hfac * asmth / sqrt(v_sum[type] / count_sum[type]);
00607 #endif
00608 
00609               if(ThisTask == 0)
00610                 printf("type=%d  dmean=%g asmth=%g minmass=%g a=%g  sqrt(<p^2>)=%g  dlogmax=%g\n",
00611                        type, dmean, asmth, min_mass[type], All.Time, sqrt(v_sum[type] / count_sum[type]), dt);
00612 
00613               if(dt < dt_displacement)
00614                 dt_displacement = dt;
00615             }
00616         }
00617 
00618       if(ThisTask == 0)
00619         printf("displacement time constraint: %g  (%g)\n", dt_displacement, All.MaxSizeTimestep);
00620     }
00621 }

Generated on Sun May 22 17:33:29 2005 for GADGET-2 by  doxygen 1.3.9.1