Commit 4e5ebffd authored by Volker Springel's avatar Volker Springel
Browse files

Merge branch 'master' of http://gitlab.mpcdf.mpg.de/vrs/gadget4

parents 7ee6f647 b613aaa9
...@@ -10,6 +10,7 @@ ...@@ -10,6 +10,7 @@
#PERIODIC # enables periodic boundary condistions #PERIODIC # enables periodic boundary condistions
#NTYPES=6 # number of particle types #NTYPES=6 # number of particle types
#RANDOMIZE_DOMAINCENTER # shifts the particle distribution randomly each step to reduce correlations of force errors in time #RANDOMIZE_DOMAINCENTER # shifts the particle distribution randomly each step to reduce correlations of force errors in time
#RANDOMIZE_DOMAINCENTER_TYPES # can be used to set a zoom region via a particle type which is then never placed across large node boundaries
#LEAN # selects a special 'lean' mode of code operation, which is optimized for aggressive memory saving #LEAN # selects a special 'lean' mode of code operation, which is optimized for aggressive memory saving
#LONG_X_BITS=2 # can be used to reduce periodic box-dimension in x-direction relative to nominal box size #LONG_X_BITS=2 # can be used to reduce periodic box-dimension in x-direction relative to nominal box size
#LONG_Y_BITS=2 # can be used to reduce periodic box-dimension in y-direction relative to nominal box size #LONG_Y_BITS=2 # can be used to reduce periodic box-dimension in y-direction relative to nominal box size
......
#defines for headers # defines for headers
VERSION_H VERSION_H
HDF5UTIL_H HDF5UTIL_H
ALLVARS_H ALLVARS_H
......
...@@ -324,6 +324,18 @@ recommended, and should have only positive effects. ...@@ -324,6 +324,18 @@ recommended, and should have only positive effects.
------- -------
**RANDOMIZE_DOMAINCENTER_TYPES** = 2
Can be set to select one or several types (this is a bitmask) which
will then be used to locate the extension of a certain region. When
the particle set is randomly translated throughout the box, the code
will then try to avoid intersecting large oct-tree node boundaries
with this region. When this option is not set explicitely but
PLACEHIGHRESREGION is active, then this is automatically done
with a default setting RANDOMIZE_DOMAINCENTER_TYPES=PLACEHIGHRESREGION.
-------
**EVALPOTENTIAL** **EVALPOTENTIAL**
When this is activated, the code also computes the gravitational When this is activated, the code also computes the gravitational
......
...@@ -147,6 +147,10 @@ ...@@ -147,6 +147,10 @@
#define HRPMGRID PMGRID #define HRPMGRID PMGRID
#endif #endif
#if !defined(RANDOMIZE_DOMAINCENTER_TYPES) && defined(PLACEHIGHRESREGION)
#define RANDOMIZE_DOMAINCENTER_TYPES PLACEHIGHRESREGION
#endif
#if defined(SUBFIND) && !defined(SELFGRAVITY) #if defined(SUBFIND) && !defined(SELFGRAVITY)
#error "Running SUBFIND without SELFGRAVITY enabled does not make sense." #error "Running SUBFIND without SELFGRAVITY enabled does not make sense."
#endif #endif
...@@ -221,6 +225,10 @@ ...@@ -221,6 +225,10 @@
#error "If PMGRID is used without PERIODIC, TREEPM_NOTIMESPLIT needs to be activated" #error "If PMGRID is used without PERIODIC, TREEPM_NOTIMESPLIT needs to be activated"
#endif #endif
#if defined(PMGRID) && defined(HIERARCHICAL_GRAVITY) && !defined(TREEPM_NOTIMESPLIT)
#error "If PMGRID is used together with HIERARCHICAL_GRAVITY, you also need to use TREEPM_NOTIMESPLIT"
#endif
#if defined(PLACEHIGHRESREGION) && !defined(RANDOMIZE_DOMAINCENTER) #if defined(PLACEHIGHRESREGION) && !defined(RANDOMIZE_DOMAINCENTER)
#error "PLACEHIGHRESREGION requires RANDOMIZE_DOMAINCENTER." #error "PLACEHIGHRESREGION requires RANDOMIZE_DOMAINCENTER."
#endif #endif
......
...@@ -153,6 +153,9 @@ class simparticles : public intposconvert, public setcomm ...@@ -153,6 +153,9 @@ class simparticles : public intposconvert, public setcomm
MyIntPosType Left[2][3]; MyIntPosType Left[2][3];
MyIntPosType OldMeshSize[2]; MyIntPosType OldMeshSize[2];
MyIntPosType ReferenceIntPos[2][3]; MyIntPosType ReferenceIntPos[2][3];
#endif
#if defined(RANDOMIZE_DOMAINCENTER_TYPES) || defined(PLACEHIGHRESREGION)
MyIntPosType PlacingMask; MyIntPosType PlacingMask;
MyIntPosType PlacingBlocksize; MyIntPosType PlacingBlocksize;
#endif #endif
...@@ -410,7 +413,6 @@ class simparticles : public intposconvert, public setcomm ...@@ -410,7 +413,6 @@ class simparticles : public intposconvert, public setcomm
int test_if_grav_timestep_is_too_large(int p, int bin); int test_if_grav_timestep_is_too_large(int p, int bin);
int get_timestep_bin(integertime ti_step); int get_timestep_bin(integertime ti_step);
private:
#ifdef ADAPTIVE_HYDRO_SOFTENING #ifdef ADAPTIVE_HYDRO_SOFTENING
int get_softeningtype_for_hydro_particle(int i) int get_softeningtype_for_hydro_particle(int i)
{ {
......
...@@ -255,7 +255,7 @@ class domain : public setcomm ...@@ -255,7 +255,7 @@ class domain : public setcomm
return a.targetindex < b.targetindex; return a.targetindex < b.targetindex;
} }
#if defined(PLACEHIGHRESREGION) || defined(RANDOMIZE_DOMAINCENTER) #if defined(RANDOMIZE_DOMAINCENTER_TYPES) || defined(RANDOMIZE_DOMAINCENTER)
MyIntPosType domainInnersize; MyIntPosType domainInnersize;
MyIntPosType domainReferenceIntPos[3]; MyIntPosType domainReferenceIntPos[3];
MySignedIntPosType domainXmintot[3], domainXmaxtot[3]; MySignedIntPosType domainXmintot[3], domainXmaxtot[3];
......
...@@ -119,13 +119,13 @@ void domain<simparticles>::do_box_wrapping(void) ...@@ -119,13 +119,13 @@ void domain<simparticles>::do_box_wrapping(void)
#ifdef RANDOMIZE_DOMAINCENTER #ifdef RANDOMIZE_DOMAINCENTER
/* determine new shift vector */ /* determine new shift vector */
#if defined(PLACEHIGHRESREGION) #if defined(RANDOMIZE_DOMAINCENTER_TYPES) || defined(PLACEHIGHRESREGION)
domain_find_type_extension(); domain_find_type_extension();
#endif #endif
if(ThisTask == 0) if(ThisTask == 0)
{ {
#if defined(PLACEHIGHRESREGION) #if defined(RANDOMIZE_DOMAINCENTER_TYPES) || defined(PLACEHIGHRESREGION)
int count = 0; int count = 0;
#endif #endif
...@@ -145,7 +145,7 @@ void domain<simparticles>::do_box_wrapping(void) ...@@ -145,7 +145,7 @@ void domain<simparticles>::do_box_wrapping(void)
Tp->CurrentShiftVector[j] += get_random_number() * pow(2.0, 32); Tp->CurrentShiftVector[j] += get_random_number() * pow(2.0, 32);
#endif #endif
#if defined(PLACEHIGHRESREGION) #if defined(RANDOMIZE_DOMAINCENTER_TYPES) || defined(PLACEHIGHRESREGION)
MyIntPosType boxoff = (Tp->CurrentShiftVector[j] & Tp->PlacingMask); MyIntPosType boxoff = (Tp->CurrentShiftVector[j] & Tp->PlacingMask);
MyIntPosType inboxoff = (Tp->CurrentShiftVector[j] - boxoff) % (Tp->PlacingBlocksize - domainInnersize); MyIntPosType inboxoff = (Tp->CurrentShiftVector[j] - boxoff) % (Tp->PlacingBlocksize - domainInnersize);
...@@ -190,7 +190,7 @@ void domain<simparticles>::do_box_wrapping(void) ...@@ -190,7 +190,7 @@ void domain<simparticles>::do_box_wrapping(void)
#endif #endif
} }
#if defined(PLACEHIGHRESREGION) #if defined(RANDOMIZE_DOMAINCENTER_TYPES) || defined(PLACEHIGHRESREGION)
template <typename partset> template <typename partset>
int domain<partset>::domain_type_extension_overlap(int j) int domain<partset>::domain_type_extension_overlap(int j)
...@@ -217,7 +217,7 @@ void domain<partset>::domain_find_type_extension(void) ...@@ -217,7 +217,7 @@ void domain<partset>::domain_find_type_extension(void)
for(int i = 0; i < Tp->NumPart; i++) for(int i = 0; i < Tp->NumPart; i++)
{ {
if(((1 << Tp->P[i].getType()) & (PLACEHIGHRESREGION))) if(((1 << Tp->P[i].getType()) & (RANDOMIZE_DOMAINCENTER_TYPES)))
{ {
for(int j = 0; j < 3; j++) for(int j = 0; j < 3; j++)
domainReferenceIntPos[j] = Tp->P[i].IntPos[j]; domainReferenceIntPos[j] = Tp->P[i].IntPos[j];
...@@ -232,8 +232,8 @@ void domain<partset>::domain_find_type_extension(void) ...@@ -232,8 +232,8 @@ void domain<partset>::domain_find_type_extension(void)
MPI_Allreduce(MPI_IN_PLACE, have_global, 1, MPI_2INT, MPI_MINLOC, Communicator); MPI_Allreduce(MPI_IN_PLACE, have_global, 1, MPI_2INT, MPI_MINLOC, Communicator);
if(have_global[0] >= NTask) if(have_global[0] >= NTask)
Terminate("have_global[0]=%d >= NTask=%d: Don't we have any particle? Note: PLACEHIGHRESREGION=%d is a bitmask", have_global[0], Terminate("have_global[0]=%d >= NTask=%d: Don't we have any particle? Note: RANDOMIZE_DOMAINCENTER_TYPES=%d is a bitmask",
NTask, PLACEHIGHRESREGION); have_global[0], NTask, RANDOMIZE_DOMAINCENTER_TYPES);
MPI_Bcast(domainReferenceIntPos, 3 * sizeof(MyIntPosType), MPI_BYTE, have_global[1], Communicator); MPI_Bcast(domainReferenceIntPos, 3 * sizeof(MyIntPosType), MPI_BYTE, have_global[1], Communicator);
...@@ -249,7 +249,7 @@ void domain<partset>::domain_find_type_extension(void) ...@@ -249,7 +249,7 @@ void domain<partset>::domain_find_type_extension(void)
for(int i = 0; i < Tp->NumPart; i++) for(int i = 0; i < Tp->NumPart; i++)
{ {
if(((1 << Tp->P[i].getType()) & (PLACEHIGHRESREGION))) if(((1 << Tp->P[i].getType()) & (RANDOMIZE_DOMAINCENTER_TYPES)))
{ {
MyIntPosType diff[3] = {Tp->P[i].IntPos[0] - domainReferenceIntPos[0], Tp->P[i].IntPos[1] - domainReferenceIntPos[1], MyIntPosType diff[3] = {Tp->P[i].IntPos[0] - domainReferenceIntPos[0], Tp->P[i].IntPos[1] - domainReferenceIntPos[1],
Tp->P[i].IntPos[2] - domainReferenceIntPos[2]}; Tp->P[i].IntPos[2] - domainReferenceIntPos[2]};
...@@ -280,10 +280,10 @@ void domain<partset>::domain_find_type_extension(void) ...@@ -280,10 +280,10 @@ void domain<partset>::domain_find_type_extension(void)
if((MyIntPosType)(domainXmaxtot[2] - domainXmintot[2]) > domainInnersize) if((MyIntPosType)(domainXmaxtot[2] - domainXmintot[2]) > domainInnersize)
domainInnersize = domainXmaxtot[2] - domainXmintot[2]; domainInnersize = domainXmaxtot[2] - domainXmintot[2];
domain_printf("DOMAIN: Shrink-wrap region size for PLACEHIGHRESREGION is %g\n", domainInnersize * Tp->FacIntToCoord); domain_printf("DOMAIN: Shrink-wrap region size for RANDOMIZE_DOMAINCENTER_TYPES is %g\n", domainInnersize * Tp->FacIntToCoord);
if(domainInnersize * Tp->FacIntToCoord >= 0.125 * All.BoxSize) if(domainInnersize * Tp->FacIntToCoord >= 0.125 * All.BoxSize)
Terminate("inappropriately big region selection for PLACEHIGHRESREGION"); Terminate("inappropriately big region selection for RANDOMIZE_DOMAINCENTER_TYPES");
/* increase the region by at least 1/8 of its size to still allow some randomness in placing the particles within the high-res node /* increase the region by at least 1/8 of its size to still allow some randomness in placing the particles within the high-res node
*/ */
......
...@@ -137,7 +137,7 @@ void snap_io::init_basic(simparticles *Sp_ptr) ...@@ -137,7 +137,7 @@ void snap_io::init_basic(simparticles *Sp_ptr)
#ifdef OUTPUT_PRESSURE #ifdef OUTPUT_PRESSURE
init_field("PRES", "Pressure", MEM_MY_FLOAT, FILE_MY_IO_FLOAT, SKIP_ON_READ, 1, A_NONE, 0, io_func_pressure, init_field("PRES", "Pressure", MEM_MY_FLOAT, FILE_MY_IO_FLOAT, SKIP_ON_READ, 1, A_NONE, 0, io_func_pressure,
GAS_ONLY, /* particle pressure */ GAS_ONLY, /* particle pressure */
1, -3., 2., -3., 1., 2., All.UnitDensity_in_cgs * All.UnitVelocity_in_cm_per_s * All.UnitVelocity_in_cm_per_s); 1, -3 * GAMMA, 2., -3., 1., 2., All.UnitDensity_in_cgs * All.UnitVelocity_in_cm_per_s * All.UnitVelocity_in_cm_per_s);
#endif #endif
#if defined(TIMEDEP_ART_VISC) && defined(OUTPUT_VISCOSITY_PARAMETER) #if defined(TIMEDEP_ART_VISC) && defined(OUTPUT_VISCOSITY_PARAMETER)
...@@ -415,13 +415,13 @@ void snap_io::read_ic(const char *fname) ...@@ -415,13 +415,13 @@ void snap_io::read_ic(const char *fname)
Sp->NumPart += add_numpart; Sp->NumPart += add_numpart;
#endif #endif
#ifdef GADGET2_HEADER #ifdef GADGET2_HEADER
#ifndef INITIAL_CONDITIONS_CONTAIN_ENTROPY #ifndef INITIAL_CONDITIONS_CONTAIN_ENTROPY
if(header.flag_entropy_instead_u) Terminate("Initial condition file contains entropy, but INITIAL_CONDITIONS_CONTAIN_ENTROPY is not set\n"); if(header.flag_entropy_instead_u)
#else Terminate("Initial condition file contains entropy, but INITIAL_CONDITIONS_CONTAIN_ENTROPY is not set\n");
if(! header.flag_entropy_instead_u)Terminate("Initial condition file contains uthermal, but INITIAL_CONDITIONS_CONTAIN_ENTROPY is set\n"); #else
if(!header.flag_entropy_instead_u)
Terminate("Initial condition file contains uthermal, but INITIAL_CONDITIONS_CONTAIN_ENTROPY is set\n");
#endif #endif
#endif #endif
...@@ -736,6 +736,11 @@ void snap_io::fill_file_header(int writeTask, int lastTask, long long *n_type, l ...@@ -736,6 +736,11 @@ void snap_io::fill_file_header(int writeTask, int lastTask, long long *n_type, l
header.Omega0 = All.Omega0; header.Omega0 = All.Omega0;
header.OmegaLambda = All.OmegaLambda; header.OmegaLambda = All.OmegaLambda;
#if !(defined(REARRANGE_OPTION) && defined(MERGERTREE))
header.HubbleParam = All.HubbleParam;
header.Hubble = All.Hubble;
#endif
#ifdef OUTPUT_IN_DOUBLEPRECISION #ifdef OUTPUT_IN_DOUBLEPRECISION
header.flag_doubleprecision = 1; header.flag_doubleprecision = 1;
#else #else
......
...@@ -79,10 +79,13 @@ class snap_io : public IO_Def ...@@ -79,10 +79,13 @@ class snap_io : public IO_Def
double BoxSize; /**< box-size of simulation in case periodic boundaries were used */ double BoxSize; /**< box-size of simulation in case periodic boundaries were used */
double Omega0; /**< matter density in units of critical density */ double Omega0; /**< matter density in units of critical density */
double OmegaLambda; /**< cosmological constant parameter */ double OmegaLambda; /**< cosmological constant parameter */
long long Ntrees; // this replaces the storage space for HubbleParam #if defined(REARRANGE_OPTION) && defined(MERGERTREE)
long long TotNtrees; // this replaces the storage space for Hubble long long Ntrees; // this replaces the storage space for HubbleParam
// double HubbleParam; /**< little 'h' to scale units systems */ long long TotNtrees; // this replaces the storage space for Hubble
// double Hubble; /**< Hubble constant in internal units */ #else
double HubbleParam; /**< little 'h' to scale units systems */
double Hubble; /**< Hubble constant in internal units */
#endif
unsigned int npartTotalHighWord[NTYPES_HEADER]; /**< High word of the total number of particles of each type */ unsigned int npartTotalHighWord[NTYPES_HEADER]; /**< High word of the total number of particles of each type */
int flag_entropy_instead_u; /**< flags that IC-file contains entropy instead of u */ int flag_entropy_instead_u; /**< flags that IC-file contains entropy instead of u */
int flag_doubleprecision; /**< flags that snapshot contains double-precision instead of single precision */ int flag_doubleprecision; /**< flags that snapshot contains double-precision instead of single precision */
......
...@@ -129,7 +129,6 @@ void sim::init(int RestartSnapNum) ...@@ -129,7 +129,6 @@ void sim::init(int RestartSnapNum)
Mem.myfree(tmp); Mem.myfree(tmp);
int count = 0; int count = 0;
for(int i = 0; i < Sp.NumPart; i++) for(int i = 0; i < Sp.NumPart; i++)
#ifdef SPLIT_PARTICLE_TYPE #ifdef SPLIT_PARTICLE_TYPE
...@@ -190,6 +189,11 @@ void sim::init(int RestartSnapNum) ...@@ -190,6 +189,11 @@ void sim::init(int RestartSnapNum)
Sp.P[j].setMass(Sp.P[j].getMass() * fac); Sp.P[j].setMass(Sp.P[j].getMass() * fac);
Sp.P[i].setMass(Sp.P[i].getMass() * (1 - fac)); Sp.P[i].setMass(Sp.P[i].getMass() * (1 - fac));
#ifdef SECOND_ORDER_LPT_ICS
Sp.P[j].OldAcc *= fac;
Sp.P[i].OldAcc *= (1 - fac);
#endif
Sp.P[j].setType(0); Sp.P[j].setType(0);
#if NSOFTCLASSES > 1 #if NSOFTCLASSES > 1
Sp.P[j].setSofteningClass(All.SofteningClassOfPartType[0]); Sp.P[j].setSofteningClass(All.SofteningClassOfPartType[0]);
...@@ -477,7 +481,7 @@ void sim::init(int RestartSnapNum) ...@@ -477,7 +481,7 @@ void sim::init(int RestartSnapNum)
*/ */
#ifdef PRESSURE_ENTROPY_SPH #ifdef PRESSURE_ENTROPY_SPH
#ifndef INITIAL_CONDITIONS_CONTAIN_ENTROPY #ifndef INITIAL_CONDITIONS_CONTAIN_ENTROPY
NgbTree.setup_entropy_to_invgamma(); NgbTree.setup_entropy_to_invgamma();
#endif #endif
#endif #endif
...@@ -485,15 +489,15 @@ void sim::init(int RestartSnapNum) ...@@ -485,15 +489,15 @@ void sim::init(int RestartSnapNum)
for(int i = 0; i < Sp.NumGas; i++) for(int i = 0; i < Sp.NumGas; i++)
{ {
#ifndef INITIAL_CONDITIONS_CONTAIN_ENTROPY #ifndef INITIAL_CONDITIONS_CONTAIN_ENTROPY
if(ThisTask == 0 && i == 0) if(ThisTask == 0 && i == 0)
printf("INIT: Converting u -> entropy\n"); printf("INIT: Converting u -> entropy\n");
#if !defined(PRESSURE_ENTROPY_SPH) && !defined(ISOTHERM_EQS) #if !defined(PRESSURE_ENTROPY_SPH) && !defined(ISOTHERM_EQS)
Sp.SphP[i].Entropy = GAMMA_MINUS1 * Sp.SphP[i].Entropy / pow(Sp.SphP[i].Density * All.cf_a3inv, GAMMA_MINUS1); Sp.SphP[i].Entropy = GAMMA_MINUS1 * Sp.SphP[i].Entropy / pow(Sp.SphP[i].Density * All.cf_a3inv, GAMMA_MINUS1);
#endif #endif
Sp.SphP[i].EntropyPred = Sp.SphP[i].Entropy; Sp.SphP[i].EntropyPred = Sp.SphP[i].Entropy;
#endif #endif
/* The predicted entropy values have been already set for all SPH formulation */ /* The predicted entropy values have been already set for all SPH formulation */
/* so it should be ok computing pressure and csound now */ /* so it should be ok computing pressure and csound now */
......
...@@ -172,11 +172,11 @@ void sim::run(void) ...@@ -172,11 +172,11 @@ void sim::run(void)
{ {
mpi_printf("\nFinal time=%g reached. Simulation ends.\n", All.TimeMax); mpi_printf("\nFinal time=%g reached. Simulation ends.\n", All.TimeMax);
if(All.Ti_lastoutput != All.Ti_Current) /* make a snapshot at the final time in case none has been produced at this time */ /* make a snapshot at the final time in case none has been produced at this time yet */
if(All.Ti_lastoutput != All.Ti_Current)
{ {
snap_io Snap(&Sp, Communicator, All.SnapFormat); /* get an I/O object */ All.Ti_nextoutput = All.Ti_Current;
/* this snapshot will be overwritten if All.TimeMax is increased and the run is continued */ create_snapshot_if_desired();
Snap.write_snapshot(All.SnapshotFileCount++, NORMAL_SNAPSHOT);
} }
break; break;
...@@ -521,6 +521,15 @@ void sim::create_snapshot_if_desired(void) ...@@ -521,6 +521,15 @@ void sim::create_snapshot_if_desired(void)
if(Sp.P[i].Ti_Current != All.Ti_Current) if(Sp.P[i].Ti_Current != All.Ti_Current)
Terminate("P[i].Ti_Current != All.Ti_Current"); Terminate("P[i].Ti_Current != All.Ti_Current");
#if defined(STARFORMATION) && defined(FOF)
// do an extra domain decomposition here to make sure that there are no new stars among the block of gas particles
NgbTree.treefree();
Domain.domain_free();
Domain.domain_decomposition(STANDARD);
NgbTree.treeallocate(Sp.NumGas, &Sp, &Domain);
NgbTree.treebuild(Sp.NumGas, NULL);
#endif
#ifndef OUTPUT_NON_SYNCHRONIZED_ALLOWED #ifndef OUTPUT_NON_SYNCHRONIZED_ALLOWED
NgbTree.treefree(); NgbTree.treefree();
Sp.TimeBinsGravity.timebins_free(); Sp.TimeBinsGravity.timebins_free();
......
...@@ -556,6 +556,9 @@ void sph::density(int *list, int ntarget) ...@@ -556,6 +556,9 @@ void sph::density(int *list, int ntarget)
All.Timebase_interval; All.Timebase_interval;
double dtime = All.cf_atime * dt / All.cf_atime_hubble_a; double dtime = All.cf_atime * dt / All.cf_atime_hubble_a;
SphP[target].set_viscosity_coefficient(dtime); SphP[target].set_viscosity_coefficient(dtime);
#endif
#ifdef ADAPTIVE_HYDRO_SOFTENING
Tp->P[target].setSofteningClass(Tp->get_softeningtype_for_hydro_particle(target));
#endif #endif
/* now check whether we had enough neighbours */ /* now check whether we had enough neighbours */
double desnumngb = All.DesNumNgb; double desnumngb = All.DesNumNgb;
...@@ -890,9 +893,9 @@ void sph::density_evaluate_kernel(pinfo &pdat) ...@@ -890,9 +893,9 @@ void sph::density_evaluate_kernel(pinfo &pdat)
} }
if(All.ComovingIntegrationOn) if(All.ComovingIntegrationOn)
vdotr2 += All.cf_atime2_hubble_a * r2; vdotr2 += All.cf_atime2_hubble_a * r2;
Vec4d mu_ij = vdotr2 /(All.cf_afac3 * All.Time * r); Vec4d mu_ij = vdotr2 / (All.cf_afac3 * All.cf_atime * r);
Vec4d cs_j(ngb0->Csnd, ngb1->Csnd, ngb2->Csnd, ngb3->Csnd); Vec4d cs_j(ngb0->Csnd, ngb1->Csnd, ngb2->Csnd, ngb3->Csnd);
Vec4d cs_sum = cs_i + cs_j; Vec4d cs_sum = cs_i + cs_j;
...@@ -903,10 +906,9 @@ void sph::density_evaluate_kernel(pinfo &pdat) ...@@ -903,10 +906,9 @@ void sph::density_evaluate_kernel(pinfo &pdat)
Vec4d decay_vel = select(decision, cs_sum, decay_vel_2); Vec4d decay_vel = select(decision, cs_sum, decay_vel_2);
Vec4d fac_decay_vel = select(decision_r_gt_0, 1, 0);
Vec4d fac_decay_vel = select(decision_r_gt_0, 1, 0); decay_vel = decay_vel * fac_decay_vel;
decay_vel = decay_vel * fac_decay_vel;
decayVel_loc = max(decayVel_loc, decay_vel); decayVel_loc = max(decayVel_loc, decay_vel);
...@@ -914,10 +916,11 @@ void sph::density_evaluate_kernel(pinfo &pdat) ...@@ -914,10 +916,11 @@ void sph::density_evaluate_kernel(pinfo &pdat)
} }
#ifdef TIMEDEP_ART_VISC #ifdef TIMEDEP_ART_VISC
for(int i = 0; i < vector_length; i++) { for(int i = 0; i < vector_length; i++)
if(decayVel_loc[i] > targetSphP->decayVel) {
targetSphP->decayVel = decayVel_loc[i]; if(decayVel_loc[i] > targetSphP->decayVel)
} targetSphP->decayVel = decayVel_loc[i];
}
#endif #endif
} }
...@@ -1025,9 +1028,9 @@ void sph::density_evaluate_kernel(pinfo &pdat) ...@@ -1025,9 +1028,9 @@ void sph::density_evaluate_kernel(pinfo &pdat)
} }
if(All.ComovingIntegrationOn) if(All.ComovingIntegrationOn)
vdotr2 += All.cf_atime2_hubble_a * r2; vdotr2 += All.cf_atime2_hubble_a * r2;
double mu_ij = vdotr2 / (All.cf_afac3 * All.Time * kernel.r); double mu_ij = vdotr2 / (All.cf_afac3 * All.cf_atime * kernel.r);
double decay_vel; double decay_vel;
if(vdotr2 < 0) if(vdotr2 < 0)
......
...@@ -151,7 +151,7 @@ class subdens_comm : public generic_comm<subdens_in, subdens_out, T_tree, T_doma ...@@ -151,7 +151,7 @@ class subdens_comm : public generic_comm<subdens_in, subdens_out, T_tree, T_doma
{ {
int p, type; int p, type;
double mass, r2; double mass, r2;
particle_data *P; typename T_partset::pdata *P;
if(no < Tree->MaxPart) /* single particle */ if(no < Tree->MaxPart) /* single particle */
{ {
......
...@@ -122,7 +122,7 @@ class nearest_comm : public generic_comm<nearest_in, nearest_out, T_tree, T_doma ...@@ -122,7 +122,7 @@ class nearest_comm : public generic_comm<nearest_in, nearest_out, T_tree, T_doma
{ {
if(no < Tree->MaxPart) /* single particle */ if(no < Tree->MaxPart) /* single particle */
{ {
particle_data *P = Tree->get_Pp(no, shmrank); auto *P = Tree->get_Pp(no, shmrank);
no = Tree->get_nextnodep(shmrank)[no]; /* note: here shmrank cannot change */ no = Tree->get_nextnodep(shmrank)[no]; /* note: here shmrank cannot change */
......
...@@ -189,7 +189,7 @@ class sngb_comm : public generic_comm<sngb_in, sngb_out, T_tree, T_domain, T_par ...@@ -189,7 +189,7 @@ class sngb_comm : public generic_comm<sngb_in, sngb_out, T_tree, T_domain, T_par
{ {
if(no < Tree->MaxPart) /* single particle */ if(no < Tree->MaxPart) /* single particle */
{ {
particle_data *P = Tree->get_Pp(no, shmrank); auto *P = Tree->get_Pp(no, shmrank);
subfind_data *PS = Tree->get_PSp(no, shmrank); subfind_data *PS = Tree->get_PSp(no, shmrank);
no = Tree->get_nextnodep(shmrank)[no]; /* note: here shmrank cannot change */ no = Tree->get_nextnodep(shmrank)[no]; /* note: here shmrank cannot change */
......
...@@ -557,6 +557,9 @@ void sim::hydro_force(int step_indicator) ...@@ -557,6 +557,9 @@ void sim::hydro_force(int step_indicator)
if(step_indicator == SECOND_HALF_STEP) if(step_indicator == SECOND_HALF_STEP)
{ {
Sp.SphP[target].EntropyPred = Sp.SphP[target].Entropy; Sp.SphP[target].EntropyPred = Sp.SphP[target].Entropy;
#ifdef PRESSURE_ENTROPY_SPH
Sp.SphP[target].EntropyToInvGammaPred = pow(Sp.SphP[target].EntropyPred, 1.0 / GAMMA);
#endif
Sp.SphP[target].set_thermodynamic_variables(); Sp.SphP[target].set_thermodynamic_variables();
Sp.SphP[target].VelPred[0] += (Sp.SphP[target].HydroAccel[0] - Old[i].HydroAccel[0]) * dt_hydrokick; Sp.SphP[target].VelPred[0] += (Sp.SphP[target].HydroAccel[0] - Old[i].HydroAccel[0]) * dt_hydrokick;
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment