Skip to content

Commit 9ffa101

Browse files
tdavidcltdavidcl
authored
Merge branch 'main' into features/clang-tidy-fix-1
2 parents fa48faa + f9423cc commit 9ffa101

5 files changed

Lines changed: 145 additions & 25 deletions

File tree

‎doc/sphinx/examples/physics/run_eos.py‎

Lines changed: 67 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,67 @@
1+
"""
2+
Barotropic EOS functions
3+
=======================================
4+
5+
"""
6+
7+
import matplotlib.pyplot as plt
8+
import numpy as np
9+
10+
import shamrock
11+
12+
cs = 190.0
13+
rho_c1 = 1.92e-13 * 1000 # g/cm^3 -> kg/m^3
14+
rho_c2 = 3.84e-8 * 1000 # g/cm^3 -> kg/m^3
15+
rho_c3 = 1.92e-3 * 1000 # g/cm^3 -> kg/m^3
16+
17+
18+
si = shamrock.UnitSystem()
19+
sicte = shamrock.Constants(si)
20+
kb = sicte.kb()
21+
print(kb)
22+
mu = 2.375
23+
mh = 1.00784 * sicte.dalton()
24+
print(mu * mh * kb)
25+
26+
rho_plot = np.logspace(-15, 5, 1000)
27+
P_plot = []
28+
cs_plot = []
29+
T_plot = []
30+
for rho in rho_plot:
31+
P, _cs, T = shamrock.phys.eos.eos_Machida06(
32+
cs=cs, rho=rho, rho_c1=rho_c1, rho_c2=rho_c2, rho_c3=rho_c3, mu=mu, mh=mh, kb=kb
33+
)
34+
P_plot.append(P)
35+
cs_plot.append(_cs)
36+
T_plot.append(T)
37+
38+
plt.figure()
39+
plt.plot(rho_plot, P_plot, label="P")
40+
plt.yscale("log")
41+
plt.xscale("log")
42+
plt.xlabel("$\\rho$ [kg.m^-3]")
43+
plt.ylabel("$P$ [Pa]")
44+
plt.axvspan(rho_c1, rho_c2, color="grey", alpha=0.5)
45+
plt.axvspan(rho_c3, rho_plot[-1], color="grey", alpha=0.5)
46+
47+
plt.figure()
48+
plt.plot(rho_plot, cs_plot, label="cs")
49+
plt.yscale("log")
50+
plt.xlabel("$\\rho$ [kg.m^-3]")
51+
plt.xscale("log")
52+
plt.ylabel("$c_s$ [m/s]")
53+
plt.axvspan(rho_c1, rho_c2, color="grey", alpha=0.5)
54+
plt.axvspan(rho_c3, rho_plot[-1], color="grey", alpha=0.5)
55+
56+
plt.figure()
57+
plt.plot(rho_plot, T_plot, label="T")
58+
plt.yscale("log")
59+
plt.xscale("log")
60+
plt.xlabel("$\\rho$ [kg.m^-3]")
61+
plt.ylabel("$T$ [K]")
62+
plt.axvspan(rho_c1, rho_c2, color="grey", alpha=0.5)
63+
plt.axvspan(rho_c3, rho_plot[-1], color="grey", alpha=0.5)
64+
65+
66+
plt.tight_layout()
67+
plt.show()

‎src/shambackends/src/SyclMpiTypes.cpp‎

Lines changed: 23 additions & 25 deletions
Original file line numberDiff line numberDiff line change
@@ -47,24 +47,17 @@ inline MPI_Datatype __tmp_mpi_type_f64_3;
4747

4848
#define __SYCL_TYPE_COMMIT_len2(base_name, src_type) \
4949
{ \
50-
base_name a; \
51-
MPI_Aint offset_##base_name = ((size_t) ((char *) &(a.x()) - (char *) &(a))); \
52-
MPICHECK(MPI_Type_create_struct( \
53-
1, &__len_vec2, &offset_##base_name, &mpi_type_##src_type, &mpi_type_##base_name)); \
50+
check_offset_validity<base_name>(); \
51+
MPICHECK(MPI_Type_contiguous(__len_vec2, mpi_type_##src_type, &mpi_type_##base_name)); \
5452
MPICHECK(MPI_Type_commit(&mpi_type_##base_name)); \
5553
shamlog_debug_mpi_ln("SyclMpiTypes", "init mpi type for : " #base_name); \
5654
}
5755

5856
#define __SYCL_TYPE_COMMIT_len3(base_name, src_type) \
5957
{ \
60-
base_name a; \
61-
MPI_Aint offset_##base_name = ((size_t) ((char *) &(a.x()) - (char *) &(a))); \
62-
MPICHECK(MPI_Type_create_struct( \
63-
1, \
64-
&__len_vec3, \
65-
&offset_##base_name, \
66-
&mpi_type_##src_type, \
67-
&__tmp_mpi_type_##base_name)); \
58+
check_offset_validity<base_name>(); \
59+
MPICHECK( \
60+
MPI_Type_contiguous(__len_vec3, mpi_type_##src_type, &__tmp_mpi_type_##base_name)); \
6861
MPICHECK(MPI_Type_create_resized( \
6962
__tmp_mpi_type_##base_name, 0, sizeof(base_name), &mpi_type_##base_name)); \
7063
MPICHECK(MPI_Type_commit(&mpi_type_##base_name)); \
@@ -73,35 +66,40 @@ inline MPI_Datatype __tmp_mpi_type_f64_3;
7366

7467
#define __SYCL_TYPE_COMMIT_len4(base_name, src_type) \
7568
{ \
76-
base_name a; \
77-
MPI_Aint offset_##base_name = ((size_t) ((char *) &(a.x()) - (char *) &(a))); \
78-
MPICHECK(MPI_Type_create_struct( \
79-
1, &__len_vec4, &offset_##base_name, &mpi_type_##src_type, &mpi_type_##base_name)); \
69+
check_offset_validity<base_name>(); \
70+
MPICHECK(MPI_Type_contiguous(__len_vec4, mpi_type_##src_type, &mpi_type_##base_name)); \
8071
MPICHECK(MPI_Type_commit(&mpi_type_##base_name)); \
8172
shamlog_debug_mpi_ln("SyclMpiTypes", "init mpi type for : " #base_name); \
8273
}
8374

8475
#define __SYCL_TYPE_COMMIT_len8(base_name, src_type) \
8576
{ \
86-
base_name a; \
87-
MPI_Aint offset_##base_name = ((size_t) ((char *) &(a.s0()) - (char *) &(a))); \
88-
MPICHECK(MPI_Type_create_struct( \
89-
1, &__len_vec8, &offset_##base_name, &mpi_type_##src_type, &mpi_type_##base_name)); \
77+
check_offset_validity<base_name>(); \
78+
MPICHECK(MPI_Type_contiguous(__len_vec8, mpi_type_##src_type, &mpi_type_##base_name)); \
9079
MPICHECK(MPI_Type_commit(&mpi_type_##base_name)); \
9180
shamlog_debug_mpi_ln("SyclMpiTypes", "init mpi type for : " #base_name); \
9281
}
9382

9483
#define __SYCL_TYPE_COMMIT_len16(base_name, src_type) \
9584
{ \
96-
base_name a; \
97-
MPI_Aint offset_##base_name = ((size_t) ((char *) &(a.s0()) - (char *) &(a))); \
98-
MPICHECK(MPI_Type_create_struct( \
99-
1, &__len_vec16, &offset_##base_name, &mpi_type_##src_type, &mpi_type_##base_name)); \
85+
check_offset_validity<base_name>(); \
86+
MPICHECK(MPI_Type_contiguous(__len_vec16, mpi_type_##src_type, &mpi_type_##base_name)); \
10087
MPICHECK(MPI_Type_commit(&mpi_type_##base_name)); \
10188
shamlog_debug_mpi_ln("SyclMpiTypes", "init mpi type for : " #base_name); \
10289
}
10390

104-
// TODO check mpi errors
91+
template<class T>
92+
void check_offset_validity() {
93+
T a{};
94+
95+
std::ptrdiff_t base = reinterpret_cast<std::ptrdiff_t>(&a);
96+
std::ptrdiff_t s0 = reinterpret_cast<std::ptrdiff_t>(&a.s0());
97+
98+
if (s0 - base != 0) {
99+
throw shambase::make_except_with_loc<std::runtime_error>(shambase::format(
100+
"Offset is not valid for type {}, base = {}, s0 = {}", typeid(T).name(), base, s0));
101+
}
102+
}
105103

106104
void create_sycl_mpi_types() {
107105

‎src/shamphys/include/shamphys/eos.hpp‎

Lines changed: 30 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -67,4 +67,34 @@ namespace shamphys {
6767
static constexpr T pressure_from_cs(T cs0sq, T rho) { return cs0sq * rho; }
6868
};
6969

70+
template<class T>
71+
struct EOS_Machida06 {
72+
73+
static constexpr T soundspeed(T P, T rho, T rho_c1, T rho_c2, T rho_c3) {
74+
const T gamma = (rho < rho_c1) ? T(1.0)
75+
: (rho < rho_c2) ? T(7.0 / 5.0)
76+
: (rho < rho_c3) ? T(1.1)
77+
: T(5.0 / 3.0);
78+
return sycl::sqrt(gamma * P / rho);
79+
}
80+
81+
static constexpr T temperature(T P, T rho, T mu, T mh, T kb) {
82+
return mu * mh * P / (rho * kb);
83+
}
84+
85+
static constexpr T pressure(T cs, T rho, T rho_c1, T rho_c2, T rho_c3) {
86+
if (rho < rho_c1) {
87+
return cs * cs * rho;
88+
} else if (rho < rho_c2) {
89+
return cs * cs * rho_c1 * sycl::pow(rho / rho_c1, 7. / 5.);
90+
} else if (rho < rho_c3) {
91+
return cs * cs * rho_c1 * sycl::pow(rho_c2 / rho_c1, 7. / 5.)
92+
* sycl::pow(rho / rho_c2, 1.1);
93+
} else {
94+
return cs * cs * rho_c1 * sycl::pow(rho_c2 / rho_c1, 7. / 5.)
95+
* sycl::pow(rho_c3 / rho_c2, 1.1) * sycl::pow(rho / rho_c3, 5. / 3.);
96+
}
97+
}
98+
};
99+
70100
} // namespace shamphys

‎src/shampylib/src/pyShamphys.cpp‎

Lines changed: 23 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -23,6 +23,7 @@
2323
#include "shamphys/Planets.hpp"
2424
#include "shamphys/SedovTaylor.hpp"
2525
#include "shamphys/SodTube.hpp"
26+
#include "shamphys/eos.hpp"
2627
#include "shamphys/fmm/GreenFuncGravCartesian.hpp"
2728
#include "shamphys/fmm/contract_grav_moment.hpp"
2829
#include "shamphys/fmm/grav_moment_offset.hpp"
@@ -409,4 +410,26 @@ Register_pymod(shamphyslibinit) {
409410
py::arg("dM"),
410411
py::arg("from"),
411412
py::arg("to"));
413+
414+
// EOS
415+
416+
py::module eos_module = shamphys_module.def_submodule("eos");
417+
418+
eos_module.def(
419+
"eos_Machida06",
420+
[](f64 cs, f64 rho, f64 rho_c1, f64 rho_c2, f64 rho_c3, f64 mu, f64 mh, f64 kb) {
421+
auto P = shamphys::EOS_Machida06<f64>::pressure(cs, rho, rho_c1, rho_c2, rho_c3);
422+
auto _cs = shamphys::EOS_Machida06<f64>::soundspeed(P, rho, rho_c1, rho_c2, rho_c3);
423+
auto T = shamphys::EOS_Machida06<f64>::temperature(P, rho, mu, mh, kb);
424+
return std::tuple<f64, f64, f64>{P, _cs, T};
425+
},
426+
py::kw_only(),
427+
py::arg("cs"),
428+
py::arg("rho"),
429+
py::arg("rho_c1"),
430+
py::arg("rho_c2"),
431+
py::arg("rho_c3"),
432+
py::arg("mu"),
433+
py::arg("mh"),
434+
py::arg("kb"));
412435
}

‎src/shamunits/include/shamunits/Constants.hpp‎

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -70,6 +70,7 @@
7070
X(jupiter_mass /***/, Uget(kg, 1)) \
7171
X(sol_mass /*******/, Uget(kg, 1)) \
7272
X(planck_mass /****/, Uget(kg, 1)) \
73+
X(dalton /*********/, Uget(kg, 1)) \
7374
/* densities */ \
7475
X(guiness_density, Uget(kg, 1) * Uget(m, -1)) \
7576
/* derived ctes */ \
@@ -155,6 +156,7 @@ namespace shamunits {
155156
static constexpr T jupiter_mass = 1.898e27; //(kg)
156157
static constexpr T sol_mass = 1.98847e30; //(kg)
157158
static constexpr T planck_mass = 2.17643424e-8; //(kg)
159+
static constexpr T dalton = 1.66053906892e-27; //(kg)
158160

159161
static constexpr T guiness_density = Conv::gcm3_to_guiness_density * 1000; //(kg.m-3)
160162
};

0 commit comments

Comments
 (0)