Skip to content
GitLab
Commit f3b44d59 authored by Matthieu Schaller's avatar Matthieu Schaller
Browse files

Merge branch 'bug_fixes' into 'master' 

Bug fixes



See merge request !202
parents f53617b2 5ca16c2f
Hide whitespace changes
Inline Side-by-side
examples/CosmoVolume/cosmoVolume.yml
View file @ f3b44d59
......@@ -23,11 +23,6 @@ Snapshots:
basename: cosmo # Common part of the name of output files
time_first: 0. # Time of the first output (in internal units)
delta_time: 0.05 # Time difference between consecutive outputs (in internal units)
UnitMass_in_cgs: 1 # Grams
UnitLength_in_cgs: 1 # Centimeters
UnitVelocity_in_cgs: 1 # Centimeters per second
UnitCurrent_in_cgs: 1 # Amperes
UnitTemp_in_cgs: 1 # Kelvin
# Parameters governing the conserved quantities statistics
Statistics:
......
examples/EAGLE_12/eagle_12.yml
View file @ f3b44d59
# Define the system of units to use internally.
InternalUnitSystem:
UnitMass_in_cgs: 1 # Grams
UnitLength_in_cgs: 1 # Centimeters
UnitVelocity_in_cgs: 1 # Centimeters per second
UnitCurrent_in_cgs: 1 # Amperes
UnitTemp_in_cgs: 1 # Kelvin
UnitMass_in_cgs: 1.989e43 # 10^10 M_sun in grams
UnitLength_in_cgs: 3.085678e24 # Mpc in centimeters
UnitVelocity_in_cgs: 1e5 # km/s in centimeters per second
UnitCurrent_in_cgs: 1 # Amperes
UnitTemp_in_cgs: 1 # Kelvin
# Parameters for the task scheduling
Scheduler:
......@@ -14,20 +14,15 @@ Scheduler:
# Parameters governing the time integration
TimeIntegration:
time_begin: 0. # The starting time of the simulation (in internal units).
time_end: 1. # The end time of the simulation (in internal units).
dt_min: 1e-7 # The minimal time-step size of the simulation (in internal units).
dt_max: 1e-2 # The maximal time-step size of the simulation (in internal units).
time_end: 1e-2 # The end time of the simulation (in internal units).
dt_min: 1e-10 # The minimal time-step size of the simulation (in internal units).
dt_max: 1e-4 # The maximal time-step size of the simulation (in internal units).
# Parameters governing the snapshots
Snapshots:
basename: eagle # Common part of the name of output files
time_first: 0. # Time of the first output (in internal units)
delta_time: 0.05 # Time difference between consecutive outputs (in internal units)
UnitMass_in_cgs: 1 # Grams
UnitLength_in_cgs: 1 # Centimeters
UnitVelocity_in_cgs: 1 # Centimeters per second
UnitCurrent_in_cgs: 1 # Amperes
UnitTemp_in_cgs: 1 # Kelvin
delta_time: 1e-3 # Time difference between consecutive outputs (in internal units)
# Parameters governing the conserved quantities statistics
Statistics:
......
examples/EAGLE_25/eagle_25.yml
View file @ f3b44d59
# Define the system of units to use internally.
InternalUnitSystem:
UnitMass_in_cgs: 1 # Grams
UnitLength_in_cgs: 1 # Centimeters
UnitVelocity_in_cgs: 1 # Centimeters per second
UnitCurrent_in_cgs: 1 # Amperes
UnitTemp_in_cgs: 1 # Kelvin
UnitMass_in_cgs: 1.989e43 # 10^10 M_sun in grams
UnitLength_in_cgs: 3.085678e24 # Mpc in centimeters
UnitVelocity_in_cgs: 1e5 # km/s in centimeters per second
UnitCurrent_in_cgs: 1 # Amperes
UnitTemp_in_cgs: 1 # Kelvin
# Parameters for the task scheduling
Scheduler:
......@@ -14,20 +14,15 @@ Scheduler:
# Parameters governing the time integration
TimeIntegration:
time_begin: 0. # The starting time of the simulation (in internal units).
time_end: 1. # The end time of the simulation (in internal units).
dt_min: 1e-7 # The minimal time-step size of the simulation (in internal units).
dt_max: 1e-2 # The maximal time-step size of the simulation (in internal units).
time_end: 1e-2 # The end time of the simulation (in internal units).
dt_min: 1e-10 # The minimal time-step size of the simulation (in internal units).
dt_max: 1e-4 # The maximal time-step size of the simulation (in internal units).
# Parameters governing the snapshots
Snapshots:
basename: eagle # Common part of the name of output files
time_first: 0. # Time of the first output (in internal units)
delta_time: 0.05 # Time difference between consecutive outputs (in internal units)
UnitMass_in_cgs: 1 # Grams
UnitLength_in_cgs: 1 # Centimeters
UnitVelocity_in_cgs: 1 # Centimeters per second
UnitCurrent_in_cgs: 1 # Amperes
UnitTemp_in_cgs: 1 # Kelvin
delta_time: 1e-3 # Time difference between consecutive outputs (in internal units)
# Parameters governing the conserved quantities statistics
Statistics:
......
examples/EAGLE_50/eagle_50.yml
View file @ f3b44d59
# Define the system of units to use internally.
InternalUnitSystem:
UnitMass_in_cgs: 1 # Grams
UnitLength_in_cgs: 1 # Centimeters
UnitVelocity_in_cgs: 1 # Centimeters per second
UnitCurrent_in_cgs: 1 # Amperes
UnitTemp_in_cgs: 1 # Kelvin
UnitMass_in_cgs: 1.989e43 # 10^10 M_sun in grams
UnitLength_in_cgs: 3.085678e24 # Mpc in centimeters
UnitVelocity_in_cgs: 1e5 # km/s in centimeters per second
UnitCurrent_in_cgs: 1 # Amperes
UnitTemp_in_cgs: 1 # Kelvin
# Parameters for the task scheduling
Scheduler:
......@@ -14,20 +14,15 @@ Scheduler:
# Parameters governing the time integration
TimeIntegration:
time_begin: 0. # The starting time of the simulation (in internal units).
time_end: 1. # The end time of the simulation (in internal units).
dt_min: 1e-7 # The minimal time-step size of the simulation (in internal units).
dt_max: 1e-2 # The maximal time-step size of the simulation (in internal units).
time_end: 1e-2 # The end time of the simulation (in internal units).
dt_min: 1e-10 # The minimal time-step size of the simulation (in internal units).
dt_max: 1e-4 # The maximal time-step size of the simulation (in internal units).
# Parameters governing the snapshots
Snapshots:
basename: eagle # Common part of the name of output files
time_first: 0. # Time of the first output (in internal units)
delta_time: 0.05 # Time difference between consecutive outputs (in internal units)
UnitMass_in_cgs: 1 # Grams
UnitLength_in_cgs: 1 # Centimeters
UnitVelocity_in_cgs: 1 # Centimeters per second
UnitCurrent_in_cgs: 1 # Amperes
UnitTemp_in_cgs: 1 # Kelvin
delta_time: 1e-3 # Time difference between consecutive outputs (in internal units)
# Parameters governing the conserved quantities statistics
Statistics:
......
examples/ExternalPointMass/externalPointMass.yml
View file @ f3b44d59
......@@ -18,11 +18,6 @@ Snapshots:
basename: pointMass # Common part of the name of output files
time_first: 0. # Time of the first output (in internal units)
delta_time: 0.02 # Time difference between consecutive outputs (in internal units)
UnitMass_in_cgs: 1.9885e33 # Grams
UnitLength_in_cgs: 3.0856776e21 # Centimeters
UnitVelocity_in_cgs: 1e5 # Centimeters per second
UnitCurrent_in_cgs: 1 # Amperes
UnitTemp_in_cgs: 1 # Kelvin
# Parameters governing the conserved quantities statistics
Statistics:
......
examples/MultiTypes/multiTypes.yml
View file @ f3b44d59
......@@ -18,11 +18,6 @@ Snapshots:
basename: multiTypes # Common part of the name of output files
time_first: 0. # Time of the first output (in internal units)
delta_time: 0.01 # Time difference between consecutive outputs (in internal units)
UnitMass_in_cgs: 1 # Grams
UnitLength_in_cgs: 1 # Centimeters
UnitVelocity_in_cgs: 1 # Centimeters per second
UnitCurrent_in_cgs: 1 # Amperes
UnitTemp_in_cgs: 1 # Kelvin
# Parameters governing the conserved quantities statistics
Statistics:
......
examples/SedovBlast/sedov.yml
View file @ f3b44d59
......@@ -18,11 +18,6 @@ Snapshots:
basename: sedov # Common part of the name of output files
time_first: 0. # Time of the first output (in internal units)
delta_time: 0.1 # Time difference between consecutive outputs (in internal units)
UnitMass_in_cgs: 1 # Grams
UnitLength_in_cgs: 1 # Centimeters
UnitVelocity_in_cgs: 1 # Centimeters per second
UnitCurrent_in_cgs: 1 # Amperes
UnitTemp_in_cgs: 1 # Kelvin
# Parameters governing the conserved quantities statistics
Statistics:
......
examples/SodShock/sodShock.yml
View file @ f3b44d59
......@@ -18,11 +18,6 @@ Snapshots:
basename: sod # Common part of the name of output files
time_first: 0. # Time of the first output (in internal units)
delta_time: 0.1 # Time difference between consecutive outputs (in internal units)
UnitMass_in_cgs: 1 # Grams
UnitLength_in_cgs: 1 # Centimeters
UnitVelocity_in_cgs: 1 # Centimeters per second
UnitCurrent_in_cgs: 1 # Amperes
UnitTemp_in_cgs: 1 # Kelvin
# Parameters governing the conserved quantities statistics
Statistics:
......
examples/UniformBox/uniformBox.yml
View file @ f3b44d59
......@@ -18,11 +18,6 @@ Snapshots:
basename: uniformBox # Common part of the name of output files
time_first: 0. # Time of the first output (in internal units)
delta_time: 0.01 # Time difference between consecutive outputs (in internal units)
UnitMass_in_cgs: 1 # Grams
UnitLength_in_cgs: 1 # Centimeters
UnitVelocity_in_cgs: 1 # Centimeters per second
UnitCurrent_in_cgs: 1 # Amperes
UnitTemp_in_cgs: 1 # Kelvin
# Parameters governing the conserved quantities statistics
Statistics:
......
examples/main.c
View file @ f3b44d59
......@@ -358,8 +358,8 @@ int main(int argc, char *argv[]) {
&periodic, &flag_entropy_ICs, myrank, nr_nodes,
MPI_COMM_WORLD, MPI_INFO_NULL, dry_run);
#else
read_ic_serial(ICfileName, &us, dim, &parts, &gparts, &Ngas, &Ngpart, &periodic,
&flag_entropy_ICs, myrank, nr_nodes, MPI_COMM_WORLD,
read_ic_serial(ICfileName, &us, dim, &parts, &gparts, &Ngas, &Ngpart,
&periodic, &flag_entropy_ICs, myrank, nr_nodes, MPI_COMM_WORLD,
MPI_INFO_NULL, dry_run);
#endif
#else
......
examples/parameter_example.yml
View file @ f3b44d59
......@@ -25,11 +25,11 @@ Snapshots:
basename: output # Common part of the name of output files
time_first: 0. # Time of the first output (in internal units)
delta_time: 0.01 # Time difference between consecutive outputs (in internal units)
UnitMass_in_cgs: 1 # Unit system for the outputs (Grams)
UnitLength_in_cgs: 1 # Unit system for the outputs (Centimeters)
UnitVelocity_in_cgs: 1 # Unit system for the outputs (Centimeters per second)
UnitCurrent_in_cgs: 1 # Unit system for the outputs (Amperes)
UnitTemp_in_cgs: 1 # Unit system for the outputs (Kelvin)
UnitMass_in_cgs: 1 # (Optional) Unit system for the outputs (Grams)
UnitLength_in_cgs: 1 # (Optional) Unit system for the outputs (Centimeters)
UnitVelocity_in_cgs: 1 # (Optional) Unit system for the outputs (Centimeters per second)
UnitCurrent_in_cgs: 1 # (Optional) Unit system for the outputs (Amperes)
UnitTemp_in_cgs: 1 # (Optional) Unit system for the outputs (Kelvin)
# Parameters governing the conserved quantities statistics
Statistics:
......
src/adiabatic_index.h
View file @ f3b44d59
......@@ -60,15 +60,12 @@ __attribute__((always_inline)) INLINE static float pow_gamma(float x) {
#if defined(HYDRO_GAMMA_5_3)
const float x2 = x * x;
const float x5 = x2 * x2 * x;
return cbrtf(x5);
const float cbrt = cbrtf(x); /* x^(1/3) */
return cbrt * cbrt * x; /* x^(5/3) */
#elif defined(HYDRO_GAMMA_4_3)
const float x2 = x * x;
const float x4 = x2 * x2;
return cbrtf(x4);
return cbrtf(x) * x; /* x^(4/3) */
#elif defined(HYDRO_GAMMA_2_1)
......@@ -93,11 +90,12 @@ __attribute__((always_inline)) INLINE static float pow_gamma_minus_one(
#if defined(HYDRO_GAMMA_5_3)
return cbrtf(x * x);
const float cbrt = cbrtf(x); /* x^(1/3) */
return cbrt * cbrt; /* x^(2/3) */
#elif defined(HYDRO_GAMMA_4_3)
return cbrtf(x);
return cbrtf(x); /* x^(1/3) */
#elif defined(HYDRO_GAMMA_2_1)
......@@ -122,11 +120,12 @@ __attribute__((always_inline)) INLINE static float pow_minus_gamma_minus_one(
#if defined(HYDRO_GAMMA_5_3)
return 1.f / cbrtf(x * x);
const float inv_cbrt = 1.f / cbrtf(x); /* x^(-1/3) */
return inv_cbrt * inv_cbrt; /* x^(-2/3) */
#elif defined(HYDRO_GAMMA_4_3)
return 1.f / cbrtf(x);
return 1.f / cbrtf(x); /* x^(-1/3) */
#elif defined(HYDRO_GAMMA_2_1)
......
src/common_io.c
View file @ f3b44d59
......@@ -120,7 +120,7 @@ size_t sizeOfType(enum DATA_TYPE type) {
* Returns an error if the type is not FLOAT or DOUBLE
*/
int isDoublePrecision(enum DATA_TYPE type) {
switch (type) {
case FLOAT:
return 0;
......
src/drift.h
View file @ f3b44d59
......@@ -37,7 +37,8 @@
* @param ti_current Integer end of time-step
*/
__attribute__((always_inline)) INLINE static void drift_gpart(
struct gpart* gp, float dt, double timeBase, int ti_old, int ti_current) {
struct gpart *restrict gp, float dt, double timeBase, int ti_old,
int ti_current) {
/* Drift... */
gp->x[0] += gp->v_full[0] * dt;
gp->x[1] += gp->v_full[1] * dt;
......@@ -60,8 +61,8 @@ __attribute__((always_inline)) INLINE static void drift_gpart(
* @param ti_current Integer end of time-step
*/
__attribute__((always_inline)) INLINE static void drift_part(
struct part* p, struct xpart* xp, float dt, double timeBase, int ti_old,
int ti_current) {
struct part *restrict p, struct xpart *restrict xp, float dt,
double timeBase, int ti_old, int ti_current) {
/* Useful quantity */
const float h_inv = 1.0f / p->h;
......
src/engine.c
View file @ f3b44d59
......@@ -2672,7 +2672,7 @@ void engine_dump_snapshot(struct engine *e) {
struct clocks_time time1, time2;
clocks_gettime(&time1);
if (e->verbose) message("writing snapshot at t=%f.", e->time);
if (e->verbose) message("writing snapshot at t=%e.", e->time);
/* Dump... */
#if defined(WITH_MPI)
......@@ -2808,7 +2808,7 @@ void engine_init(struct engine *e, struct space *s,
e->ti_nextSnapshot = 0;
parser_get_param_string(params, "Snapshots:basename", e->snapshotBaseName);
e->snapshotUnits = malloc(sizeof(struct UnitSystem));
units_init(e->snapshotUnits, params, "Snapshots");
units_init_default(e->snapshotUnits, params, "Snapshots", internal_units);
e->dt_min = parser_get_param_double(params, "TimeIntegration:dt_min");
e->dt_max = parser_get_param_double(params, "TimeIntegration:dt_max");
e->file_stats = NULL;
......@@ -3228,6 +3228,6 @@ void engine_compute_next_snapshot_time(struct engine *e) {
const float next_snapshot_time =
e->ti_nextSnapshot * e->timeBase + e->timeBegin;
if (e->verbose)
message("Next output time set to t=%f.", next_snapshot_time);
message("Next output time set to t=%e.", next_snapshot_time);
}
}
src/hydro/Default/hydro.h
View file @ f3b44d59
......@@ -28,8 +28,8 @@
*
*/
__attribute__((always_inline)) INLINE static float hydro_compute_timestep(
const struct part* p, const struct xpart* xp,
const struct hydro_props* hydro_properties) {
const struct part *restrict p, const struct xpart *restrict xp,
const struct hydro_props *restrict hydro_properties) {
const float CFL_condition = hydro_properties->CFL_condition;
......@@ -55,7 +55,7 @@ __attribute__((always_inline)) INLINE static float hydro_compute_timestep(
* @param xp The extended particle data to act upon
*/
__attribute__((always_inline)) INLINE static void hydro_first_init_part(
struct part* p, struct xpart* xp) {}
struct part *restrict p, struct xpart *restrict xp) {}
/**
* @brief Prepares a particle for the density calculation.
......@@ -66,15 +66,15 @@ __attribute__((always_inline)) INLINE static void hydro_first_init_part(
* @param p The particle to act upon
*/
__attribute__((always_inline)) INLINE static void hydro_init_part(
struct part* p) {
struct part *restrict p) {
p->density.wcount = 0.f;
p->density.wcount_dh = 0.f;
p->rho = 0.f;
p->rho_dh = 0.f;
p->density.div_v = 0.f;
p->density.curl_v[0] = 0.f;
p->density.curl_v[1] = 0.f;
p->density.curl_v[2] = 0.f;
p->density.rot_v[0] = 0.f;
p->density.rot_v[1] = 0.f;
p->density.rot_v[2] = 0.f;
}
/**
......@@ -87,7 +87,7 @@ __attribute__((always_inline)) INLINE static void hydro_init_part(
* @param time The current time
*/
__attribute__((always_inline)) INLINE static void hydro_end_density(
struct part* p, float time) {
struct part *restrict p, float time) {
/* Some smoothing length multiples. */
const float h = p->h;
......@@ -110,6 +110,14 @@ __attribute__((always_inline)) INLINE static void hydro_end_density(
/* Compute the derivative term */
p->rho_dh = 1.f / (1.f + 0.33333333f * p->h * p->rho_dh * irho);
/* Finish calculation of the velocity curl components */
p->density.rot_v[0] *= ih4 * irho;
p->density.rot_v[1] *= ih4 * irho;
p->density.rot_v[2] *= ih4 * irho;
/* Finish calculation of the velocity divergence */
p->density.div_v *= ih4 * irho;
}
/**
......@@ -122,7 +130,8 @@ __attribute__((always_inline)) INLINE static void hydro_end_density(
* @param time The current time
*/
__attribute__((always_inline)) INLINE static void hydro_prepare_force(
struct part* p, struct xpart* xp, int ti_current, double timeBase) {
struct part *restrict p, struct xpart *restrict xp, int ti_current,
double timeBase) {
/* Some smoothing length multiples. */
const float h = p->h;
......@@ -132,10 +141,10 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force(
/* Pre-compute some stuff for the balsara switch. */
const float normDiv_v = fabs(p->density.div_v / p->rho * ih4);
const float normCurl_v = sqrtf(p->density.curl_v[0] * p->density.curl_v[0] +
p->density.curl_v[1] * p->density.curl_v[1] +
p->density.curl_v[2] * p->density.curl_v[2]) /
p->rho * ih4;
const float normRot_v = sqrtf(p->density.rot_v[0] * p->density.rot_v[0] +
p->density.rot_v[1] * p->density.rot_v[1] +
p->density.rot_v[2] * p->density.rot_v[2]) /
p->rho * ih4;
/* Compute this particle's sound speed. */
const float u = p->u;
......@@ -146,7 +155,7 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force(
p->force.POrho2 = u * hydro_gamma_minus_one / (p->rho * xp->omega);
/* Balsara switch */
p->force.balsara = normDiv_v / (normDiv_v + normCurl_v + 0.0001f * fc * ih);
p->force.balsara = normDiv_v / (normDiv_v + normRot_v + 0.0001f * fc * ih);
/* Viscosity parameter decay time */
const float tau = h / (2.f * const_viscosity_length * p->force.c);
......@@ -171,7 +180,7 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force(
* @param p The particle to act upon
*/
__attribute__((always_inline)) INLINE static void hydro_reset_acceleration(
struct part* p) {
struct part *restrict p) {
/* Reset the acceleration. */
p->a_hydro[0] = 0.0f;
......@@ -194,7 +203,8 @@ __attribute__((always_inline)) INLINE static void hydro_reset_acceleration(
* @param timeBase The minimal time-step size
*/
__attribute__((always_inline)) INLINE static void hydro_predict_extra(
struct part* p, struct xpart* xp, int t0, int t1, double timeBase) {
struct part *restrict p, struct xpart *restrict xp, int t0, int t1,
double timeBase) {
float u, w;
const float dt = (t1 - t0) * timeBase;
......@@ -218,7 +228,7 @@ __attribute__((always_inline)) INLINE static void hydro_predict_extra(
* @param p The particle to act upon
*/
__attribute__((always_inline)) INLINE static void hydro_end_force(
struct part* p) {}
struct part *restrict p) {}
/**
* @brief Kick the additional variables
......@@ -229,17 +239,18 @@ __attribute__((always_inline)) INLINE static void hydro_end_force(
* @param half_dt The half time-step for this kick
*/
__attribute__((always_inline)) INLINE static void hydro_kick_extra(
struct part* p, struct xpart* xp, float dt, float half_dt) {}
struct part *restrict p, struct xpart *restrict xp, float dt,
float half_dt) {}
/**
* @brief Converts hydro quantity of a particle
* @brief Converts hydro quantity of a particle at the start of a run
*
* Requires the density to be known
*
* @param p The particle to act upon
*/
__attribute__((always_inline)) INLINE static void hydro_convert_quantities(
struct part* p) {}
struct part *restrict p) {}
/**
* @brief Returns the internal energy of a particle
......@@ -248,7 +259,7 @@ __attribute__((always_inline)) INLINE static void hydro_convert_quantities(
* @param dt Time since the last kick
*/
__attribute__((always_inline)) INLINE static float hydro_get_internal_energy(
const struct part* p, float dt) {
const struct part *restrict p, float dt) {
return p->u;
}
src/hydro/Default/hydro_iact.h
View file @ f3b44d59
......@@ -79,8 +79,8 @@ __attribute__((always_inline)) INLINE static void runner_iact_density(
pi->density.wcount += wi;
pi->density.wcount_dh -= xi * wi_dx;
pi->density.div_v += mj * dvdr * wi_dx;
for (k = 0; k < 3; k++) pi->density.curl_v[k] += mj * curlvr[k] * wi_dx;
pi->density.div_v -= mj * dvdr * wi_dx;
for (k = 0; k < 3; k++) pi->density.rot_v[k] += mj * curlvr[k] * wi_dx;
/* Compute density of pj. */
h_inv = 1.0 / hj;
......@@ -92,8 +92,8 @@ __attribute__((always_inline)) INLINE static void runner_iact_density(
pj->density.wcount += wj;
pj->density.wcount_dh -= xj * wj_dx;
pj->density.div_v += mi * dvdr * wj_dx;
for (k = 0; k < 3; k++) pj->density.curl_v[k] += mi * curlvr[k] * wj_dx;
pj->density.div_v -= mi * dvdr * wj_dx;
for (k = 0; k < 3; k++) pj->density.rot_v[k] += mi * curlvr[k] * wj_dx;
}
/**
......@@ -111,7 +111,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_vec_density(
vector dx[3], dv[3];
vector vi[3], vj[3];
vector dvdr, div_vi, div_vj;
vector curlvr[3], curl_vi[3], curl_vj[3];
vector curlvr[3], rot_vi[3], rot_vj[3];
int k, j;
#if VEC_SIZE == 8
......@@ -178,28 +178,28 @@ __attribute__((always_inline)) INLINE static void runner_iact_vec_density(
wcounti.v = wi.v;
wcounti_dh.v = xi.v * wi_dx.v;
div_vi.v = mj.v * dvdr.v * wi_dx.v;
for (k = 0; k < 3; k++) curl_vi[k].v = mj.v * curlvr[k].v * wi_dx.v;
for (k = 0; k < 3; k++) rot_vi[k].v = mj.v * curlvr[k].v * wi_dx.v;
rhoj.v = mi.v * wj.v;
rhoj_dh.v = mi.v * (vec_set1(3.0f) * wj.v + xj.v * wj_dx.v);
wcountj.v = wj.v;
wcountj_dh.v = xj.v * wj_dx.v;
div_vj.v = mi.v * dvdr.v * wj_dx.v;
for (k = 0; k < 3; k++) curl_vj[k].v = mi.v * curlvr[k].v * wj_dx.v;
for (k = 0; k < 3; k++) rot_vj[k].v = mi.v * curlvr[k].v * wj_dx.v;
for (k = 0; k < VEC_SIZE; k++) {
pi[k]->rho += rhoi.f[k];
pi[k]->rho_dh -= rhoi_dh.f[k];
pi[k]->density.wcount += wcounti.f[k];
pi[k]->density.wcount_dh -= wcounti_dh.f[k];
pi[k]->density.div_v += div_vi.f[k];
for (j = 0; j < 3; j++) pi[k]->density.curl_v[j] += curl_vi[j].f[k];
pi[k]->density.div_v -= div_vi.f[k];
for (j = 0; j < 3; j++) pi[k]->density.rot_v[j] += rot_vi[j].f[k];
pj[k]->rho += rhoj.f[k];
pj[k]->rho_dh -= rhoj_dh.f[k];
pj[k]->density.wcount += wcountj.f[k];
pj[k]->density.wcount_dh -= wcountj_dh.f[k];
pj[k]->density.div_v += div_vj.f[k];
for (j = 0; j < 3; j++) pj[k]->density.curl_v[j] += curl_vj[j].f[k];
pj[k]->density.div_v -= div_vj.f[k];
for (j = 0; j < 3; j++) pj[k]->density.rot_v[j] += rot_vj[j].f[k];
}
#else
......@@ -254,8 +254,8 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_density(
pi->density.wcount += wi;
pi->density.wcount_dh -= xi * wi_dx;
pi->density.div_v += mj * dvdr * wi_dx;
for (k = 0; k < 3; k++) pi->density.curl_v[k] += mj * curlvr[k] * wi_dx;
pi->density.div_v -= mj * dvdr * wi_dx;
for (k = 0; k < 3; k++) pi->density.rot_v[k] += mj * curlvr[k] * wi_dx;
}
/**
......@@ -273,7 +273,7 @@ runner_iact_nonsym_vec_density(float *R2, float *Dx, float *Hi, float *Hj,
vector dx[3], dv[3];
vector vi[3], vj[3];
vector dvdr;
vector curlvr[3], curl_vi[3];
vector curlvr[3], rot_vi[3];
int k, j;
#if VEC_SIZE == 8
......@@ -331,15 +331,15 @@ runner_iact_nonsym_vec_density(float *R2, float *Dx, float *Hi, float *Hj,
wcounti.v = wi.v;
wcounti_dh.v = xi.v * wi_dx.v;
div_vi.v = mj.v * dvdr.v * wi_dx.v;
for (k = 0; k < 3; k++) curl_vi[k].v = mj.v * curlvr[k].v * wi_dx.v;
for (k = 0; k < 3; k++) rot_vi[k].v = mj.v * curlvr[k].v * wi_dx.v;
for (k = 0; k < VEC_SIZE; k++) {
pi[k]->rho += rhoi.f[k];
pi[k]->rho_dh -= rhoi_dh.f[k];
pi[k]->density.wcount += wcounti.f[k];
pi[k]->density.wcount_dh -= wcounti_dh.f[k];
pi[k]->density.div_v += div_vi.f[k];
for (j = 0; j < 3; j++) pi[k]->density.curl_v[j] += curl_vi[j].f[k];
pi[k]->density.div_v -= div_vi.f[k];
for (j = 0; j < 3; j++) pi[k]->density.rot_v[j] += rot_vi[j].f[k];
}
#else
......
src/hydro/Default/hydro_part.h
View file @ f3b44d59
......@@ -80,7 +80,7 @@ struct part {
float wcount_dh;
/* Particle velocity curl. */
float curl_v[3];
float rot_v[3];