From d7ca62d0d92d5d756468cdeb4b87055fce73e685 Mon Sep 17 00:00:00 2001
From: Antonio Ragagnin <antonio.ragagnin@inaf.it>
Date: Sat, 14 Dec 2024 13:28:14 +0000
Subject: [PATCH] Add new file
---
run_pm_dmo_NFW_fixed_timestep.md | 66 ++++++++++++++++++++++++++++++++
1 file changed, 66 insertions(+)
create mode 100644 run_pm_dmo_NFW_fixed_timestep.md
diff --git a/run_pm_dmo_NFW_fixed_timestep.md b/run_pm_dmo_NFW_fixed_timestep.md
new file mode 100644
index 0000000..7ecfa06
--- /dev/null
+++ b/run_pm_dmo_NFW_fixed_timestep.md
@@ -0,0 +1,66 @@
+To run this module you need access to hotwheels PM, integrate, and IO components.
+
+```python
+import numpy as np, os, matplotlib as plt
+from hotwheels_core import *
+from hotwheels_pm import *
+from hotwheels_integrate import *
+from hotwheels_io import *
+
+# if G=43007.1, masses are in 1e10Msun/h, and radii are in ckpc/h, then this
+# constant converts from Gyr to internal units. Not much to comment about it
+gyr_to_cu = 3.086e+16 / (1e9 * 3600 * 24 * 365)
+#call MPI_Init()
+mpi = hwc.MPI().init()
+#configure my malloc with 2GB
+mym = MyMalloc(alloc_bytes=int(2e9))
+#configure the default particle SoA with 1e5 particles
+p = SoA(maxpart=int(1e5), mem=mym)
+#add the particle soa to the multi-type particle SoAs
+soas = SoAs(p, mem=mym)
+# configure the timestep class to go from 0 to 1 Gyr
+# note: G=43007.1 in units of length=kpc, velocity=km/s, mass = 1e10Msun
+ts = integrate.FixedTimeStep(soas, G=43007.1, t_from=0., t_to=1. * gyr_to_cu, MPI=mpi)
+ic = NFWIC(rs=100., rho0=1e-6, rs_factor=10.) #init the build of a NFW profile with R<10rs, (note that mass is in units of 1e10Msun)
+# init a refined PM grids with 7 stacked PLACEHIRESREGION at smaller and smaller scales
+# note that PM needs the time-step integrator class TS to attach its DriftTable and kick callbacks
+pm = SuperHiResPM(soas=soas, mem=mym, TS=ts, MPI=mpi, pmgrid=128, grids=8, dt_displacement_factor=0.25)
+#configure to compile all modules in the current folder
+build = make.Build('./', mpi, pm, ts, mym, *soas.values())
+#configure generate SoA headers
+headers = OnTheFly(build.build_name, *build.components, generate_user_c=True)
+
+if mpi.rank == 0: #master rank compile and build headers
+ headers.write()
+ build.compile()
+
+with (utils.Panic(Build=build) as panic, #attach panic handler to C calls
+ utils.Timer(Build=build) as timer, #attach timer handler to C calls
+ build.enter(debug=mpi.rank == 0), #parse the compiled object files
+ mpi.enter(pm), #attach MPI init info to PM module
+ mym.enter(*build.components), #actually allocates the 2GB of ram
+ p, #allocate the particle data stracture in the MyMalloc area
+ ic.enter(p, mpi.ranks, p.get_maxpart()), #sample the NFW in the particle SoA `p` fields
+ pm, #call pm_init() the PM regions
+ ts #compute DriftTables if necessary
+ ):
+
+ pm.compute_accelerations() #first acc computation
+
+ while ts.time < ts.time_end: #main run.c loop
+
+ ts.find_timesteps() #integrator will find timesteps
+ #integrator will call halfstep kick, including PM registered drift and kick callbacks
+ ts.do_first_halfstep_kick()
+ ts.drift() #integretor will drift
+ pm.compute_accelerations() #accelerations
+ #integrator will call halfstep kick, including PM registered drift and kick callbacks
+ ts.do_second_halfstep_kick()
+ if mpi.rank == 0 and steps % 10 == 0: #sometimes, master rank will do a plot
+ f, ax = plt.subplots(1)
+ ax.hist2d(p['pos'][:, 0], p['pos'][:, 1], bins=128)
+ ax.set_aspect('equal')
+ pat = os.path.join(os.getenv('HW_BUILD', '.'), 'snap%d_rank%d.png' % (steps, mpi.rank))
+ f.savefig(pat, bbox_inches='tight', dpi=200)
+ plt.close(f)
+print('simualtion finished')
\ No newline at end of file
--
GitLab