diff --git a/Template-Config.sh b/Template-Config.sh
index 8091e3cb2030930ca59df77a5a042e85c87a5809..49e155bcb3e683dff0a2d8a96a7921f388d41778 100644
--- a/Template-Config.sh
+++ b/Template-Config.sh
@@ -10,6 +10,7 @@
#PERIODIC # enables periodic boundary condistions
#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_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
#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
diff --git a/defines_extra b/defines_extra
index fb52addc8147aacd4380ab1fcd489de9f6c8df83..efebea56f7b44cb52644382096e4b1d3bc8d45a4 100644
--- a/defines_extra
+++ b/defines_extra
@@ -1,4 +1,4 @@
-#defines for headers
+# defines for headers
VERSION_H
HDF5UTIL_H
ALLVARS_H
diff --git a/documentation/04_config-options.md b/documentation/04_config-options.md
index cf80f934b9f357ae0226c9d69861ae91bd4449b1..0c9bd761861e3dc1d798c72ce93ee38e385515a6 100644
--- a/documentation/04_config-options.md
+++ b/documentation/04_config-options.md
@@ -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**
When this is activated, the code also computes the gravitational
diff --git a/src/data/constants.h b/src/data/constants.h
index 2fc5303aa9552c88ea297548995f86e350e69262..a66aea27b283b0722eb3d8e221886ea174023cb8 100644
--- a/src/data/constants.h
+++ b/src/data/constants.h
@@ -147,6 +147,10 @@
#define HRPMGRID PMGRID
#endif
+#if !defined(RANDOMIZE_DOMAINCENTER_TYPES) && defined(PLACEHIGHRESREGION)
+#define RANDOMIZE_DOMAINCENTER_TYPES PLACEHIGHRESREGION
+#endif
+
#if defined(SUBFIND) && !defined(SELFGRAVITY)
#error "Running SUBFIND without SELFGRAVITY enabled does not make sense."
#endif
@@ -221,6 +225,10 @@
#error "If PMGRID is used without PERIODIC, TREEPM_NOTIMESPLIT needs to be activated"
#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)
#error "PLACEHIGHRESREGION requires RANDOMIZE_DOMAINCENTER."
#endif
diff --git a/src/data/simparticles.h b/src/data/simparticles.h
index 6b4f2279bdbfa4b79baff00af79a2f985d21ce68..dbf89705d5b5c2b33b0e387c63c779ccb7e09a65 100644
--- a/src/data/simparticles.h
+++ b/src/data/simparticles.h
@@ -153,6 +153,9 @@ class simparticles : public intposconvert, public setcomm
MyIntPosType Left[2][3];
MyIntPosType OldMeshSize[2];
MyIntPosType ReferenceIntPos[2][3];
+#endif
+
+#if defined(RANDOMIZE_DOMAINCENTER_TYPES) || defined(PLACEHIGHRESREGION)
MyIntPosType PlacingMask;
MyIntPosType PlacingBlocksize;
#endif
@@ -410,7 +413,6 @@ class simparticles : public intposconvert, public setcomm
int test_if_grav_timestep_is_too_large(int p, int bin);
int get_timestep_bin(integertime ti_step);
- private:
#ifdef ADAPTIVE_HYDRO_SOFTENING
int get_softeningtype_for_hydro_particle(int i)
{
diff --git a/src/domain/domain.h b/src/domain/domain.h
index 58d42041cad5112fb46031cd9dbf4b432bce7279..1ac2308a3876d37aa44baa876d82a6a7f7e2dda1 100644
--- a/src/domain/domain.h
+++ b/src/domain/domain.h
@@ -255,7 +255,7 @@ class domain : public setcomm
return a.targetindex < b.targetindex;
}
-#if defined(PLACEHIGHRESREGION) || defined(RANDOMIZE_DOMAINCENTER)
+#if defined(RANDOMIZE_DOMAINCENTER_TYPES) || defined(RANDOMIZE_DOMAINCENTER)
MyIntPosType domainInnersize;
MyIntPosType domainReferenceIntPos[3];
MySignedIntPosType domainXmintot[3], domainXmaxtot[3];
diff --git a/src/domain/domain_box.cc b/src/domain/domain_box.cc
index 2c1520c11647ad6b52bab7681ebf355a1d7b7ffc..10bc0b008ceb64de12be2089d0deb7f93de79af4 100644
--- a/src/domain/domain_box.cc
+++ b/src/domain/domain_box.cc
@@ -119,13 +119,13 @@ void domain<simparticles>::do_box_wrapping(void)
#ifdef RANDOMIZE_DOMAINCENTER
/* determine new shift vector */
-#if defined(PLACEHIGHRESREGION)
+#if defined(RANDOMIZE_DOMAINCENTER_TYPES) || defined(PLACEHIGHRESREGION)
domain_find_type_extension();
#endif
if(ThisTask == 0)
{
-#if defined(PLACEHIGHRESREGION)
+#if defined(RANDOMIZE_DOMAINCENTER_TYPES) || defined(PLACEHIGHRESREGION)
int count = 0;
#endif
@@ -145,7 +145,7 @@ void domain<simparticles>::do_box_wrapping(void)
Tp->CurrentShiftVector[j] += get_random_number() * pow(2.0, 32);
#endif
-#if defined(PLACEHIGHRESREGION)
+#if defined(RANDOMIZE_DOMAINCENTER_TYPES) || defined(PLACEHIGHRESREGION)
MyIntPosType boxoff = (Tp->CurrentShiftVector[j] & Tp->PlacingMask);
MyIntPosType inboxoff = (Tp->CurrentShiftVector[j] - boxoff) % (Tp->PlacingBlocksize - domainInnersize);
@@ -190,7 +190,7 @@ void domain<simparticles>::do_box_wrapping(void)
#endif
}
-#if defined(PLACEHIGHRESREGION)
+#if defined(RANDOMIZE_DOMAINCENTER_TYPES) || defined(PLACEHIGHRESREGION)
template <typename partset>
int domain<partset>::domain_type_extension_overlap(int j)
@@ -217,7 +217,7 @@ void domain<partset>::domain_find_type_extension(void)
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++)
domainReferenceIntPos[j] = Tp->P[i].IntPos[j];
@@ -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);
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],
- NTask, PLACEHIGHRESREGION);
+ Terminate("have_global[0]=%d >= NTask=%d: Don't we have any particle? Note: RANDOMIZE_DOMAINCENTER_TYPES=%d is a bitmask",
+ have_global[0], NTask, RANDOMIZE_DOMAINCENTER_TYPES);
MPI_Bcast(domainReferenceIntPos, 3 * sizeof(MyIntPosType), MPI_BYTE, have_global[1], Communicator);
@@ -249,7 +249,7 @@ void domain<partset>::domain_find_type_extension(void)
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],
Tp->P[i].IntPos[2] - domainReferenceIntPos[2]};
@@ -280,10 +280,10 @@ void domain<partset>::domain_find_type_extension(void)
if((MyIntPosType)(domainXmaxtot[2] - domainXmintot[2]) > domainInnersize)
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)
- 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
*/
diff --git a/src/io/snap_io.cc b/src/io/snap_io.cc
index be030887546353e72acaab7d408985b39c7f4de7..81af319f29c43aa999feb780e3fb5c193bd2d9cc 100644
--- a/src/io/snap_io.cc
+++ b/src/io/snap_io.cc
@@ -137,7 +137,7 @@ void snap_io::init_basic(simparticles *Sp_ptr)
#ifdef OUTPUT_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 */
- 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
#if defined(TIMEDEP_ART_VISC) && defined(OUTPUT_VISCOSITY_PARAMETER)
@@ -415,13 +415,13 @@ void snap_io::read_ic(const char *fname)
Sp->NumPart += add_numpart;
#endif
-
-
#ifdef GADGET2_HEADER
#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");
-#else
- if(! header.flag_entropy_instead_u)Terminate("Initial condition file contains uthermal, but INITIAL_CONDITIONS_CONTAIN_ENTROPY is set\n");
+ if(header.flag_entropy_instead_u)
+ Terminate("Initial condition file contains entropy, but INITIAL_CONDITIONS_CONTAIN_ENTROPY is not 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
@@ -736,6 +736,11 @@ void snap_io::fill_file_header(int writeTask, int lastTask, long long *n_type, l
header.Omega0 = All.Omega0;
header.OmegaLambda = All.OmegaLambda;
+#if !(defined(REARRANGE_OPTION) && defined(MERGERTREE))
+ header.HubbleParam = All.HubbleParam;
+ header.Hubble = All.Hubble;
+#endif
+
#ifdef OUTPUT_IN_DOUBLEPRECISION
header.flag_doubleprecision = 1;
#else
diff --git a/src/io/snap_io.h b/src/io/snap_io.h
index 371474a2e6b710e7ed0dd57253249ea352752f89..ca9014108aee31d73fd90c2fd6439991f864e562 100644
--- a/src/io/snap_io.h
+++ b/src/io/snap_io.h
@@ -79,10 +79,13 @@ class snap_io : public IO_Def
double BoxSize; /**< box-size of simulation in case periodic boundaries were used */
double Omega0; /**< matter density in units of critical density */
double OmegaLambda; /**< cosmological constant parameter */
- long long Ntrees; // this replaces the storage space for HubbleParam
- long long TotNtrees; // this replaces the storage space for Hubble
- // double HubbleParam; /**< little 'h' to scale units systems */
- // double Hubble; /**< Hubble constant in internal units */
+#if defined(REARRANGE_OPTION) && defined(MERGERTREE)
+ long long Ntrees; // this replaces the storage space for HubbleParam
+ long long TotNtrees; // this replaces the storage space for Hubble
+#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 */
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 */
diff --git a/src/main/init.cc b/src/main/init.cc
index 198f9fef983e9733052e6b210de04adca52ddeec..0a026f660b92e6208d2462876d38e1f3522e835e 100644
--- a/src/main/init.cc
+++ b/src/main/init.cc
@@ -129,7 +129,6 @@ void sim::init(int RestartSnapNum)
Mem.myfree(tmp);
-
int count = 0;
for(int i = 0; i < Sp.NumPart; i++)
#ifdef SPLIT_PARTICLE_TYPE
@@ -190,6 +189,11 @@ void sim::init(int RestartSnapNum)
Sp.P[j].setMass(Sp.P[j].getMass() * 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);
#if NSOFTCLASSES > 1
Sp.P[j].setSofteningClass(All.SofteningClassOfPartType[0]);
@@ -477,7 +481,7 @@ void sim::init(int RestartSnapNum)
*/
#ifdef PRESSURE_ENTROPY_SPH
#ifndef INITIAL_CONDITIONS_CONTAIN_ENTROPY
- NgbTree.setup_entropy_to_invgamma();
+ NgbTree.setup_entropy_to_invgamma();
#endif
#endif
@@ -485,15 +489,15 @@ void sim::init(int RestartSnapNum)
for(int i = 0; i < Sp.NumGas; i++)
{
#ifndef INITIAL_CONDITIONS_CONTAIN_ENTROPY
-
- if(ThisTask == 0 && i == 0)
- printf("INIT: Converting u -> entropy\n");
+
+ if(ThisTask == 0 && i == 0)
+ printf("INIT: Converting u -> entropy\n");
#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
- Sp.SphP[i].EntropyPred = Sp.SphP[i].Entropy;
-
+ Sp.SphP[i].EntropyPred = Sp.SphP[i].Entropy;
+
#endif
/* The predicted entropy values have been already set for all SPH formulation */
/* so it should be ok computing pressure and csound now */
diff --git a/src/main/run.cc b/src/main/run.cc
index ee509728e2096f56e203d0fd69297da8f42f3cf0..70ebf393936167f5ccfa60080619cc97c5cdc4a6 100644
--- a/src/main/run.cc
+++ b/src/main/run.cc
@@ -172,11 +172,11 @@ void sim::run(void)
{
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 */
- /* this snapshot will be overwritten if All.TimeMax is increased and the run is continued */
- Snap.write_snapshot(All.SnapshotFileCount++, NORMAL_SNAPSHOT);
+ All.Ti_nextoutput = All.Ti_Current;
+ create_snapshot_if_desired();
}
break;
@@ -521,6 +521,15 @@ void sim::create_snapshot_if_desired(void)
if(Sp.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
NgbTree.treefree();
Sp.TimeBinsGravity.timebins_free();
diff --git a/src/sph/density.cc b/src/sph/density.cc
index 35859a3a23c1a3389bf8bd91a112b66101d9a571..d152b6779009830362b434bc16c8f77a06c12aea 100644
--- a/src/sph/density.cc
+++ b/src/sph/density.cc
@@ -556,6 +556,9 @@ void sph::density(int *list, int ntarget)
All.Timebase_interval;
double dtime = All.cf_atime * dt / All.cf_atime_hubble_a;
SphP[target].set_viscosity_coefficient(dtime);
+#endif
+#ifdef ADAPTIVE_HYDRO_SOFTENING
+ Tp->P[target].setSofteningClass(Tp->get_softeningtype_for_hydro_particle(target));
#endif
/* now check whether we had enough neighbours */
double desnumngb = All.DesNumNgb;
@@ -890,9 +893,9 @@ void sph::density_evaluate_kernel(pinfo &pdat)
}
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_sum = cs_i + cs_j;
@@ -903,10 +906,9 @@ void sph::density_evaluate_kernel(pinfo &pdat)
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);
@@ -914,10 +916,11 @@ void sph::density_evaluate_kernel(pinfo &pdat)
}
#ifdef TIMEDEP_ART_VISC
- for(int i = 0; i < vector_length; i++) {
- if(decayVel_loc[i] > targetSphP->decayVel)
- targetSphP->decayVel = decayVel_loc[i];
- }
+ for(int i = 0; i < vector_length; i++)
+ {
+ if(decayVel_loc[i] > targetSphP->decayVel)
+ targetSphP->decayVel = decayVel_loc[i];
+ }
#endif
}
@@ -1025,9 +1028,9 @@ void sph::density_evaluate_kernel(pinfo &pdat)
}
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;
if(vdotr2 < 0)
diff --git a/src/subfind/subfind_density.cc b/src/subfind/subfind_density.cc
index e148ddda5f90e410b4b90a68ceacb1e15fda58ef..4bee2aa8649526515dd88636be84c977b7cfa7be 100644
--- a/src/subfind/subfind_density.cc
+++ b/src/subfind/subfind_density.cc
@@ -151,7 +151,7 @@ class subdens_comm : public generic_comm<subdens_in, subdens_out, T_tree, T_doma
{
int p, type;
double mass, r2;
- particle_data *P;
+ typename T_partset::pdata *P;
if(no < Tree->MaxPart) /* single particle */
{
diff --git a/src/subfind/subfind_findlinkngb.cc b/src/subfind/subfind_findlinkngb.cc
index 6e0edb0e8f18fe0b9486347ba89ff9d2df7dccb2..69bcc845569d6b2395e75f30db72a3a5def003db 100644
--- a/src/subfind/subfind_findlinkngb.cc
+++ b/src/subfind/subfind_findlinkngb.cc
@@ -122,7 +122,7 @@ class nearest_comm : public generic_comm<nearest_in, nearest_out, T_tree, T_doma
{
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 */
diff --git a/src/subfind/subfind_nearesttwo.cc b/src/subfind/subfind_nearesttwo.cc
index fb3cef0bcca198e5c5f0f96417c9d01d8a9439fc..074a0346807e1ef4d5b0f0f7d58b5621884add8d 100644
--- a/src/subfind/subfind_nearesttwo.cc
+++ b/src/subfind/subfind_nearesttwo.cc
@@ -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 */
{
- particle_data *P = Tree->get_Pp(no, shmrank);
+ auto *P = Tree->get_Pp(no, shmrank);
subfind_data *PS = Tree->get_PSp(no, shmrank);
no = Tree->get_nextnodep(shmrank)[no]; /* note: here shmrank cannot change */
diff --git a/src/time_integration/kicks.cc b/src/time_integration/kicks.cc
index bc2697b0d3366abfefeb44a674557fdcf522a741..1d552ae959efad629c2f5178e941caf61873a140 100644
--- a/src/time_integration/kicks.cc
+++ b/src/time_integration/kicks.cc
@@ -557,6 +557,9 @@ void sim::hydro_force(int step_indicator)
if(step_indicator == SECOND_HALF_STEP)
{
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].VelPred[0] += (Sp.SphP[target].HydroAccel[0] - Old[i].HydroAccel[0]) * dt_hydrokick;