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

src/shammodels/sph/CMakeLists.txt

@@ -34,6 +34,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

src/shammodels/sph/include/shammodels/sph/Solver.hpp

@@ -12,7 +12,7 @@
 /**
  * @file Solver.hpp
  * @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
  */
 
@@ -186,6 +186,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();
 
@@ -204,6 +207,7 @@ namespace shammodels::sph {
         void prepare_corrector();
         /// @brief Updates time derivatives and applies external forces
         void update_derivs(Tscal dt_hydro);
+
         /**
          * @brief
          *

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

@@ -538,7 +538,21 @@ struct shammodels::sph::SolverConfig {
         mhd_config.set(v);
     }
 
-    inline void set_NonIdealMHD(typename MHDConfig::NonIdealMHD v) { mhd_config.set(v); }
+    inline void set_NonIdealMHD(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.");
+        mhd_config.set(v);
+    }
 
     //////////////////////////////////////////////////////////////////////////////////////////////
     // MHD Config (END)
@@ -973,6 +987,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") {