Skip to content
Snippets Groups Projects
Commit 7394924f authored by Matthieu Schaller's avatar Matthieu Schaller
Browse files

Updated unit tests for the new definitions.

parent 4f5948af
No related branches found
No related tags found
1 merge request!516More cosmology work
......@@ -3,8 +3,9 @@ Scheduler:
# Parameters for the self-gravity scheme
Gravity:
eta: 0.025 # Constant dimensionless multiplier for time integration.
theta: 0.7 # Opening angle (Multipole acceptance criterion)
epsilon: 0.00001 # Softening length (in internal units).
a_smooth: 0.
r_cut: 0.
eta: 0.025 # Constant dimensionless multiplier for time integration.
theta: 0.7 # Opening angle (Multipole acceptance criterion)
comoving_softening: 0.00001 # Comoving softening length (in internal units).
max_physical_softening: 0.00001 # Physical softening length (in internal units).
a_smooth: 0.
r_cut: 0.
......@@ -57,15 +57,18 @@ int main() {
struct swift_params *params = malloc(sizeof(struct swift_params));
parser_read_file("fft_params.yml", params);
struct cosmology cosmo;
cosmology_init_no_cosmo(&cosmo);
/* Initialise the gravity properties */
struct gravity_props gravity_properties;
gravity_props_init(&gravity_properties, params);
gravity_props_init(&gravity_properties, params, &cosmo);
/* Build the infrastructure */
struct space space;
double dim[3] = {1., 1., 1.};
space_init(&space, params, dim, NULL, gparts, NULL, 0, nr_gparts, 0, 1, 1, 1,
0, 0);
space_init(&space, params, &cosmo, dim, NULL, gparts, NULL, 0, nr_gparts, 0,
1, 1, 0, 1, 0, 0);
struct engine engine;
engine.s = &space;
......
......@@ -106,6 +106,7 @@ int main() {
props.a_smooth = 1.25;
props.r_cut_min = 0.;
props.theta_crit2 = 0.;
props.epsilon_cur = eps;
e.gravity_properties = &props;
struct runner r;
......@@ -178,7 +179,6 @@ int main() {
gp->x[1] = 0.5;
gp->x[2] = 0.5;
gp->mass = 0.;
gp->epsilon = eps;
gp->time_bin = 1;
gp->type = swift_type_dark_matter;
gp->id_or_neg_offset = n + 1;
......@@ -196,7 +196,6 @@ int main() {
ci.gparts[0].x[1] = 0.5;
ci.gparts[0].x[2] = 0.5;
ci.gparts[0].mass = 1.;
ci.gparts[0].epsilon = eps;
ci.gparts[0].time_bin = 1;
ci.gparts[0].type = swift_type_dark_matter;
ci.gparts[0].id_or_neg_offset = 1;
......@@ -211,11 +210,12 @@ int main() {
for (int n = 0; n < num_tests; ++n) {
const struct gpart *gp = &cj.gparts[n];
const struct gpart *gp2 = &ci.gparts[0];
const double epsilon = gravity_get_softening(gp, &props);
double pot_true =
potential(ci.gparts[0].mass, gp->x[0] - gp2->x[0], gp->epsilon, rlr);
potential(ci.gparts[0].mass, gp->x[0] - gp2->x[0], epsilon, rlr);
double acc_true =
acceleration(ci.gparts[0].mass, gp->x[0] - gp2->x[0], gp->epsilon, rlr);
acceleration(ci.gparts[0].mass, gp->x[0] - gp2->x[0], epsilon, rlr);
message("x=%e f=%e f_true=%e pot=%e pot_true=%e", gp->x[0] - gp2->x[0],
gp->a_grav[0], acc_true, gp->potential, pot_true);
......@@ -252,12 +252,12 @@ int main() {
for (int n = 0; n < num_tests; ++n) {
const struct gpart *gp = &cj.gparts[n];
const struct gravity_tensors *mpole = ci.multipole;
const double epsilon = gravity_get_softening(gp, &props);
double pot_true = potential(mpole->m_pole.M_000, gp->x[0] - mpole->CoM[0],
gp->epsilon, rlr * FLT_MAX);
double acc_true =
acceleration(mpole->m_pole.M_000, gp->x[0] - mpole->CoM[0], gp->epsilon,
rlr * FLT_MAX);
epsilon, rlr * FLT_MAX);
double acc_true = acceleration(
mpole->m_pole.M_000, gp->x[0] - mpole->CoM[0], epsilon, rlr * FLT_MAX);
message("x=%e f=%e f_true=%e pot=%e pot_true=%e", gp->x[0] - mpole->CoM[0],
gp->a_grav[0], acc_true, gp->potential, pot_true);
......@@ -304,7 +304,6 @@ int main() {
ci.gparts[n].mass = 1. / 8.;
ci.gparts[n].epsilon = eps;
ci.gparts[n].time_bin = 1;
ci.gparts[n].type = swift_type_dark_matter;
ci.gparts[n].id_or_neg_offset = 1;
......@@ -333,18 +332,19 @@ int main() {
for (int i = 0; i < 8; ++i) {
const struct gpart *gp2 = &ci.gparts[i];
const double epsilon = gravity_get_softening(gp, &props);
const double dx[3] = {gp2->x[0] - gp->x[0], gp2->x[1] - gp->x[1],
gp2->x[2] - gp->x[2]};
const double d = sqrt(dx[0] * dx[0] + dx[1] * dx[1] + dx[2] * dx[2]);
pot_true += potential(gp2->mass, d, gp->epsilon, rlr * FLT_MAX);
pot_true += potential(gp2->mass, d, epsilon, rlr * FLT_MAX);
acc_true[0] -=
acceleration(gp2->mass, d, gp->epsilon, rlr * FLT_MAX) * dx[0] / d;
acceleration(gp2->mass, d, epsilon, rlr * FLT_MAX) * dx[0] / d;
acc_true[1] -=
acceleration(gp2->mass, d, gp->epsilon, rlr * FLT_MAX) * dx[1] / d;
acceleration(gp2->mass, d, epsilon, rlr * FLT_MAX) * dx[1] / d;
acc_true[2] -=
acceleration(gp2->mass, d, gp->epsilon, rlr * FLT_MAX) * dx[2] / d;
acceleration(gp2->mass, d, epsilon, rlr * FLT_MAX) * dx[2] / d;
}
message("x=%e f=%e f_true=%e pot=%e pot_true=%e %e %e",
......
......@@ -103,6 +103,7 @@ int main() {
struct gravity_props props;
props.a_smooth = 1.25;
props.epsilon_cur = eps;
e.gravity_properties = &props;
struct runner r;
......@@ -139,7 +140,6 @@ int main() {
c.gparts[0].x[1] = 0.5;
c.gparts[0].x[2] = 0.5;
c.gparts[0].mass = 1.;
c.gparts[0].epsilon = eps;
c.gparts[0].time_bin = 1;
c.gparts[0].type = swift_type_dark_matter;
c.gparts[0].id_or_neg_offset = 1;
......@@ -156,7 +156,6 @@ int main() {
gp->x[1] = 0.5;
gp->x[2] = 0.5;
gp->mass = 0.;
gp->epsilon = eps;
gp->time_bin = 1;
gp->type = swift_type_dark_matter;
gp->id_or_neg_offset = n + 1;
......@@ -172,9 +171,10 @@ int main() {
for (int n = 1; n < num_tests + 1; ++n) {
const struct gpart *gp = &c.gparts[n];
double pot_true = potential(c.gparts[0].mass, gp->x[0], gp->epsilon, rlr);
double acc_true =
acceleration(c.gparts[0].mass, gp->x[0], gp->epsilon, rlr);
const double epsilon = gravity_get_softening(gp, &props);
double pot_true = potential(c.gparts[0].mass, gp->x[0], epsilon, rlr);
double acc_true = acceleration(c.gparts[0].mass, gp->x[0], epsilon, rlr);
// message("x=%e f=%e f_true=%e", gp->x[0], gp->a_grav[0], acc_true);
......
......@@ -48,7 +48,8 @@ int main() {
/* Read data */
read_ic_single("input.hdf5", &us, dim, &parts, &gparts, &sparts, &Ngas,
&Ngpart, &Nspart, &periodic, &flag_entropy_ICs, 1, 1, 0, 1, 0);
&Ngpart, &Nspart, &periodic, &flag_entropy_ICs, 1, 1, 0, 0, 1.,
1, 0);
/* Check global properties read are correct */
assert(dim[0] == boxSize);
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment