Update run_pm_dmo_NFW_fixed_timestep.md

run_pm_dmo_NFW_fixed_timestep.md

-To run this module you need access to hotwheels PM, integrate, and IO components.
+In this example we will use create IC that sample a **NFW profile** and evolve it for 1Gyr with gravity computed using several layers of a **refined particle-mesh (PM)**.
+
+![density map and acceleration of a NFW mock profile with 7 stacked PLACEHIGHRESREGION PM components](NFW_PM_fixed_timestep.png)
+
+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).
+
+
+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.
 
 ```python
 import numpy as np, os, matplotlib as plt
@@ -7,9 +15,14 @@ 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)
+#
+# step 1: config of components
+#
+# at this stage we do not allocate any resource, we just
+# pass the right config parameters to the constructors
+# in order to compile the underlying C libraries.
+#
+
 #call MPI_Init()
 mpi = hwc.MPI().init()
 #configure my malloc with 2GB
@@ -20,6 +33,9 @@ p = SoA(maxpart=int(1e5), mem=mym)
 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
+# 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)
 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
@@ -30,10 +46,20 @@ build = make.Build('./', mpi, pm, ts, mym, *soas.values())
 #configure generate SoA headers
 headers = OnTheFly(build.build_name, *build.components, generate_user_c=True)
 
+#
+# step 2: build SoAs headers and compile C files
+#
+
+
 if mpi.rank == 0: #master rank compile and build headers
     headers.write()
     build.compile()
 
+
+#
+# step 3: allocate resources (e.g. MPI, MyMalloc, and allocations within it)
+#
+
 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
@@ -47,6 +73,12 @@ with (utils.Panic(Build=build) as panic, #attach panic handler to C calls
 
         pm.compute_accelerations() #first acc computation
 
+
+        #
+        # step 4: the actual run
+        #
+
+
         while ts.time < ts.time_end: #main run.c loop
 
             ts.find_timesteps() #integrator will find timesteps
@@ -56,6 +88,7 @@ with (utils.Panic(Build=build) as panic, #attach panic handler to C calls
             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)
@@ -63,4 +96,6 @@ with (utils.Panic(Build=build) as panic, #attach panic handler to C calls
                 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