Merge branch 'main' into features/clang-tidy-fix-3

Author: tdavidcl

doc/sphinx/examples/physics/run_eos.py

@@ -0,0 +1,67 @@
+"""
+Barotropic EOS functions
+=======================================
+
+"""
+
+import matplotlib.pyplot as plt
+import numpy as np
+
+import shamrock
+
+cs = 190.0
+rho_c1 = 1.92e-13 * 1000  # g/cm^3 -> kg/m^3
+rho_c2 = 3.84e-8 * 1000  # g/cm^3 -> kg/m^3
+rho_c3 = 1.92e-3 * 1000  # g/cm^3 -> kg/m^3
+
+
+si = shamrock.UnitSystem()
+sicte = shamrock.Constants(si)
+kb = sicte.kb()
+print(kb)
+mu = 2.375
+mh = 1.00784 * sicte.dalton()
+print(mu * mh * kb)
+
+rho_plot = np.logspace(-15, 5, 1000)
+P_plot = []
+cs_plot = []
+T_plot = []
+for rho in rho_plot:
+    P, _cs, T = shamrock.phys.eos.eos_Machida06(
+        cs=cs, rho=rho, rho_c1=rho_c1, rho_c2=rho_c2, rho_c3=rho_c3, mu=mu, mh=mh, kb=kb
+    )
+    P_plot.append(P)
+    cs_plot.append(_cs)
+    T_plot.append(T)
+
+plt.figure()
+plt.plot(rho_plot, P_plot, label="P")
+plt.yscale("log")
+plt.xscale("log")
+plt.xlabel("$\\rho$ [kg.m^-3]")
+plt.ylabel("$P$ [Pa]")
+plt.axvspan(rho_c1, rho_c2, color="grey", alpha=0.5)
+plt.axvspan(rho_c3, rho_plot[-1], color="grey", alpha=0.5)
+
+plt.figure()
+plt.plot(rho_plot, cs_plot, label="cs")
+plt.yscale("log")
+plt.xlabel("$\\rho$ [kg.m^-3]")
+plt.xscale("log")
+plt.ylabel("$c_s$ [m/s]")
+plt.axvspan(rho_c1, rho_c2, color="grey", alpha=0.5)
+plt.axvspan(rho_c3, rho_plot[-1], color="grey", alpha=0.5)
+
+plt.figure()
+plt.plot(rho_plot, T_plot, label="T")
+plt.yscale("log")
+plt.xscale("log")
+plt.xlabel("$\\rho$ [kg.m^-3]")
+plt.ylabel("$T$ [K]")
+plt.axvspan(rho_c1, rho_c2, color="grey", alpha=0.5)
+plt.axvspan(rho_c3, rho_plot[-1], color="grey", alpha=0.5)
+
+
+plt.tight_layout()
+plt.show()

src/shamphys/include/shamphys/eos.hpp

@@ -67,4 +67,34 @@ namespace shamphys {
         static constexpr T pressure_from_cs(T cs0sq, T rho) { return cs0sq * rho; }
     };
 
+    template<class T>
+    struct EOS_Machida06 {
+
+        static constexpr T soundspeed(T P, T rho, T rho_c1, T rho_c2, T rho_c3) {
+            const T gamma = (rho < rho_c1)   ? T(1.0)
+                            : (rho < rho_c2) ? T(7.0 / 5.0)
+                            : (rho < rho_c3) ? T(1.1)
+                                             : T(5.0 / 3.0);
+            return sycl::sqrt(gamma * P / rho);
+        }
+
+        static constexpr T temperature(T P, T rho, T mu, T mh, T kb) {
+            return mu * mh * P / (rho * kb);
+        }
+
+        static constexpr T pressure(T cs, T rho, T rho_c1, T rho_c2, T rho_c3) {
+            if (rho < rho_c1) {
+                return cs * cs * rho;
+            } else if (rho < rho_c2) {
+                return cs * cs * rho_c1 * sycl::pow(rho / rho_c1, 7. / 5.);
+            } else if (rho < rho_c3) {
+                return cs * cs * rho_c1 * sycl::pow(rho_c2 / rho_c1, 7. / 5.)
+                       * sycl::pow(rho / rho_c2, 1.1);
+            } else {
+                return cs * cs * rho_c1 * sycl::pow(rho_c2 / rho_c1, 7. / 5.)
+                       * sycl::pow(rho_c3 / rho_c2, 1.1) * sycl::pow(rho / rho_c3, 5. / 3.);
+            }
+        }
+    };
+
 } // namespace shamphys

src/shampylib/src/pyShamphys.cpp

@@ -23,6 +23,7 @@
 #include "shamphys/Planets.hpp"
 #include "shamphys/SedovTaylor.hpp"
 #include "shamphys/SodTube.hpp"
+#include "shamphys/eos.hpp"
 #include "shamphys/fmm/GreenFuncGravCartesian.hpp"
 #include "shamphys/fmm/contract_grav_moment.hpp"
 #include "shamphys/fmm/grav_moment_offset.hpp"
@@ -409,4 +410,26 @@ Register_pymod(shamphyslibinit) {
         py::arg("dM"),
         py::arg("from"),
         py::arg("to"));
+
+    // EOS
+
+    py::module eos_module = shamphys_module.def_submodule("eos");
+
+    eos_module.def(
+        "eos_Machida06",
+        [](f64 cs, f64 rho, f64 rho_c1, f64 rho_c2, f64 rho_c3, f64 mu, f64 mh, f64 kb) {
+            auto P   = shamphys::EOS_Machida06<f64>::pressure(cs, rho, rho_c1, rho_c2, rho_c3);
+            auto _cs = shamphys::EOS_Machida06<f64>::soundspeed(P, rho, rho_c1, rho_c2, rho_c3);
+            auto T   = shamphys::EOS_Machida06<f64>::temperature(P, rho, mu, mh, kb);
+            return std::tuple<f64, f64, f64>{P, _cs, T};
+        },
+        py::kw_only(),
+        py::arg("cs"),
+        py::arg("rho"),
+        py::arg("rho_c1"),
+        py::arg("rho_c2"),
+        py::arg("rho_c3"),
+        py::arg("mu"),
+        py::arg("mh"),
+        py::arg("kb"));
 }

src/shamunits/include/shamunits/Constants.hpp

@@ -70,6 +70,7 @@
     X(jupiter_mass /***/, Uget(kg, 1))                                                             \
     X(sol_mass /*******/, Uget(kg, 1))                                                             \
     X(planck_mass /****/, Uget(kg, 1))                                                             \
+    X(dalton /*********/, Uget(kg, 1))                                                             \
     /* densities */                                                                                \
     X(guiness_density, Uget(kg, 1) * Uget(m, -1))                                                  \
     /* derived ctes  */                                                                            \
@@ -155,6 +156,7 @@ namespace shamunits {
             static constexpr T jupiter_mass  = 1.898e27;                               //(kg)
             static constexpr T sol_mass      = 1.98847e30;                             //(kg)
             static constexpr T planck_mass   = 2.17643424e-8;                          //(kg)
+            static constexpr T dalton        = 1.66053906892e-27;                      //(kg)
 
             static constexpr T guiness_density = Conv::gcm3_to_guiness_density * 1000; //(kg.m-3)
         };