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

init.c

Go to the documentation of this file.
00001 #include <stdlib.h>
00002 #include <string.h>
00003 #include <math.h>
00004 #include <mpi.h>
00005 
00006 #include "allvars.h"
00007 #include "proto.h"
00008 
00009 
00020 void init(void)
00021 {
00022   int i, j;
00023   double a3;
00024 
00025   All.Time = All.TimeBegin;
00026 
00027   switch (All.ICFormat)
00028     {
00029     case 1:
00030 #if (MAKEGLASS > 1)
00031       seed_glass();
00032 #else
00033       read_ic(All.InitCondFile);
00034 #endif
00035       break;
00036     case 2:
00037     case 3:
00038       read_ic(All.InitCondFile);
00039       break;
00040     default:
00041       if(ThisTask == 0)
00042         printf("ICFormat=%d not supported.\n", All.ICFormat);
00043       endrun(0);
00044     }
00045 
00046   All.Time = All.TimeBegin;
00047   All.Ti_Current = 0;
00048 
00049   if(All.ComovingIntegrationOn)
00050     {
00051       All.Timebase_interval = (log(All.TimeMax) - log(All.TimeBegin)) / TIMEBASE;
00052       a3 = All.Time * All.Time * All.Time;
00053     }
00054   else
00055     {
00056       All.Timebase_interval = (All.TimeMax - All.TimeBegin) / TIMEBASE;
00057       a3 = 1;
00058     }
00059 
00060   set_softenings();
00061 
00062   All.NumCurrentTiStep = 0;     /* setup some counters */
00063   All.SnapshotFileCount = 0;
00064   if(RestartFlag == 2)
00065     All.SnapshotFileCount = atoi(All.InitCondFile + strlen(All.InitCondFile) - 3) + 1;
00066 
00067   All.TotNumOfForces = 0;
00068   All.NumForcesSinceLastDomainDecomp = 0;
00069 
00070   if(All.ComovingIntegrationOn)
00071     if(All.PeriodicBoundariesOn == 1)
00072       check_omega();
00073 
00074   All.TimeLastStatistics = All.TimeBegin - All.TimeBetStatistics;
00075 
00076   if(All.ComovingIntegrationOn) /*  change to new velocity variable */
00077     {
00078       for(i = 0; i < NumPart; i++)
00079         for(j = 0; j < 3; j++)
00080           P[i].Vel[j] *= sqrt(All.Time) * All.Time;
00081     }
00082 
00083   for(i = 0; i < NumPart; i++)  /*  start-up initialization */
00084     {
00085       for(j = 0; j < 3; j++)
00086         P[i].GravAccel[j] = 0;
00087 #ifdef PMGRID
00088       for(j = 0; j < 3; j++)
00089         P[i].GravPM[j] = 0;
00090 #endif
00091       P[i].Ti_endstep = 0;
00092       P[i].Ti_begstep = 0;
00093 
00094       P[i].OldAcc = 0;
00095       P[i].GravCost = 1;
00096       P[i].Potential = 0;
00097     }
00098 
00099 #ifdef PMGRID
00100   All.PM_Ti_endstep = All.PM_Ti_begstep = 0;
00101 #endif
00102 
00103 #ifdef FLEXSTEPS
00104   All.PresentMinStep = TIMEBASE;
00105 #endif
00106 
00107   for(i = 0; i < N_gas; i++)    /* initialize sph_properties */
00108     {
00109       for(j = 0; j < 3; j++)
00110         {
00111           SphP[i].VelPred[j] = P[i].Vel[j];
00112           SphP[i].HydroAccel[j] = 0;
00113         }
00114 
00115       SphP[i].DtEntropy = 0;
00116 
00117       if(RestartFlag == 0)
00118         {
00119           SphP[i].Hsml = 0;
00120           SphP[i].Density = -1;
00121         }
00122     }
00123 
00124   ngb_treeallocate(MAX_NGB);
00125 
00126   force_treeallocate(All.TreeAllocFactor * All.MaxPart, All.MaxPart);
00127 
00128   All.NumForcesSinceLastDomainDecomp = 1 + All.TotNumPart * All.TreeDomainUpdateFrequency;
00129 
00130   Flag_FullStep = 1;            /* to ensure that Peano-Hilber order is done */
00131 
00132   domain_Decomposition();       /* do initial domain decomposition (gives equal numbers of particles) */
00133 
00134   ngb_treebuild();              /* will build tree */
00135 
00136   setup_smoothinglengths();
00137 
00138   TreeReconstructFlag = 1;
00139 
00140   /* at this point, the entropy variable normally contains the 
00141    * internal energy, read in from the initial conditions file, unless the file
00142    * explicitly signals that the initial conditions contain the entropy directly. 
00143    * Once the density has been computed, we can convert thermal energy to entropy.
00144    */
00145 #ifndef ISOTHERM_EQS
00146   if(header.flag_entropy_instead_u == 0)
00147     for(i = 0; i < N_gas; i++)
00148       SphP[i].Entropy = GAMMA_MINUS1 * SphP[i].Entropy / pow(SphP[i].Density / a3, GAMMA_MINUS1);
00149 #endif
00150 }
00151 
00152 
00156 void check_omega(void)
00157 {
00158   double mass = 0, masstot, omega;
00159   int i;
00160 
00161   for(i = 0; i < NumPart; i++)
00162     mass += P[i].Mass;
00163 
00164   MPI_Allreduce(&mass, &masstot, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
00165 
00166   omega =
00167     masstot / (All.BoxSize * All.BoxSize * All.BoxSize) / (3 * All.Hubble * All.Hubble / (8 * M_PI * All.G));
00168 
00169   if(fabs(omega - All.Omega0) > 1.0e-3)
00170     {
00171       if(ThisTask == 0)
00172         {
00173           printf("\n\nI've found something odd!\n");
00174           printf
00175             ("The mass content accounts only for Omega=%g,\nbut you specified Omega=%g in the parameterfile.\n",
00176              omega, All.Omega0);
00177           printf("\nI better stop.\n");
00178 
00179           fflush(stdout);
00180         }
00181       endrun(1);
00182     }
00183 }
00184 
00185 
00186 
00193 void setup_smoothinglengths(void)
00194 {
00195   int i, no, p;
00196 
00197   if(RestartFlag == 0)
00198     {
00199 
00200       for(i = 0; i < N_gas; i++)
00201         {
00202           no = Father[i];
00203 
00204           while(10 * All.DesNumNgb * P[i].Mass > Nodes[no].u.d.mass)
00205             {
00206               p = Nodes[no].u.d.father;
00207 
00208               if(p < 0)
00209                 break;
00210 
00211               no = p;
00212             }
00213 #ifndef TWODIMS
00214           SphP[i].Hsml =
00215             pow(3.0 / (4 * M_PI) * All.DesNumNgb * P[i].Mass / Nodes[no].u.d.mass, 1.0 / 3) * Nodes[no].len;
00216 #else
00217           SphP[i].Hsml =
00218             pow(1.0 / (M_PI) * All.DesNumNgb * P[i].Mass / Nodes[no].u.d.mass, 1.0 / 2) * Nodes[no].len;
00219 #endif
00220         }
00221     }
00222 
00223   density();
00224 }
00225 
00226 
00230 #if (MAKEGLASS > 1)
00231 void seed_glass(void)
00232 {
00233   int i, k, n_for_this_task;
00234   double Range[3], LowerBound[3];
00235   double drandom, partmass;
00236   long long IDstart;
00237 
00238   All.TotNumPart = MAKEGLASS;
00239   partmass = All.Omega0 * (3 * All.Hubble * All.Hubble / (8 * M_PI * All.G))
00240     * (All.BoxSize * All.BoxSize * All.BoxSize) / All.TotNumPart;
00241 
00242   All.MaxPart = All.PartAllocFactor * (All.TotNumPart / NTask); /* sets the maximum number of particles that may */
00243 
00244   allocate_memory();
00245 
00246   header.npartTotal[1] = All.TotNumPart;
00247   header.mass[1] = partmass;
00248 
00249   if(ThisTask == 0)
00250     {
00251       printf("\nGlass initialising\nPartMass= %g\n", partmass);
00252       printf("TotNumPart= %d%09d\n\n",
00253              (int) (All.TotNumPart / 1000000000), (int) (All.TotNumPart % 1000000000));
00254     }
00255 
00256   /* set the number of particles assigned locally to this task */
00257   n_for_this_task = All.TotNumPart / NTask;
00258 
00259   if(ThisTask == NTask - 1)
00260     n_for_this_task = All.TotNumPart - (NTask - 1) * n_for_this_task;
00261 
00262   NumPart = 0;
00263   IDstart = 1 + (All.TotNumPart / NTask) * ThisTask;
00264 
00265   /* split the temporal domain into Ntask slabs in z-direction */
00266 
00267   Range[0] = Range[1] = All.BoxSize;
00268   Range[2] = All.BoxSize / NTask;
00269   LowerBound[0] = LowerBound[1] = 0;
00270   LowerBound[2] = ThisTask * Range[2];
00271 
00272   srand48(ThisTask);
00273 
00274   for(i = 0; i < n_for_this_task; i++)
00275     {
00276       for(k = 0; k < 3; k++)
00277         {
00278           drandom = drand48();
00279 
00280           P[i].Pos[k] = LowerBound[k] + Range[k] * drandom;
00281           P[i].Vel[k] = 0;
00282         }
00283 
00284       P[i].Mass = partmass;
00285       P[i].Type = 1;
00286       P[i].ID = IDstart + i;
00287 
00288       NumPart++;
00289     }
00290 }
00291 #endif

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