Skip to content
Snippets Groups Projects
Commit de96fa54 authored by Bert Vandenbroucke's avatar Bert Vandenbroucke
Browse files

Shadowswift 1D and 2D (hydro-only) now work again!

parent 0c73ae40
Branches
Tags
1 merge request!575Shadowswift
Show whitespace changes
Inline Side-by-side
src/hydro/Shadowswift/hydro.h
+ 69
− 66
View file @ de96fa54
...@@ -43,15 +43,23 @@ __attribute__((always_inline)) INLINE static float hydro_compute_timestep( ...@@ -43,15 +43,23 @@ __attribute__((always_inline)) INLINE static float hydro_compute_timestep(
const float CFL_condition = hydro_properties->CFL_condition; const float CFL_condition = hydro_properties->CFL_condition;
float R = get_radius_dimension_sphere(p->cell.volume); float vrel[3];
vrel[0] = p->primitives.v[0] - xp->v_full[0];
vrel[1] = p->primitives.v[1] - xp->v_full[1];
vrel[2] = p->primitives.v[2] - xp->v_full[2];
float vmax =
sqrtf(vrel[0] * vrel[0] + vrel[1] * vrel[1] + vrel[2] * vrel[2]) +
sqrtf(hydro_gamma * p->primitives.P / p->primitives.rho);
vmax = max(vmax, p->timestepvars.vmax);
if (p->timestepvars.vmax == 0.) { const float psize =
/* vmax can be zero in vacuum. We force the time step to become the maximal cosmo->a *
time step in this case */ powf(p->cell.volume / hydro_dimension_unit_sphere, hydro_dimension_inv);
return FLT_MAX; float dt = FLT_MAX;
} else { if (vmax > 0.) {
return CFL_condition * R / fabsf(p->timestepvars.vmax); dt = psize / vmax;
} }
return CFL_condition * dt;
} }
/** /**
...@@ -104,7 +112,7 @@ __attribute__((always_inline)) INLINE static void hydro_first_init_part( ...@@ -104,7 +112,7 @@ __attribute__((always_inline)) INLINE static void hydro_first_init_part(
p->conserved.momentum[2] = mass * p->primitives.v[2]; p->conserved.momentum[2] = mass * p->primitives.v[2];
#ifdef EOS_ISOTHERMAL_GAS #ifdef EOS_ISOTHERMAL_GAS
p->conserved.energy = mass * const_isothermal_internal_energy; p->conserved.energy = mass * gas_internal_energy_from_entropy(0.f, 0.f);
#else #else
p->conserved.energy *= mass; p->conserved.energy *= mass;
#endif #endif
...@@ -119,12 +127,21 @@ __attribute__((always_inline)) INLINE static void hydro_first_init_part( ...@@ -119,12 +127,21 @@ __attribute__((always_inline)) INLINE static void hydro_first_init_part(
p->v[0] = 0.; p->v[0] = 0.;
p->v[1] = 0.; p->v[1] = 0.;
p->v[2] = 0.; p->v[2] = 0.;
#else
p->v[0] = p->primitives.v[0];
p->v[1] = p->primitives.v[1];
p->v[2] = p->primitives.v[2];
#endif #endif
/* set the initial velocity of the cells */ /* set the initial velocity of the cells */
xp->v_full[0] = p->v[0]; xp->v_full[0] = p->v[0];
xp->v_full[1] = p->v[1]; xp->v_full[1] = p->v[1];
xp->v_full[2] = p->v[2]; xp->v_full[2] = p->v[2];
/* ignore accelerations present in the initial condition */
p->a_hydro[0] = 0.0f;
p->a_hydro[1] = 0.0f;
p->a_hydro[2] = 0.0f;
} }
/** /**
...@@ -290,6 +307,12 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force( ...@@ -290,6 +307,12 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force(
p->force.v_full[0] = xp->v_full[0]; p->force.v_full[0] = xp->v_full[0];
p->force.v_full[1] = xp->v_full[1]; p->force.v_full[1] = xp->v_full[1];
p->force.v_full[2] = xp->v_full[2]; p->force.v_full[2] = xp->v_full[2];
p->conserved.flux.mass = 0.0f;
p->conserved.flux.momentum[0] = 0.0f;
p->conserved.flux.momentum[1] = 0.0f;
p->conserved.flux.momentum[2] = 0.0f;
p->conserved.flux.energy = 0.0f;
} }
/** /**
...@@ -423,53 +446,36 @@ __attribute__((always_inline)) INLINE static void hydro_kick_extra( ...@@ -423,53 +446,36 @@ __attribute__((always_inline)) INLINE static void hydro_kick_extra(
struct part* p, struct xpart* xp, float dt, const struct cosmology* cosmo, struct part* p, struct xpart* xp, float dt, const struct cosmology* cosmo,
const struct hydro_props* hydro_props) { const struct hydro_props* hydro_props) {
float vcell[3];
/* Update the conserved variables. We do this here and not in the kick, /* Update the conserved variables. We do this here and not in the kick,
since we need the updated variables below. */ since we need the updated variables below. */
p->conserved.mass += p->conserved.flux.mass; p->conserved.mass += p->conserved.flux.mass * dt;
p->conserved.momentum[0] += p->conserved.flux.momentum[0]; p->conserved.momentum[0] += p->conserved.flux.momentum[0] * dt;
p->conserved.momentum[1] += p->conserved.flux.momentum[1]; p->conserved.momentum[1] += p->conserved.flux.momentum[1] * dt;
p->conserved.momentum[2] += p->conserved.flux.momentum[2]; p->conserved.momentum[2] += p->conserved.flux.momentum[2] * dt;
p->conserved.energy += p->conserved.flux.energy;
#ifdef EOS_ISOTHERMAL_GAS #ifdef EOS_ISOTHERMAL_GAS
/* reset the thermal energy */ /* reset the thermal energy */
p->conserved.energy = p->conserved.mass * const_isothermal_internal_energy; p->conserved.energy =
p->conserved.mass * gas_internal_energy_from_entropy(0.f, 0.f);
#ifdef SHADOWFAX_TOTAL_ENERGY #else
p->conserved.energy += 0.5f * (p->conserved.momentum[0] * p->primitives.v[0] + p->conserved.energy += p->conserved.flux.energy * dt;
p->conserved.momentum[1] * p->primitives.v[1] +
p->conserved.momentum[2] * p->primitives.v[2]);
#endif #endif
#endif #if defined(SHADOWFAX_FIX_CELLS)
p->v[0] = 0.0f;
p->v[1] = 0.0f;
p->v[2] = 0.0f;
#else
if (p->conserved.mass > 0.0f && p->primitives.rho > 0.0f) {
/* reset fluxes */ const float inverse_mass = 1.f / p->conserved.mass;
/* we can only do this here, since we need to keep the fluxes for inactive
particles */
p->conserved.flux.mass = 0.0f;
p->conserved.flux.momentum[0] = 0.0f;
p->conserved.flux.momentum[1] = 0.0f;
p->conserved.flux.momentum[2] = 0.0f;
p->conserved.flux.energy = 0.0f;
if (p->conserved.mass > 0.) { /* Normal case: set particle velocity to fluid velocity. */
/* We want the cell velocity to be as close as possible to the fluid p->v[0] = p->conserved.momentum[0] * inverse_mass;
velocity */ p->v[1] = p->conserved.momentum[1] * inverse_mass;
vcell[0] = p->conserved.momentum[0] / p->conserved.mass; p->v[2] = p->conserved.momentum[2] * inverse_mass;
vcell[1] = p->conserved.momentum[1] / p->conserved.mass;
vcell[2] = p->conserved.momentum[2] / p->conserved.mass;
} else {
vcell[0] = 0.;
vcell[1] = 0.;
vcell[2] = 0.;
}
#ifdef SHADOWFAX_STEER_CELL_MOTION #ifdef SHADOWFAX_STEER_CELL_MOTION
/* To prevent stupid things like cell crossovers or generators that move
outside their cell, we steer the motion of the cell somewhat */
if (p->primitives.rho) {
float centroid[3], d[3]; float centroid[3], d[3];
float volume, csnd, R, vfac, fac, dnrm; float volume, csnd, R, vfac, fac, dnrm;
voronoi_get_centroid(&p->cell, centroid); voronoi_get_centroid(&p->cell, centroid);
...@@ -487,32 +493,29 @@ __attribute__((always_inline)) INLINE static void hydro_kick_extra( ...@@ -487,32 +493,29 @@ __attribute__((always_inline)) INLINE static void hydro_kick_extra(
} else { } else {
vfac = csnd / dnrm; vfac = csnd / dnrm;
} }
} else { p->v[0] += vfac * d[0];
vfac = 0.0f; p->v[1] += vfac * d[1];
} p->v[2] += vfac * d[2];
vcell[0] += vfac * d[0];
vcell[1] += vfac * d[1];
vcell[2] += vfac * d[2];
} }
#endif #endif
#if defined(SHADOWFAX_FIX_CELLS) } else {
xp->v_full[0] = 0.;
xp->v_full[1] = 0.;
xp->v_full[2] = 0.;
p->v[0] = 0.; p->v[0] = 0.;
p->v[1] = 0.; p->v[1] = 0.;
p->v[2] = 0.; p->v[2] = 0.;
#else }
xp->v_full[0] = vcell[0];
xp->v_full[1] = vcell[1];
xp->v_full[2] = vcell[2];
p->v[0] = xp->v_full[0];
p->v[1] = xp->v_full[1];
p->v[2] = xp->v_full[2];
#endif #endif
/* Now make sure all velocity variables are up to date. */
xp->v_full[0] = p->v[0];
xp->v_full[1] = p->v[1];
xp->v_full[2] = p->v[2];
if (p->gpart) {
p->gpart->v_full[0] = p->v[0];
p->gpart->v_full[1] = p->v[1];
p->gpart->v_full[2] = p->v[2];
}
} }
/** /**
... ...
......
src/hydro/Shadowswift/hydro_gradients.h
+ 1
− 54
View file @ de96fa54
...@@ -86,7 +86,7 @@ __attribute__((always_inline)) INLINE static void hydro_gradients_finalize( ...@@ -86,7 +86,7 @@ __attribute__((always_inline)) INLINE static void hydro_gradients_finalize(
*/ */
__attribute__((always_inline)) INLINE static void hydro_gradients_predict( __attribute__((always_inline)) INLINE static void hydro_gradients_predict(
struct part* pi, struct part* pj, float hi, float hj, const float* dx, struct part* pi, struct part* pj, float hi, float hj, const float* dx,
float r, float* xij_i, float* Wi, float* Wj, float mindt) { float r, float* xij_i, float* Wi, float* Wj) {
float dWi[5], dWj[5]; float dWi[5], dWj[5];
float xij_j[3]; float xij_j[3];
...@@ -132,59 +132,6 @@ __attribute__((always_inline)) INLINE static void hydro_gradients_predict( ...@@ -132,59 +132,6 @@ __attribute__((always_inline)) INLINE static void hydro_gradients_predict(
hydro_slope_limit_face(Wi, Wj, dWi, dWj, xij_i, xij_j, r); hydro_slope_limit_face(Wi, Wj, dWi, dWj, xij_i, xij_j, r);
/* time */
dWi[0] -= 0.5 * mindt * (Wi[1] * pi->primitives.gradients.rho[0] +
Wi[2] * pi->primitives.gradients.rho[1] +
Wi[3] * pi->primitives.gradients.rho[2] +
Wi[0] * (pi->primitives.gradients.v[0][0] +
pi->primitives.gradients.v[1][1] +
pi->primitives.gradients.v[2][2]));
dWi[1] -= 0.5 * mindt * (Wi[1] * pi->primitives.gradients.v[0][0] +
Wi[2] * pi->primitives.gradients.v[0][1] +
Wi[3] * pi->primitives.gradients.v[0][2] +
pi->primitives.gradients.P[0] / Wi[0]);
dWi[2] -= 0.5 * mindt * (Wi[1] * pi->primitives.gradients.v[1][0] +
Wi[2] * pi->primitives.gradients.v[1][1] +
Wi[3] * pi->primitives.gradients.v[1][2] +
pi->primitives.gradients.P[1] / Wi[0]);
dWi[3] -= 0.5 * mindt * (Wi[1] * pi->primitives.gradients.v[2][0] +
Wi[2] * pi->primitives.gradients.v[2][1] +
Wi[3] * pi->primitives.gradients.v[2][2] +
pi->primitives.gradients.P[2] / Wi[0]);
dWi[4] -=
0.5 * mindt * (Wi[1] * pi->primitives.gradients.P[0] +
Wi[2] * pi->primitives.gradients.P[1] +
Wi[3] * pi->primitives.gradients.P[2] +
hydro_gamma * Wi[4] * (pi->primitives.gradients.v[0][0] +
pi->primitives.gradients.v[1][1] +
pi->primitives.gradients.v[2][2]));
dWj[0] -= 0.5 * mindt * (Wj[1] * pj->primitives.gradients.rho[0] +
Wj[2] * pj->primitives.gradients.rho[1] +
Wj[3] * pj->primitives.gradients.rho[2] +
Wj[0] * (pj->primitives.gradients.v[0][0] +
pj->primitives.gradients.v[1][1] +
pj->primitives.gradients.v[2][2]));
dWj[1] -= 0.5 * mindt * (Wj[1] * pj->primitives.gradients.v[0][0] +
Wj[2] * pj->primitives.gradients.v[0][1] +
Wj[3] * pj->primitives.gradients.v[0][2] +
pj->primitives.gradients.P[0] / Wj[0]);
dWj[2] -= 0.5 * mindt * (Wj[1] * pj->primitives.gradients.v[1][0] +
Wj[2] * pj->primitives.gradients.v[1][1] +
Wj[3] * pj->primitives.gradients.v[1][2] +
pj->primitives.gradients.P[1] / Wj[0]);
dWj[3] -= 0.5 * mindt * (Wj[1] * pj->primitives.gradients.v[2][0] +
Wj[2] * pj->primitives.gradients.v[2][1] +
Wj[3] * pj->primitives.gradients.v[2][2] +
pj->primitives.gradients.P[2] / Wj[0]);
dWj[4] -=
0.5 * mindt * (Wj[1] * pj->primitives.gradients.P[0] +
Wj[2] * pj->primitives.gradients.P[1] +
Wj[3] * pj->primitives.gradients.P[2] +
hydro_gamma * Wj[4] * (pj->primitives.gradients.v[0][0] +
pj->primitives.gradients.v[1][1] +
pj->primitives.gradients.v[2][2]));
Wi[0] += dWi[0]; Wi[0] += dWi[0];
Wi[1] += dWi[1]; Wi[1] += dWi[1];
Wi[2] += dWi[2]; Wi[2] += dWi[2];
... ...
......
src/hydro/Shadowswift/hydro_iact.h
+ 19
− 28
View file @ de96fa54
...@@ -143,7 +143,6 @@ __attribute__((always_inline)) INLINE static void runner_iact_fluxes_common( ...@@ -143,7 +143,6 @@ __attribute__((always_inline)) INLINE static void runner_iact_fluxes_common(
float vmax, dvdotdx; float vmax, dvdotdx;
float vi[3], vj[3], vij[3]; float vi[3], vj[3], vij[3];
float Wi[5], Wj[5]; float Wi[5], Wj[5];
float dti, dtj, mindt;
float n_unit[3]; float n_unit[3];
A = voronoi_get_face(&pi->cell, pj->id, xij_i); A = voronoi_get_face(&pi->cell, pj->id, xij_i);
...@@ -168,9 +167,6 @@ __attribute__((always_inline)) INLINE static void runner_iact_fluxes_common( ...@@ -168,9 +167,6 @@ __attribute__((always_inline)) INLINE static void runner_iact_fluxes_common(
Wj[3] = pj->primitives.v[2]; Wj[3] = pj->primitives.v[2];
Wj[4] = pj->primitives.P; Wj[4] = pj->primitives.P;
dti = pi->force.dt;
dtj = pj->force.dt;
/* calculate the maximal signal velocity */ /* calculate the maximal signal velocity */
vmax = 0.0f; vmax = 0.0f;
if (Wi[0] > 0.) { if (Wi[0] > 0.) {
...@@ -192,10 +188,6 @@ __attribute__((always_inline)) INLINE static void runner_iact_fluxes_common( ...@@ -192,10 +188,6 @@ __attribute__((always_inline)) INLINE static void runner_iact_fluxes_common(
pj->timestepvars.vmax = fmaxf(pj->timestepvars.vmax, vmax); pj->timestepvars.vmax = fmaxf(pj->timestepvars.vmax, vmax);
} }
/* The flux will be exchanged using the smallest time step of the two
* particles */
mindt = fminf(dti, dtj);
/* compute the normal vector of the interface */ /* compute the normal vector of the interface */
for (k = 0; k < 3; ++k) { for (k = 0; k < 3; ++k) {
n_unit[k] = -dx[k] / r; n_unit[k] = -dx[k] / r;
...@@ -219,13 +211,12 @@ __attribute__((always_inline)) INLINE static void runner_iact_fluxes_common( ...@@ -219,13 +211,12 @@ __attribute__((always_inline)) INLINE static void runner_iact_fluxes_common(
Wj[2] -= vij[1]; Wj[2] -= vij[1];
Wj[3] -= vij[2]; Wj[3] -= vij[2];
hydro_gradients_predict(pi, pj, hi, hj, dx, r, xij_i, Wi, Wj, mindt); hydro_gradients_predict(pi, pj, hi, hj, dx, r, xij_i, Wi, Wj);
/* we don't need to rotate, we can use the unit vector in the Riemann problem /* we don't need to rotate, we can use the unit vector in the Riemann problem
* itself (see GIZMO) */ * itself (see GIZMO) */
if (Wi[0] < 0.0f || Wj[0] < 0.0f || Wi[4] < 0.0f || Wj[4] < 0.0f) { if (Wi[0] < 0.0f || Wj[0] < 0.0f || Wi[4] < 0.0f || Wj[4] < 0.0f) {
printf("mindt: %g\n", mindt);
printf("WL: %g %g %g %g %g\n", pi->primitives.rho, pi->primitives.v[0], printf("WL: %g %g %g %g %g\n", pi->primitives.rho, pi->primitives.v[0],
pi->primitives.v[1], pi->primitives.v[2], pi->primitives.P); pi->primitives.v[1], pi->primitives.v[2], pi->primitives.P);
#ifdef USE_GRADIENTS #ifdef USE_GRADIENTS
...@@ -266,20 +257,20 @@ __attribute__((always_inline)) INLINE static void runner_iact_fluxes_common( ...@@ -266,20 +257,20 @@ __attribute__((always_inline)) INLINE static void runner_iact_fluxes_common(
/* Update conserved variables */ /* Update conserved variables */
/* eqn. (16) */ /* eqn. (16) */
pi->conserved.flux.mass -= mindt * A * totflux[0]; pi->conserved.flux.mass -= A * totflux[0];
pi->conserved.flux.momentum[0] -= mindt * A * totflux[1]; pi->conserved.flux.momentum[0] -= A * totflux[1];
pi->conserved.flux.momentum[1] -= mindt * A * totflux[2]; pi->conserved.flux.momentum[1] -= A * totflux[2];
pi->conserved.flux.momentum[2] -= mindt * A * totflux[3]; pi->conserved.flux.momentum[2] -= A * totflux[3];
pi->conserved.flux.energy -= mindt * A * totflux[4]; pi->conserved.flux.energy -= A * totflux[4];
#ifndef SHADOWFAX_TOTAL_ENERGY #ifndef SHADOWFAX_TOTAL_ENERGY
float ekin = 0.5f * (pi->primitives.v[0] * pi->primitives.v[0] + float ekin = 0.5f * (pi->primitives.v[0] * pi->primitives.v[0] +
pi->primitives.v[1] * pi->primitives.v[1] + pi->primitives.v[1] * pi->primitives.v[1] +
pi->primitives.v[2] * pi->primitives.v[2]); pi->primitives.v[2] * pi->primitives.v[2]);
pi->conserved.flux.energy += mindt * A * totflux[1] * pi->primitives.v[0]; pi->conserved.flux.energy += A * totflux[1] * pi->primitives.v[0];
pi->conserved.flux.energy += mindt * A * totflux[2] * pi->primitives.v[1]; pi->conserved.flux.energy += A * totflux[2] * pi->primitives.v[1];
pi->conserved.flux.energy += mindt * A * totflux[3] * pi->primitives.v[2]; pi->conserved.flux.energy += A * totflux[3] * pi->primitives.v[2];
pi->conserved.flux.energy -= mindt * A * totflux[0] * ekin; pi->conserved.flux.energy -= A * totflux[0] * ekin;
#endif #endif
/* here is how it works: /* here is how it works:
...@@ -295,20 +286,20 @@ __attribute__((always_inline)) INLINE static void runner_iact_fluxes_common( ...@@ -295,20 +286,20 @@ __attribute__((always_inline)) INLINE static void runner_iact_fluxes_common(
==> we update particle j if (MODE IS 1) OR (j IS INACTIVE) ==> we update particle j if (MODE IS 1) OR (j IS INACTIVE)
*/ */
if (mode == 1 || pj->force.active == 0) { if (mode == 1 || pj->force.active == 0) {
pj->conserved.flux.mass += mindt * A * totflux[0]; pj->conserved.flux.mass += A * totflux[0];
pj->conserved.flux.momentum[0] += mindt * A * totflux[1]; pj->conserved.flux.momentum[0] += A * totflux[1];
pj->conserved.flux.momentum[1] += mindt * A * totflux[2]; pj->conserved.flux.momentum[1] += A * totflux[2];
pj->conserved.flux.momentum[2] += mindt * A * totflux[3]; pj->conserved.flux.momentum[2] += A * totflux[3];
pj->conserved.flux.energy += mindt * A * totflux[4]; pj->conserved.flux.energy += A * totflux[4];
#ifndef SHADOWFAX_TOTAL_ENERGY #ifndef SHADOWFAX_TOTAL_ENERGY
ekin = 0.5f * (pj->primitives.v[0] * pj->primitives.v[0] + ekin = 0.5f * (pj->primitives.v[0] * pj->primitives.v[0] +
pj->primitives.v[1] * pj->primitives.v[1] + pj->primitives.v[1] * pj->primitives.v[1] +
pj->primitives.v[2] * pj->primitives.v[2]); pj->primitives.v[2] * pj->primitives.v[2]);
pj->conserved.flux.energy -= mindt * A * totflux[1] * pj->primitives.v[0]; pj->conserved.flux.energy -= A * totflux[1] * pj->primitives.v[0];
pj->conserved.flux.energy -= mindt * A * totflux[2] * pj->primitives.v[1]; pj->conserved.flux.energy -= A * totflux[2] * pj->primitives.v[1];
pj->conserved.flux.energy -= mindt * A * totflux[3] * pj->primitives.v[2]; pj->conserved.flux.energy -= A * totflux[3] * pj->primitives.v[2];
pj->conserved.flux.energy += mindt * A * totflux[0] * ekin; pj->conserved.flux.energy += A * totflux[0] * ekin;
#endif #endif
} }
} }
... ...
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please to comment