diff --git a/.github/workflows/shamrock-acpp-phys-test.yml b/.github/workflows/shamrock-acpp-phys-test.yml
index 9fcfaf2758..30a1e0c45c 100644
--- a/.github/workflows/shamrock-acpp-phys-test.yml
+++ b/.github/workflows/shamrock-acpp-phys-test.yml
@@ -126,6 +126,8 @@ jobs:
             testfile: regression_godunov_soundwave_3d.py
           - worldsize: 1
             testfile: run_sg_compare_error_sph.py
+          - worldsize: 1
+            testfile: run_mhd_wavedamping.py
           - worldsize: 1
             testfile: run_compare_shamrock_ph_disc.py
           - worldsize: 1
diff --git a/examples/sph/run_mhd_wavedamping.py b/examples/sph/run_mhd_wavedamping.py
new file mode 100644
index 0000000000..0113f3a3b0
--- /dev/null
+++ b/examples/sph/run_mhd_wavedamping.py
@@ -0,0 +1,283 @@
+"""
+Wave damping test in SPH with non-ideal MHD
+===========================================
+
+This example runs the wave damping test, aimed at evaluating the behaviour of the ambipolar diffusion term.
+The RMS of the magnetic field components is
+monitored and compared to the analytical dispersion relation.
+"""
+
+# sphinx_gallery_multi_image = "single"
+
+import os
+
+import matplotlib.pyplot as plt
+import numpy as np
+
+import shamrock
+
+shamrock.enable_experimental_features()
+
+# Initialize shamrock (if not already done by the executable)
+if not shamrock.sys.is_initialized():
+    shamrock.change_loglevel(1)
+    shamrock.sys.init("0:0")
+
+# Use shamrock plotting style
+shamrock.matplotlib.set_shamrock_mpl_style()
+
+# %%
+# Define physical parameters and unit system
+# ------------------------------------------
+# We use a unit system where mu_0 = 1 in code units (defined via UnitSystem).
+# The code unit system is set to have unit_length = 1, unit_mass = 1.2566370621219e-06
+# so that mu_0 = 1 exactly.
+
+Lx = 1.0  # box length
+dr = 1 / 128  # particle spacing
+rho0 = 1.0  # initial density
+Bx0 = 1.0  # background field in x
+C_ADc = 0.01  # ambipolar diffusion coefficient (Phantom convention)
+cs = 1.0  # isothermal sound speed
+t_target = 5.0  # total simulation time
+dt_dump = 0.01  # dump interval
+
+# Unit system and constants
+codeu = shamrock.UnitSystem(
+    unit_time=1.0,
+    unit_length=1.0,
+    unit_mass=1.2566370621219e-06,
+)
+ucte = shamrock.Constants(codeu)
+mu_0 = ucte.mu_0()  # = 1 in these units
+c = ucte.c()  # speed of light (used only for unit conversion)
+
+# Derived quantities
+vA = Bx0 / np.sqrt(rho0)  # Alfven speed
+etaAD_cgs = C_ADc * vA * vA  # ambipolar diffusivity (cgs-like)
+etaAD_si = etaAD_cgs * (4 * np.pi / c)  # SI conversion (not used in code)
+
+print(f"mu_0 = {mu_0}, vA = {vA:.3f}, etaAD = {etaAD_cgs:.3e}")
+
+# %%
+# Create context and SPH model
+# ----------------------------
+ctx = shamrock.Context()
+ctx.pdata_layout_new()
+model = shamrock.get_Model_SPH(context=ctx, vector_type="f64_3", sph_kernel="C4")
+
+# %%
+# Set up simulation configuration
+# -------------------------------
+cfg = model.gen_default_config()
+cfg.set_units(codeu)
+cfg.set_artif_viscosity_None()  # no artificial viscosity
+cfg.set_NonIdealMHD(
+    sigma_mhd=0, sigma_u=0, etaO=0, etaH=0, etaAD=etaAD_cgs, alpha_B=0, alpha_AV=0, beta_AV=0
+)
+cfg.set_boundary_periodic()  # periodic boundaries in all directions
+cfg.set_eos_isothermal(cs)  # isothermal equation of state
+cfg.print_status()
+model.set_solver_config(cfg)
+
+# Initialize scheduler and particle container
+model.init_scheduler(int(1e6), 1)
+
+# %%
+# Generate particle distribution in an FCC lattice
+# ------------------------------------------------
+bmin = (-Lx / 2.0, -np.sqrt(3) / 4.0 * Lx, -np.sqrt(6) / 4.0 * Lx)
+bmax = (Lx / 2.0, np.sqrt(3) / 4.0 * Lx, np.sqrt(6) / 4.0 * Lx)
+
+# Adjust box to exactly fit the lattice
+bmin, bmax = model.get_ideal_fcc_box(dr, bmin, bmax)
+xm, ym, zm = bmin
+xM, yM, zM = bmax
+Lx_actual = xM - xm
+
+model.resize_simulation_box(bmin, bmax)
+model.add_cube_fcc_3d(dr, bmin, bmax)
+
+vol_b = (xM - xm) * (yM - ym) * (zM - zm)
+totmass = rho0 * vol_b
+pmass = model.total_mass_to_part_mass(totmass)
+model.set_particle_mass(pmass)
+
+# %%
+# Set initial conditions
+# ----------------------
+# Background magnetic field in x, and a sinusoidal velocity perturbation in z.
+k = 2 * np.pi / Lx_actual
+v0 = 0.01 * vA
+
+
+def B_func(r):
+    return (Bx0, 0.0, 0.0)
+
+
+def vel_func(r):
+    x, y, z = r
+    vz = v0 * np.sin(k * (x - xm))
+    return (0.0, 0.0, vz)
+
+
+def u_func(r):
+    return 0.0
+
+
+model.set_field_value_lambda_f64_3("B/rho", B_func)
+model.set_field_value_lambda_f64_3("vxyz", vel_func)
+model.set_field_value_lambda_f64("uint", u_func)
+
+# %%
+# Set CFL parameters
+# ------------------
+model.set_cfl_cour(0.3)
+model.set_cfl_force(0.25)
+
+# %%
+# Prepare storage for time series and output directory
+# ----------------------------------------------------
+times = []
+Brmsx = []
+Brmsy = []
+Brmsz = []
+
+dump_folder = "_wave_dump"
+if shamrock.sys.world_rank() == 0:
+    os.makedirs(dump_folder, exist_ok=True)
+
+# %%
+# Time loop with data collection
+# ------------------------------
+t_sum = 0.0
+i_dump = 0
+next_dt_target = t_sum + dt_dump
+
+while next_dt_target <= t_target + 1e-12:
+    # Evolve until next dump time
+    model.evolve_until(next_dt_target)
+    t_now = model.get_time()
+
+    # Collect particle data from the context
+    data = ctx.collect_data()
+    h_arr = data["hpart"]
+    hfac = model.get_hfact()
+    rho = pmass * (hfac / h_arr) * (hfac / h_arr) * (hfac / h_arr)
+    Bx = data["B/rho"][:, 0] * rho
+    By = data["B/rho"][:, 1] * rho
+    Bz = data["B/rho"][:, 2] * rho
+
+    # Compute RMS values
+    rms_x = np.sqrt(np.mean(Bx**2))
+    rms_y = np.sqrt(np.mean(By**2))
+    rms_z = np.sqrt(np.mean(Bz**2))
+
+    times.append(t_now)
+    Brmsx.append(rms_x)
+    Brmsy.append(rms_y)
+    Brmsz.append(rms_z)
+
+    # (Optional) write VTK for visualisation
+    model.do_vtk_dump(os.path.join(dump_folder, f"wave_{i_dump:04d}.vtk"), True)
+
+    print(f"t = {t_now:.3f}, rms Bz = {rms_z:.5f}")
+
+    i_dump += 1
+    next_dt_target += dt_dump
+
+# Convert to numpy arrays
+times = np.array(times)
+Brmsx = np.array(Brmsx)
+Brmsy = np.array(Brmsy)
+Brmsz = np.array(Brmsz)
+
+# %%
+# Analytical solution
+# -------------------
+# For the damping of a sinusoidal Alfven wave with ambipolar diffusion,
+# the dispersion relation gives:
+#   omega = omega_R + i omega_I
+# with omega_I = - (k^2 etaAD)/2   (damping rate)
+# and   omega_R = 0.5 * sqrt( - (k^2 etaAD)^2 - 4 (k vA)^2 )
+# The z-component of B (the perturbed component) evolves as:
+#   Bz(t) = Bz(0) * |sin(omega_R t)| * exp(omega_I t)
+# where Bz(0) = (rho0 * v0 * Bx0) / (vA * sqrt(2))
+
+Lx = Lx_actual
+# Phantom's exact dispersion relation
+Bx0 = 1.0
+rho0 = 1.0
+C_ADc = 0.01  # ion-neutral coupling, same as Phantom
+vA = Bx0 / np.sqrt(rho0)  # no mu_0 since mu_0=1 in your code units
+v0 = 0.01 * vA
+k = 2 * np.pi / Lx_actual  # = 2*pi since Lx=1
+
+etaAD_cgs = C_ADc * vA * vA
+etaAD_si = etaAD_cgs * 4 * np.pi / c
+etaAD = etaAD_cgs
+
+quadb = (k) ** 2 * etaAD
+quadc = -((k * vA) ** 2)
+omegaI = -0.5 * quadb  # negative = damping
+omegaR = 0.5 * np.sqrt(-(quadb**2) - 4 * quadc)
+
+h0 = v0 * Bx0 / (vA * np.sqrt(2.0))  #  (4 * np.pi) *
+
+print(f"omegaR = {omegaR:.4f}, omegaI = {omegaI:.4f}, h0 = {h0:.4f}")
+
+time_th = np.linspace(0, t_target, 1000)
+theory = h0 * np.abs(np.sin(omegaR * time_th)) * np.exp(omegaI * time_th)
+
+# %%
+# L2 error against the analytical solution
+# Evaluated at the simulation dump times (dt=0.01, same interval as Phantom's, see Phantom paper section 5.7.1).
+theory_at_times = h0 * np.abs(np.sin(omegaR * times)) * np.exp(omegaI * times)
+l2_error = np.sqrt(np.mean((Brmsz - theory_at_times) ** 2))
+print(f"L2 error (rms Bz vs theory, dt={dt_dump}) = {l2_error:.3e}")
+
+np.savez(
+    os.path.join(dump_folder, "wave_damping_data.npz"),
+    times=times,
+    Brmsx=Brmsx,
+    Brmsy=Brmsy,
+    Brmsz=Brmsz,
+    theory_at_times=theory_at_times,
+    l2_error=l2_error,
+    dr=dr,
+    omegaR=omegaR,
+    omegaI=omegaI,
+    h0=h0,
+)
+
+# %%
+# Plot results
+
+fig, axs = plt.subplots(1, 4, figsize=(12, 8))
+
+axs[0].plot(times, Brmsz, "b-", linewidth=2)
+axs[0].plot(time_th, theory, "r--", linewidth=2)
+axs[0].set_xlabel("Time (s)")
+axs[0].set_ylabel("rms Bz")
+axs[0].grid(alpha=0.3)
+
+axs[1].plot(times, Brmsy, "g-", linewidth=2)
+axs[1].set_xlabel("Time (s)")
+axs[1].set_ylabel("rms By")
+axs[1].grid(alpha=0.3)
+
+axs[2].plot(times, Brmsx, "m-", linewidth=2)
+axs[2].set_xlabel("Time (s)")
+axs[2].set_ylabel("rms Bx")
+axs[2].grid(alpha=0.3)
+
+axs[3].plot(time_th, theory, "r--", linewidth=2)
+axs[3].set_xlabel("Time (s)")
+axs[3].set_ylabel("rms Bz theory")
+axs[3].grid(alpha=0.3)
+
+plt.tight_layout()
+plt.savefig(os.path.join(dump_folder, "wave_damping_analysis.png"), dpi=150)
+plt.show()
+
+print("Analysis completed. Results saved in", dump_folder)
diff --git a/examples/tests_ci/run_mhd_wavedamping.py b/examples/tests_ci/run_mhd_wavedamping.py
new file mode 100644
index 0000000000..1983b06330
--- /dev/null
+++ b/examples/tests_ci/run_mhd_wavedamping.py
@@ -0,0 +1,186 @@
+"""
+Wave damping test for the non-ideal MHD ambipolar diffusion
+================================================================
+
+(Choi et al. 2009, Wurster, Price & Ayliffe 2014, Wurster, Price & Bate 2016 section 4.1) at the same
+Resolution nx=128 as in the Phantom paper section 5.7.1.
+"""
+
+import numpy as np
+
+import shamrock
+
+shamrock.enable_experimental_features()
+
+# If we use the shamrock executable to run this script instead of the python interpreter,
+# we should not initialize the system as the shamrock executable needs to handle specific MPI logic
+if not shamrock.sys.is_initialized():
+    shamrock.change_loglevel(1)
+    shamrock.sys.init("0:0")
+
+# %%
+# Parameters of the test (matches examples/sph/run_mhd_wavedamping.py / Phantom nx=128)
+
+L2_ERROR_THRESHOLD = 7.5e-5
+
+Lx = 1.0  # box length
+dr = 1 / 128  # particle spacing (nx=128, matching Phantom's benchmark resolution)
+rho0 = 1.0  # initial density
+Bx0 = 1.0  # background field in x
+C_ADc = 0.01  # ambipolar diffusion coefficient (Phantom convention)
+cs = 1.0  # isothermal sound speed
+t_target = 5.0  # total simulation time
+dt_dump = 0.01  # dump interval (matches Phantom's L2-error sampling interval)
+
+# Unit system chosen so that mu_0 = 1 exactly in code units
+codeu = shamrock.UnitSystem(
+    unit_time=1.0,
+    unit_length=1.0,
+    unit_mass=1.2566370621219e-06,
+)
+ucte = shamrock.Constants(codeu)
+mu_0 = ucte.mu_0()  # = 1 in these units
+
+vA = Bx0 / np.sqrt(rho0)  # Alfven speed
+etaAD = C_ADc * vA * vA  # ambipolar diffusivity
+
+if shamrock.sys.world_rank() == 0:
+    print(f"mu_0 = {mu_0}, vA = {vA:.3f}, etaAD = {etaAD:.3e}")
+
+# %%
+# Create context and SPH model
+
+ctx = shamrock.Context()
+ctx.pdata_layout_new()
+model = shamrock.get_Model_SPH(context=ctx, vector_type="f64_3", sph_kernel="C4")
+
+# %%
+# Set up simulation configuration
+# All artificial dissipation terms are turned off (alpha_B, alpha_AV, beta_AV, sigma_mhd),
+# matching Phantom's stated methodology for this test.
+
+cfg = model.gen_default_config()
+cfg.set_units(codeu)
+cfg.set_artif_viscosity_None()
+cfg.set_NonIdealMHD(
+    sigma_mhd=0, sigma_u=0, etaO=0, etaH=0, etaAD=etaAD, alpha_B=0, alpha_AV=0, beta_AV=0
+)
+cfg.set_boundary_periodic()
+cfg.set_eos_isothermal(cs)
+cfg.print_status()
+model.set_solver_config(cfg)
+
+model.init_scheduler(int(1e6), 1)
+
+# %%
+# Generate particle distribution in an FCC lattice
+
+bmin = (-Lx / 2.0, -np.sqrt(3) / 4.0 * Lx, -np.sqrt(6) / 4.0 * Lx)
+bmax = (Lx / 2.0, np.sqrt(3) / 4.0 * Lx, np.sqrt(6) / 4.0 * Lx)
+
+bmin, bmax = model.get_ideal_fcc_box(dr, bmin, bmax)
+xm, ym, zm = bmin
+xM, yM, zM = bmax
+Lx_actual = xM - xm
+
+model.resize_simulation_box(bmin, bmax)
+model.add_cube_fcc_3d(dr, bmin, bmax)
+
+vol_b = (xM - xm) * (yM - ym) * (zM - zm)
+totmass = rho0 * vol_b
+pmass = model.total_mass_to_part_mass(totmass)
+model.set_particle_mass(pmass)
+
+# %%
+# Set initial conditions: background field in x, sinusoidal velocity perturbation in z.
+
+k = 2 * np.pi / Lx_actual
+v0 = 0.01 * vA
+
+
+def B_func(r):
+    return (Bx0, 0.0, 0.0)
+
+
+def vel_func(r):
+    x, y, z = r
+    vz = v0 * np.sin(k * (x - xm))
+    return (0.0, 0.0, vz)
+
+
+def u_func(r):
+    return 0.0
+
+
+model.set_field_value_lambda_f64_3("B/rho", B_func)
+model.set_field_value_lambda_f64_3("vxyz", vel_func)
+model.set_field_value_lambda_f64("uint", u_func)
+
+model.set_cfl_cour(0.3)
+model.set_cfl_force(0.25)
+
+# %%
+# Time loop with data collection
+
+times = []
+Brmsz = []
+
+t_sum = 0.0
+next_dt_target = t_sum + dt_dump
+
+while next_dt_target <= t_target + 1e-12:
+    model.evolve_until(next_dt_target)
+    t_now = model.get_time()
+
+    data = ctx.collect_data()
+    h_arr = data["hpart"]
+    hfac = model.get_hfact()
+    rho = pmass * (hfac / h_arr) * (hfac / h_arr) * (hfac / h_arr)
+    Bz = data["B/rho"][:, 2] * rho
+
+    rms_z = np.sqrt(np.mean(Bz**2))
+
+    times.append(t_now)
+    Brmsz.append(rms_z)
+
+    if shamrock.sys.world_rank() == 0:
+        print(f"t = {t_now:.3f}, rms Bz = {rms_z:.5f}")
+
+    next_dt_target += dt_dump
+
+times = np.array(times)
+Brmsz = np.array(Brmsz)
+
+# %%
+# Analytical solution (damped Alfven wave dispersion relation, see
+# examples/sph/run_mhd_wavedamping.py for the full derivation)
+
+k = 2 * np.pi / Lx_actual
+quadb = k**2 * etaAD
+quadc = -((k * vA) ** 2)
+omegaI = -0.5 * quadb
+omegaR = 0.5 * np.sqrt(-(quadb**2) - 4 * quadc)
+h0 = v0 * Bx0 / (vA * np.sqrt(2.0))
+
+if shamrock.sys.world_rank() == 0:
+    print(f"omegaR = {omegaR:.4f}, omegaI = {omegaI:.4f}, h0 = {h0:.4f}")
+
+theory_at_times = h0 * np.abs(np.sin(omegaR * times)) * np.exp(omegaI * times)
+l2_error = np.sqrt(np.mean((Brmsz - theory_at_times) ** 2))
+
+if shamrock.sys.world_rank() == 0:
+    print(f"L2 error (rms Bz vs theory) = {l2_error:.3e} (threshold = {L2_ERROR_THRESHOLD:.3e})")
+
+# %%
+# Check the L2 error against Phantom's published benchmark threshold
+
+to_raise = []
+
+if l2_error > L2_ERROR_THRESHOLD:
+    to_raise.append(
+        f"L2 error of rms Bz vs analytical solution is out of tolerance: "
+        f"{l2_error:.3e} > {L2_ERROR_THRESHOLD:.3e}"
+    )
+
+for to_raise_item in to_raise:
+    raise ValueError(to_raise_item)
diff --git a/src/shammodels/sph/CMakeLists.txt b/src/shammodels/sph/CMakeLists.txt
index 30838aadcb..6c4c5d2d52 100644
--- a/src/shammodels/sph/CMakeLists.txt
+++ b/src/shammodels/sph/CMakeLists.txt
@@ -37,6 +37,7 @@ set(Sources
     src/modules/UpdateDerivs.cpp
     src/modules/ComputeLoadBalanceValue.cpp
     src/modules/ComputeOmega.cpp
+    src/modules/ComputeJ.cpp
     src/modules/ComputeLuminosity.cpp
     src/modules/NeighbourCache.cpp
     src/modules/ParticleReordering.cpp
diff --git a/src/shammodels/sph/include/shammodels/sph/Solver.hpp b/src/shammodels/sph/include/shammodels/sph/Solver.hpp
index 1abecabcdf..582e896688 100644
--- a/src/shammodels/sph/include/shammodels/sph/Solver.hpp
+++ b/src/shammodels/sph/include/shammodels/sph/Solver.hpp
@@ -13,7 +13,7 @@
  * @file Solver.hpp
  * @author Léodasce Sewanou (leodasce.sewanou@ens-lyon.fr)
  * @author Timothée David--Cléris (tim.shamrock@proton.me)
- * @author Yona Lapeyre (yona.lapeyre@ens-lyon.fr) --no git blame--
+ * @author Yona Lapeyre (yona.lapeyre@ens-lyon.fr)
  * @brief
  */
 
@@ -154,7 +154,12 @@ namespace shammodels::sph {
         static constexpr u32 dim = shambase::VectorProperties<Tvec>::dimension;
         using Kernel             = SPHKernel<Tscal>;
 
-        using Config = SolverConfig<Tvec, SPHKernel>;
+        using Config  = SolverConfig<Tvec, SPHKernel>;
+        using Cfg_MHD = typename Config::MHDConfig;
+
+        using NoneMHD     = typename Cfg_MHD::None;
+        using IdealMHD    = typename Cfg_MHD::IdealMhdConstrainedHyperPara;
+        using NonIdealMHD = typename Cfg_MHD::NonIdealMHD;
 
         using u_morton = typename Config::u_morton;
 
@@ -323,6 +328,9 @@ namespace shammodels::sph {
         /// @brief Updates artificial viscosity coefficients for shock capturing
         void update_artificial_viscosity(Tscal dt);
 
+        /// @brief Updates the magnetic current field (for NIMHD)
+        void update_j();
+
         /// @brief Initializes data layout for ghost particle fields
         void init_ghost_layout();
 
@@ -341,6 +349,7 @@ namespace shammodels::sph {
         void prepare_corrector();
         /// @brief Updates time derivatives and applies external forces
         void update_derivs(Tscal dt_hydro);
+
         /**
          * @brief
          *
diff --git a/src/shammodels/sph/include/shammodels/sph/SolverConfig.hpp b/src/shammodels/sph/include/shammodels/sph/SolverConfig.hpp
index f3461f93d0..09dfa38bd3 100644
--- a/src/shammodels/sph/include/shammodels/sph/SolverConfig.hpp
+++ b/src/shammodels/sph/include/shammodels/sph/SolverConfig.hpp
@@ -71,12 +71,18 @@ namespace shammodels::sph {
         /**
          * @brief The CFL condition for the courant factor
          */
-        Tscal cfl_cour;
+        Tscal cfl_cour = 0.3;
 
         /**
          * @brief The CFL condition for the force
          */
-        Tscal cfl_force;
+        Tscal cfl_force = 0.25;
+
+        Tscal _pi = shambase::constants::pi<Tscal>;
+        /**
+         * @brief The CFL condition for the force
+         */
+        Tscal cfl_NIMHD = 1. / (2 * _pi); // as in phantom
 
         /**
          * @brief The CFL multiplier stiffness
@@ -552,7 +558,7 @@ struct shammodels::sph::SolverConfig {
     /// The radius of the sph kernel
     static constexpr Tscal Rkern = Kernel::Rkern;
 
-    Tscal gpart_mass; ///< The mass of each gas particle
+    Tscal gpart_mass{0}; ///< The mass of each gas particle (must be set before use)
 
     bool track_particles_id = false;
 
@@ -651,11 +657,27 @@ struct shammodels::sph::SolverConfig {
     }
 
     /// Enable the ideal MHD hydro solver
-    inline void set_IdealMHD(typename MHDConfig::IdealMHD_constrained_hyper_para v) {
+    inline void set_ideal_mhd(typename MHDConfig::IdealMhdConstrainedHyperPara v) {
         mhd_config.set(v);
     }
 
-    inline void set_NonIdealMHD(typename MHDConfig::NonIdealMHD v) { mhd_config.set(v); }
+    inline void set_non_ideal_mhd(typename MHDConfig::NonIdealMHD v) {
+        logger::raw_ln("$DANGER$DANGER$DANGER$DANGER$DANGER$DANGER$DANGER$DANGER$");
+        logger::raw_ln(" ______   _______  __    _  _______  _______  ______  ");
+        logger::raw_ln("|      | |   _   ||  |  | ||       ||       ||    _ | ");
+        logger::raw_ln("|  _    ||  |_|  ||   |_| ||    ___||    ___||   | ||");
+        logger::raw_ln("| | |   ||       ||       ||   | __ |   |___ |   |_||_");
+        logger::raw_ln("| |_|   ||       ||  _    ||   ||  ||    ___||    __  |");
+        logger::raw_ln("|       ||   _   || | |   ||   |_| ||   |___ |   |  | |");
+        logger::raw_ln("|______| |__| |__||_|  |__||_______||_______||___|  |_|");
+        logger::raw_ln("$DANGER$DANGER$DANGER$DANGER$DANGER$DANGER$DANGER$DANGER$");
+        logger::raw_ln("The Non-ideal MHD solver is UNDER DEVELOPMENT.");
+        logger::raw_ln("It is. NOT. FULLY. TESTED. YET.");
+        logger::raw_ln("Use at your own risk.");
+        shamrock::experimental_feature_check(
+            "Non-ideal MHD is experimental, please enable experimental features to use it");
+        mhd_config.set(v);
+    }
 
     //////////////////////////////////////////////////////////////////////////////////////////////
     // MHD Config (END)
@@ -1121,20 +1143,23 @@ struct shammodels::sph::SolverConfig {
         return artif_viscosity.has_field_soundspeed() || is_eos_locally_isothermal();
     }
 
+    /// @brief Whether the solver is set for non ideal MHD
+    inline bool do_nimhd() { return mhd_config.do_nimhd(); }
+
     /// @brief Whether the solver has a field for B_on_rho
-    inline bool has_field_B_on_rho() { return mhd_config.has_B_field() && (dim == 3); }
+    inline bool has_field_b_on_rho() { return mhd_config.has_b_field() && (dim == 3); }
 
     /// @brief Whether the solver has a field for psi_on_ch
     inline bool has_field_psi_on_ch() { return mhd_config.has_psi_field(); }
 
     /// @brief Whether the solver has a field for divB
-    inline bool has_field_divB() { return mhd_config.has_divB_field(); }
+    inline bool has_field_div_b() { return mhd_config.has_div_b_field(); }
 
     /// @brief Whether the solver has a field for curlB
-    inline bool has_field_curlB() { return mhd_config.has_curlB_field() && (dim == 3); }
+    inline bool has_field_curl_b() { return mhd_config.has_curl_b_field() && (dim == 3); }
 
     /// @brief Whether the solver has a field for dt divB
-    inline bool has_field_dtdivB() { return mhd_config.has_dtdivB_field(); }
+    inline bool has_field_dtdiv_b() { return mhd_config.has_dtdiv_b_field(); }
 
     /// @brief Whether to store luminosity
     bool compute_luminosity = false;
@@ -1159,6 +1184,7 @@ struct shammodels::sph::SolverConfig {
 
     inline void check_config() {
         dust_config.check_config();
+        mhd_config.check_config();
 
         if (track_particles_id && false /*particle injection when added*/) {
             shamrock::experimental_feature_check(
@@ -1173,6 +1199,11 @@ struct shammodels::sph::SolverConfig {
             shamrock::experimental_feature_check(
                 "Self gravity is experimental, please enable experimental features to use it");
         }
+
+        if (mhd_config.do_nimhd()) {
+            shamrock::experimental_feature_check(
+                "Non-ideal MHD is experimental, please enable experimental features to use it");
+        }
     }
 
     void set_layout(shamrock::patch::PatchDataLayerLayout &pdl);
diff --git a/src/shammodels/sph/include/shammodels/sph/config/MHDConfig.hpp b/src/shammodels/sph/include/shammodels/sph/config/MHDConfig.hpp
index 33daf37ee0..599a30be59 100644
--- a/src/shammodels/sph/include/shammodels/sph/config/MHDConfig.hpp
+++ b/src/shammodels/sph/include/shammodels/sph/config/MHDConfig.hpp
@@ -18,6 +18,7 @@
  */
 
 #include "shambackends/vec.hpp"
+#include "shamrock/experimental_features.hpp"
 #include "shamsys/legacy/log.hpp"
 #include <nlohmann/json.hpp>
 #include <variant>
@@ -36,74 +37,104 @@ struct shammodels::sph::MHDConfig {
 
     struct None {};
 
-    struct IdealMHD_constrained_hyper_para {
+    struct IdealMhdConstrainedHyperPara {
         Tscal sigma_mhd = 0.1;
         Tscal alpha_u   = 1.;
+        Tscal alpha_B   = 1.;
+        Tscal alpha_AV  = 1.;
+        Tscal beta_AV   = 1.;
     };
 
     struct NonIdealMHD {
         Tscal sigma_mhd = 0.1;
         Tscal alpha_u   = 1.;
+        Tscal alpha_B   = 1.;
+        Tscal alpha_AV  = 1.;
+        Tscal beta_AV   = 1.;
+        Tscal etaO      = 1.;
+        Tscal etaH      = 1.;
+        Tscal etaAD     = 1.;
     };
 
     // how to set a new state of a variant as a dummy:
     // a) do everything right
     // b) forget to add the state to the variant
     //-> question your life choices
-    using Variant = std::variant<None, IdealMHD_constrained_hyper_para, NonIdealMHD>;
+    using Variant = std::variant<None, IdealMhdConstrainedHyperPara, NonIdealMHD>;
 
-    Variant config = None{};
+    Variant configMHD = None{};
 
-    void set(Variant v) { config = v; }
+    void set(Variant v) { configMHD = v; }
 
-    inline bool has_B_field() {
-        bool is_B = bool(std::get_if<IdealMHD_constrained_hyper_para>(&config))
-                    || bool(std::get_if<NonIdealMHD>(&config));
+    inline bool do_nimhd() {
+        bool is_NIMHD = bool(std::get_if<NonIdealMHD>(&configMHD));
+        return is_NIMHD;
+    }
+
+    inline bool has_b_field() {
+        bool is_B = bool(std::get_if<IdealMhdConstrainedHyperPara>(&configMHD))
+                    || bool(std::get_if<NonIdealMHD>(&configMHD));
         return is_B;
     }
 
     inline bool has_psi_field() {
-        bool is_psi = bool(std::get_if<IdealMHD_constrained_hyper_para>(&config))
-                      || bool(std::get_if<NonIdealMHD>(&config));
+        bool is_psi = bool(std::get_if<IdealMhdConstrainedHyperPara>(&configMHD))
+                      || bool(std::get_if<NonIdealMHD>(&configMHD));
         return is_psi;
     }
 
-    inline bool has_divB_field() {
-        bool is_divB = bool(std::get_if<IdealMHD_constrained_hyper_para>(&config));
+    inline bool has_div_b_field() {
+        bool is_divB = bool(std::get_if<IdealMhdConstrainedHyperPara>(&configMHD));
         return is_divB;
     }
 
-    inline bool has_curlB_field() {
-        bool is_curlB = bool(std::get_if<NonIdealMHD>(&config));
+    inline bool has_curl_b_field() {
+        bool is_curlB = bool(std::get_if<NonIdealMHD>(&configMHD));
         return is_curlB;
     }
 
-    inline bool has_dtdivB_field() {
-        bool is_dtdivB = bool(std::get_if<NonIdealMHD>(&config));
+    inline bool has_dtdiv_b_field() {
+        bool is_dtdivB = bool(std::get_if<NonIdealMHD>(&configMHD));
         return is_dtdivB;
     }
 
     inline void print_status() {
-        logger::raw_ln("--- MHD config");
+        logger::raw_ln("--- MHD configMHD");
 
-        if (None *v = std::get_if<None>(&config)) {
+        if (None *v = std::get_if<None>(&configMHD)) {
             logger::raw_ln("  Config MHD Type : None (No MHD)");
         } else if (
-            IdealMHD_constrained_hyper_para *v
-            = std::get_if<IdealMHD_constrained_hyper_para>(&config)) {
+            IdealMhdConstrainedHyperPara *v
+            = std::get_if<IdealMhdConstrainedHyperPara>(&configMHD)) {
             logger::raw_ln("  Config MHD  : Ideal MHD, constrained hyperbolic/parabolic treatment");
             logger::raw_ln("  sigma_mhd  =", v->sigma_mhd);
-        } else if (NonIdealMHD *v = std::get_if<NonIdealMHD>(&config)) {
+            logger::raw_ln("  alpha_B    =", v->alpha_B);
+            logger::raw_ln("  alpha_AV   =", v->alpha_AV);
+            logger::raw_ln("  beta_AV    =", v->beta_AV);
+        } else if (NonIdealMHD *v = std::get_if<NonIdealMHD>(&configMHD)) {
             logger::raw_ln("  Config MHD Type : Non Ideal MHD");
             logger::raw_ln("  sigma_mhd   =", v->sigma_mhd);
+            logger::raw_ln("  alpha_B     =", v->alpha_B);
+            logger::raw_ln("  alpha_AV    =", v->alpha_AV);
+            logger::raw_ln("  beta_AV     =", v->beta_AV);
         } else {
             shambase::throw_unimplemented();
         }
 
-        logger::raw_ln("--- MHD config (deduced)");
+        logger::raw_ln("--- MHD configMHD (deduced)");
 
         logger::raw_ln("-------------");
     }
+
+    inline void check_config() {
+
+        if (do_nimhd()) {
+
+            if (!shamrock::are_experimental_features_allowed()) {
+                shambase::throw_with_loc<std::runtime_error>("Non Ideal MHD is experimental");
+            }
+        }
+    }
 };
 
 namespace shammodels::sph {
@@ -119,26 +150,35 @@ namespace shammodels::sph {
         using T = MHDConfig<Tvec>;
 
         using None        = typename T::None;
-        using IMHD        = typename T::IdealMHD_constrained_hyper_para;
+        using IMHD        = typename T::IdealMhdConstrainedHyperPara;
         using NonIdealMHD = typename T::NonIdealMHD;
 
-        // Write the config type into the JSON object
-        if (const None *v = std::get_if<None>(&p.config)) {
+        // Write the configMHD type into the JSON object
+        if (const None *v = std::get_if<None>(&p.configMHD)) {
             j = {
                 {"mhd_type", "none"},
             };
-        } else if (const IMHD *v = std::get_if<IMHD>(&p.config)) {
+        } else if (const IMHD *v = std::get_if<IMHD>(&p.configMHD)) {
             j = {
                 {"mhd_type", "ideal_mhd_constrained_hyper_para"},
                 {"sigma_mhd", v->sigma_mhd},
                 {"alpha_u", v->alpha_u},
+                {"alpha_B", v->alpha_B},
+                {"alpha_AV", v->alpha_AV},
+                {"beta_AV", v->beta_AV},
             };
-        } else if (const NonIdealMHD *v = std::get_if<NonIdealMHD>(&p.config)) {
+        } else if (const NonIdealMHD *v = std::get_if<NonIdealMHD>(&p.configMHD)) {
             // Write the shear base, direction, and speed into the JSON object
             j = {
                 {"mhd_type", "non_ideal_mhd"},
                 {"sigma_mhd", v->sigma_mhd},
                 {"alpha_u", v->alpha_u},
+                {"alpha_B", v->alpha_B},
+                {"alpha_AV", v->alpha_AV},
+                {"beta_AV", v->beta_AV},
+                {"etaO", v->etaO},
+                {"etaH", v->etaH},
+                {"etaAD", v->etaAD},
             };
         } else {
             shambase::throw_unimplemented();
@@ -162,15 +202,15 @@ namespace shammodels::sph {
             shambase::throw_with_loc<std::runtime_error>("no field mhd_type is found in this json");
         }
 
-        // Read the config type from the JSON object
+        // Read the configMHD type from the JSON object
         std::string mhd_type;
         j.at("mhd_type").get_to(mhd_type);
 
         using None        = typename T::None;
-        using IMHD        = typename T::IdealMHD_constrained_hyper_para;
+        using IMHD        = typename T::IdealMhdConstrainedHyperPara;
         using NonIdealMHD = typename T::NonIdealMHD;
 
-        // Set the BCConfig based on the config type
+        // Set the BCConfig based on the configMHD type
         if (mhd_type == "none") {
             p.set(None{});
         } else if (mhd_type == "ideal_mhd_constrained_hyper_para") {
@@ -178,12 +218,21 @@ namespace shammodels::sph {
                 IMHD{
                     j.at("sigma_mhd").get<Tscal>(),
                     j.at("alpha_u").get<Tscal>(),
+                    j.at("alpha_B").get<Tscal>(),
+                    j.at("alpha_AV").get<Tscal>(),
+                    j.at("beta_AV").get<Tscal>(),
                 });
         } else if (mhd_type == "non_ideal_mhd") {
             p.set(
                 NonIdealMHD{
                     j.at("sigma_mhd").get<Tscal>(),
                     j.at("alpha_u").get<Tscal>(),
+                    j.at("etaO").get<Tscal>(),
+                    j.at("etaH").get<Tscal>(),
+                    j.at("etaAD").get<Tscal>(),
+                    j.at("alpha_B").get<Tscal>(),
+                    j.at("alpha_AV").get<Tscal>(),
+                    j.at("beta_AV").get<Tscal>(),
                 });
         } else {
             shambase::throw_unimplemented("wtf !");
diff --git a/src/shammodels/sph/include/shammodels/sph/math/mhd.hpp b/src/shammodels/sph/include/shammodels/sph/math/mhd.hpp
index a45f62a4ae..6f50e8cc02 100644
--- a/src/shammodels/sph/include/shammodels/sph/math/mhd.hpp
+++ b/src/shammodels/sph/include/shammodels/sph/math/mhd.hpp
@@ -32,9 +32,98 @@ namespace shamrock::sph::mhd {
 
     enum MHDType { Ideal = 0, NonIdeal = 1 };
 
+    template<class Tvec, class Tscal, MHDType mhd_mode = NonIdeal>
+    inline Tvec mag_current_j_sum(
+        Tscal m_b, Tvec B_a, Tvec B_b, Tvec nabla_Wab_ha, Tscal sub_fact_a, Tscal mu_0) {
+
+        // J = curl(B)/mu_0 (mu_0 explicit, SI/Heaviside-Lorentz-like convention, not
+        // Gaussian-cgs 4*pi/c)
+
+        return m_b * sham::inv_sat_zero(sub_fact_a) * sycl::cross(B_a - B_b, nabla_Wab_ha) / mu_0;
+        // return {0., 0., 0.};
+    }
+
+    template<class Tvec, class Tscal, MHDType mhd_mode = NonIdeal>
+    inline Tvec wurster_d(Tvec B, Tvec J, Tscal etaO, Tscal etaH, Tscal etaAD, Tscal mu_0) {
+
+        Tvec Bhat  = B * sham::inv_sat_zero(sycl::length(B));
+        Tvec curlB = mu_0 * J; // diffusivities in L^2/T in any unit system
+        Tvec D     = etaO * curlB + etaH * sycl::cross(curlB, Bhat)
+                     - etaAD * sycl::cross(sycl::cross(curlB, Bhat), Bhat);
+
+        return D;
+    }
+
+    template<class Tvec, class Tscal, MHDType mhd_mode = NonIdeal>
+    inline Tscal u_ni_heating(
+        Tvec B, Tvec J, Tscal rho, Tscal etaO, Tscal etaH, Tscal etaAD, Tscal mu_0) {
+
+        // return sycl::dot(D, J) * sham::inv_sat_zero(rho);
+        // Tscal BdB       = sycl::dot(B, B);
+        // Tscal JdJ       = sycl::dot(J, J);
+        // Tscal BdJ       = sycl::dot(B, J);
+        // Tscal BdJBdJhat = sham::inv_sat_zero(BdB) * BdJ * BdJ;
+
+        // return (etaO * JdJ + etaAD * (JdJ - BdJBdJhat)) * sham::inv_sat_zero(rho); @ to check
+        Tvec D = wurster_d(B, J, etaO, etaH, etaAD, mu_0);
+        return sycl::dot(D, J) * sham::inv_sat_zero(rho);
+    }
+
+    template<class Tvec, class Tscal, MHDType mhd_mode = NonIdeal>
+    inline Tvec b_ni_terms(
+        Tvec D_a,
+        Tvec D_b,
+        Tscal m_b,
+        Tscal rho_a_sq,
+        Tscal rho_b_sq,
+        Tscal omega_a,
+        Tscal omega_b,
+        Tvec nabla_Wab_ha,
+        Tvec nabla_Wab_hb) {
+
+        Tscal sub_fact_a = rho_a_sq * omega_a;
+        Tscal sub_fact_b = rho_b_sq * omega_b;
+
+        Tvec acc_a = sham::inv_sat_zero(sub_fact_a) * (sycl::cross(D_a, nabla_Wab_ha));
+        Tvec acc_b = sham::inv_sat_zero(sub_fact_b) * (sycl::cross(D_b, nabla_Wab_hb));
+        return m_b * (acc_a + acc_b);
+    }
+
+    // not using Whurster D, developping with J. Equivalent to b_ni_terms
+    template<class Tvec, class Tscal, MHDType mhd_mode = NonIdeal>
+    inline Tvec b_ni_ad(
+        Tscal eta_AD,
+        Tvec J_a,
+        Tvec J_b,
+        Tscal m_b,
+        Tscal rho_a_sq,
+        Tscal rho_b_sq,
+        Tvec B_a,
+        Tvec B_b,
+        Tscal omega_a,
+        Tscal omega_b,
+        Tvec nabla_Wab_ha,
+        Tvec nabla_Wab_hb) {
+
+        Tscal sub_fact_a = rho_a_sq * omega_a;
+        Tscal sub_fact_b = rho_b_sq * omega_b;
+
+        Tvec Bhat_a = B_a * sham::inv_sat_zero(sycl::length(B_a));
+        Tvec Bhat_b = B_b * sham::inv_sat_zero(sycl::length(B_b));
+        Tvec JcB_a  = sycl::cross(J_a, Bhat_a);
+        Tvec JcB_b  = sycl::cross(J_b, Bhat_b);
+
+        Tvec JcBcB_a = sycl::cross(JcB_a, Bhat_a);
+        Tvec JcBcB_b = sycl::cross(JcB_b, Bhat_b);
+
+        Tvec acc_a = sham::inv_sat_zero(sub_fact_a) * eta_AD * sycl::cross(JcBcB_a, nabla_Wab_ha);
+        Tvec acc_b = sham::inv_sat_zero(sub_fact_b) * eta_AD * sycl::cross(JcBcB_b, nabla_Wab_hb);
+        return -m_b * (acc_a + acc_b);
+    }
+
     // mag tension form the Tricco 2023 formula
     template<class Tvec, class Tscal>
-    inline Tvec B_dot_grad_W(
+    inline Tvec b_dot_grad_w(
         Tscal m_b,
         Tscal rho_a_sq,
         Tscal rho_b_sq,
@@ -59,7 +148,7 @@ namespace shamrock::sph::mhd {
     }
 
     // from the Phantom paper formula
-    template<class Tvec, class Tscal, MHDType MHD_mode = Ideal>
+    template<class Tvec, class Tscal, MHDType mhd_mode = Ideal>
     inline Tvec mag_tension(
         Tscal m_b,
         Tvec B_a,
@@ -88,8 +177,8 @@ namespace shamrock::sph::mhd {
         return magnetic_tension_term;
     }
 
-    template<class Tvec, class Tscal, MHDType MHD_mode = Ideal>
-    inline Tvec fdivB(
+    template<class Tvec, class Tscal, MHDType mhd_mode = Ideal>
+    inline Tvec fdiv_b(
         Tscal m_b,
         Tvec B_a,
         Tvec B_b,
@@ -130,8 +219,8 @@ namespace shamrock::sph::mhd {
         return artres;
     }
 
-    template<class Tvec, class Tscal, MHDType MHD_mode = Ideal>
-    inline Tscal dB_on_rho_induction_term(
+    template<class Tvec, class Tscal, MHDType mhd_mode = Ideal>
+    inline Tscal d_b_on_rho_induction_term(
         Tscal m_b, Tscal rho_a_sq, Tvec B_a, Tscal omega_a, Tvec nabla_Wab_ha) {
 
         Tscal sub_fact_a = rho_a_sq * omega_a;
@@ -141,8 +230,8 @@ namespace shamrock::sph::mhd {
         return induction_term_no_vab;
     }
 
-    template<class Tvec, class Tscal, MHDType MHD_mode = Ideal>
-    inline Tvec dB_on_rho_psi_term(
+    template<class Tvec, class Tscal, MHDType mhd_mode = Ideal>
+    inline Tvec d_b_on_rho_psi_term(
         Tscal m_b,
         Tscal rho_a_sq,
         Tscal rho_b_sq,
@@ -159,12 +248,12 @@ namespace shamrock::sph::mhd {
         Tvec psisubterm_a = ((psi_a) *sham::inv_sat_zero(sub_fact_a)) * nabla_Wab_ha;
         Tvec psisubterm_b = ((psi_b) *sham::inv_sat_zero(sub_fact_b)) * nabla_Wab_hb;
 
-        Tvec psiterm = -m_b * (psisubterm_a + psisubterm_a);
+        Tvec psiterm = -m_b * (psisubterm_a + psisubterm_b);
 
         return psiterm;
     }
 
-    template<class Tvec, class Tscal, MHDType MHD_mode = Ideal>
+    template<class Tvec, class Tscal, MHDType mhd_mode = Ideal>
     inline Tscal dpsi_on_ch_parabolic_propag(
         Tscal m_b, Tscal rho_a, Tvec B_a, Tvec B_b, Tscal omega_a, Tvec nabla_Wab_ha, Tscal ch_a) {
 
@@ -179,7 +268,7 @@ namespace shamrock::sph::mhd {
         return parabolic_propag;
     }
 
-    template<class Tvec, class Tscal, MHDType MHD_mode = Ideal>
+    template<class Tvec, class Tscal, MHDType mhd_mode = Ideal>
     inline Tscal dpsi_on_ch_parabolic_diff(
         Tscal m_b,
         Tscal rho_a,
@@ -197,7 +286,7 @@ namespace shamrock::sph::mhd {
         return parabolic_diff;
     }
 
-    template<class Kernel, class Tvec, class Tscal, MHDType MHD_mode = Ideal>
+    template<class Kernel, class Tvec, class Tscal, MHDType mhd_mode = Ideal>
     inline void add_to_derivs_spmhd(
         Tscal pmass,
         Tvec dr,
@@ -223,16 +312,26 @@ namespace shamrock::sph::mhd {
         Tscal h_b,
 
         Tscal alpha_u,
+        Tscal alpha_B,
+        Tscal alpha_AV,
+        Tscal beta_AV,
 
         Tvec B_a,
         Tvec B_b,
 
+        Tvec J_a,
+        Tvec J_b,
+
         Tscal psi_a,
         Tscal psi_b,
 
         Tscal mu_0,
         Tscal sigma_mhd,
 
+        Tscal etaO,
+        Tscal etaH,
+        Tscal etaAD,
+
         Tvec &dv_dt,
         Tscal &du_dt,
         Tvec &dB_on_rho_dt,
@@ -260,9 +359,9 @@ namespace shamrock::sph::mhd {
 
         Tscal vsig_u = shamrock::sph::vsig_u(P_a, P_b, rho_a, rho_b);
         Tscal vsig_a = shamphys::MHD_physics<Tvec, Tscal>::vsig_MHD(
-            v_ab, r_ab_unit, cs_a, B_a, rho_a, mu_0, 1., 1.);
+            v_ab, r_ab_unit, cs_a, B_a, rho_a, mu_0, alpha_AV, beta_AV);
         Tscal vsig_b = shamphys::MHD_physics<Tvec, Tscal>::vsig_MHD(
-            v_ab, r_ab_unit, cs_a, B_b, rho_b, mu_0, 1., 1.);
+            v_ab, r_ab_unit, cs_b, B_b, rho_b, mu_0, alpha_AV, beta_AV);
 
         Tscal dWab_a = Fab_a;
         Tscal dWab_b = Fab_b;
@@ -290,7 +389,7 @@ namespace shamrock::sph::mhd {
         // update_derivs)
 
         // dv/dt terms
-        sum_fdivB += fdivB(
+        sum_fdivB += fdiv_b(
             pmass, B_a, B_b, r_ab_unit * dWab_a, r_ab_unit * dWab_b, sub_fact_a, sub_fact_b, mu_0);
 
         Tvec gas_pressure_pishock = sph::sph_pressure_symetric(
@@ -304,7 +403,7 @@ namespace shamrock::sph::mhd {
             r_ab_unit * dWab_a,
             r_ab_unit * dWab_b);
 
-        sum_mag_tension += -B_dot_grad_W(
+        sum_mag_tension += -b_dot_grad_w(
             pmass,
             rho_a_sq,
             rho_b * rho_b,
@@ -352,8 +451,18 @@ namespace shamrock::sph::mhd {
             dWab_a * omega_a_rho_a_inv,
             dWab_b / (rho_b * omega_b));
 
-        du_dt += lambda_artes(
-            pmass, rho_a_sq, rho_b * rho_b, vsig_B, B_a, B_b, omega_a, omega_b, Fab_a, Fab_b);
+        du_dt += alpha_B
+                 * lambda_artes(
+                     pmass,
+                     rho_a_sq,
+                     rho_b * rho_b,
+                     vsig_B,
+                     B_a,
+                     B_b,
+                     omega_a,
+                     omega_b,
+                     Fab_a,
+                     Fab_b);
 
         // end du/dt terms
 
@@ -364,10 +473,11 @@ namespace shamrock::sph::mhd {
         Tvec dB_on_rho_dissipation_term
             = 0.5 * pmass * (rho_diss_term_a + rho_diss_term_b) * (B_a - B_b) * vsig_B;
 
-        dB_on_rho_dt
-            += v_ab * dB_on_rho_induction_term(pmass, rho_a_sq, B_a, omega_a, r_ab_unit * dWab_b);
+        dB_on_rho_dt += v_ab
+                        * d_b_on_rho_induction_term(
+                            pmass, rho_a_sq, B_a, omega_a, r_ab_unit * dWab_a); // @@@ dWab_b ?
 
-        dB_on_rho_dt += dB_on_rho_psi_term(
+        dB_on_rho_dt += d_b_on_rho_psi_term(
             pmass,
             rho_a_sq,
             rho_b * rho_b,
@@ -378,7 +488,7 @@ namespace shamrock::sph::mhd {
             r_ab_unit * dWab_a,
             r_ab_unit * dWab_b);
 
-        dB_on_rho_dt += dB_on_rho_dissipation_term;
+        dB_on_rho_dt += alpha_B * dB_on_rho_dissipation_term;
 
         // end d(B/rho)/dt terms
 
@@ -401,6 +511,46 @@ namespace shamrock::sph::mhd {
 
         // for conservative checks
         drho_dt += (1. / omega_a) * pmass * sycl::dot(v_ab, r_ab_unit * dWab_a);
+
+        // Non-ideal MHD terms
+        if constexpr (mhd_mode == NonIdeal) {
+
+            Tvec D_a = wurster_d<Tvec, Tscal, mhd_mode>(B_a, J_a, etaO, etaH, etaAD, mu_0);
+            Tvec D_b = wurster_d<Tvec, Tscal, mhd_mode>(B_b, J_b, etaO, etaH, etaAD, mu_0);
+
+            Tvec B_NI = b_ni_terms<Tvec, Tscal, mhd_mode>(
+                D_a,
+                D_b,
+                pmass,
+                rho_a_sq,
+                rho_b * rho_b,
+                omega_a,
+                omega_b,
+                r_ab_unit * dWab_a,
+                r_ab_unit * dWab_b);
+
+            dB_on_rho_dt += B_NI;
+
+            // Tvec B_NI_ADterm = b_ni_ad<Tvec, Tscal, mhd_mode>(
+            //     etaAD,
+            //     J_a,
+            //     J_b,
+            //     pmass,
+            //     rho_a_sq,
+            //     rho_b * rho_b,
+            //     B_a,
+            //     B_b,
+            //     omega_a,
+            //     omega_b,
+            //     r_ab_unit * dWab_a,
+            //     r_ab_unit * dWab_b);
+
+            // dB_on_rho_dt += B_NI_ADterm;
+
+            // Tscal u_NI = u_ni_heating<Tvec, Tscal, mhd_mode>(B_a, J_a, rho_a, etaO, etaAD) * 0.5
+            //              + u_ni_heating<Tvec, Tscal, mhd_mode>(B_b, J_b, rho_b, etaO, etaAD) *
+            //              0.5;
+        }
     }
 
 } // namespace shamrock::sph::mhd
diff --git a/src/shammodels/sph/include/shammodels/sph/modules/ComputeCFLNIMHD.hpp b/src/shammodels/sph/include/shammodels/sph/modules/ComputeCFLNIMHD.hpp
new file mode 100644
index 0000000000..e193fa66ae
--- /dev/null
+++ b/src/shammodels/sph/include/shammodels/sph/modules/ComputeCFLNIMHD.hpp
@@ -0,0 +1,76 @@
+// -------------------------------------------------------//
+//
+// SHAMROCK code for hydrodynamics
+// Copyright (c) 2021-2026 Timothée David--Cléris <tim.shamrock@proton.me>
+// SPDX-License-Identifier: CeCILL Free Software License Agreement v2.1
+// Shamrock is licensed under the CeCILL 2.1 License, see LICENSE for more information
+//
+// -------------------------------------------------------//
+
+#pragma once
+
+/**
+ * @file ComputeCFLNIMHD.hpp
+ * @author Yona Lapeyre (yona.lapeyre@ens-lyon.fr)
+ * @brief
+ *
+ */
+
+#include "shambackends/kernel_call_distrib.hpp"
+#include "shamrock/solvergraph/IFieldSpan.hpp"
+#include "shamrock/solvergraph/Indexes.hpp"
+#include "shamsolvergraph/edge/IDataEdge.hpp"
+#include "shamsolvergraph/node/INode.hpp"
+#include "shamsys/NodeInstance.hpp"
+
+#define NODE_EDGES(X_RO, X_RW)                                                                     \
+    X_RO(shamrock::solvergraph::Indexes<u32>, part_counts)                                         \
+    X_RO(shamrock::solvergraph::IDataEdge<Tscal>, C_nimhd)                                         \
+    X_RO(shamrock::solvergraph::IDataEdge<Tscal>, eta_o)                                           \
+    X_RO(shamrock::solvergraph::IDataEdge<Tscal>, eta_ad)                                          \
+    X_RO(shamrock::solvergraph::IDataEdge<Tscal>, eta_h)                                           \
+    X_RO(shamrock::solvergraph::IFieldSpan<Tscal>, hpart)                                          \
+    X_RW(shamrock::solvergraph::IFieldSpan<Tscal>, cfl_dt)
+
+template<class Tvec>
+class ComputeCFLNIMHD : public shamrock::solvergraph::INode {
+
+    using Tscal = shambase::VecComponent<Tvec>;
+
+    public:
+    ComputeCFLNIMHD() {}
+
+    EXPAND_NODE_EDGES(NODE_EDGES)
+
+    inline void _impl_evaluate_internal() {
+        auto edges = get_edges();
+
+        auto dev_sched = shamsys::instance::get_compute_scheduler_ptr();
+
+        Tscal C_nimhd = edges.C_nimhd.data;
+        Tscal eta_o   = edges.eta_o.data;
+        Tscal eta_ad  = edges.eta_ad.data;
+        Tscal eta_h   = edges.eta_h.data;
+
+        sham::distributed_data_kernel_call(
+            dev_sched,
+            sham::DDMultiRef{edges.hpart.get_spans()},
+            sham::DDMultiRef{edges.cfl_dt.get_spans()},
+            edges.part_counts.indexes,
+            [C_nimhd, eta_o, eta_ad, eta_h](u32 id_a, const Tscal *hpart, Tscal *cfl_dt) {
+                Tscal h_a     = hpart[id_a];
+                Tscal max_eta = sycl::max(
+                    sycl::max(sycl::fabs(eta_o), sycl::fabs(eta_ad)), sycl::fabs(eta_h));
+
+                Tscal dt_nimhd = C_nimhd * h_a * h_a / max_eta;
+
+                cfl_dt[id_a] = sycl::min(cfl_dt[id_a], dt_nimhd);
+            });
+    }
+
+    inline virtual std::string _impl_get_label() const { return "ComputeCFLNIMHD"; };
+
+    inline virtual std::string _impl_get_tex() const { return "C_{NIMHD}"; };
+};
+
+#undef NODE_EDGES
diff --git a/src/shammodels/sph/include/shammodels/sph/modules/ComputeJ.hpp b/src/shammodels/sph/include/shammodels/sph/modules/ComputeJ.hpp
new file mode 100644
index 0000000000..60cf6b81d6
--- /dev/null
+++ b/src/shammodels/sph/include/shammodels/sph/modules/ComputeJ.hpp
@@ -0,0 +1,79 @@
+// -------------------------------------------------------//
+//
+// SHAMROCK code for hydrodynamics
+// Copyright (c) 2021-2026 Timothée David--Cléris <tim.shamrock@proton.me>
+// SPDX-License-Identifier: CeCILL Free Software License Agreement v2.1
+// Shamrock is licensed under the CeCILL 2.1 License, see LICENSE for more information
+//
+// -------------------------------------------------------//
+
+#pragma once
+
+/**
+ * @file ComputeJ.hpp
+ * @author Timothée David--Cléris (tim.shamrock@proton.me)
+ * @brief
+ *
+ */
+
+#include "shambackends/typeAliasVec.hpp"
+#include "shambackends/vec.hpp"
+#include "shammodels/sph/SolverConfig.hpp"
+#include "shammodels/sph/modules/SolverStorage.hpp"
+#include "shamrock/scheduler/ShamrockCtx.hpp"
+
+namespace shammodels::sph::modules {
+
+    template<class Tvec, template<class> class SPHKernel>
+    class NodeComputeJ : public shamrock::solvergraph::INode {
+
+        using Tscal = shambase::VecComponent<Tvec>;
+
+        Tscal part_mass;
+        Tscal mu_0;
+
+        public:
+        NodeComputeJ(Tscal part_mass, Tscal mu_0) : part_mass(part_mass), mu_0(mu_0) {}
+
+        struct Edges {
+            const shamrock::solvergraph::Indexes<u32> &part_counts;
+            const shammodels::sph::solvergraph::NeighCache &neigh_cache;
+            const shamrock::solvergraph::IFieldSpan<Tvec> &xyz;
+            const shamrock::solvergraph::IFieldSpan<Tscal> &hpart;
+            const shamrock::solvergraph::IFieldSpan<Tscal> &omega;
+            const shamrock::solvergraph::IFieldSpan<Tvec> &B_on_rho;
+            shamrock::solvergraph::IFieldSpan<Tvec> &J;
+        };
+
+        inline void set_edges(
+            std::shared_ptr<shamrock::solvergraph::Indexes<u32>> part_counts,
+            std::shared_ptr<shammodels::sph::solvergraph::NeighCache> neigh_cache,
+            std::shared_ptr<shamrock::solvergraph::IFieldSpan<Tvec>> xyz,
+            std::shared_ptr<shamrock::solvergraph::IFieldSpan<Tscal>> hpart,
+            std::shared_ptr<shamrock::solvergraph::IFieldSpan<Tscal>> omega,
+            std::shared_ptr<shamrock::solvergraph::IFieldSpan<Tvec>> B_on_rho,
+            std::shared_ptr<shamrock::solvergraph::IFieldSpan<Tvec>> J) {
+            __internal_set_ro_edges({part_counts, neigh_cache, xyz, hpart, omega, B_on_rho});
+            __internal_set_rw_edges({J});
+        }
+
+        inline Edges get_edges() {
+            return Edges{
+                get_ro_edge<shamrock::solvergraph::Indexes<u32>>(0),
+                get_ro_edge<shammodels::sph::solvergraph::NeighCache>(1),
+                get_ro_edge<shamrock::solvergraph::IFieldSpan<Tvec>>(2),
+                get_ro_edge<shamrock::solvergraph::IFieldSpan<Tscal>>(3),
+                get_ro_edge<shamrock::solvergraph::IFieldSpan<Tscal>>(4),
+                get_ro_edge<shamrock::solvergraph::IFieldSpan<Tvec>>(5),
+                get_rw_edge<shamrock::solvergraph::IFieldSpan<Tvec>>(0),
+            };
+        }
+
+        void _impl_evaluate_internal();
+
+        inline virtual std::string _impl_get_label() const { return "ComputeJ"; };
+
+        virtual std::string _impl_get_tex() const;
+    };
+
+} // namespace shammodels::sph::modules
diff --git a/src/shammodels/sph/include/shammodels/sph/modules/SolverStorage.hpp b/src/shammodels/sph/include/shammodels/sph/modules/SolverStorage.hpp
index c6acb50948..0af2f1db81 100644
--- a/src/shammodels/sph/include/shammodels/sph/modules/SolverStorage.hpp
+++ b/src/shammodels/sph/include/shammodels/sph/modules/SolverStorage.hpp
@@ -104,6 +104,10 @@ namespace shammodels::sph {
         Component<shamrock::ComputeField<Tvec>> old_dB_on_rho;
         Component<shamrock::ComputeField<Tscal>> old_dpsi_on_ch;
 
+        Component<shambase::DistributedData<PatchDataField<Tvec>>> MagCurrentJ_ghost;
+        std::shared_ptr<shamrock::solvergraph::Field<Tvec>> MagCurrentJ;
+        std::shared_ptr<shamrock::solvergraph::ExchangeGhostField<Tvec>> exchange_gz_J;
+
         Component<shamrock::ComputeField<Tscal>> old_dtepsilon;
         Component<shamrock::ComputeField<Tvec>> old_dtdeltav;
         Component<shamrock::ComputeField<Tscal>> old_ds_j_dt;
diff --git a/src/shammodels/sph/include/shammodels/sph/modules/UpdateDerivs.hpp b/src/shammodels/sph/include/shammodels/sph/modules/UpdateDerivs.hpp
index 8f45371438..6e66e95cae 100644
--- a/src/shammodels/sph/include/shammodels/sph/modules/UpdateDerivs.hpp
+++ b/src/shammodels/sph/include/shammodels/sph/modules/UpdateDerivs.hpp
@@ -11,6 +11,7 @@
 
 /**
  * @file UpdateDerivs.hpp
+ * @author Yona Lapeyre (yona.lapeyre@ens-lyon.fr)
  * @author Timothée David--Cléris (tim.shamrock@proton.me)
  * @brief
  *
@@ -19,6 +20,7 @@
 #include "shambackends/typeAliasVec.hpp"
 #include "shambackends/vec.hpp"
 #include "shammodels/sph/SolverConfig.hpp"
+#include "shammodels/sph/math/mhd.hpp"
 #include "shammodels/sph/modules/SolverStorage.hpp"
 #include "shamrock/scheduler/ShamrockCtx.hpp"
 
@@ -67,10 +69,28 @@ namespace shammodels::sph::modules {
         using Cfg_MHD = typename Config::MHDConfig;
 
         using NoneMHD     = typename Cfg_MHD::None;
-        using IdealMHD    = typename Cfg_MHD::IdealMHD_constrained_hyper_para;
+        using IdealMHD    = typename Cfg_MHD::IdealMhdConstrainedHyperPara;
         using NonIdealMHD = typename Cfg_MHD::NonIdealMHD;
 
-        void update_derivs_MHD(IdealMHD cfg);
+        template<shamrock::sph::mhd::MHDType mhd_mode>
+        void compute_j(Tscal mu_0);
+
+        // void update_derivs_mhd(Cfg_MHD cfg);
+        //  One templated implementation, specialised per MHDType at the call sites below.
+        template<shamrock::sph::mhd::MHDType mhd_mode>
+        void update_derivs_mhd_impl(
+            Tscal sigma_mhd,
+            Tscal alpha_u,
+            Tscal alpha_B,
+            Tscal alpha_AV,
+            Tscal beta_AV,
+            Tscal etaO,
+            Tscal etaH,
+            Tscal etaAD);
+
+        // Thin wrappers that unpack the variant and forward to the template above.
+        void update_derivs_mhd(IdealMHD cfg);
+        void update_derivs_mhd(NonIdealMHD cfg);
     };
 
 } // namespace shammodels::sph::modules
diff --git a/src/shammodels/sph/src/Solver.cpp b/src/shammodels/sph/src/Solver.cpp
index dc5d2db72c..35b5e86b01 100644
--- a/src/shammodels/sph/src/Solver.cpp
+++ b/src/shammodels/sph/src/Solver.cpp
@@ -53,8 +53,10 @@
 #include "shammodels/sph/modules/ComputeCFLDust1Fluid.hpp"
 #include "shammodels/sph/modules/ComputeCFLDustDrift.hpp"
 #include "shammodels/sph/modules/ComputeCFLForce.hpp"
+#include "shammodels/sph/modules/ComputeCFLNIMHD.hpp"
 #include "shammodels/sph/modules/ComputeCFLSinkSink.hpp"
 #include "shammodels/sph/modules/ComputeEos.hpp"
+#include "shammodels/sph/modules/ComputeJ.hpp"
 #include "shammodels/sph/modules/ComputeLoadBalanceValue.hpp"
 #include "shammodels/sph/modules/ComputeLuminosity.hpp"
 #include "shammodels/sph/modules/ComputeNeighStats.hpp"
@@ -267,7 +269,7 @@ void shammodels::sph::Solver<Tvec, Kern>::init_solver_graph() {
     auto &sync_data = sched.synchronized_data;
 
     shamrock::patch::PatchDataLayerLayout &pdl = scheduler().pdl_old();
-    bool has_B_field                           = solver_config.has_field_B_on_rho();
+    bool has_b_field                           = solver_config.has_field_b_on_rho();
     bool has_psi_field                         = solver_config.has_field_psi_on_ch();
     bool has_epsilon_field                     = solver_config.dust_config.has_epsilon_field();
     bool has_deltav_field                      = solver_config.dust_config.has_deltav_field();
@@ -291,7 +293,7 @@ void shammodels::sph::Solver<Tvec, Kern>::init_solver_graph() {
     solver_graph.register_edge("duint", FieldRefs<Tscal>("duint", "du_{\\rm int}"));
     solver_graph.register_edge("hpart", FieldRefs<Tscal>("hpart", "h_{\\rm part}"));
 
-    if (has_B_field) {
+    if (has_b_field) {
         solver_graph.register_edge("B/rho", FieldRefs<Tvec>("B/rho", "B_{\\rho}"));
         solver_graph.register_edge("dB/rho", FieldRefs<Tvec>("dB/rho", "dB_{\\rho}"));
     }
@@ -416,7 +418,7 @@ void shammodels::sph::Solver<Tvec, Kern>::init_solver_graph() {
             attach_field_sequence.push_back(attach_hpart);
         }
 
-        if (has_B_field) {
+        if (has_b_field) {
             auto attach_B_on_rho = solver_graph.register_node(
                 "attach_B_on_rho", GetFieldRefFromLayer<Tvec>(pdl, "B/rho"));
             shambase::get_check_ref(attach_B_on_rho)
@@ -426,7 +428,7 @@ void shammodels::sph::Solver<Tvec, Kern>::init_solver_graph() {
             attach_field_sequence.push_back(attach_B_on_rho);
         }
 
-        if (has_B_field) {
+        if (has_b_field) {
             auto attach_dB_on_rho = solver_graph.register_node(
                 "attach_dB_on_rho", GetFieldRefFromLayer<Tvec>(pdl, "dB/rho"));
             shambase::get_check_ref(attach_dB_on_rho)
@@ -553,7 +555,7 @@ void shammodels::sph::Solver<Tvec, Kern>::init_solver_graph() {
                 half_step_sequence.push_back(half_step_uint);
             }
 
-            if (has_B_field) {
+            if (has_b_field) {
                 auto half_step_B_on_rho = solver_graph.register_node(
                     prefix + "_B_on_rho", shammodels::common::modules::ForwardEuler<Tvec>{});
                 shambase::get_check_ref(half_step_B_on_rho)
@@ -775,6 +777,9 @@ void shammodels::sph::Solver<Tvec, Kern>::init_solver_graph() {
 
     storage.omega = std::make_shared<shamrock::solvergraph::Field<Tscal>>(1, "omega", "\\Omega");
 
+    storage.MagCurrentJ
+        = std::make_shared<shamrock::solvergraph::Field<Tvec>>(1, "MagCurrentJ", "\\mathbf{J}");
+
     if (solver_config.has_field_alphaAV()) {
         storage.alpha_av_updated = std::make_shared<shamrock::solvergraph::Field<Tscal>>(
             1, "alpha_av_updated", "\\alpha_{\\rm AV}");
@@ -791,6 +796,8 @@ void shammodels::sph::Solver<Tvec, Kern>::init_solver_graph() {
     storage.exchange_gz_positions
         = std::make_shared<shamrock::solvergraph::ExchangeGhostLayer>(storage.xyzh_ghost_layout);
 
+    storage.exchange_gz_J = std::make_shared<shamrock::solvergraph::ExchangeGhostField<Tvec>>();
+
     ////////////////////////////////////////////////////////////////////////////////////////
     // sink accretion
     ////////////////////////////////////////////////////////////////////////////////////////
@@ -1928,9 +1935,9 @@ void shammodels::sph::Solver<Tvec, Kern>::communicate_merge_ghosts_fields() {
     bool has_alphaAV_field    = solver_config.has_field_alphaAV();
     bool has_soundspeed_field = solver_config.ghost_has_soundspeed();
 
-    bool has_B_field       = solver_config.has_field_B_on_rho();
+    bool has_b_field       = solver_config.has_field_b_on_rho();
     bool has_psi_field     = solver_config.has_field_psi_on_ch();
-    bool has_curlB_field   = solver_config.has_field_curlB();
+    bool has_curl_b_field  = solver_config.has_field_curl_b();
     bool has_epsilon_field = solver_config.dust_config.has_epsilon_field();
     bool has_deltav_field  = solver_config.dust_config.has_deltav_field();
     bool has_s_j_field     = solver_config.dust_config.has_s_j_field();
@@ -1946,11 +1953,11 @@ void shammodels::sph::Solver<Tvec, Kern>::communicate_merge_ghosts_fields() {
     const u32 ialpha_AV   = (has_alphaAV_field) ? pdl.get_field_idx<Tscal>("alpha_AV") : 0;
     const u32 isoundspeed = (has_soundspeed_field) ? pdl.get_field_idx<Tscal>("soundspeed") : 0;
 
-    const u32 iB_on_rho   = (has_B_field) ? pdl.get_field_idx<Tvec>("B/rho") : 0;
-    const u32 idB_on_rho  = (has_B_field) ? pdl.get_field_idx<Tvec>("dB/rho") : 0;
+    const u32 iB_on_rho   = (has_b_field) ? pdl.get_field_idx<Tvec>("B/rho") : 0;
+    const u32 idB_on_rho  = (has_b_field) ? pdl.get_field_idx<Tvec>("dB/rho") : 0;
     const u32 ipsi_on_ch  = (has_psi_field) ? pdl.get_field_idx<Tscal>("psi/ch") : 0;
     const u32 idpsi_on_ch = (has_psi_field) ? pdl.get_field_idx<Tscal>("dpsi/ch") : 0;
-    const u32 icurlB      = (has_curlB_field) ? pdl.get_field_idx<Tvec>("curlB") : 0;
+    const u32 icurlB      = (has_curl_b_field) ? pdl.get_field_idx<Tvec>("curlB") : 0;
 
     bool do_MHD_debug       = solver_config.do_MHD_debug();
     const u32 imag_pressure = (do_MHD_debug) ? pdl.get_field_idx<Tvec>("mag_pressure") : -1;
@@ -1979,9 +1986,9 @@ void shammodels::sph::Solver<Tvec, Kern>::communicate_merge_ghosts_fields() {
     const u32 isoundspeed_interf
         = (has_soundspeed_field) ? ghost_layout.get_field_idx<Tscal>("soundspeed") : 0;
 
-    const u32 iB_interf     = (has_B_field) ? ghost_layout.get_field_idx<Tvec>("B/rho") : 0;
+    const u32 iB_interf     = (has_b_field) ? ghost_layout.get_field_idx<Tvec>("B/rho") : 0;
     const u32 ipsi_interf   = (has_psi_field) ? ghost_layout.get_field_idx<Tscal>("psi/ch") : 0;
-    const u32 icurlB_interf = (has_curlB_field) ? ghost_layout.get_field_idx<Tvec>("curlB") : 0;
+    const u32 icurlB_interf = (has_curl_b_field) ? ghost_layout.get_field_idx<Tvec>("curlB") : 0;
 
     const u32 iepsilon_interf
         = (has_epsilon_field) ? ghost_layout.get_field_idx<Tscal>("epsilon") : 0;
@@ -2035,7 +2042,7 @@ void shammodels::sph::Solver<Tvec, Kern>::communicate_merge_ghosts_fields() {
                     .append_subset_to(buf_idx, cnt, pdat.get_field<Tscal>(isoundspeed_interf));
             }
 
-            if (has_B_field) {
+            if (has_b_field) {
                 sender_patch.get_field<Tvec>(iB_on_rho).append_subset_to(
                     buf_idx, cnt, pdat.get_field<Tvec>(iB_interf));
             }
@@ -2045,7 +2052,7 @@ void shammodels::sph::Solver<Tvec, Kern>::communicate_merge_ghosts_fields() {
                     .append_subset_to(buf_idx, cnt, pdat.get_field<Tscal>(ipsi_interf));
             }
 
-            if (has_curlB_field) {
+            if (has_curl_b_field) {
                 sender_patch.get_field<Tvec>(icurlB).append_subset_to(
                     buf_idx, cnt, pdat.get_field<Tvec>(icurlB_interf));
             }
@@ -2118,7 +2125,7 @@ void shammodels::sph::Solver<Tvec, Kern>::communicate_merge_ghosts_fields() {
                         .insert(pdat.get_field<Tscal>(isoundspeed));
                 }
 
-                if (has_B_field) {
+                if (has_b_field) {
                     pdat_new.get_field<Tvec>(iB_interf).insert(pdat.get_field<Tvec>(iB_on_rho));
                 }
 
@@ -2127,7 +2134,7 @@ void shammodels::sph::Solver<Tvec, Kern>::communicate_merge_ghosts_fields() {
                         .insert(pdat.get_field<Tscal>(ipsi_on_ch));
                 }
 
-                if (has_curlB_field) {
+                if (has_curl_b_field) {
                     pdat_new.get_field<Tvec>(icurlB_interf).insert(pdat.get_field<Tvec>(icurlB));
                 }
 
@@ -2172,6 +2179,63 @@ void shammodels::sph::Solver<Tvec, Kern>::update_artificial_viscosity(Tscal dt)
         .update_artificial_viscosity(dt);
 }
 
+template<class Tvec, template<class> class Kern>
+void shammodels::sph::Solver<Tvec, Kern>::update_j() {
+
+    using namespace shamrock::patch;
+    PatchDataLayerLayout &pdl = scheduler().pdl_old();
+
+    const u32 iB_on_rho = pdl.get_field_idx<Tvec>("B/rho");
+    std::shared_ptr<shamrock::solvergraph::FieldRefs<Tvec>> B_on_rho_edge
+        = std::make_shared<shamrock::solvergraph::FieldRefs<Tvec>>("", "");
+
+    shamrock::patch::PatchDataLayerLayout &ghost_layout
+        = shambase::get_check_ref(storage.ghost_layout.get());
+    u32 iB_on_rho_interf = ghost_layout.get_field_idx<Tvec>("B/rho");
+    u32 ihpart_interf    = ghost_layout.get_field_idx<Tscal>("hpart");
+
+    shamrock::solvergraph::DDPatchDataFieldRef<Tvec> B_on_rho_refs = {};
+    scheduler().for_each_patchdata_nonempty([&](const Patch p, PatchDataLayer &pdat) {
+        auto &field = storage.merged_patchdata_ghost.get()
+                          .get(p.id_patch)
+                          .template get_field<Tvec>(iB_on_rho_interf);
+        B_on_rho_refs.add_obj(p.id_patch, std::ref(field));
+    });
+
+    B_on_rho_edge->set_refs(B_on_rho_refs);
+
+    // Use the "hpart" field of merged_patchdata_ghost (refreshed every corrector iteration by
+    // communicate_merge_ghosts_fields(), same as update_derivs)
+    std::shared_ptr<shamrock::solvergraph::FieldRefs<Tscal>> hpart_edge
+        = std::make_shared<shamrock::solvergraph::FieldRefs<Tscal>>("", "");
+
+    shamrock::solvergraph::DDPatchDataFieldRef<Tscal> hpart_refs = {};
+    scheduler().for_each_patchdata_nonempty([&](const Patch p, PatchDataLayer &pdat) {
+        auto &field = storage.merged_patchdata_ghost.get()
+                          .get(p.id_patch)
+                          .template get_field<Tscal>(ihpart_interf);
+        hpart_refs.add_obj(p.id_patch, std::ref(field));
+    });
+
+    hpart_edge->set_refs(hpart_refs);
+
+    Tscal const mu_0 = solver_config.get_constant_mu_0();
+
+    shambase::get_check_ref(storage.MagCurrentJ);
+    // use MagCurrenJ: on active particles (no gz)
+    modules::NodeComputeJ<Tvec, Kern> computeJ{solver_config.gpart_mass, mu_0};
+    computeJ.set_edges(
+        storage.part_counts,
+        storage.neigh_cache,
+        storage.positions_with_ghosts,
+        hpart_edge,
+        storage.omega,
+        B_on_rho_edge,
+        storage.MagCurrentJ);
+
+    computeJ.evaluate();
+}
+
 ////////////////////////////////////////////////////////////////////////////////////////////////////
 // end artificial viscosity section ////////////////////////////////////////////////////////////////
 ////////////////////////////////////////////////////////////////////////////////////////////////////
@@ -2198,7 +2262,7 @@ void shammodels::sph::Solver<Tvec, Kern>::prepare_corrector() {
     shamrock::SchedulerUtility utility(scheduler());
     PatchDataLayerLayout &pdl = scheduler().pdl_old();
 
-    bool has_B_field       = solver_config.has_field_B_on_rho();
+    bool has_b_field       = solver_config.has_field_b_on_rho();
     bool has_psi_field     = solver_config.has_field_psi_on_ch();
     bool has_epsilon_field = solver_config.dust_config.has_epsilon_field();
     bool has_deltav_field  = solver_config.dust_config.has_deltav_field();
@@ -2206,14 +2270,14 @@ void shammodels::sph::Solver<Tvec, Kern>::prepare_corrector() {
 
     const u32 iduint      = pdl.get_field_idx<Tscal>("duint");
     const u32 iaxyz       = pdl.get_field_idx<Tvec>("axyz");
-    const u32 idB_on_rho  = (has_B_field) ? pdl.get_field_idx<Tvec>("dB/rho") : 0;
+    const u32 idB_on_rho  = (has_b_field) ? pdl.get_field_idx<Tvec>("dB/rho") : 0;
     const u32 idpsi_on_ch = (has_psi_field) ? pdl.get_field_idx<Tscal>("dpsi/ch") : 0;
 
     shamlog_debug_ln("sph::BasicGas", "save old fields");
     storage.old_axyz.set(utility.save_field<Tvec>(iaxyz, "axyz_old"));
     storage.old_duint.set(utility.save_field<Tscal>(iduint, "duint_old"));
 
-    if (has_B_field) {
+    if (has_b_field) {
         storage.old_dB_on_rho.set(utility.save_field<Tvec>(idB_on_rho, "dB/rho_old"));
     }
     if (has_psi_field) {
@@ -2505,12 +2569,14 @@ shammodels::sph::TimestepLog shammodels::sph::Solver<Tvec, Kern>::evolve_once()
     using namespace shamrock;
     using namespace shamrock::patch;
 
-    bool has_B_field       = solver_config.has_field_B_on_rho();
+    bool has_b_field       = solver_config.has_field_b_on_rho();
     bool has_psi_field     = solver_config.has_field_psi_on_ch();
     bool has_epsilon_field = solver_config.dust_config.has_epsilon_field();
     bool has_deltav_field  = solver_config.dust_config.has_deltav_field();
     bool has_s_j_field     = solver_config.dust_config.has_s_j_field();
 
+    bool do_nimhd = solver_config.do_nimhd();
+
     PatchDataLayerLayout &pdl = scheduler().pdl_old();
 
     const u32 ixyz        = pdl.get_field_idx<Tvec>("xyz");
@@ -2519,8 +2585,8 @@ shammodels::sph::TimestepLog shammodels::sph::Solver<Tvec, Kern>::evolve_once()
     const u32 iuint       = pdl.get_field_idx<Tscal>("uint");
     const u32 iduint      = pdl.get_field_idx<Tscal>("duint");
     const u32 ihpart      = pdl.get_field_idx<Tscal>("hpart");
-    const u32 iB_on_rho   = (has_B_field) ? pdl.get_field_idx<Tvec>("B/rho") : 0;
-    const u32 idB_on_rho  = (has_B_field) ? pdl.get_field_idx<Tvec>("dB/rho") : 0;
+    const u32 iB_on_rho   = (has_b_field) ? pdl.get_field_idx<Tvec>("B/rho") : 0;
+    const u32 idB_on_rho  = (has_b_field) ? pdl.get_field_idx<Tvec>("dB/rho") : 0;
     const u32 ipsi_on_ch  = (has_psi_field) ? pdl.get_field_idx<Tscal>("psi/ch") : 0;
     const u32 idpsi_on_ch = (has_psi_field) ? pdl.get_field_idx<Tscal>("dpsi/ch") : 0;
     const u32 iepsilon    = (has_epsilon_field) ? pdl.get_field_idx<Tscal>("epsilon") : 0;
@@ -2593,7 +2659,7 @@ shammodels::sph::TimestepLog shammodels::sph::Solver<Tvec, Kern>::evolve_once()
     u32 iuint_interf       = ghost_layout.get_field_idx<Tscal>("uint");
     u32 ivxyz_interf       = ghost_layout.get_field_idx<Tvec>("vxyz");
     u32 iomega_interf      = ghost_layout.get_field_idx<Tscal>("omega");
-    u32 iB_on_rho_interf   = (has_B_field) ? ghost_layout.get_field_idx<Tvec>("B/rho") : 0;
+    u32 iB_on_rho_interf   = (has_b_field) ? ghost_layout.get_field_idx<Tvec>("B/rho") : 0;
     u32 ipsi_on_rho_interf = (has_psi_field) ? ghost_layout.get_field_idx<Tscal>("psi/ch") : 0;
 
     using RTreeField = RadixTreeField<Tscal>;
@@ -2670,12 +2736,12 @@ shammodels::sph::TimestepLog shammodels::sph::Solver<Tvec, Kern>::evolve_once()
             }
         }
 
-        // if (solver_config.has_field_divB()) {
+        // if (solver_config.has_field_div_b()) {
         //     sph::modules::DiffOperatorsB<Tvec, Kern>(context, solver_config, storage)
         //         .update_divB();
         // }
 
-        // if (solver_config.has_field_curlB()) {
+        // if (solver_config.has_field_curl_b()) {
         //     sph::modules::DiffOperatorsB<Tvec, Kern>(context, solver_config, storage)
         //         .update_curlB();
         // }
@@ -2726,6 +2792,67 @@ shammodels::sph::TimestepLog shammodels::sph::Solver<Tvec, Kern>::evolve_once()
             storage.alpha_av_ghost.set(std::move(merged_field));
         }
 
+        if (do_nimhd) {
+
+            // communicate needed fields (B,b, hb): done just before
+
+            // @@@ is this the correct hpart ? the one updated bu sph_prestep ?
+            shambase::get_check_ref(storage.hpart_with_ghosts)
+                .set_refs(storage.merged_xyzh.get()
+                              .template map<std::reference_wrapper<PatchDataField<Tscal>>>(
+                                  [&](u64 id, shamrock::patch::PatchDataLayer &mpdat) {
+                                      return std::ref(mpdat.get_field<Tscal>(
+                                          1)); // hpart is at index 1 in merged_xyzh
+                                  }));
+
+            // compute J field
+            update_j();
+
+            // communicate J field
+            shamrock::solvergraph::Field<Tvec> &comp_field_send
+                = shambase::get_check_ref(storage.MagCurrentJ);
+
+            using InterfaceBuildInfos =
+                typename sph::BasicSPHGhostHandler<Tvec>::InterfaceBuildInfos;
+
+            shambase::Timer time_interf;
+            time_interf.start();
+
+            auto field_interf = ghost_handle.template build_interface_native<PatchDataField<Tvec>>(
+                storage.ghost_patch_cache.get(),
+                [&](u64 sender,
+                    u64 /*receiver*/,
+                    InterfaceBuildInfos binfo,
+                    sham::DeviceBuffer<u32> &buf_idx,
+                    u32 cnt) -> PatchDataField<Tvec> {
+                    PatchDataField<Tvec> &sender_field = comp_field_send.get_field(sender);
+
+                    return sender_field.make_new_from_subset(buf_idx, cnt);
+                });
+
+            shambase::DistributedDataShared<PatchDataField<Tvec>> interf_pdat
+                = ghost_handle.communicate_pdatfield(
+                    std::move(field_interf), 1, storage.exchange_gz_J);
+
+            shambase::DistributedData<PatchDataField<Tvec>> merged_field
+                = ghost_handle.template merge_native<PatchDataField<Tvec>, PatchDataField<Tvec>>(
+                    std::move(interf_pdat),
+                    [&](const shamrock::patch::Patch p, shamrock::patch::PatchDataLayer &pdat) {
+                        PatchDataField<Tvec> &receiver_field
+                            = comp_field_send.get_field(p.id_patch);
+                        return receiver_field.duplicate();
+                    },
+                    [](PatchDataField<Tvec> &mpdat, PatchDataField<Tvec> &pdat_interf) {
+                        mpdat.insert(pdat_interf);
+                    });
+
+            time_interf.stop();
+            storage.timings_details.interface += time_interf.elapsed_sec();
+
+            // we get J with ghosts !
+            storage.MagCurrentJ_ghost.set(std::move(merged_field));
+        }
+
         // compute pressure
         compute_eos_fields();
 
@@ -2969,13 +3096,13 @@ shammodels::sph::TimestepLog shammodels::sph::Solver<Tvec, Kern>::evolve_once()
         utility.fields_leapfrog_corrector<Tscal>(
             iuint, iduint, storage.old_duint.get(), uepsilon_u_sq, dt / 2);
 
-        if (solver_config.has_field_B_on_rho()) {
+        if (solver_config.has_field_b_on_rho()) {
             ComputeField<Tscal> BOR_epsilon_BOR_sq
                 = utility.make_compute_field<Tscal>("B/rho epsilon_B/rho^2", 1);
             utility.fields_leapfrog_corrector<Tvec>(
                 iB_on_rho, idB_on_rho, storage.old_dB_on_rho.get(), BOR_epsilon_BOR_sq, dt / 2);
         }
-        if (solver_config.has_field_B_on_rho()) {
+        if (solver_config.has_field_b_on_rho()) {
             ComputeField<Tscal> POC_epsilon_POC_sq
                 = utility.make_compute_field<Tscal>("psi/ch epsilon_psi/ch^2", 1);
             utility.fields_leapfrog_corrector<Tscal>(
@@ -3039,10 +3166,10 @@ shammodels::sph::TimestepLog shammodels::sph::Solver<Tvec, Kern>::evolve_once()
 
         storage.old_axyz.reset();
         storage.old_duint.reset();
-        if (solver_config.has_field_B_on_rho()) {
+        if (solver_config.has_field_b_on_rho()) {
             storage.old_dB_on_rho.reset();
         }
-        if (solver_config.has_field_B_on_rho()) {
+        if (solver_config.has_field_b_on_rho()) {
             storage.old_dpsi_on_ch.reset();
         }
 
@@ -3131,6 +3258,36 @@ shammodels::sph::TimestepLog shammodels::sph::Solver<Tvec, Kern>::evolve_once()
                 });
             }
 
+            if (do_nimhd) {
+
+                const u32 iJ = pdl.get_field_idx<Tvec>("J");
+                shamrock::solvergraph::Field<Tvec> &MagCurrentJ
+                    = shambase::get_check_ref(storage.MagCurrentJ);
+
+                scheduler().for_each_patchdata_nonempty([&](Patch cur_p, PatchDataLayer &pdat) {
+                    sham::DeviceBuffer<Tvec> &buf_J = pdat.get_field<Tvec>(iJ).get_buf();
+
+                    sham::DeviceBuffer<Tvec> &buf_MagCurrentJ
+                        = MagCurrentJ.get_field(cur_p.id_patch).get_buf();
+
+                    auto &q = shamsys::instance::get_compute_scheduler().get_queue();
+                    sham::EventList depends_list;
+
+                    auto J           = buf_J.get_write_access(depends_list);
+                    auto MagCurrentJ = buf_MagCurrentJ.get_read_access(depends_list);
+
+                    auto e = q.submit(depends_list, [&](sycl::handler &cgh) {
+                        shambase::parallel_for(
+                            cgh, pdat.get_obj_cnt(), "write back J", [=](i32 id_a) {
+                                J[id_a] = MagCurrentJ[id_a];
+                            });
+                    });
+
+                    buf_J.complete_event_state(e);
+                    buf_MagCurrentJ.complete_event_state(e);
+                });
+            }
+
             shamlog_debug_ln("BasicGas", "computing next CFL");
 
             // Update element counts
@@ -3397,6 +3554,40 @@ shammodels::sph::TimestepLog shammodels::sph::Solver<Tvec, Kern>::evolve_once()
             compute_cfl_force->set_edges(
                 storage.part_counts, C_force_edge, hpart_refs, axyz_refs, cfl_dt);
 
+            std::shared_ptr<ComputeCFLNIMHD<Tvec>> compute_cfl_NIMHD;
+            if (do_nimhd) {
+                compute_cfl_NIMHD = std::make_shared<ComputeCFLNIMHD<Tvec>>();
+
+                Tscal C_NIMHD   = solver_config.cfl_config.cfl_NIMHD * get_cfl_multipler();
+                Cfg_MHD cfg_mhd = solver_config.mhd_config;
+                auto *nimhd     = std::get_if<typename Cfg_MHD::NonIdealMHD>(&cfg_mhd.configMHD);
+
+                Tscal eta_AD = nimhd->etaAD;
+                Tscal eta_O  = nimhd->etaO;
+                Tscal eta_H  = nimhd->etaH;
+                std::shared_ptr<shamrock::solvergraph::IDataEdge<Tscal>> C_NIMHD_edge
+                    = shamrock::solvergraph::IDataEdge<Tscal>::make_shared("C_NIMHD", "C_{NIMHD}");
+                C_NIMHD_edge->data = C_NIMHD;
+                std::shared_ptr<shamrock::solvergraph::IDataEdge<Tscal>> eta_O_edge
+                    = shamrock::solvergraph::IDataEdge<Tscal>::make_shared("eta_O", "eta_{O}");
+                eta_O_edge->data = eta_O;
+                std::shared_ptr<shamrock::solvergraph::IDataEdge<Tscal>> eta_AD_edge
+                    = shamrock::solvergraph::IDataEdge<Tscal>::make_shared("eta_AD", "eta_{AD}");
+                eta_AD_edge->data = eta_AD;
+                std::shared_ptr<shamrock::solvergraph::IDataEdge<Tscal>> eta_H_edge
+                    = shamrock::solvergraph::IDataEdge<Tscal>::make_shared("eta_H", "eta_{H}");
+                eta_H_edge->data = eta_H;
+
+                compute_cfl_NIMHD->set_edges(
+                    storage.part_counts,
+                    C_NIMHD_edge,
+                    eta_O_edge,
+                    eta_AD_edge,
+                    eta_H_edge,
+                    hpart_refs,
+                    cfl_dt);
+            }
+
             std::shared_ptr<ComputeCFLDivBCleaning<Tscal>> compute_cfl_divB_cleaning;
             if (has_psi_field) {
                 compute_cfl_divB_cleaning = std::make_shared<ComputeCFLDivBCleaning<Tscal>>();
@@ -3511,6 +3702,11 @@ shammodels::sph::TimestepLog shammodels::sph::Solver<Tvec, Kern>::evolve_once()
             compute_cfl_force->evaluate();
             save_cfl_detail("force");
 
+            if (do_nimhd) {
+                compute_cfl_NIMHD->evaluate();
+                save_cfl_detail("NIMHD");
+            }
+
             if (has_psi_field) {
                 compute_cfl_divB_cleaning->evaluate();
                 save_cfl_detail("divB_cleaning");
@@ -3636,6 +3832,11 @@ shammodels::sph::TimestepLog shammodels::sph::Solver<Tvec, Kern>::evolve_once()
         if (solver_config.has_field_alphaAV()) {
             storage.alpha_av_ghost.reset();
         }
+
+        if (do_nimhd) {
+            storage.MagCurrentJ_ghost.reset();
+        }
+
     } while (need_rerun_corrector);
 
     reset_merge_ghosts_fields();
diff --git a/src/shammodels/sph/src/SolverConfig.cpp b/src/shammodels/sph/src/SolverConfig.cpp
index db62f1639a..cc4d13df33 100644
--- a/src/shammodels/sph/src/SolverConfig.cpp
+++ b/src/shammodels/sph/src/SolverConfig.cpp
@@ -61,7 +61,7 @@ namespace shammodels::sph {
             pdl.add_field<Tscal>("soundspeed", 1);
         }
 
-        if (has_field_B_on_rho()) {
+        if (has_field_b_on_rho()) {
 
             pdl.add_field<Tvec>("B/rho", 1);
             pdl.add_field<Tvec>("dB/rho", 1);
@@ -72,14 +72,18 @@ namespace shammodels::sph {
             pdl.add_field<Tscal>("psi/ch", 1);
             pdl.add_field<Tscal>("dpsi/ch", 1);
         }
-        if (has_field_divB()) {
+        if (has_field_div_b()) {
             pdl.add_field<Tscal>("divB", 1);
         }
 
-        if (has_field_curlB()) {
+        if (has_field_curl_b()) {
             pdl.add_field<Tvec>("curlB", 1);
         }
 
+        if (do_nimhd()) {
+            pdl.add_field<Tvec>("J", 1);
+        }
+
         if (dust_config.has_epsilon_field()) {
             u32 ndust = dust_config.get_dust_nvar();
             pdl.add_field<Tscal>("epsilon", ndust);
@@ -139,7 +143,7 @@ namespace shammodels::sph {
             ghost_layout.add_field<Tscal>("soundspeed", 1);
         }
 
-        if (has_field_B_on_rho()) {
+        if (has_field_b_on_rho()) {
             ghost_layout.add_field<Tvec>("B/rho", 1);
         }
 
@@ -147,7 +151,7 @@ namespace shammodels::sph {
             ghost_layout.add_field<Tscal>("psi/ch", 1);
         }
 
-        if (has_field_curlB()) {
+        if (has_field_curl_b()) {
             ghost_layout.add_field<Tvec>("curlB", 1);
         }
 
diff --git a/src/shammodels/sph/src/modules/ComputeJ.cpp b/src/shammodels/sph/src/modules/ComputeJ.cpp
new file mode 100644
index 0000000000..799f5c3164
--- /dev/null
+++ b/src/shammodels/sph/src/modules/ComputeJ.cpp
@@ -0,0 +1,116 @@
+// -------------------------------------------------------//
+//
+// SHAMROCK code for hydrodynamics
+// Copyright (c) 2021-2026 Timothée David--Cléris <tim.shamrock@proton.me>
+// SPDX-License-Identifier: CeCILL Free Software License Agreement v2.1
+// Shamrock is licensed under the CeCILL 2.1 License, see LICENSE for more information
+//
+// -------------------------------------------------------//
+
+/**
+ * @file ComputeJ.cpp
+ * @author Yona Lapeyre (yona.lapeyre@ens-lyon.fr)
+ * @brief
+ *
+ */
+
+#include "shambase/stacktrace.hpp"
+#include "shambackends/kernel_call_distrib.hpp"
+#include "shammodels/sph/SPHUtilities.hpp"
+#include "shammodels/sph/math/mhd.hpp"
+#include "shammodels/sph/modules/ComputeJ.hpp"
+#include "shamrock/scheduler/SchedulerUtility.hpp"
+#include "shamrock/solvergraph/IFieldSpan.hpp"
+#include <fmt/ranges.h>
+
+template<class Tvec, template<class> class SPHKernel>
+void shammodels::sph::modules::NodeComputeJ<Tvec, SPHKernel>::_impl_evaluate_internal() {
+
+    __shamrock_stack_entry();
+    auto edges = get_edges();
+
+    auto dev_sched = shamsys::instance::get_compute_scheduler_ptr();
+    edges.J.ensure_sizes(edges.part_counts.indexes);
+
+    sham::distributed_data_kernel_call(
+        dev_sched,
+        sham::DDMultiRef{
+            edges.xyz.get_spans(),
+            edges.hpart.get_spans(),
+            edges.neigh_cache.neigh_cache,
+            edges.omega.get_spans(),
+            edges.B_on_rho.get_spans()},
+        sham::DDMultiRef{edges.J.get_spans()},
+        edges.part_counts.indexes,
+        [part_mass = this->part_mass, mu_0 = this->mu_0](
+            u32 id_a,
+            const Tvec *r,
+            const Tscal *hpart,
+            const auto ploop_ptrs,
+            const Tscal *omega,
+            const Tvec *B_on_rho,
+            Tvec *J) {
+            shamrock::tree::ObjectCacheIterator particle_looper(ploop_ptrs);
+
+            using namespace shamrock::sph;
+            using namespace shamrock::sph::mhd;
+
+            Tvec xyz_a = r[id_a]; // could be recovered from lambda
+
+            Tscal h_a = hpart[id_a];
+
+            Tscal rho_a    = rho_h(part_mass, h_a, SPHKernel<Tscal>::hfactd);
+            Tscal rho_a_sq = rho_a * rho_a;
+
+            Tvec B_a         = B_on_rho[id_a] * rho_a;
+            Tscal omega_a    = omega[id_a];
+            Tscal sub_fact_a = rho_a * omega_a;
+
+            Tscal part_omega_sum = 0;
+            Tvec J_sum{0, 0, 0};
+
+            constexpr Tscal rker2 = SPHKernel<Tscal>::Rkern * SPHKernel<Tscal>::Rkern;
+
+            particle_looper.for_each_object(id_a, [&](u32 id_b) {
+                Tvec dr    = xyz_a - r[id_b];
+                Tscal rab2 = sycl::dot(dr, dr);
+                Tscal h_b  = hpart[id_b];
+
+                if (rab2 > h_a * h_a * rker2 && rab2 > h_b * h_b * rker2) {
+                    return;
+                }
+
+                Tscal rab   = sycl::sqrt(rab2);
+                Tscal rho_b = rho_h(part_mass, h_b, SPHKernel<Tscal>::hfactd);
+                // if (id_a == 0) {
+                //     logger::raw_ln("@@@@@@ h_b", h_b);
+                //     logger::raw_ln("idb", id_b);
+                //     logger::raw_ln("@@@@@@ xyz_b", r[id_b]);
+                // }
+                Tvec B_b = B_on_rho[id_b] * rho_b;
+
+                Tscal Fab_a       = SPHKernel<Tscal>::dW_3d(rab, h_a);
+                Tvec r_ab_unit    = dr * sham::inv_sat_positive(rab);
+                Tvec nabla_Wab_ha = r_ab_unit * Fab_a;
+
+                J_sum += shamrock::sph::mhd::mag_current_j_sum(
+                    part_mass, B_a, B_b, nabla_Wab_ha, sub_fact_a, mu_0);
+            });
+
+            J[id_a] = J_sum; //* 4 * _pi / c;
+        });
+}
+
+template<class Tvec, template<class> class SPHKernel>
+std::string shammodels::sph::modules::NodeComputeJ<Tvec, SPHKernel>::_impl_get_tex() const {
+    return "TODO";
+}
+
+using namespace shammath;
+template class shammodels::sph::modules::NodeComputeJ<f64_3, M4>;
+template class shammodels::sph::modules::NodeComputeJ<f64_3, M6>;
+template class shammodels::sph::modules::NodeComputeJ<f64_3, M8>;
+
+template class shammodels::sph::modules::NodeComputeJ<f64_3, C2>;
+template class shammodels::sph::modules::NodeComputeJ<f64_3, C4>;
+template class shammodels::sph::modules::NodeComputeJ<f64_3, C6>;
diff --git a/src/shammodels/sph/src/modules/ConservativeCheck.cpp b/src/shammodels/sph/src/modules/ConservativeCheck.cpp
index 70c8a44915..2cd582d15a 100644
--- a/src/shammodels/sph/src/modules/ConservativeCheck.cpp
+++ b/src/shammodels/sph/src/modules/ConservativeCheck.cpp
@@ -45,10 +45,10 @@ void shammodels::sph::modules::ConservativeCheck<Tvec, SPHKernel>::check_conserv
     const u32 iduint    = pdl.get_field_idx<Tscal>("duint");
     const u32 ihpart    = pdl.get_field_idx<Tscal>("hpart");
 
-    bool has_B_field     = solver_config.has_field_B_on_rho();
-    const u32 iB_on_rho  = (has_B_field) ? pdl.get_field_idx<Tvec>("B/rho") : -1;
-    const u32 idB_on_rho = (has_B_field) ? pdl.get_field_idx<Tvec>("dB/rho") : -1;
-    const u32 idrho_dt   = (has_B_field) ? pdl.get_field_idx<Tscal>("drho/dt") : -1;
+    bool has_b_field     = solver_config.has_field_b_on_rho();
+    const u32 iB_on_rho  = (has_b_field) ? pdl.get_field_idx<Tvec>("B/rho") : -1;
+    const u32 idB_on_rho = (has_b_field) ? pdl.get_field_idx<Tvec>("dB/rho") : -1;
+    const u32 idrho_dt   = (has_b_field) ? pdl.get_field_idx<Tscal>("drho/dt") : -1;
 
     std::string cv_checks = "conservation infos :\n";
 
@@ -133,7 +133,7 @@ void shammodels::sph::modules::ConservativeCheck<Tvec, SPHKernel>::check_conserv
                 de[item] = pmass * (sycl::dot(v[item], a[item]) + du[item]);
             });
 
-        if (has_B_field) {
+        if (has_b_field) {
             PatchDataField<Tvec> &field_B_on_rho  = pdat.get_field<Tvec>(iB_on_rho);
             PatchDataField<Tvec> &field_dB_on_rho = pdat.get_field<Tvec>(idB_on_rho);
             PatchDataField<Tscal> &field_drho_dt  = pdat.get_field<Tscal>(idrho_dt);
diff --git a/src/shammodels/sph/src/modules/UpdateDerivs.cpp b/src/shammodels/sph/src/modules/UpdateDerivs.cpp
index e80f7a2919..585289af05 100644
--- a/src/shammodels/sph/src/modules/UpdateDerivs.cpp
+++ b/src/shammodels/sph/src/modules/UpdateDerivs.cpp
@@ -41,6 +41,7 @@
 #include "shammodels/sph/modules/UpdateDerivs.hpp"
 #include "shamphys/mhd.hpp"
 #include "shamrock/patch/PatchDataFieldSpan.hpp"
+#include "shamrock/scheduler/SchedulerUtility.hpp"
 #include "shamrock/solvergraph/FieldRefs.hpp"
 #include "shamrock/solvergraph/IFieldSpan.hpp"
 #include "shamrock/solvergraph/Indexes.hpp"
@@ -63,11 +64,11 @@ void shammodels::sph::modules::UpdateDerivs<Tvec, SPHKernel>::update_derivs(Tsca
         update_derivs_cd10(*v);
     } else if (ConstantDisc *v = std::get_if<ConstantDisc>(&cfg_av.config)) {
         update_derivs_disc_visco(*v);
-    } else if (IdealMHD *v = std::get_if<IdealMHD>(&cfg_mhd.config)) {
-        update_derivs_MHD(*v);
-    } else if (NonIdealMHD *v = std::get_if<NonIdealMHD>(&cfg_mhd.config)) {
-        shambase::throw_unimplemented();
-    } else if (NoneMHD *v = std::get_if<NoneMHD>(&cfg_mhd.config)) {
+    } else if (IdealMHD *v = std::get_if<IdealMHD>(&cfg_mhd.configMHD)) {
+        update_derivs_mhd(*v);
+    } else if (NonIdealMHD *v = std::get_if<NonIdealMHD>(&cfg_mhd.configMHD)) {
+        update_derivs_mhd(*v);
+    } else if (NoneMHD *v = std::get_if<NoneMHD>(&cfg_mhd.configMHD)) {
         shambase::throw_unimplemented();
     } else if (None *v = std::get_if<None>(&cfg_av.config)) {
         shambase::throw_unimplemented();
@@ -779,7 +780,43 @@ void shammodels::sph::modules::UpdateDerivs<Tvec, SPHKernel>::update_derivs_disc
 }
 
 template<class Tvec, template<class> class SPHKernel>
-void shammodels::sph::modules::UpdateDerivs<Tvec, SPHKernel>::update_derivs_MHD(IdealMHD cfg) {
+void shammodels::sph::modules::UpdateDerivs<Tvec, SPHKernel>::update_derivs_mhd(IdealMHD cfg) {
+    update_derivs_mhd_impl<shamrock::sph::mhd::MHDType::Ideal>(
+        cfg.sigma_mhd,
+        cfg.alpha_u,
+        cfg.alpha_B,
+        cfg.alpha_AV,
+        cfg.beta_AV,
+        /*etaO=*/Tscal(0),
+        /*etaH=*/Tscal(0),
+        /*etaAD=*/Tscal(0));
+}
+
+template<class Tvec, template<class> class SPHKernel>
+void shammodels::sph::modules::UpdateDerivs<Tvec, SPHKernel>::update_derivs_mhd(NonIdealMHD cfg) {
+    update_derivs_mhd_impl<shamrock::sph::mhd::MHDType::NonIdeal>(
+        cfg.sigma_mhd,
+        cfg.alpha_u,
+        cfg.alpha_B,
+        cfg.alpha_AV,
+        cfg.beta_AV,
+        cfg.etaO,
+        cfg.etaH,
+        cfg.etaAD);
+}
+
+template<class Tvec, template<class> class SPHKernel>
+template<shamrock::sph::mhd::MHDType mhd_mode>
+void shammodels::sph::modules::UpdateDerivs<Tvec, SPHKernel>::update_derivs_mhd_impl(
+    Tscal sigma_mhd,
+    Tscal alpha_u,
+    Tscal alpha_B,
+    Tscal alpha_AV,
+    Tscal beta_AV,
+    Tscal etaO,
+    Tscal etaH,
+    Tscal etaAD) {
+
     StackEntry stack_loc{};
 
     using namespace shamrock;
@@ -809,7 +846,6 @@ void shammodels::sph::modules::UpdateDerivs<Tvec, SPHKernel>::update_derivs_MHD(
     const u32 ipsi_cons     = (do_MHD_debug) ? pdl.get_field_idx<Tscal>("psi_cons") : -1;
     const u32 iu_mhd        = (do_MHD_debug) ? pdl.get_field_idx<Tscal>("u_mhd") : -1;
 
-    // Tscal mu_0 = 1.;
     Tscal const mu_0 = solver_config.get_constant_mu_0();
 
     shamrock::patch::PatchDataLayerLayout &ghost_layout
@@ -821,7 +857,7 @@ void shammodels::sph::modules::UpdateDerivs<Tvec, SPHKernel>::update_derivs_MHD(
     u32 iB_on_rho_interf  = ghost_layout.get_field_idx<Tvec>("B/rho");
     u32 ipsi_on_ch_interf = ghost_layout.get_field_idx<Tscal>("psi/ch");
 
-    // logger::raw_ln("charged the ghost fields.");
+    bool do_nimhd = solver_config.do_nimhd();
 
     auto &merged_xyzh                                 = storage.merged_xyzh.get();
     shamrock::solvergraph::Field<Tscal> &omega        = shambase::get_check_ref(storage.omega);
@@ -846,22 +882,18 @@ void shammodels::sph::modules::UpdateDerivs<Tvec, SPHKernel>::update_derivs_MHD(
         sham::DeviceBuffer<Tvec> &buf_dB_on_rho   = pdat.get_field_buf_ref<Tvec>(idB_on_rho);
         sham::DeviceBuffer<Tscal> &buf_dpsi_on_ch = pdat.get_field_buf_ref<Tscal>(idpsi_on_ch);
         sham::DeviceBuffer<Tscal> &buf_drho_dt    = pdat.get_field_buf_ref<Tscal>(idrho_dt);
-        // logger::raw_ln("charged dB dpsi");
 
         sham::DeviceBuffer<Tvec> &buf_B_on_rho = mpdat.get_field_buf_ref<Tvec>(iB_on_rho_interf);
         sham::DeviceBuffer<Tscal> &buf_psi_on_ch
             = mpdat.get_field_buf_ref<Tscal>(ipsi_on_ch_interf);
 
-        // logger::raw_ln("charged B psi");
-        //  ADD curlBBBBBBBBB
-
-        sycl::range range_npart{pdat.get_obj_cnt()};
+        bool do_nimhd = solver_config.do_nimhd();
+        sham::DeviceBuffer<Tvec> *buf_J
+            = (do_nimhd) ? &storage.MagCurrentJ_ghost.get().get(cur_p.id_patch).get_buf() : nullptr;
 
         tree::ObjectCache &pcache
             = shambase::get_check_ref(storage.neigh_cache).get_cache(cur_p.id_patch);
 
-        /////////////////////////////////////////////
-
         sham::DeviceQueue &q = shamsys::instance::get_compute_scheduler().get_queue();
         sham::EventList depends_list;
 
@@ -879,6 +911,7 @@ void shammodels::sph::modules::UpdateDerivs<Tvec, SPHKernel>::update_derivs_MHD(
         auto dB_on_rho  = buf_dB_on_rho.get_write_access(depends_list);
         auto dpsi_on_ch = buf_dpsi_on_ch.get_write_access(depends_list);
         auto drho_dt    = buf_drho_dt.get_write_access(depends_list);
+        auto J_field    = (do_nimhd) ? buf_J->get_read_access(depends_list) : nullptr;
 
         Tvec *mag_pressure
             = (do_MHD_debug)
@@ -896,7 +929,6 @@ void shammodels::sph::modules::UpdateDerivs<Tvec, SPHKernel>::update_derivs_MHD(
             = (do_MHD_debug)
                   ? pdat.get_field_buf_ref<Tvec>(itensile_corr).get_write_access(depends_list)
                   : nullptr;
-
         Tscal *psi_propag
             = (do_MHD_debug)
                   ? pdat.get_field_buf_ref<Tscal>(ipsi_propag).get_write_access(depends_list)
@@ -909,7 +941,6 @@ void shammodels::sph::modules::UpdateDerivs<Tvec, SPHKernel>::update_derivs_MHD(
             = (do_MHD_debug)
                   ? pdat.get_field_buf_ref<Tscal>(ipsi_cons).get_write_access(depends_list)
                   : nullptr;
-
         Tscal *u_mhd = (do_MHD_debug)
                            ? pdat.get_field_buf_ref<Tscal>(iu_mhd).get_write_access(depends_list)
                            : nullptr;
@@ -918,12 +949,24 @@ void shammodels::sph::modules::UpdateDerivs<Tvec, SPHKernel>::update_derivs_MHD(
 
         auto e = q.submit(depends_list, [&](sycl::handler &cgh) {
             const Tscal pmass     = solver_config.gpart_mass;
-            const Tscal sigma_mhd = cfg.sigma_mhd;
-            const Tscal alpha_u   = cfg.alpha_u;
+            const Tscal _sigma    = sigma_mhd;
+            const Tscal _alpha_u  = alpha_u;
+            const Tscal _alpha_B  = alpha_B;
+            const Tscal _alpha_AV = alpha_AV;
+            const Tscal _beta_AV  = beta_AV;
+            const Tscal _etaO     = etaO;
+            const Tscal _etaH     = etaH;
+            const Tscal _etaAD    = etaAD;
 
             shamlog_debug_ln("@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@", "");
-            shamlog_debug_sycl_ln("deriv kernel", "sigma_mhd  :", sigma_mhd);
-            shamlog_debug_sycl_ln("deriv kernel", "alpha_u  :", alpha_u);
+            shamlog_debug_sycl_ln("deriv kernel", "sigma_mhd  :", _sigma);
+            shamlog_debug_sycl_ln("deriv kernel", "alpha_u    :", _alpha_u);
+            shamlog_debug_sycl_ln("deriv kernel", "alpha_B    :", _alpha_B);
+            shamlog_debug_sycl_ln("deriv kernel", "alpha_AV   :", _alpha_AV);
+            shamlog_debug_sycl_ln("deriv kernel", "beta_AV    :", _beta_AV);
+            shamlog_debug_sycl_ln("deriv kernel", "etaO       :", _etaO);
+            shamlog_debug_sycl_ln("deriv kernel", "etaH       :", _etaH);
+            shamlog_debug_sycl_ln("deriv kernel", "etaAD      :", _etaAD);
             shamlog_debug_ln("@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@", "");
 
             tree::ObjectCacheIterator particle_looper(ploop_ptrs);
@@ -935,16 +978,15 @@ void shammodels::sph::modules::UpdateDerivs<Tvec, SPHKernel>::update_derivs_MHD(
 
                 using namespace shamrock::sph;
 
-                Tvec sum_axyz  = {0, 0, 0};
-                Tscal sum_du_a = 0;
+                Tscal h_a     = hpart[id_a];
+                Tvec xyz_a    = xyz[id_a];
+                Tvec vxyz_a   = vxyz[id_a];
+                Tscal P_a     = pressure[id_a];
+                Tscal cs_a    = cs[id_a];
+                Tscal omega_a = omega[id_a];
+                Tscal u_a     = u[id_a];
 
-                Tscal h_a       = hpart[id_a];
-                Tvec xyz_a      = xyz[id_a];
-                Tvec vxyz_a     = vxyz[id_a];
-                Tscal P_a       = pressure[id_a];
-                Tscal cs_a      = cs[id_a];
-                Tscal omega_a   = omega[id_a];
-                const Tscal u_a = u[id_a];
+                Tvec J_a = (do_nimhd) ? J_field[id_a] : Tvec{0, 0, 0};
 
                 Tscal rho_a     = rho_h(pmass, h_a, Kernel::hfactd);
                 Tscal rho_a_sq  = rho_a * rho_a;
@@ -955,7 +997,7 @@ void shammodels::sph::modules::UpdateDerivs<Tvec, SPHKernel>::update_derivs_MHD(
                 Tscal v_shock_a  = sycl::sqrt(cs_a * cs_a + v_alfven_a * v_alfven_a);
                 Tscal psi_a      = psi_on_ch[id_a] * v_shock_a;
 
-                Tscal omega_a_rho_a_inv = 1 / (omega_a * rho_a);
+                Tscal omega_a_rho_a_inv = 1. / (omega_a * rho_a);
 
                 Tvec force_pressure{0, 0, 0};
                 Tscal tmpdU_pressure = 0;
@@ -967,15 +1009,12 @@ void shammodels::sph::modules::UpdateDerivs<Tvec, SPHKernel>::update_derivs_MHD(
                 Tvec mag_tension_term{0, 0, 0};
                 Tvec gas_pressure_term{0, 0, 0};
                 Tvec tensile_corr_term{0, 0, 0};
-
                 Tscal psi_propag_term = 0;
                 Tscal psi_diff_term   = 0;
                 Tscal psi_cons_term   = 0;
-
-                Tscal u_mhd_term = 0;
+                Tscal u_mhd_term      = 0;
 
                 particle_looper.for_each_object(id_a, [&](u32 id_b) {
-                    // compute only omega_a
                     Tvec dr    = xyz_a - xyz[id_b];
                     Tscal rab2 = sycl::dot(dr, dr);
                     Tscal h_b  = hpart[id_b];
@@ -984,26 +1023,25 @@ void shammodels::sph::modules::UpdateDerivs<Tvec, SPHKernel>::update_derivs_MHD(
                         return;
                     }
 
-                    Tvec vxyz_b     = vxyz[id_b];
-                    const Tscal u_b = u[id_b];
-                    Tscal P_b       = pressure[id_b];
-                    Tscal omega_b   = omega[id_b];
-                    Tscal cs_b      = cs[id_b];
+                    Tvec vxyz_b   = vxyz[id_b];
+                    Tscal u_b     = u[id_b];
+                    Tscal P_b     = pressure[id_b];
+                    Tscal omega_b = omega[id_b];
+                    Tscal cs_b    = cs[id_b];
+                    Tscal rab     = sycl::sqrt(rab2);
 
-                    Tscal rab = sycl::sqrt(rab2);
+                    Tvec J_b = (do_nimhd) ? J_field[id_b] : Tvec{0, 0, 0};
 
                     Tscal rho_b      = rho_h(pmass, h_b, Kernel::hfactd);
                     Tvec B_b         = B_on_rho[id_b] * rho_b;
                     Tscal v_alfven_b = sycl::sqrt(sycl::dot(B_b, B_b) / (mu_0 * rho_b));
                     Tscal v_shock_b  = sycl::sqrt(cs_b * cs_b + v_alfven_b * v_alfven_b);
                     Tscal psi_b      = psi_on_ch[id_b] * v_shock_b;
-                    // const Tscal alpha_a = alpha_AV;
-                    // const Tscal alpha_b = alpha_AV;
+
                     Tscal Fab_a = Kernel::dW_3d(rab, h_a);
                     Tscal Fab_b = Kernel::dW_3d(rab, h_b);
 
-                    // Tscal sigma_mhd = 0.3;
-                    shamrock::sph::mhd::add_to_derivs_spmhd<Kernel, Tvec, Tscal>(
+                    shamrock::sph::mhd::add_to_derivs_spmhd<Kernel, Tvec, Tscal, mhd_mode>(
                         pmass,
                         dr,
                         rab,
@@ -1026,18 +1064,21 @@ void shammodels::sph::modules::UpdateDerivs<Tvec, SPHKernel>::update_derivs_MHD(
                         cs_b,
                         h_a,
                         h_b,
-
-                        alpha_u,
-
+                        _alpha_u,
+                        _alpha_B,
+                        _alpha_AV,
+                        _beta_AV,
                         B_a,
                         B_b,
-
+                        J_a,
+                        J_b,
                         psi_a,
                         psi_b,
-
                         mu_0,
-                        sigma_mhd,
-
+                        _sigma,
+                        _etaO,
+                        _etaH,
+                        _etaAD,
                         force_pressure,
                         tmpdU_pressure,
                         magnetic_eq,
@@ -1047,7 +1088,6 @@ void shammodels::sph::modules::UpdateDerivs<Tvec, SPHKernel>::update_derivs_MHD(
                         mag_tension_term,
                         gas_pressure_term,
                         tensile_corr_term,
-
                         psi_propag_term,
                         psi_diff_term,
                         psi_cons_term,
@@ -1060,17 +1100,22 @@ void shammodels::sph::modules::UpdateDerivs<Tvec, SPHKernel>::update_derivs_MHD(
                 dpsi_on_ch[id_a] = psi_eq - psi_a / h_a;
                 drho_dt[id_a]    = drho_eq;
 
+                if (do_nimhd) {
+                    // only add once per particle
+                    Tscal u_NI = shamrock::sph::mhd::u_ni_heating<Tvec, Tscal, mhd_mode>(
+                        B_a, J_a, rho_a, etaO, etaH, etaAD, mu_0);
+                    du[id_a] += u_NI;
+                }
+
                 if (do_MHD_debug) {
                     mag_pressure[id_a] = mag_pressure_term;
                     mag_tension[id_a]  = mag_tension_term;
                     gas_pressure[id_a] = gas_pressure_term;
                     tensile_corr[id_a] = tensile_corr_term;
-
-                    psi_propag[id_a] = psi_propag_term;
-                    psi_diff[id_a]   = psi_diff_term;
-                    psi_cons[id_a]   = -psi_a / h_a;
-
-                    u_mhd[id_a] = u_mhd_term;
+                    psi_propag[id_a]   = psi_propag_term;
+                    psi_diff[id_a]     = psi_diff_term;
+                    psi_cons[id_a]     = -psi_a / h_a;
+                    u_mhd[id_a]        = u_mhd_term;
                 }
             });
         });
@@ -1081,7 +1126,6 @@ void shammodels::sph::modules::UpdateDerivs<Tvec, SPHKernel>::update_derivs_MHD(
         buf_vxyz.complete_event_state(e);
         buf_hpart.complete_event_state(e);
         buf_omega.complete_event_state(e);
-        buf_uint.complete_event_state(e);
         buf_pressure.complete_event_state(e);
         buf_cs.complete_event_state(e);
         buf_B_on_rho.complete_event_state(e);
@@ -1090,16 +1134,20 @@ void shammodels::sph::modules::UpdateDerivs<Tvec, SPHKernel>::update_derivs_MHD(
         buf_dpsi_on_ch.complete_event_state(e);
         buf_drho_dt.complete_event_state(e);
 
+        if (do_nimhd) {
+            buf_J->complete_event_state(e);
+        }
+
+        buf_uint.complete_event_state(e);
+
         if (do_MHD_debug) {
             pdat.get_field_buf_ref<Tvec>(imag_pressure).complete_event_state(e);
             pdat.get_field_buf_ref<Tvec>(imag_tension).complete_event_state(e);
             pdat.get_field_buf_ref<Tvec>(igas_pressure).complete_event_state(e);
             pdat.get_field_buf_ref<Tvec>(itensile_corr).complete_event_state(e);
-
             pdat.get_field_buf_ref<Tscal>(ipsi_propag).complete_event_state(e);
             pdat.get_field_buf_ref<Tscal>(ipsi_diff).complete_event_state(e);
             pdat.get_field_buf_ref<Tscal>(ipsi_cons).complete_event_state(e);
-
             pdat.get_field_buf_ref<Tscal>(iu_mhd).complete_event_state(e);
         }
 
@@ -1107,6 +1155,8 @@ void shammodels::sph::modules::UpdateDerivs<Tvec, SPHKernel>::update_derivs_MHD(
         resulting_events.add_event(e);
         pcache.complete_event_state(resulting_events);
     });
+
+    // storage.MagCurrentJ.reset();
 }
 
 template<class Tvec, template<class> class SPHKernel>
diff --git a/src/shammodels/sph/src/pySPHModel.cpp b/src/shammodels/sph/src/pySPHModel.cpp
index cf3dcceb28..04bfebf4b4 100644
--- a/src/shammodels/sph/src/pySPHModel.cpp
+++ b/src/shammodels/sph/src/pySPHModel.cpp
@@ -219,13 +219,44 @@ void add_instance(py::module &m, std::string name_config, std::string name_model
             py::arg("alpha_u"),
             py::arg("beta_AV"))
         .def(
-            "set_IdealMHD",
-            [](TConfig &self, Tscal sigma_mhd, Tscal sigma_u) {
-                self.set_IdealMHD({sigma_mhd, sigma_u});
+            "set_ideal_mhd",
+            [](TConfig &self,
+               Tscal sigma_mhd,
+               Tscal sigma_u,
+               Tscal alpha_B,
+               Tscal alpha_AV,
+               Tscal beta_AV) {
+                self.set_ideal_mhd({sigma_mhd, sigma_u, alpha_B, alpha_AV, beta_AV});
+            },
+            py::kw_only(),
+            py::arg("sigma_mhd"),
+            py::arg("sigma_u"),
+            py::arg("alpha_B")  = 1.0,
+            py::arg("alpha_AV") = 1.0,
+            py::arg("beta_AV")  = 1.0)
+        .def(
+            "set_non_ideal_mhd",
+            [](TConfig &self,
+               Tscal sigma_mhd,
+               Tscal sigma_u,
+               Tscal etaO,
+               Tscal etaH,
+               Tscal etaAD,
+               Tscal alpha_B,
+               Tscal alpha_AV,
+               Tscal beta_AV) {
+                self.set_non_ideal_mhd(
+                    {sigma_mhd, sigma_u, alpha_B, alpha_AV, beta_AV, etaO, etaH, etaAD});
             },
             py::kw_only(),
             py::arg("sigma_mhd"),
-            py::arg("sigma_u"))
+            py::arg("sigma_u"),
+            py::arg("etaO"),
+            py::arg("etaH"),
+            py::arg("etaAD"),
+            py::arg("alpha_B")  = 1.0,
+            py::arg("alpha_AV") = 1.0,
+            py::arg("beta_AV")  = 1.0)
         .def(
             "set_self_gravity_none",
             [](TConfig &self) {