[SPH][NIMHD] Add Non-Ideal MHD terms

src/shammodels/sph/include/shammodels/sph/SolverConfig.hpp

@@ -970,6 +970,9 @@ 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); }
 

src/shammodels/sph/include/shammodels/sph/config/MHDConfig.hpp

@@ -44,6 +44,9 @@ struct shammodels::sph::MHDConfig {
     struct NonIdealMHD {
         Tscal sigma_mhd = 0.1;
         Tscal alpha_u   = 1.;
+        Tscal etaO      = 1.;
+        Tscal etaH      = 1.;
+        Tscal etaAD     = 1.;
     };
 
     // how to set a new state of a variant as a dummy:
@@ -52,55 +55,60 @@ struct shammodels::sph::MHDConfig {
     //-> question your life choices
     using Variant = std::variant<None, IdealMHD_constrained_hyper_para, NonIdealMHD>;
 
-    Variant config = None{};
+    Variant configMHD = None{};
 
-    void set(Variant v) { config = v; }
+    void set(Variant v) { configMHD = v; }
+
+    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<IdealMHD_constrained_hyper_para>(&config))
-                    || bool(std::get_if<NonIdealMHD>(&config));
+        bool is_B = bool(std::get_if<IdealMHD_constrained_hyper_para>(&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<IdealMHD_constrained_hyper_para>(&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));
+        bool is_divB = bool(std::get_if<IdealMHD_constrained_hyper_para>(&configMHD));
         return is_divB;
     }
 
     inline bool has_curlB_field() {
-        bool is_curlB = bool(std::get_if<NonIdealMHD>(&config));
+        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));
+        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)) {
+            = std::get_if<IdealMHD_constrained_hyper_para>(&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)) {
+        } 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);
         } else {
             shambase::throw_unimplemented();
         }
 
-        logger::raw_ln("--- MHD config (deduced)");
+        logger::raw_ln("--- MHD configMHD (deduced)");
 
         logger::raw_ln("-------------");
     }
@@ -122,18 +130,18 @@ namespace shammodels::sph {
         using IMHD        = typename T::IdealMHD_constrained_hyper_para;
         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},
             };
-        } 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"},
@@ -162,15 +170,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 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") {

src/shammodels/sph/include/shammodels/sph/math/mhd.hpp

@@ -24,6 +24,7 @@
 #include "shammodels/sph/SolverConfig.hpp"
 #include "shammodels/sph/math/forces.hpp"
 #include "shammodels/sph/math/q_ab.hpp"
+#include "shammodels/sph/modules/SolverStorage.hpp"
 #include "shamphys/mhd.hpp"
 #include "shamunits/Constants.hpp"
 #include <tuple>
@@ -32,6 +33,49 @@ namespace shamrock::sph::mhd {
 
     enum MHDType { Ideal = 0, NonIdeal = 1 };
 
+    template<class Tvec, class Tscal, MHDType MHD_mode = NonIdeal>
+    inline void MagCurrentJ_sum(
+        Tscal m_b, Tvec B_a, Tvec B_b, Tvec nabla_Wab_ha, Tscal sub_fact_a, Tscal mu_0, Tvec &J_a) {
+
+        J_a += m_b * sham::inv_sat_zero(sub_fact_a) * sycl::cross(B_a - B_b, nabla_Wab_ha) / mu_0;
+    }
+
+    template<class Tvec, class Tscal, MHDType MHD_mode = NonIdeal>
+    inline Tvec WursterD(Tvec B, Tvec J, Tscal etaO, Tscal etaH, Tscal etaAD) {
+
+        Tvec Bhat = B / sham::inv_sat_zero(sycl::length(B));
+        Tvec D    = etaO * J + etaH * sycl::cross(J, Bhat)
+                    + etaAD * sycl::cross(sycl::cross(J, Bhat), Bhat);
+
+        return D;
+    }
+
+    template<class Tvec, class Tscal, MHDType MHD_mode = NonIdeal>
+    inline Tscal u_NI_heating(Tvec D, Tvec J, Tscal rho) {
+
+        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);
+    }
+
     // mag tension form the Tricco 2023 formula
     template<class Tvec, class Tscal>
     inline Tvec B_dot_grad_W(
@@ -159,7 +203,7 @@ 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;
     }
@@ -227,12 +271,19 @@ namespace shamrock::sph::mhd {
         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,
@@ -364,8 +415,9 @@ 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
+                        * dB_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(
             pmass,
@@ -401,6 +453,37 @@ 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) {
+            // logger::raw_ln("############# NON IDEAL MHD #############");
+
+            // Tvec J_a = MagCurrentJ<Tvec, Tscal, MHD_mode>(
+            //     pmass, B_a, B_b, r_ab_unit * dWab_a, sub_fact_a, mu_0);
+            // Tvec J_b = MagCurrentJ<Tvec, Tscal, MHD_mode>(
+            //     pmass, B_a, B_b, r_ab_unit * dWab_b, sub_fact_b, mu_0);
+
+            Tvec D_a = WursterD<Tvec, Tscal, MHD_mode>(B_a, J_a, etaO, etaH, etaAD);
+            Tvec D_b = WursterD<Tvec, Tscal, MHD_mode>(B_b, J_b, etaO, etaH, etaAD);
+
+            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;
+
+            Tscal u_NI = u_NI_heating<Tvec, Tscal, MHD_mode>(D_a, J_a, rho_a) * 0.5
+                         + u_NI_heating<Tvec, Tscal, MHD_mode>(D_b, J_b, rho_b) * 0.5;
+
+            du_dt += u_NI;
+        }
     }
 
 } // namespace shamrock::sph::mhd

src/shammodels/sph/include/shammodels/sph/modules/SolverStorage.hpp

@@ -105,6 +105,8 @@ namespace shammodels::sph {
         Component<shamrock::ComputeField<Tvec>> old_dB_on_rho;
         Component<shamrock::ComputeField<Tscal>> old_dpsi_on_ch;
 
+        Component<shamrock::ComputeField<Tvec>> MagCurrentJ;
+
         Component<shamrock::ComputeField<Tscal>> old_dtepsilon;
         Component<shamrock::ComputeField<Tvec>> old_dtdeltav;
         Component<shamrock::ComputeField<Tscal>> old_ds_j_dt;

src/shammodels/sph/include/shammodels/sph/modules/UpdateDerivs.hpp

@@ -19,6 +19,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"
 
@@ -70,7 +71,18 @@ namespace shammodels::sph::modules {
         using IdealMHD    = typename Cfg_MHD::IdealMHD_constrained_hyper_para;
         using NonIdealMHD = typename Cfg_MHD::NonIdealMHD;
 
+        // 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 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);
+
+        template<shamrock::sph::mhd::MHDType MHD_mode>
+        void compute_J(Tscal mu_0);
     };
 
 } // namespace shammodels::sph::modules