-
Notifications
You must be signed in to change notification settings - Fork 66
Commit
This commit does not belong to any branch on this repository, and may belong to a fork outside of the repository.
Merge pull request #133 from PhoebeSandhaus/main
Implementation of Gas Damping Forces
- Loading branch information
Showing
8 changed files
with
613 additions
and
3 deletions.
There are no files selected for viewing
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -0,0 +1,48 @@ | ||
| export OPENGL=0 | ||
|
|
||
| ifndef REB_DIR | ||
| ifneq ($(wildcard ../../../rebound/.*),) # Check for REBOUND in default location | ||
| REB_DIR=../../../rebound | ||
| endif | ||
| ifneq ($(wildcard ../../../../rebound/.*),) # Check for REBOUNDx being inside REBOUND directory | ||
| REB_DIR=../../../ | ||
| endif | ||
| endif | ||
| ifndef REB_DIR # REBOUND is not in default location and REB_DIR is not set | ||
| $(error REBOUNDx not in the same directory as REBOUND. To use a custom location, you Must set the REB_DIR environment variable for the path to your rebound directory, e.g., export REB_DIR=/Users/dtamayo/rebound. See reboundx.readthedocs.org) | ||
| endif | ||
| PROBLEMDIR=$(shell basename `dirname \`pwd\``)"/"$(shell basename `pwd`) | ||
|
|
||
| include $(REB_DIR)/src/Makefile.defs | ||
|
|
||
| REBX_DIR=../../ | ||
|
|
||
| all: librebound.so libreboundx.so | ||
| @echo "" | ||
| @echo "Compiling problem file ..." | ||
| $(CC) -I$(REBX_DIR)/src/ -I$(REB_DIR)/src/ -Wl,-rpath,./ $(OPT) $(PREDEF) problem.c -L. -lreboundx -lrebound $(LIB) -o rebound | ||
| @echo "" | ||
| @echo "Problem file compiled successfully." | ||
|
|
||
| librebound.so: | ||
| @echo "Compiling shared library librebound.so ..." | ||
| $(MAKE) -C $(REB_DIR)/src/ | ||
| @echo "Creating link for shared library librebound.so ..." | ||
| @-rm -f librebound.so | ||
| @ln -s $(REB_DIR)/src/librebound.so . | ||
|
|
||
| libreboundx.so: | ||
| @echo "Compiling shared library libreboundx.so ..." | ||
| $(MAKE) -C $(REBX_DIR)/src/ | ||
| @-rm -f libreboundx.so | ||
| @ln -s $(REBX_DIR)/src/libreboundx.so . | ||
|
|
||
| clean: | ||
| @echo "Cleaning up shared library librebound.so ..." | ||
| @-rm -f librebound.so | ||
| $(MAKE) -C $(REB_DIR)/src/ clean | ||
| @echo "Cleaning up shared library libreboundx.so ..." | ||
| @-rm -f libreboundx.so | ||
| $(MAKE) -C $(REBX_DIR)/src/ clean | ||
| @echo "Cleaning up local directory ..." | ||
| @-rm -vf rebound |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -0,0 +1,74 @@ | ||
| /** | ||
| * Gas Damping Timescale Example | ||
| * | ||
| * This example shows how to add gas damping to a REBOUND simulation. | ||
| * | ||
| */ | ||
| #include <stdio.h> | ||
| #include <stdlib.h> | ||
| #include <unistd.h> | ||
| #include <math.h> | ||
| #include <string.h> | ||
| #include "rebound.h" | ||
| #include "reboundx.h" | ||
|
|
||
| void heartbeat(struct reb_simulation* sim); | ||
|
|
||
| int main(int argc, char* argv[]){ | ||
|
|
||
| /* start the rebound simulation here */ | ||
| struct reb_simulation* sim = reb_simulation_create(); | ||
|
|
||
| // set units to use throughout the simulation | ||
| sim->G = 4*M_PI*M_PI; // Gravitational constant in AU^3/M_sun/yr^2 | ||
|
|
||
| sim->heartbeat = heartbeat; | ||
|
|
||
| // add the host star | ||
| struct reb_particle star = {0}; | ||
| star.m = 1.; | ||
| reb_simulation_add(sim, star); | ||
|
|
||
| double m = 3.e-6; // roughly 1 Earth Mass in Solar Masses | ||
| double a = 0.1; | ||
| double e = 0.05; | ||
| double inc = 5.*M_PI/180.; // 5 degrees in radians | ||
| double Omega = 0.; | ||
| double omega = 0.; | ||
| double f = 0.; | ||
|
|
||
| struct reb_particle planet = reb_particle_from_orbit(sim->G, star, m, a, e, inc, Omega, omega, f); | ||
| reb_simulation_add(sim,planet); | ||
|
|
||
| // move simulation to center-of-mass frame | ||
| reb_simulation_move_to_com(sim); | ||
|
|
||
| // add in reboundx and load/add new force | ||
| struct rebx_extras* rebx = rebx_attach(sim); | ||
| struct rebx_force* gd = rebx_load_force(rebx, "gas_damping_timescale"); | ||
| rebx_add_force(rebx, gd); | ||
|
|
||
| // Set the total simulation time | ||
| double tmax = 5.e2; // 500 years | ||
|
|
||
| // Set parameters for simulation | ||
| rebx_set_param_double(rebx, &gd->ap, "cs_coeff", 0.272); // set speed sound coefficient to 0.272 AU^(3/4) yr^-1 | ||
| rebx_set_param_double(rebx, &gd->ap, "tau_coeff", 0.003); // set timescale coefficient to 0.003 yr AU^-2 | ||
|
|
||
| // Set parameter for particle | ||
| rebx_set_param_double(rebx, &sim->particles[1].ap, "d_factor", 5.); // set depletion factor to 5 | ||
|
|
||
| // integrate simulation | ||
| reb_simulation_integrate(sim, tmax); | ||
| rebx_free(rebx); // free all the memory allocated by reboundx | ||
| reb_simulation_free(sim); // free all the memory allocated by rebound | ||
| } | ||
|
|
||
| void heartbeat(struct reb_simulation* sim){ | ||
| // output a e i of the planet | ||
| if(reb_simulation_output_check(sim, 50.)){ | ||
| const struct reb_particle star = sim->particles[0]; | ||
| const struct reb_orbit orbit = reb_orbit_from_particle(sim->G, sim->particles[1], star); // calculate orbit of planet | ||
| printf("%f\t%f\t%f\t%f\n",sim->t, orbit.a, orbit.e, orbit.inc); | ||
| } | ||
| } |
Large diffs are not rendered by default.
Oops, something went wrong.
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -0,0 +1,141 @@ | ||
| /** | ||
| * @file gas_damping_timescale.c | ||
| * @brief Update orbits with prescribed timescales by directly changing orbital elements after each timestep | ||
| * @author Phoebe Sandhaus <[email protected]> | ||
| * | ||
| * @section LICENSE | ||
| * Copyright (c) 2015 Dan Tamayo, Hanno Rein | ||
| * | ||
| * This file is part of reboundx. | ||
| * | ||
| * reboundx is free software: you can redistribute it and/or modify | ||
| * it under the terms of the GNU General Public License as published by | ||
| * the Free Software Foundation, either version 3 of the License, or | ||
| * (at your option) any later version. | ||
| * | ||
| * reboundx is distributed in the hope that it will be useful, | ||
| * but WITHOUT ANY WARRANTY; without even the implied warranty of | ||
| * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the | ||
| * GNU General Public License for more details. | ||
| * | ||
| * You should have received a copy of the GNU General Public License | ||
| * along with rebound. If not, see <http://www.gnu.org/licenses/>. | ||
| * | ||
| * The section after the dollar signs gets built into the documentation by a script. All lines must start with space * space like below. | ||
| * Tables always must be preceded and followed by a blank line. See http://docutils.sourceforge.net/docs/user/rst/quickstart.html for a primer on rst. | ||
| * $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ | ||
| * | ||
| * $Gas Effects$ // Effect category (must be the first non-blank line after dollar signs and between dollar signs to be detected by script). | ||
| * | ||
| * ======================= =============================================== | ||
| * Authors Phoebe Sandhaus | ||
| * Implementation Paper `Sandhaus et al. in prep` | ||
| * Based on `Dawson et al. 2016 <https://ui.adsabs.harvard.edu/abs/2016ApJ...822...54D/abstract>; Kominami & Ida 2002 <https://ui.adsabs.harvard.edu/abs/2002Icar..157...43K/abstract>`_ | ||
| * C Example :ref:`c_example_gas_damping_timescale` | ||
| * Python Example `GasDampingTimescale.ipynb <https://github.com/PhoebeSandhaus/reboundx_gas_damping/tree/main/ipython_examples/GasDampingTimescale.ipynb>`_ | ||
| * ======================= =============================================== | ||
| * | ||
| * This updates particles' positions and velocities between timesteps by first calculating a damping timescale for each individual particle, and then applying the timescale to damp both the eccentricity and inclination of the particle. Note: The timescale of damping should be much greater than a particle's orbital period. The damping force should also be small as compared to the gravitational forces on the particle. | ||
| * | ||
| * **Effect Parameters** | ||
| * | ||
| * ============================ =========== ================================================================== | ||
| * Field (C type) Required Description | ||
| * ============================ =========== ================================================================== | ||
| * None - - | ||
| * ============================ =========== ================================================================== | ||
| * | ||
| * **Particle Parameters** | ||
| * | ||
| * ============================ =========== ============================================================================== | ||
| * Field (C type) Required Description | ||
| * ============================ =========== ============================================================================== | ||
| * d_factor (double) Yes Depletion factor d in Equation 16 from Dawson et al. 2016; d=1 corresponds roughly to the full minimum mass solar nebula with Sigma_gas (surface gas density) = 1700 g cm^-2 at 1 AU [for d=1]; d>1 corresponds to a more depleted nebula | ||
| * cs_coeff (double) Yes Sound speed coefficient; Changing the value will change assumed units; Example: If you are using the units: AU, M_sun, yr --> cs_coeff = 0.272 # AU^(3/4) yr^-1 | ||
| * tau_coeff (double) Yes Timescale coefficient; Changing the value will change assumed units; Example: If you are using the units: AU, M_sun, yr --> tau_coeff = 0.003 # yr AU^-2 | ||
| * ============================ =========== ============================================================================== | ||
| * | ||
| */ | ||
|
|
||
| #include <stdio.h> | ||
| #include <stdlib.h> | ||
| #include <math.h> | ||
| #include "rebound.h" | ||
| #include "reboundx.h" | ||
| #include "rebxtools.h" | ||
|
|
||
| static struct reb_vec3d rebx_calculate_gas_damping_timescale(struct reb_simulation* const sim, struct rebx_force* const force, struct reb_particle* planet, struct reb_particle* star){ | ||
| struct rebx_extras* const rebx = sim->extras; | ||
| struct reb_orbit o = reb_orbit_from_particle(sim->G, *planet, *star); | ||
|
|
||
| const double* const d_factor = rebx_get_param(rebx, planet->ap, "d_factor"); | ||
| const double* const cs_coeff = rebx_get_param(rebx, force->ap, "cs_coeff"); | ||
| const double* const tau_coeff = rebx_get_param(rebx, force->ap, "tau_coeff"); | ||
|
|
||
| struct reb_vec3d a = {0}; | ||
|
|
||
| if (d_factor == NULL || cs_coeff == NULL || tau_coeff == NULL){ | ||
| rebx_error(rebx, "Need to set d_factor, cs_coeff, tau_coeff parameters. See examples in documentation.\n"); | ||
| return a; | ||
| } | ||
|
|
||
| // initialize positions and velocities | ||
| const double dvx = planet->vx - star->vx; | ||
| const double dvy = planet->vy - star->vy; | ||
| const double dvz = planet->vz - star->vz; | ||
| const double dx = planet->x-star->x; | ||
| const double dy = planet->y-star->y; | ||
| const double dz = planet->z-star->z; | ||
| const double r2 = dx*dx + dy*dy + dz*dz; | ||
|
|
||
| // initial semimajor axis, eccentricity, and inclination | ||
| const double a0 = o.a; | ||
| const double e0 = o.e; | ||
| const double inc0 = o.inc; | ||
| const double starMass = star->m; | ||
| const double planetMass = planet->m; | ||
|
|
||
| // eccentricity and inclination timescales from Dawson+16 Eqn 16 | ||
| double coeff; | ||
|
|
||
| double vk = sqrt(sim->G*starMass/a0); | ||
| double v = sqrt(e0*e0+inc0*inc0)*vk; | ||
| double cs = *cs_coeff/sqrt(sqrt(a0)); | ||
| double v_over_cs = v/cs; | ||
|
|
||
| if (v <= cs){ | ||
| coeff = 1.; | ||
| } | ||
| else { | ||
| if (inc0 < cs/vk) { | ||
| coeff = v_over_cs*v_over_cs*v_over_cs; | ||
| } | ||
| else { | ||
| coeff = v_over_cs*v_over_cs*v_over_cs*v_over_cs; | ||
| } | ||
| } | ||
|
|
||
| double tau_e = -*tau_coeff*(*d_factor)*a0*a0*(starMass/planetMass)*coeff; | ||
| double tau_inc = 2.*tau_e; // from Kominami & Ida 2002 [Eqs. 2.9 and 2.10] | ||
|
|
||
|
|
||
| if (tau_e < INFINITY || tau_inc < INFINITY){ | ||
| const double vdotr = dx*dvx + dy*dvy + dz*dvz; | ||
| const double prefac = 2*vdotr/r2/tau_e; | ||
| a.x += prefac*dx; | ||
| a.y += prefac*dy; | ||
| a.z += prefac*dz + 2.*dvz/tau_inc; | ||
| } | ||
| return a; | ||
| } | ||
|
|
||
| void rebx_gas_damping_timescale(struct reb_simulation* const sim, struct rebx_force* const force, struct reb_particle* const particles, const int N){ | ||
| int* ptr = rebx_get_param(sim->extras, force->ap, "coordinates"); | ||
| enum REBX_COORDINATES coordinates = REBX_COORDINATES_JACOBI; // Default | ||
| if (ptr != NULL){ | ||
| coordinates = *ptr; | ||
| } | ||
| const int back_reactions_inclusive = 1; | ||
| const char* reference_name = "primary"; | ||
| rebx_com_force(sim, force, coordinates, back_reactions_inclusive, reference_name, rebx_calculate_gas_damping_timescale, particles, N); | ||
| } |