Skip to content

[SPH] enable Machida 2006 EOS #1444

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Draft
tdavidcl wants to merge 2 commits into
base: main
Choose a base branch
from
Draft

[SPH] enable Machida 2006 EOS #1444

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
86f6972
[SPH] enable Machida 2006 EOS
tdavidcl Dec 16, 2025
70a0edc
Merge branch 'main' into patch-2025-12-16-13-39
tdavidcl Jan 5, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
22 changes: 21 additions & 1 deletion src/shammodels/common/include/shammodels/common/EOSConfig.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -66,14 +66,18 @@ namespace shammodels {
using LocallyIsothermalFA2014
= shamphys::EOS_Config_LocallyIsothermalDisc_Farris2014<Tscal>;

/// Machida 2006 equation of state configuration
using Machida06 = shamphys::EOS_Config_Machida06<Tscal>;

/// Variant type to store the EOS configuration
using Variant = std::variant<
Isothermal,
Adiabatic,
Polytropic,
LocallyIsothermal,
LocallyIsothermalLP07,
LocallyIsothermalFA2014>;
LocallyIsothermalFA2014,
Machida06>;

/// Current EOS configuration
Variant config = Adiabatic{};
Expand Down Expand Up @@ -123,6 +127,22 @@ namespace shammodels {
config = LocallyIsothermalFA2014{h_over_r};
}

/**
* @brief Set the EOS configuration to a Machida 2006 equation of state
*
* @param rho_c1 Critical density 1
* @param rho_c2 Critical density 2
* @param rho_c3 Critical density 3
* @param cs Sound speed
* @param mu Mean molecular weight
* @param mh Mass of hydrogen
* @param kb Boltzmann constant
*/
inline void set_machida06(
Tscal rho_c1, Tscal rho_c2, Tscal rho_c3, Tscal cs, Tscal mu, Tscal mh, Tscal kb) {
config = Machida06{rho_c1, rho_c2, rho_c3, cs, mu, mh, kb};
}

/**
* @brief Print current status of the EOSConfig
*/
Expand Down
24 changes: 22 additions & 2 deletions src/shammodels/common/src/EOSConfig.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -61,7 +61,7 @@ namespace shammodels {
using LocIsoT = typename EOSConfig<Tvec>::LocallyIsothermal;
using LocIsoTLP07 = typename EOSConfig<Tvec>::LocallyIsothermalLP07;
using LocIsoTFA2014 = typename EOSConfig<Tvec>::LocallyIsothermalFA2014;

using Machida06 = typename EOSConfig<Tvec>::Machida06;
if (const Isothermal *eos_config = std::get_if<Isothermal>(&p.config)) {
j = json{{"Tvec", type_id}, {"eos_type", "isothermal"}, {"cs", eos_config->cs}};
} else if (const Adiabatic *eos_config = std::get_if<Adiabatic>(&p.config)) {
Expand All @@ -86,6 +86,17 @@ namespace shammodels {
{"Tvec", type_id},
{"eos_type", "locally_isothermal_fa2014"},
{"h_over_r", eos_config->h_over_r}};
} else if (const Machida06 *eos_config = std::get_if<Machida06>(&p.config)) {
j = json{
{"Tvec", type_id},
{"eos_type", "machida06"},
{"rho_c1", eos_config->rho_c1},
{"rho_c2", eos_config->rho_c2},
{"rho_c3", eos_config->rho_c3},
{"cs", eos_config->cs},
{"mu", eos_config->mu},
{"mh", eos_config->mh},
{"kb", eos_config->kb}};
} else {
shambase::throw_unimplemented(); // should never be reached
}
Expand Down Expand Up @@ -138,7 +149,7 @@ namespace shammodels {
using LocIsoT = typename EOSConfig<Tvec>::LocallyIsothermal;
using LocIsoTLP07 = typename EOSConfig<Tvec>::LocallyIsothermalLP07;
using LocIsoTFA2014 = typename EOSConfig<Tvec>::LocallyIsothermalFA2014;

using Machida06 = typename EOSConfig<Tvec>::Machida06;
if (eos_type == "isothermal") {
p.config = Isothermal{j.at("cs").get<Tscal>()};
} else if (eos_type == "adiabatic") {
Expand All @@ -152,6 +163,15 @@ namespace shammodels {
j.at("cs0").get<Tscal>(), j.at("q").get<Tscal>(), j.at("r0").get<Tscal>()};
} else if (eos_type == "locally_isothermal_fa2014") {
p.config = LocIsoTFA2014{j.at("h_over_r").get<Tscal>()};
} else if (eos_type == "machida06") {
p.config = Machida06{
j.at("rho_c1").get<Tscal>(),
j.at("rho_c2").get<Tscal>(),
j.at("rho_c3").get<Tscal>(),
j.at("cs").get<Tscal>(),
j.at("mu").get<Tscal>(),
j.at("mh").get<Tscal>(),
j.at("kb").get<Tscal>()};
} else {
shambase::throw_unimplemented("wtf !");
}
Expand Down
16 changes: 16 additions & 0 deletions src/shammodels/sph/include/shammodels/sph/SolverConfig.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -558,6 +558,22 @@ struct shammodels::sph::SolverConfig {
eos_config.set_locally_isothermalFA2014(h_over_r);
}

/**
* @brief Set the EOS configuration to a Machida 2006 equation of state
*
* @param rho_c1 Critical density 1
* @param rho_c2 Critical density 2
* @param rho_c3 Critical density 3
* @param cs Sound speed
* @param mu Mean molecular weight
* @param mh Mass of hydrogen
* @param kb Boltzmann constant
*/
inline void set_eos_machida06(
Tscal rho_c1, Tscal rho_c2, Tscal rho_c3, Tscal cs, Tscal mu, Tscal mh, Tscal kb) {
eos_config.set_machida06(rho_c1, rho_c2, rho_c3, cs, mu, mh, kb);
}

//////////////////////////////////////////////////////////////////////////////////////////////
// EOS Config (END)
//////////////////////////////////////////////////////////////////////////////////////////////
Expand Down
53 changes: 53 additions & 0 deletions src/shammodels/sph/src/modules/ComputeEos.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -108,6 +108,7 @@ void shammodels::sph::modules::ComputeEos<Tvec, SPHKernel>::compute_eos_internal
using SolverEOS_LocallyIsothermal = typename SolverConfigEOS::LocallyIsothermal;
using SolverEOS_LocallyIsothermalLP07 = typename SolverConfigEOS::LocallyIsothermalLP07;
using SolverEOS_LocallyIsothermalFA2014 = typename SolverConfigEOS::LocallyIsothermalFA2014;
using SolverEOS_Machida06 = typename SolverConfigEOS::Machida06;

sham::DeviceQueue &q = shamsys::instance::get_compute_scheduler().get_queue();

Expand Down Expand Up @@ -389,6 +390,58 @@ void shammodels::sph::modules::ComputeEos<Tvec, SPHKernel>::compute_eos_internal
buf_xyz.complete_event_state(e);
});

} else if (
SolverEOS_Machida06 *eos_config
= std::get_if<SolverEOS_Machida06>(&solver_config.eos_config.config)) {

using EOS = shamphys::EOS_Machida06<Tscal>;

storage.merged_patchdata_ghost.get().for_each([&](u64 id, PatchDataLayer &mpdat) {
auto &mfield = storage.merged_xyzh.get().get(id);

sham::DeviceBuffer<Tvec> &buf_xyz = mfield.template get_field_buf_ref<Tvec>(0);

sham::DeviceBuffer<Tscal> &buf_P
= shambase::get_check_ref(storage.pressure).get_field(id).get_buf();
sham::DeviceBuffer<Tscal> &buf_cs
= shambase::get_check_ref(storage.soundspeed).get_field(id).get_buf();
sham::DeviceBuffer<Tscal> &buf_uint = mpdat.get_field_buf_ref<Tscal>(iuint_interf);
auto rho_getter = rho_getter_gen(mpdat);

Tscal rho_c1 = eos_config->rho_c1;
Tscal rho_c2 = eos_config->rho_c2;
Tscal rho_c3 = eos_config->rho_c3;
Tscal cs = eos_config->cs;
Tscal mu = eos_config->mu;
Tscal mh = eos_config->mh;
Tscal kb = eos_config->kb;

u32 total_elements
= shambase::get_check_ref(storage.part_counts_with_ghost).indexes.get(id);

sham::kernel_call(
q,
sham::MultiRef{rho_getter, buf_uint, buf_xyz},
sham::MultiRef{buf_P, buf_cs},
total_elements,
[rho_c1, rho_c2, rho_c3, cs0 = cs, mu, mh, kb](
u32 i,
auto rho,
const Tscal *__restrict U,
const Tvec *__restrict xyz,
Tscal *__restrict P,
Tscal *__restrict cs) {
Comment on lines +427 to +433
Copy link

@gemini-code-assist gemini-code-assist bot Dec 16, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

medium

The lambda captures mu, mh, and kb, but these variables are not used within the lambda's body. This can be misleading and adds unnecessary overhead. Additionally, the parameters U and xyz are also unused. It's good practice to indicate this by removing their names. Please apply the suggested changes to improve code clarity.

                [rho_c1, rho_c2, rho_c3, cs0 = cs](
                    u32 i,
                    auto rho,
                    const Tscal *__restrict /*U*/,
                    const Tvec *__restrict /*xyz*/,
                    Tscal *__restrict P,
                    Tscal *__restrict cs) {

using namespace shamrock::sph;

Tscal rho_a = rho(i);

Tscal P_a = EOS::pressure(cs0, rho_a, rho_c1, rho_c2, rho_c3);
Tscal cs_a = EOS::soundspeed(P_a, rho_a, rho_c1, rho_c2, rho_c3);

P[i] = P_a;
cs[i] = cs_a;
});
});
} else {
shambase::throw_unimplemented();
}
Expand Down
20 changes: 20 additions & 0 deletions src/shammodels/sph/src/pySPHModel.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -95,6 +95,26 @@ void add_instance(py::module &m, std::string name_config, std::string name_model
},
py::kw_only(),
py::arg("h_over_r"))
.def(
"set_eos_machida06",
[](TConfig &self,
Tscal rho_c1,
Tscal rho_c2,
Tscal rho_c3,
Tscal cs,
Tscal mu,
Tscal mh,
Tscal kb) {
self.set_eos_machida06(rho_c1, rho_c2, rho_c3, cs, mu, mh, kb);
},
py::kw_only(),
py::arg("rho_c1"),
py::arg("rho_c2"),
py::arg("rho_c3"),
py::arg("cs"),
py::arg("mu"),
py::arg("mh"),
py::arg("kb"))
.def("set_artif_viscosity_None", &TConfig::set_artif_viscosity_None)
.def(
"to_json",
Expand Down
37 changes: 37 additions & 0 deletions src/shamphys/include/shamphys/eos_config.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -187,4 +187,41 @@ namespace shamphys {
return (lhs.h_over_r == rhs.h_over_r);
}

/**
* @brief Configuration struct for the Machida 2006 equation of state
*
* @tparam Tscal Scalar type
*/
template<class Tscal>
struct EOS_Config_Machida06 {
Tscal rho_c1;
Tscal rho_c2;
Tscal rho_c3;
Tscal cs;
Tscal mu;
Tscal mh;
Tscal kb;
};

/**
* @brief Equal operator for the EOS_Config_Machida06 struct
*
* @tparam Tscal Scalar type
* @param lhs First EOS_Config_Machida06 struct to compare
* @param rhs Second EOS_Config_Machida06 struct to compare
*
* This function checks if two EOS_Config_Machida06 structs are equal by comparing their rho_c1,
* rho_c2, rho_c3, cs, mu, mh, and kb values.
*
* @return true if the two structs have the same rho_c1, rho_c2, rho_c3, cs, mu, mh, and kb
* values, false otherwise
*/
template<class Tscal>
inline bool operator==(
const EOS_Config_Machida06<Tscal> &lhs, const EOS_Config_Machida06<Tscal> &rhs) {
return (lhs.rho_c1 == rhs.rho_c1) && (lhs.rho_c2 == rhs.rho_c2)
&& (lhs.rho_c3 == rhs.rho_c3) && (lhs.cs == rhs.cs) && (lhs.mu == rhs.mu)
&& (lhs.mh == rhs.mh) && (lhs.kb == rhs.kb);
}

} // namespace shamphys