diff --git a/examples/ExternalGravity/externalGravity.yml b/examples/ExternalGravity/externalGravity.yml index 857bba0068df71431f77e91683e29731a22d9e99..959d7c1c2885a60ac3c589ecb0084f04d18e9cba 100644 --- a/examples/ExternalGravity/externalGravity.yml +++ b/examples/ExternalGravity/externalGravity.yml @@ -1,7 +1,7 @@ # Define the system of units to use internally. UnitSystem: - UnitMass_in_cgs: 1.898e33 # Grams + UnitMass_in_cgs: 1.989e33 # Grams UnitLength_in_cgs: 3.0856776e21 # Centimeters UnitVelocity_in_cgs: 1e5 # Centimeters per second UnitCurrent_in_cgs: 1 # Amperes @@ -46,3 +46,11 @@ DomainDecomposition: initial_grid_z: 10 repartition_type: b # The re-decomposition strategy ("n", "b", "v", "e" or "x"). See documentation for details. + +# External potential parameters +PointMass: + position_x: 50. # location of external point mass in internal units + position_y: 50. + position_z: 50. + mass: 1e10 # mass of external point mass in internal units + diff --git a/examples/ExternalGravity/makeIC.py b/examples/ExternalGravity/makeIC.py index 1bb4c69b75d5d16a943f9770c41b4f7980f68ee7..e82d970eba30f11ce509367fb9f1ab9dc82af50d 100644 --- a/examples/ExternalGravity/makeIC.py +++ b/examples/ExternalGravity/makeIC.py @@ -38,6 +38,10 @@ const_unit_length_in_cgs = (1000*PARSEC_IN_CGS) const_unit_mass_in_cgs = (SOLAR_MASS_IN_CGS) const_unit_velocity_in_cgs = (1e5) +print "UnitMass_in_cgs: ", const_unit_mass_in_cgs +print "UnitLength_in_cgs: ", const_unit_length_in_cgs +print "UnitVelocity_in_cgs: ", const_unit_velocity_in_cgs + # derived units const_unit_time_in_cgs = (const_unit_length_in_cgs / const_unit_velocity_in_cgs) const_G = ((NEWTON_GRAVITY_CGS*const_unit_mass_in_cgs*const_unit_time_in_cgs*const_unit_time_in_cgs/(const_unit_length_in_cgs*const_unit_length_in_cgs*const_unit_length_in_cgs))) diff --git a/examples/main.c b/examples/main.c index 3e4101c2d14e0f73ced84415f251dfb102194f6b..4631a09aa115285ecfe56a9f819d3fc7baf59095 100644 --- a/examples/main.c +++ b/examples/main.c @@ -264,7 +264,7 @@ int main(int argc, char *argv[]) { struct external_potential potential; units_init(&us, params); initPhysicalConstants(&us, &prog_const); - initPotentialProperties(&us, &potential); + initPotentialProperties(params, &us, &potential); if (myrank == 0) { message("Unit system: U_M = %e g.", us.UnitMass_in_cgs); message("Unit system: U_L = %e cm.", us.UnitLength_in_cgs); diff --git a/src/potentials.c b/src/potentials.c index 856bdac3dcac00b060f928b667787bfe7ba8fd31..772289e7b599276b4922ce12ec64629cfe8b9cc5 100644 --- a/src/potentials.c +++ b/src/potentials.c @@ -31,17 +31,23 @@ * @param us The current internal system of units * @param potential The external potential properties to initialize */ -void initPotentialProperties(struct UnitSystem* us, +void initPotentialProperties(const struct swift_params * parameter_file, + struct UnitSystem* us, struct external_potential* potential) { message(" %e\t %e", PARSEC_IN_CGS, units_conversion_factor(us, UNIT_CONV_LENGTH)); + + potential->point_mass.x = parser_get_param_double(parameter_file, "PointMass:position_x"); + potential->point_mass.y = parser_get_param_double(parameter_file, "PointMass:position_y"); + potential->point_mass.z = parser_get_param_double(parameter_file, "PointMass:position_z"); + potential->point_mass.mass = parser_get_param_double(parameter_file, "PointMass:mass"); - potential->point_mass.x = - 50000 * PARSEC_IN_CGS / units_conversion_factor(us, UNIT_CONV_LENGTH); - potential->point_mass.y = - 50000 * PARSEC_IN_CGS / units_conversion_factor(us, UNIT_CONV_LENGTH); - potential->point_mass.z = - 50000 * PARSEC_IN_CGS / units_conversion_factor(us, UNIT_CONV_LENGTH); - potential->point_mass.mass = - 1e10 * SOLAR_MASS_IN_CGS / units_conversion_factor(us, UNIT_CONV_MASS); + /* potential->point_mass.x = */ + /* 50000 * PARSEC_IN_CGS / units_conversion_factor(us, UNIT_CONV_LENGTH); */ + /* potential->point_mass.y = */ + /* 50000 * PARSEC_IN_CGS / units_conversion_factor(us, UNIT_CONV_LENGTH); */ + /* potential->point_mass.z = */ + /* 50000 * PARSEC_IN_CGS / units_conversion_factor(us, UNIT_CONV_LENGTH); */ + /* potential->point_mass.mass = */ + /* 1e10 * SOLAR_MASS_IN_CGS / units_conversion_factor(us, UNIT_CONV_MASS); */ } diff --git a/src/potentials.h b/src/potentials.h index 2c75669539a5a52083077e7941baae6c2f44cf91..318daeb8df8ed9d7d6e337f5a1f1d12c8da9fce5 100644 --- a/src/potentials.h +++ b/src/potentials.h @@ -115,7 +115,8 @@ __attribute__((always_inline)) INLINE static void external_gravity_pointmass(con * @param us The current internal system of units * @param potential The external potential properties to initialize */ -void initPotentialProperties(struct UnitSystem* us, +void initPotentialProperties(const struct swift_params * paramter_file, + struct UnitSystem* us, struct external_potential* potential); #endif /* SWIFT_POTENTIALS_H */ diff --git a/src/runner.c b/src/runner.c index 8ac4bc7538354288229dc461346aac3a18eedb94..af1a489b074821fe51ee30b464d24551ca2051c2 100644 --- a/src/runner.c +++ b/src/runner.c @@ -176,7 +176,8 @@ void runner_dograv_external(struct runner *r, struct cell *c) { // g->x[1], vx, vy); message("%f\t %f %f %f %f %f %f %f %f %f %f %f %f %f\n", r->e->time, g->tx, g->tv, dt, v2, fg, fga, dr, E, L[2], g->x[0], g->x[1], vx, vy); - // message(" G=%e M=%e\n", const_G, External_Potential_Mass); + message(" G=%e M=%e\n", r->e->physical_constants->newton_gravity, r->e->potential->point_mass.mass); + exit(-1); } } }