Skip to content
Snippets Groups Projects
Commit 4644157a authored by Antonio Ragagnin's avatar Antonio Ragagnin :speech_balloon:
Browse files

Update run_pm_dmo_NFW_fixed_timestep.md

parent 15c6caa3
Branches
No related tags found
No related merge requests found
...@@ -5,15 +5,22 @@ In this example we will use create IC that sample a **NFW profile** and evolve i ...@@ -5,15 +5,22 @@ In this example we will use create IC that sample a **NFW profile** and evolve i
Here is a first benefit of a modular code: since in `hotwheels `PM is a self-contained module, we can instantiate it an arbitrary number of times. So one can stack seven [PLACEHIGHRESREGION](https://wwwmpa.mpa-garching.mpg.de/gadget4/03_simtypes/) Here is a first benefit of a modular code: since in `hotwheels `PM is a self-contained module, we can instantiate it an arbitrary number of times. So one can stack seven [PLACEHIGHRESREGION](https://wwwmpa.mpa-garching.mpg.de/gadget4/03_simtypes/)
on smaller and smaller regions (a sort of refined mesh) on top of a sampled NFW halo and use PM-only to **get accurate force down to a kpc** (see image). on smaller and smaller regions (a sort of refined mesh) on top of a sampled NFW halo and use PM-only to **get accurate force down to a kpc** (see image).
Install the PM module (will install core,IO, and timestep as dependencies):
Note that to run this module you need access to hotwheels **core, IO, PM,** and **integrate** components. Note that `hotwheels` do not provide parameter or config files. It is up to the user to initalise its sub-library components and connect them. ```bash
pip install 'git+ssh://git@git.ia2.inaf.it/hotwheels/PM.git@v0.0.0alpha'
```
And you can run this code:
```python ```python
import numpy as np, os, matplotlib as plt import numpy as np, os, matplotlib as plt
from hotwheels_core import * from hotwheels.utils import *
from hotwheels_pm import * from hotwheels.wrap import *
from hotwheels_integrate import * from hotwheels.soas import *
from hotwheels_io import * from hotwheels.PM import *
from hotwheels.integrate import *
from hotwheels.io import *
# #
# Step 1: Configure components # Step 1: Configure components
...@@ -21,14 +28,14 @@ from hotwheels_io import * ...@@ -21,14 +28,14 @@ from hotwheels_io import *
# Configurations are passed to constructors to compile the underlying C libraries. # Configurations are passed to constructors to compile the underlying C libraries.
# #
mpi = hwc.MPI().init() # Initialize MPI mpi = MPI().init() # Initialize MPI
mym = MyMalloc(alloc_bytes=int(2e9)) # Configure memory allocator with 2GB mym = MyMalloc(alloc_bytes=int(2e9)) # Configure memory allocator with 2GB
p = SoA(maxpart=int(1e5), mem=mym) # Configure P to hold 1e5 particles p = SoA(maxpart=int(1e5), mem=mym) # Configure P to hold 1e5 particles
soas = SoAs(p, mem=mym) # Add P to a multi-type SoA container soas = SoAs(p) # Add P to a multi-type SoA container
# Set up a fixed time-step integrator from 0 to 1 Gyr # Set up a fixed time-step integrator from 0 to 1 Gyr
# Conversion factor for Gyr to internal units # Conversion factor for Gyr to internal units
gyr_to_cu = 3.086e+16 / (1e9 * 3600 * 24 * 365) gyr_to_cu = 3.086e+16 / (1e9 * 3600 * 24 * 365)
ts = integrate.FixedTimeStep( ts = FixedTimeStep(
soas, soas,
G=43007.1, # Gravitational constant in specific units G=43007.1, # Gravitational constant in specific units
t_from=0., t_from=0.,
...@@ -36,7 +43,7 @@ ts = integrate.FixedTimeStep( ...@@ -36,7 +43,7 @@ ts = integrate.FixedTimeStep(
MPI=mpi MPI=mpi
) )
# Initialize a NFW profile with scale radius `rs=100` and density `rho0=1e-6` # Initialize a NFW profile with scale radius `rs=100` and density `rho0=1e-6`
ic = NFWIC(rs=100., rho0=1e-6, rs_factor=10.) ic = NFWIC(r_s=100., rho_0=1e-6, r_max_f=10.)
# Configure a refined PM grid with 7 stacked high-resolution regions # Configure a refined PM grid with 7 stacked high-resolution regions
pm = SuperHiResPM( #wrapper to the PM C library pm = SuperHiResPM( #wrapper to the PM C library
soas=soas, soas=soas,
...@@ -59,13 +66,13 @@ if mpi.rank == 0: # Master rank handles compilation ...@@ -59,13 +66,13 @@ if mpi.rank == 0: # Master rank handles compilation
# #
with ( with (
utils.Panic(Build=build) as panic, # Attach panic handler Panic(Build=build) as panic, # Attach panic handler
utils.Timer(Build=build) as timer, # Attach timer handler Timer(Build=build) as timer, # Attach timer handler
build.enter(debug=mpi.rank == 0), # Parse compiled objects build.enter(debug=mpi.rank == 0), # Parse compiled objects
mpi.enter(pm), # Initialize MPI in the PM module mpi.enter(pm), # Initialize MPI in the PM module
mym.enter(*build.components), # Allocate 2GB memory mym.enter(*build.components), # Allocate 2GB memory
p, # Allocate particle data structure in MyMalloc p, # Allocate particle data structure in MyMalloc
ic.enter(p, mpi.ranks, p.get_maxpart()), # Sample NFW profile ic.enter(p, mpi.ranks, p.get_maxpart(), ts.G), # Sample NFW profile
pm, # Initialize PM and compute first accelerations pm, # Initialize PM and compute first accelerations
ts # Compute DriftTables if needed ts # Compute DriftTables if needed
): ):
...@@ -73,7 +80,6 @@ with ( ...@@ -73,7 +80,6 @@ with (
# #
# Step 3: Main simulation loop # Step 3: Main simulation loop
# #
while ts.time < ts.time_end: while ts.time < ts.time_end:
ts.find_timesteps() # Determine timesteps ts.find_timesteps() # Determine timesteps
ts.do_first_halfstep_kick() # First kick (includes drift/kick callbacks) ts.do_first_halfstep_kick() # First kick (includes drift/kick callbacks)
...@@ -90,3 +96,4 @@ with ( ...@@ -90,3 +96,4 @@ with (
plt.close(fig) plt.close(fig)
print('Simulation finished') print('Simulation finished')
```
\ No newline at end of file
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment