Commit 1aeb605a authored by Matthieu Schaller's avatar Matthieu Schaller
Browse files

Tests now also succesful for DEFAULT_SPH

parent 3fef1086
......@@ -56,8 +56,8 @@
/* SPH variant to use */
//#define MINIMAL_SPH
#define GADGET2_SPH
//#define DEFAULT_SPH
//#define GADGET2_SPH
#define DEFAULT_SPH
/* External potential properties */
#define EXTERNAL_POTENTIAL_POINTMASS
......
......@@ -20,6 +20,8 @@
#include "adiabatic_index.h"
#include "approx_math.h"
#include <float.h>
/**
* @brief Computes the hydro time-step of a given particle
*
......@@ -145,7 +147,7 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force(
/* Compute this particle's sound speed. */
const float u = p->u;
const float fc = p->force.c = sqrtf(hydro_gamma * hydro_gamma_minus_one * u);
const float fc = p->force.soundspeed = sqrtf(hydro_gamma * hydro_gamma_minus_one * u);
/* Compute the P/Omega/rho2. */
xp->omega = 1.0f + 0.3333333333f * h * p->rho_dh / p->rho;
......@@ -155,7 +157,7 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force(
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);
const float tau = h / (2.f * const_viscosity_length * p->force.soundspeed);
/* Viscosity source term */
const float S = fmaxf(-normDiv_v, 0.f);
......@@ -262,3 +264,39 @@ __attribute__((always_inline)) INLINE static float hydro_get_internal_energy(
return p->u;
}
/**
* @brief Returns the pressure of a particle
*
* @param p The particle of interest
* @param dt Time since the last kick
*/
__attribute__((always_inline)) INLINE static float hydro_get_pressure(
const struct part *restrict p, float dt) {
return p->force.POrho2 * p->rho * p->rho / p->rho_dh;
}
/**
* @brief Returns the entropy of a particle
*
* @param p The particle of interest
* @param dt Time since the last kick
*/
__attribute__((always_inline)) INLINE static float hydro_get_entropy(
const struct part *restrict p, float dt) {
return hydro_gamma_minus_one * p->u * pow_minus_gamma_minus_one(p->rho);
}
/**
* @brief Returns the sound speed of a particle
*
* @param p The particle of interest
* @param dt Time since the last kick
*/
__attribute__((always_inline)) INLINE static float hydro_get_soundspeed(
const struct part *restrict p, float dt) {
return p->force.soundspeed;
}
......@@ -398,7 +398,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_force(
omega_ij = fminf(dvdr, 0.f);
/* Compute signal velocity */
v_sig = pi->force.c + pj->force.c - 2.0f * omega_ij;
v_sig = pi->force.soundspeed + pj->force.soundspeed - 2.0f * omega_ij;
/* Compute viscosity parameter */
alpha_ij = -0.5f * (pi->alpha + pj->alpha);
......@@ -496,11 +496,11 @@ __attribute__((always_inline)) INLINE static void runner_iact_vec_force(
pju.v = vec_set(pj[0]->u, pj[1]->u, pj[2]->u, pj[3]->u, pj[4]->u, pj[5]->u,
pj[6]->u, pj[7]->u);
ci.v =
vec_set(pi[0]->force.c, pi[1]->force.c, pi[2]->force.c, pi[3]->force.c,
pi[4]->force.c, pi[5]->force.c, pi[6]->force.c, pi[7]->force.c);
vec_set(pi[0]->force.soundspeed, pi[1]->force.soundspeed, pi[2]->force.soundspeed, pi[3]->force.soundspeed,
pi[4]->force.soundspeed, pi[5]->force.soundspeed, pi[6]->force.soundspeed, pi[7]->force.soundspeed);
cj.v =
vec_set(pj[0]->force.c, pj[1]->force.c, pj[2]->force.c, pj[3]->force.c,
pj[4]->force.c, pj[5]->force.c, pj[6]->force.c, pj[7]->force.c);
vec_set(pj[0]->force.soundspeed, pj[1]->force.soundspeed, pj[2]->force.soundspeed, pj[3]->force.soundspeed,
pj[4]->force.soundspeed, pj[5]->force.soundspeed, pj[6]->force.soundspeed, pj[7]->force.soundspeed);
vi_sig.v = vec_set(pi[0]->force.v_sig, pi[1]->force.v_sig, pi[2]->force.v_sig,
pi[3]->force.v_sig, pi[4]->force.v_sig, pi[5]->force.v_sig,
pi[6]->force.v_sig, pi[7]->force.v_sig);
......@@ -539,9 +539,9 @@ __attribute__((always_inline)) INLINE static void runner_iact_vec_force(
piu.v = vec_set(pi[0]->u, pi[1]->u, pi[2]->u, pi[3]->u);
pju.v = vec_set(pj[0]->u, pj[1]->u, pj[2]->u, pj[3]->u);
ci.v =
vec_set(pi[0]->force.c, pi[1]->force.c, pi[2]->force.c, pi[3]->force.c);
vec_set(pi[0]->force.soundspeed, pi[1]->force.soundspeed, pi[2]->force.soundspeed, pi[3]->force.soundspeed);
cj.v =
vec_set(pj[0]->force.c, pj[1]->force.c, pj[2]->force.c, pj[3]->force.c);
vec_set(pj[0]->force.soundspeed, pj[1]->force.soundspeed, pj[2]->force.soundspeed, pj[3]->force.soundspeed);
vi_sig.v = vec_set(pi[0]->force.v_sig, pi[1]->force.v_sig, pi[2]->force.v_sig,
pi[3]->force.v_sig);
vj_sig.v = vec_set(pj[0]->force.v_sig, pj[1]->force.v_sig, pj[2]->force.v_sig,
......@@ -711,7 +711,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
omega_ij = fminf(dvdr, 0.f);
/* Compute signal velocity */
v_sig = pi->force.c + pj->force.c - 2.0f * omega_ij;
v_sig = pi->force.soundspeed + pj->force.soundspeed - 2.0f * omega_ij;
/* Compute viscosity parameter */
alpha_ij = -0.5f * (pi->alpha + pj->alpha);
......@@ -801,11 +801,11 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_vec_force(
pju.v = vec_set(pj[0]->u, pj[1]->u, pj[2]->u, pj[3]->u, pj[4]->u, pj[5]->u,
pj[6]->u, pj[7]->u);
ci.v =
vec_set(pi[0]->force.c, pi[1]->force.c, pi[2]->force.c, pi[3]->force.c,
pi[4]->force.c, pi[5]->force.c, pi[6]->force.c, pi[7]->force.c);
vec_set(pi[0]->force.soundspeed, pi[1]->force.soundspeed, pi[2]->force.soundspeed, pi[3]->force.soundspeed,
pi[4]->force.soundspeed, pi[5]->force.soundspeed, pi[6]->force.soundspeed, pi[7]->force.soundspeed);
cj.v =
vec_set(pj[0]->force.c, pj[1]->force.c, pj[2]->force.c, pj[3]->force.c,
pj[4]->force.c, pj[5]->force.c, pj[6]->force.c, pj[7]->force.c);
vec_set(pj[0]->force.soundspeed, pj[1]->force.soundspeed, pj[2]->force.soundspeed, pj[3]->force.soundspeed,
pj[4]->force.soundspeed, pj[5]->force.soundspeed, pj[6]->force.soundspeed, pj[7]->force.soundspeed);
vi_sig.v = vec_set(pi[0]->force.v_sig, pi[1]->force.v_sig, pi[2]->force.v_sig,
pi[3]->force.v_sig, pi[4]->force.v_sig, pi[5]->force.v_sig,
pi[6]->force.v_sig, pi[7]->force.v_sig);
......@@ -843,9 +843,9 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_vec_force(
piu.v = vec_set(pi[0]->u, pi[1]->u, pi[2]->u, pi[3]->u);
pju.v = vec_set(pj[0]->u, pj[1]->u, pj[2]->u, pj[3]->u);
ci.v =
vec_set(pi[0]->force.c, pi[1]->force.c, pi[2]->force.c, pi[3]->force.c);
vec_set(pi[0]->force.soundspeed, pi[1]->force.soundspeed, pi[2]->force.soundspeed, pi[3]->force.soundspeed);
cj.v =
vec_set(pj[0]->force.c, pj[1]->force.c, pj[2]->force.c, pj[3]->force.c);
vec_set(pj[0]->force.soundspeed, pj[1]->force.soundspeed, pj[2]->force.soundspeed, pj[3]->force.soundspeed);
vi_sig.v = vec_set(pi[0]->force.v_sig, pi[1]->force.v_sig, pi[2]->force.v_sig,
pi[3]->force.v_sig);
vj_sig.v = vec_set(pj[0]->force.v_sig, pj[1]->force.v_sig, pj[2]->force.v_sig,
......
......@@ -99,7 +99,7 @@ struct part {
float v_sig;
/* Sound speed */
float c;
float soundspeed;
/* Change in smoothing length over time. */
float h_dt;
......
......@@ -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(DEFAULT_SPH)
part->u = pressure / (hydro_gamma_minus_one * density);
#else
error("Need to define pressure here !");
#endif
......@@ -188,6 +190,8 @@ void get_solution(const struct cell *main_cell, struct solution_part *solution,
* separation.
* @param density The density of the fluid.
* @param partId The running counter of IDs.
* @param vel The type of velocity field.
* @param press The type of pressure field.
*/
struct cell *make_cell(size_t n, const double offset[3], double size, double h,
double density, long long *partId,
......@@ -302,7 +306,14 @@ void dump_particle_fields(char *fileName, struct cell *main_cell,
main_cell->parts[pid].a_hydro[0], main_cell->parts[pid].a_hydro[1],
main_cell->parts[pid].a_hydro[2], main_cell->parts[pid].force.h_dt,
main_cell->parts[pid].force.v_sig,
main_cell->parts[pid].entropy_dt, 0.f);
#if defined(GADGET2_SPH)
main_cell->parts[pid].entropy_dt, 0.f
#elif defined(DEFAULT_SPH)
0.f, main_cell->parts[pid].force.u_dt
#else
0.f, 0.f
#endif
);
}
if (with_solution) {
......
......@@ -79,16 +79,8 @@ void dump_indv_particle_fields(char *fileName, struct part *p) {
"%13e %13e %13e\n",
p->id, p->x[0], p->x[1], p->x[2], p->v[0], p->v[1], p->v[2], p->rho,
p->rho_dh, p->density.wcount, p->density.wcount_dh,
#if defined(GADGET2_SPH)
p->density.div_v, p->density.rot_v[0], p->density.rot_v[1],
p->density.rot_v[2]
#elif defined(DEFAULT_SPH)
p->density.div_v, p->density.curl_v[0], p->density.curl_v[1],
p->density.curl_v[2]
#else
0., 0., 0., 0.
#endif
);
p->density.rot_v[2]);
fclose(file);
}
......
......@@ -61,9 +61,9 @@ struct cell *make_cell(size_t N, float cellSize, int offset[3], int id_offset) {
offset[2] * cellSize + z * cellSize / N + cellSize / (2 * N);
part->h = h;
part->id = x * N * N + y * N + z + id_offset;
++part;
part->ti_begin = 0;
part->ti_end = 1;
++part;
}
}
}
......@@ -193,16 +193,20 @@ int main() {
runner_do_ghost(&r, ci);
message("h=%f rho=%f N_ngb=%f", p->h, p->rho, p->density.wcount);
message("c=%f", p->force.c);
message("c=%f", p->force.soundspeed);
runner_doself2_force(&r, ci);
runner_do_kick(&r, ci, 1);
message("ti_end=%d", p->ti_end);
free(ci->parts);
free(ci->xparts);
for (int j = 0; j < 27; ++j) {
free(cells[j]->parts);
free(cells[j]->xparts);
free(cells[j]->sort);
free(cells[j]);
}
return 0;
#endif
}
......@@ -81,9 +81,9 @@ int main(int argc, char *argv[]) {
p2.rho = 1.0f;
p2.mass = 9.7059e-4;
p2.h = 0.222871287 / 2;
p1.force.c = 0.0040824829f;
p1.force.soundspeed = 0.0040824829f;
p1.force.balsara = 0.0f;
p2.force.c = 58.8972740361f;
p2.force.soundspeed = 58.8972740361f;
p2.force.balsara = 0.0f;
p1.u = 1.e-5 / (hydro_gamma_minus_one * p1.rho);
p2.u = 1.e-5 / (hydro_gamma_minus_one * p2.rho) + 100.0f / (33 * p2.mass);
......
......@@ -76,16 +76,8 @@ void dump_indv_particle_fields(char *fileName, struct part *p) {
"%13e %13e %13e\n",
p->id, p->x[0], p->x[1], p->x[2], p->v[0], p->v[1], p->v[2], p->rho,
p->rho_dh, p->density.wcount, p->density.wcount_dh,
#if defined(GADGET2_SPH)
p->density.div_v, p->density.rot_v[0], p->density.rot_v[1],
p->density.rot_v[2]
#elif defined(DEFAULT_SPH)
p->density.div_v, p->density.curl_v[0], p->density.curl_v[1],
p->density.curl_v[2]
#else
0., 0., 0., 0.
#endif
);
p->density.rot_v[2]);
fclose(file);
}
......
Supports Markdown
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