diff --git a/src/adiabatic_index.h b/src/adiabatic_index.h index 31a4f56699b46b086d1e54cd2fd9643eee00c686..59937db3c8d09fc7e6e969de0aee60f01a2e5a2c 100644 --- a/src/adiabatic_index.h +++ b/src/adiabatic_index.h @@ -432,40 +432,4 @@ __attribute__((always_inline)) INLINE static float pow_one_over_gamma(float x) { #endif } -/** - * @brief Return the argument to the power one minus two over the adiabatic - * index - * - * Computes \f$x^{1 - \frac{2}{\gamma}}\f$. - * - * @param x Argument - * @return Argument to the power one minus two over the adiabatic index - */ -__attribute__((always_inline)) INLINE static float pow_one_minus_two_over_gamma( - float x) { - -#if defined(HYDRO_GAMMA_5_3) - - return powf(x, -0.2f); - -#elif defined(HYDRO_GAMMA_7_5) - - return powf(x, -0.428571429f); - -#elif defined(HYDRO_GAMMA_4_3) - - return 1.f / sqrtf(x); - -#elif defined(HYDRO_GAMMA_2_1) - - return 1.f; - -#else - - error("The adiabatic index is not defined !"); - return 0.f; - -#endif -} - #endif /* SWIFT_ADIABATIC_INDEX_H */ diff --git a/src/const.h b/src/const.h index 4bdcc3aff1cbe134c1c7c4235dcda864136a6cea..afad73211546cf81cbb0e6a1179569e588952fec 100644 --- a/src/const.h +++ b/src/const.h @@ -40,9 +40,9 @@ #define const_isothermal_internal_energy 20.2615290634f /* Dimensionality of the problem */ -//#define HYDRO_DIMENSION_3D +#define HYDRO_DIMENSION_3D //#define HYDRO_DIMENSION_2D -#define HYDRO_DIMENSION_1D +//#define HYDRO_DIMENSION_1D /* Hydrodynamical adiabatic index. */ #define HYDRO_GAMMA_5_3 @@ -64,8 +64,8 @@ /* SPH variant to use */ //#define MINIMAL_SPH -//#define GADGET2_SPH -#define HOPKINS_PE_SPH +#define GADGET2_SPH +//#define HOPKINS_PE_SPH //#define DEFAULT_SPH //#define GIZMO_SPH diff --git a/src/hydro/PressureEntropy/hydro_io.h b/src/hydro/PressureEntropy/hydro_io.h index 3263b6b7a9f3a37ac9ba36c2aad53685795246fc..9914a656466f3f0d0a5eeb79b511706d7068ffc6 100644 --- a/src/hydro/PressureEntropy/hydro_io.h +++ b/src/hydro/PressureEntropy/hydro_io.h @@ -56,9 +56,6 @@ void hydro_read_particles(struct part* parts, struct io_props* list, parts, mass); list[3] = io_make_input_field("SmoothingLength", FLOAT, 1, COMPULSORY, UNIT_CONV_LENGTH, parts, h); - /* list[4] = io_make_input_field("InternalEnergy", FLOAT, 1, COMPULSORY, */ - /* UNIT_CONV_ENERGY_PER_UNIT_MASS, parts, - * entropy); */ list[4] = io_make_input_field("InternalEnergy", FLOAT, 1, COMPULSORY, UNIT_CONV_ENTROPY_PER_UNIT_MASS, parts, entropy); @@ -127,8 +124,11 @@ void writeSPHflavour(hid_t h_grpsph) { /* Viscosity and thermal conduction */ /* Nothing in this minimal model... */ writeAttribute_s(h_grpsph, "Thermal Conductivity Model", "No treatment"); - writeAttribute_s(h_grpsph, "Viscosity Model", - "Minimal treatment as in Monaghan (1992)"); + writeAttribute_s( + h_grpsph, "Viscosity Model", + "as in Springel (2005), i.e. Monaghan (1992) with Balsara (1995) switch"); + writeAttribute_f(h_grpsph, "Viscosity alpha", const_viscosity_alpha); + writeAttribute_f(h_grpsph, "Viscosity beta", 3.f); /* Time integration properties */ writeAttribute_f(h_grpsph, "Maximal Delta u change over dt", diff --git a/tests/test125cells.c b/tests/test125cells.c index e666658f43de135e3e72521b52f2a688c596a6f6..d423efba3fda764310721586e14192afdee18a85 100644 --- a/tests/test125cells.c +++ b/tests/test125cells.c @@ -95,6 +95,8 @@ void set_energy_state(struct part *part, enum pressure_field press, float size, #if defined(GADGET2_SPH) part->entropy = pressure / pow_gamma(density); +#elif defined(HOPKINS_PE_SPH) + part->entropy = pressure / pow_gamma(density); #elif defined(DEFAULT_SPH) part->u = pressure / (hydro_gamma_minus_one * density); #elif defined(MINIMAL_SPH) diff --git a/tests/test27cells.c b/tests/test27cells.c index 1a1ab88748d922b3e7fbb30a73a10809dca10863..b4c79462153e31d272915808ed5b2fdd086dd5d2 100644 --- a/tests/test27cells.c +++ b/tests/test27cells.c @@ -104,11 +104,18 @@ struct cell *make_cell(size_t n, double *offset, double size, double h, } part->h = size * h / (float)n; part->id = ++(*partId); + #ifdef GIZMO_SPH part->conserved.mass = density * volume / count; #else part->mass = density * volume / count; #endif + +#if defined(HOPKINS_PE_SPH) + part->entropy = 1.f; + part->entropy_one_over_gamma = 1.f; +#endif + part->ti_begin = 0; part->ti_end = 1; ++part; @@ -192,17 +199,14 @@ void dump_particle_fields(char *fileName, struct cell *main_cell, hydro_get_density(&main_cell->parts[pid]), #if defined(GIZMO_SPH) 0.f, +#elif defined(HOPKINS_PE_SPH) + main_cell->parts[pid].density.rho_dh, #else main_cell->parts[pid].rho_dh, #endif main_cell->parts[pid].density.wcount, main_cell->parts[pid].density.wcount_dh, -#if defined(GADGET2_SPH) - main_cell->parts[pid].density.div_v, - main_cell->parts[pid].density.rot_v[0], - main_cell->parts[pid].density.rot_v[1], - main_cell->parts[pid].density.rot_v[2] -#elif defined(DEFAULT_SPH) +#if defined(GADGET2_SPH) || defined(DEFAULT_SPH) || defined(HOPKINS_PE_SPH) main_cell->parts[pid].density.div_v, main_cell->parts[pid].density.rot_v[0], main_cell->parts[pid].density.rot_v[1], @@ -234,14 +238,13 @@ void dump_particle_fields(char *fileName, struct cell *main_cell, cj->parts[pjd].v[2], hydro_get_density(&cj->parts[pjd]), #if defined(GIZMO_SPH) 0.f, +#elif defined(HOPKINS_PE_SPH) + main_cell->parts[pjd].density.rho_dh, #else main_cell->parts[pjd].rho_dh, #endif cj->parts[pjd].density.wcount, cj->parts[pjd].density.wcount_dh, -#if defined(GADGET2_SPH) - cj->parts[pjd].density.div_v, cj->parts[pjd].density.rot_v[0], - cj->parts[pjd].density.rot_v[1], cj->parts[pjd].density.rot_v[2] -#elif defined(DEFAULT_SPH) +#if defined(GADGET2_SPH) || defined(DEFAULT_SPH) || defined(HOPKINS_PE_SPH) cj->parts[pjd].density.div_v, cj->parts[pjd].density.rot_v[0], cj->parts[pjd].density.rot_v[1], cj->parts[pjd].density.rot_v[2] #else diff --git a/tests/testPair.c b/tests/testPair.c index efa1e628c2d57bf7922be8affdd5436ebca2f9cf..d19d852a857e00665b82c604f31986aebcca4d6c 100644 --- a/tests/testPair.c +++ b/tests/testPair.c @@ -138,11 +138,13 @@ void dump_particle_fields(char *fileName, struct cell *ci, struct cell *cj) { ci->parts[pid].v[2], hydro_get_density(&ci->parts[pid]), #if defined(GIZMO_SPH) 0.f, +#elif defined(HOPKINS_PE_SPH) + ci->parts[pid].density.rho_dh, #else - cj->parts[pid].rho_dh, + ci->parts[pid].rho_dh, #endif ci->parts[pid].density.wcount, ci->parts[pid].density.wcount_dh, -#if defined(GADGET2_SPH) || defined(DEFAULT_SPH) +#if defined(GADGET2_SPH) || defined(DEFAULT_SPH) || defined(HOPKINS_PE_SPH) ci->parts[pid].density.div_v, ci->parts[pid].density.rot_v[0], ci->parts[pid].density.rot_v[1], ci->parts[pid].density.rot_v[2] #else @@ -162,11 +164,13 @@ void dump_particle_fields(char *fileName, struct cell *ci, struct cell *cj) { cj->parts[pjd].v[2], hydro_get_density(&cj->parts[pjd]), #if defined(GIZMO_SPH) 0.f, +#elif defined(HOPKINS_PE_SPH) + cj->parts[pjd].density.rho_dh, #else cj->parts[pjd].rho_dh, #endif cj->parts[pjd].density.wcount, cj->parts[pjd].density.wcount_dh, -#if defined(GADGET2_SPH) || defined(DEFAULT_SPH) +#if defined(GADGET2_SPH) || defined(DEFAULT_SPH) || defined(HOPKINS_PE_SPH) cj->parts[pjd].density.div_v, cj->parts[pjd].density.rot_v[0], cj->parts[pjd].density.rot_v[1], cj->parts[pjd].density.rot_v[2] #else