[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,51 @@ 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,
+        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 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 +205,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 +273,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 +417,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 +455,39 @@ 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,
+                B_a,
+                B_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

src/shammodels/sph/src/Solver.cpp

@@ -1321,6 +1321,7 @@ void shammodels::sph::Solver<Tvec, Kern>::communicate_merge_ghosts_fields() {
 
     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 iJ          = (has_B_field) ? pdl.get_field_idx<Tvec>("J") : 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;
@@ -1345,6 +1346,7 @@ void shammodels::sph::Solver<Tvec, Kern>::communicate_merge_ghosts_fields() {
     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 iJ_interf     = ghost_layout.get_field_idx<Tvec>("J");
 
     const u32 iaxyz_interf
         = (solver_config.has_axyz_in_ghost()) ? ghost_layout.get_field_idx<Tvec>("axyz") : 0;
@@ -1411,6 +1413,8 @@ void shammodels::sph::Solver<Tvec, Kern>::communicate_merge_ghosts_fields() {
             if (has_B_field) {
                 sender_patch.get_field<Tvec>(iB_on_rho).append_subset_to(
                     buf_idx, cnt, pdat.get_field<Tvec>(iB_interf));
+                sender_patch.get_field<Tvec>(iJ).append_subset_to(
+                    buf_idx, cnt, pdat.get_field<Tvec>(iJ_interf));
             }
 
             if (has_psi_field) {
@@ -1493,6 +1497,7 @@ void shammodels::sph::Solver<Tvec, Kern>::communicate_merge_ghosts_fields() {
 
                 if (has_B_field) {
                     pdat_new.get_field<Tvec>(iB_interf).insert(pdat.get_field<Tvec>(iB_on_rho));
+                    pdat_new.get_field<Tvec>(iJ_interf).insert(pdat.get_field<Tvec>(iJ));
                 }
 
                 if (has_psi_field) {
@@ -1677,6 +1682,8 @@ shammodels::sph::TimestepLog shammodels::sph::Solver<Tvec, Kern>::evolve_once()
     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");
@@ -2574,6 +2581,25 @@ shammodels::sph::TimestepLog shammodels::sph::Solver<Tvec, Kern>::evolve_once()
                     vclean_buf.complete_event_state(e);
                 };
 
+                if (do_NIMHD) {
+
+                    auto e = q.submit(depends_list, [&](sycl::handler &cgh) {
+                        auto *nimhd = std::get_if<typename MHDConfig<Tvec>::NonIdealMHD>(
+                            &solver_config.mhd_config.configMHD);
+                        Tscal etaAD = nimhd ? nimhd->etaAD : Tscal(0);
+
+                        cgh.parallel_for(
+                            sycl::range<1>{pdat.get_obj_cnt()}, [=](sycl::item<1> item) {
+                                Tscal h_a = hpart[item];
+
+                                Tscal dt_AD = (etaAD > 0) ? (h_a * h_a / etaAD)
+                                                          : shambase::get_infty<Tscal>();
+
+                                cfl_dt[item] = sycl::min(cfl_dt[item], dt_AD);
+                            });
+                    });
+                };
+
                 buf_hpart.complete_event_state(e);
                 buf_axyz.complete_event_state(e);
                 vsig_buf.complete_event_state(e);

src/shammodels/sph/src/SolverConfig.cpp

@@ -66,6 +66,7 @@ namespace shammodels::sph {
             pdl.add_field<Tvec>("B/rho", 1);
             pdl.add_field<Tvec>("dB/rho", 1);
             pdl.add_field<Tscal>("drho/dt", 1);
+            pdl.add_field<Tvec>("J", 1);
         }
 
         if (has_field_psi_on_ch()) {
@@ -138,6 +139,7 @@ namespace shammodels::sph {
 
         if (has_field_B_on_rho()) {
             ghost_layout.add_field<Tvec>("B/rho", 1);
+            ghost_layout.add_field<Tvec>("J", 1); // @@@ only for ni mhd
         }
 
         if (has_field_psi_on_ch()) {