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
from hotwheels_core import *
from hotwheels_pm import *
from hotwheels_integrate import *
from hotwheels_io import *

#
# Step 1: Configure components
# This stage configures components without allocating resources.
# Configurations are passed to constructors to compile the underlying C libraries.
#

mpi = hwc.MPI().init() # Initialize MPI
mym = MyMalloc(alloc_bytes=int(2e9)) # Configure memory allocator with 2GB
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
# Set up a fixed time-step integrator from 0 to 1 Gyr
# Conversion factor for Gyr to internal units
gyr_to_cu = 3.086e+16 / (1e9 * 3600 * 24 * 365)
ts = integrate.FixedTimeStep(
    soas, 
    G=43007.1,  # Gravitational constant in specific units
    t_from=0., 
    t_to=1. * gyr_to_cu, 
    MPI=mpi
)
# Initialize a NFW profile with scale radius `rs=100` and density `rho0=1e-6`
ic = NFWIC(rs=100., rho0=1e-6, rs_factor=10.)
# Configure a refined PM grid with 7 stacked high-resolution regions
pm = SuperHiResPM( #wrapper to the PM C library
    soas=soas, 
    mem=mym, 
    TS=ts, #will use it to attach gravkick callback
    MPI=mpi, 
    pmgrid=128, 
    grids=8, # number of grids to instantiate
    dt_displacement_factor=0.25 #factor for DtDisplacement
)
build = make.Build('./', mpi, pm, ts, mym, *soas.values()) # Compile all modules in the current directory
headers = OnTheFly(build.build_name, *build.components, generate_user_c=True) # Generate SoA headers

if mpi.rank == 0:  # Master rank handles compilation
    headers.write()
    build.compile()

#
# Step 2: Allocate resources
#

with (
    utils.Panic(Build=build) as panic,  # Attach panic handler
    utils.Timer(Build=build) as timer,  # Attach timer handler
    build.enter(debug=mpi.rank == 0),  # Parse compiled objects
    mpi.enter(pm),  # Initialize MPI in the PM module
    mym.enter(*build.components),  # Allocate 2GB memory
    p,  # Allocate particle data structure in MyMalloc
    ic.enter(p, mpi.ranks, p.get_maxpart()),  # Sample NFW profile
    pm,  # Initialize PM and compute first accelerations
    ts  # Compute DriftTables if needed
):

    #
    # Step 3: Main simulation loop
    #

    while ts.time < ts.time_end:
        ts.find_timesteps()  # Determine timesteps
        ts.do_first_halfstep_kick()  # First kick (includes drift/kick callbacks)
        ts.drift()  # Update particle positions
        pm.compute_accelerations()  # Recompute accelerations
        ts.do_second_halfstep_kick()  # Second kick

        # Occasionally, generate plots on the master rank
        if mpi.rank == 0 and ts.steps % 10 == 0:
            fig, ax = plt.subplots(1)
            ax.hist2d(p['pos'][:, 0], p['pos'][:, 1], bins=128)
            ax.set_aspect('equal')
            fig.savefig(f'snap{ts.steps}_rank{mpi.rank}.png', bbox_inches='tight', dpi=200)
            plt.close(fig)

print('Simulation finished')