Commit 765a58c6 authored by Matthieu Schaller's avatar Matthieu Schaller
Browse files

Potential forces don't contain Newton's constant any more.

parent af3f90f0
......@@ -725,7 +725,8 @@ int cell_are_neighbours(const struct cell *restrict ci,
#endif
/* Maximum allowed distance */
const double min_dist = 1.2 * ci->width[0]; /* 1.2 accounts for rounding errors */
const double min_dist =
1.2 * ci->width[0]; /* 1.2 accounts for rounding errors */
/* (Manhattan) Distance between the cells */
for (int k = 0; k < 3; k++) {
......
......@@ -47,7 +47,7 @@ void potential_init(const struct swift_params* parameter_file,
potential->point_mass.mass =
parser_get_param_double(parameter_file, "PointMass:mass");
potential->point_mass.timestep_mult =
parser_get_param_double(parameter_file, "PointMass:timestep_mult");
parser_get_param_float(parameter_file, "PointMass:timestep_mult");
#endif /* EXTERNAL_POTENTIAL_POINTMASS */
......@@ -61,7 +61,7 @@ void potential_init(const struct swift_params* parameter_file,
parser_get_param_double(parameter_file, "IsothermalPotential:position_z");
potential->isothermal_potential.vrot =
parser_get_param_double(parameter_file, "IsothermalPotential:vrot");
potential->isothermal_potential.timestep_mult = parser_get_param_double(
potential->isothermal_potential.timestep_mult = parser_get_param_float(
parameter_file, "IsothermalPotential:timestep_mult");
#endif /* EXTERNAL_POTENTIAL_ISOTHERMALPOTENTIAL */
......
......@@ -42,7 +42,7 @@ struct external_potential {
struct {
double x, y, z;
double mass;
double timestep_mult;
float timestep_mult;
} point_mass;
#endif
......@@ -50,7 +50,7 @@ struct external_potential {
struct {
double x, y, z;
double vrot;
double timestep_mult;
float timestep_mult;
} isothermal_potential;
#endif
};
......@@ -95,6 +95,8 @@ external_gravity_isothermalpotential_timestep(
* @brief Computes the gravitational acceleration of a particle due to a point
* mass
*
* Note that the accelerations are multiplied by Newton's G constant later on.
*
* @param potential The #external_potential used in the run.
* @param phys_const The physical constants in internal units.
* @param g Pointer to the g-particle data.
......@@ -104,15 +106,19 @@ external_gravity_isothermalpotential(const struct external_potential* potential,
const struct phys_const* const phys_const,
struct gpart* g) {
const float G_newton = phys_const->const_newton_G;
const float dx = g->x[0] - potential->isothermal_potential.x;
const float dy = g->x[1] - potential->isothermal_potential.y;
const float dz = g->x[2] - potential->isothermal_potential.z;
const float rinv2 = 1.f / (dx * dx + dy * dy + dz * dz);
const double vrot = potential->isothermal_potential.vrot;
g->a_grav[0] += -vrot * vrot * rinv2 * dx;
g->a_grav[1] += -vrot * vrot * rinv2 * dy;
g->a_grav[2] += -vrot * vrot * rinv2 * dz;
const double term = -vrot * vrot * rinv2 / G_newton;
g->a_grav[0] += term * dx;
g->a_grav[1] += term * dy;
g->a_grav[2] += term * dz;
// error(" %f %f %f %f", vrot, rinv2, dx, g->a_grav[0]);
}
......
......@@ -145,7 +145,8 @@ __attribute__((always_inline)) INLINE static void runner_dopair_grav_pp(
#ifdef SWIFT_DEBUG_CHECKS
if (ci->width[0] != cj->width[0]) // MATTHIEU sanity check
error("Non matching cell sizes !! h_i=%f h_j=%f", ci->width[0], cj->width[0]);
error("Non matching cell sizes !! h_i=%f h_j=%f", ci->width[0],
cj->width[0]);
#endif
/* Anything to do here? */
......@@ -253,8 +254,8 @@ __attribute__((always_inline)) INLINE static void runner_doself_grav_pp(
if (gp->id_or_neg_offset == ICHECK)
message("id=%lld loc=[ %f %f %f ] size= %f count= %d",
gp->id_or_neg_offset, c->loc[0], c->loc[1], c->loc[2], c->width[0],
c->gcount);
gp->id_or_neg_offset, c->loc[0], c->loc[1], c->loc[2],
c->width[0], c->gcount);
}
#endif
......@@ -325,7 +326,8 @@ static void runner_dopair_grav(struct runner *r, struct cell *ci,
/* Bad stuff will happen if cell sizes are different */
if (ci->width[0] != cj->width[0])
error("Non matching cell sizes !! h_i=%f h_j=%f", ci->width[0], cj->width[0]);
error("Non matching cell sizes !! h_i=%f h_j=%f", ci->width[0],
cj->width[0]);
/* Sanity check */
if (ci == cj)
......@@ -336,10 +338,11 @@ static void runner_dopair_grav(struct runner *r, struct cell *ci,
/* Are the cells direct neighbours? */
if (!cell_are_neighbours(ci, cj))
error(
"Non-neighbouring cells ! ci->x=[%f %f %f] ci->width=%f cj->loc=[%f %f %f] "
"Non-neighbouring cells ! ci->x=[%f %f %f] ci->width=%f cj->loc=[%f %f "
"%f] "
"cj->width=%f",
ci->loc[0], ci->loc[1], ci->loc[2], ci->width[0], cj->loc[0], cj->loc[1],
cj->loc[2], cj->width[0]);
ci->loc[0], ci->loc[1], ci->loc[2], ci->width[0], cj->loc[0],
cj->loc[1], cj->loc[2], cj->width[0]);
#endif
......@@ -488,7 +491,8 @@ static void runner_do_grav_mm(struct runner *r, struct cell *ci, int timer) {
struct cell *cells = e->s->cells;
const int nr_cells = e->s->nr_cells;
const int ti_current = e->ti_current;
const double max_d = const_gravity_a_smooth * const_gravity_r_cut * ci->width[0];
const double max_d =
const_gravity_a_smooth * const_gravity_r_cut * ci->width[0];
const double max_d2 = max_d * max_d;
const double pos_i[3] = {ci->loc[0], ci->loc[1], ci->loc[2]};
......
......@@ -134,14 +134,14 @@ int main() {
memcpy(forcegrid, workspace + (sendmin - meshmin[0]) * dimy * dimz,
(sendmax - sendmin + 1) * dimy * dimz * sizeof(double));
/* for (size_t i = 0; i < fft_size; ++i) */
/* if (forcegrid[i] != workspace[i]) error("wrong"); */
/* Prepare the density grid */
double* rhogrid = fftw_malloc(fft_size * sizeof(double));
bzero(rhogrid, fft_size * sizeof(double));
/* Now get the density */
for (size_t slab_x = recvmin; slab_x <= recvmax; slab_x++) {
......@@ -169,18 +169,16 @@ int main() {
/* } */
/* } */
/* FFT of the density field */
fftw_complex* fftgrid = fftw_malloc(fft_size * sizeof(fftw_complex));
fftw_plan plan_forward = fftw_plan_dft_r2c_3d(PMGRID, PMGRID, PMGRID, rhogrid, fftgrid, FFTW_ESTIMATE);
fftw_plan plan_forward = fftw_plan_dft_r2c_3d(PMGRID, PMGRID, PMGRID, rhogrid,
fftgrid, FFTW_ESTIMATE);
fftw_execute(plan_forward);
for(size_t i = 0; i < 640; i++) {
for (size_t i = 0; i < 640; i++) {
message("workspace[%zd]= %f", i, fftgrid[i][0]);
}
/* Clean-up */
fftw_destroy_plan(plan_forward);
fftw_free(forcegrid);
......
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