-
Notifications
You must be signed in to change notification settings - Fork 13
Stellar oscillations
As a standard test of the numerical method for simulating self-gravitating fluids, we present the evolution of a stable isolated star in equilibrium. This test checks consistency and conservation properties for the coupled system of hydrodynamics and gravity.
For initial data, we first solve the Lane-Emden equation. This results
in a white dwarf with mass 0.2 solar mass and radius 4790 km. We provide
radial density profile here.
If you are interested in solving yourself, please
check wdID directoy in tools/starID_tools directory.
Using this profile, we set up particle configuration on a perturbed icosahedral lattice.
To do that, create initial_n32.par using below parameters
# Sperical particle distribution: initial data parameter file
#
initial_data_prefix = "wd_initial_n32"
# geometry: spherical domain with radius R = 1
domain_type = 1 # 0:box, 1:sphere
density_profile = "from file"
input_density_file = "density_profile_nr8000.dat"
sphere_radius = 4.79283e+08
# icosahedra lattice with small perturbations
lattice_nx = 50 # particle lattice dimension
lattice_type = 3 # 0:rectangular, 1:hcp, 2:fcc, 3:icosahedral
lattice_perturbation_amplitude = 0.09 # in units of sm. length
# equation of state type and parameters
eos_type = "polytropic"
#poly_gamma = 1.66667 # polytropic index
poly_gamma = 0.99
# density and pressure for relaxation stage
rho_initial = 5.2e+6
#pressure_initial = 1.56085e+23
pressure_initial = 1e+28
# since we only need spherical distribution of particles,
# set Sedov energy to zero
sedov_blast_energy = 0.0
# use a good kernel
sph_kernel = "Wendland C6"
sph_eta = 1.5
You can also directly get this parameter file here. Using this, we can generate initial data via below command
mpirun -np 1 ./sedov_3d_generator initial_n32.par
Note that you can find sedov_3d_generator in your ${PATH_TO_FLECSPH}/flecsph/build/app/id_generators or ${CMAKE_PRE_PATH}/bin based on your compilation procedure. After successful generation, you will have wd_initial_n32.h5part file.
After generating initial configuration, we relax this particle configuration. First, create relaxation_n32.par using below parameters
# Sperical particle distribution: relaxation phase
# Use ./hydro_3d <this_file>.par to evolve into relaxed state
# initial data
initial_data_prefix = "wd_initial_n32"
initial_iteration = 0
# geometry
domain_type = 1 # 0:box, 1:sphere
external_force_type = "spherical density support"
density_profile = "from file"
input_density_file = "density_profile_nr8000.dat"
sphere_radius = 4.79283e+08
# density and pressure for relaxation stage
rho_initial = 5.2e+6
##pressure_initial = 1.56085e+23 # actual pressure in WD
pressure_initial = 1e+28 # artificially increased for faster relaxation
# evolution
final_iteration = 1000
relaxation_steps = 1000
relaxation_beta = 1e2
relaxation_gamma = 0.0
initial_dt = 1.e-9
timestep_cfl_factor = 0.1
out_screen_every = 1
out_scalar_every = 1
out_diagnostic_every = 10
out_h5data_every = 100
output_h5data_prefix = "wd_relaxation_n32"
thermokinetic_formulation = yes
adaptive_timestep = yes
sph_variable_h = yes
evolve_internal_energy = no
# equation of state type and parameters
eos_type = "polytropic"
##poly_gamma = 1.66667 # actual polytropic index for WD
poly_gamma = 0.99 # artificially lowered for more uniform relaxation
sph_kernel = "Wendland C6"
sph_eta = 1.5
or get it from here. After that, using below command for relaxation step
mpirun -np 1 ./hydro_3d relaxation_n32.par
Note that you can find hydro_3d in your ${PATH_TO_FLECSPH}/flecsph/build/app/drivers or ${CMAKE_PRE_PATH}/bin based on your compilation procedure. After successful generation, you will have wd_relaxation_n32.h5part file.
After that, we modify the file produced in previous step by overwriting particles pressure
and internal energies to correspond self-consistent white dwarf model. To do that, create modify_n32.par
using below parameters
#
# Overwrite quantities before relaxation step
# Usage: ./sedov_3d_generator <this_file>.par
# WARNING: overwrites initial_data_prefix !!
#
modify_initial_data = yes
initial_data_prefix = "wd_modified_n32"
initial_iteration = 1000
# geometry: spherical domain with radius R = 1
domain_type = 1 # 0:box, 1:sphere
sphere_radius = 4.79283e+08
# equation of state type and parameters
eos_type = "polytropic"
poly_gamma = 1.66667 # polytropic index
# reset thermodynamical quantities
density_profile = "from file"
input_density_file = "density_profile_nr8000.dat"
external_force_type = "spherical density support"
rho_initial = 5.2e+6
pressure_initial = 1.56085e+23
sedov_blast_energy = 0.0
sedov_blast_radius = 0.1 # whatever
# good kernel
sph_kernel = "Wendland C6"
sph_eta = 1.5
or you can get it from here. Then, follow the below commands
cp wd_relaxation_n32.h5part wd_modified_n32.h5part
mpirun -np 1 ./sedov_3d_generator modify_n32.par
This will overwrite wd_modified_n32.h5part
Finally, we evolve white dwarf data from previous steps with self-gravity. To do that, create
evolve_n32.par using below parameters
#
# Test of Gravity
#
# initial data
initial_data_prefix = "wd_modified_n32"
eos_type = "polytropic"
poly_gamma = 1.66667 # polytropic index
rho_initial = 5.2e+6
pressure_initial = 1.56085e+23
sphere_radius = 4.79283e+08
domain_type = 1 # 0:box, 1:sphere
# gravity related parameters:
enable_fmm = yes
fmm_macangle = 0.3
#fmm_max_cell_mass = 0.1
# evolution parameters:
sph_kernel = "Wendland C6"
sph_eta = 1.6
initial_dt = 1.e-12
timestep_cfl_factor = 0.5
final_iteration = 10000
out_screen_every = 1
out_scalar_every = 1
out_h5data_every = 100
out_diagnostic_every = 10
output_h5data_prefix = "wd_evolution_n32"
thermokinetic_formulation = yes
adaptive_timestep = yes
sph_variable_h = yes
evolve_internal_energy = yes
relaxation_repulsion_gamma = 0.0
gravitational_constant = 6.67383e-8
or you can get it from here. Then, follow the below commands
mpirun -np 1 ./newtonian_3d evolve_n32.par
Note that you can find newtonian_3d in your ${PATH_TO_FLECSPH}/flecsph/build/app/drivers or ${CMAKE_PRE_PATH}/bin based on your compilation procedure. After successful generation, you will have wd_evolution_n32.h5part file.
The purpose of this test is checking conservation property for coupled system of hydrodynamics and gravity. For gravity computation, we use fast multipole method (FMM). We compare conservation of energy, linear momentum, and angular momentum, computed using the FMM approximation, and using exact N-body computation.
You can obtain both FMM and exact N-body results by changing fmm_macangle parameter in evolve_n32.par.
If you set fmm_macangle=0.0 this will provide exact N-body result. If you set fmm_macangle > 0.0, you will
use FMM approximation to compute gravitational force.
After successful evolution simulation, you will get scalar_reductions.dat file that contains information of
energy, linear and angular momentum during iterations. You can compare FMM results with exact N-body results by varying fmm_macangle.
There are several different ways to plot the results. We encourage user to use whatever user's convenient way but please contact [email protected] if you would like to get plotting script and example result data from us.