Skip to content
Open
Show file tree
Hide file tree
Changes from all commits
Commits
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
12 changes: 6 additions & 6 deletions .github/workflows/shamrock-acpp-phys-test.yml
Original file line number Diff line number Diff line change
Expand Up @@ -196,12 +196,12 @@ jobs:
artifact_name: gallery_sph_pairing_instab
- testfile: "examples/sph/run_cubic_reorganisation.py"
artifact_name: gallery_sph_cubic_reorg
- testfile: "examples/sph/run_dustydiffuse_tvi.py"
artifact_name: gallery_sph_dustydiffuse_tvi
- testfile: "examples/sph/run_dustywave_tvi.py"
artifact_name: gallery_sph_dustywave_tvi
- testfile: "examples/sph/run_dustysettle_tvi.py"
artifact_name: gallery_sph_dustysettle_tvi
- testfile: "examples/sph/run_dustydiffuse_tva.py"
artifact_name: gallery_sph_dustydiffuse_tva
- testfile: "examples/sph/run_dustywave_tva.py"
artifact_name: gallery_sph_dustywave_tva
- testfile: "examples/sph/run_dustysettle_tva.py"
artifact_name: gallery_sph_dustysettle_tva

# Tests CI
- testfile: "examples/tests_ci/run_compare_shamrock_ph_disc.py"
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -20,7 +20,7 @@
#
# with :math:`\rho_0 = 1`, :math:`\epsilon_0 = 0.1`, :math:`r_c = 0.25`.
#
# Then we use the dust TVI solver but force :math:`d \mathbf{v} / dt = 0` and :math:`d u / dt = 0`.
# Then we use the dust TVA solver but force :math:`d \mathbf{v} / dt = 0` and :math:`d u / dt = 0`.
# In that context the epsilon equation becomes:
#
# .. math::
Expand Down Expand Up @@ -129,7 +129,7 @@ def func_s(r):
cfg.set_artif_viscosity_VaryingCD10(
alpha_min=0.0, alpha_max=1, sigma_decay=0.1, alpha_u=1, beta_AV=2
)
cfg.set_dust_mode_monofluid_tvi(nvar=1, pure_diffusion_mode=True)
cfg.set_dust_mode_monofluid_tva(nvar=1, pure_diffusion_mode=True)
cfg.set_dust_drag_constant([ts])
cfg.set_boundary_periodic()
cfg.set_eos_isothermal(cs_g)
Expand Down Expand Up @@ -256,7 +256,7 @@ def analytic_dsdt(t):
axs[1].set_xlim(0, 0.5)
axs[1].set_ylim(-0.16, 0.4)
plt.tight_layout()
plt.savefig(f"_to_trash/dump_dustydiffuse_tvi_{t:.2f}.png")
plt.savefig(f"_to_trash/dump_dustydiffuse_tva_{t:.2f}.png")
plt.close()

####################################################
Expand All @@ -278,13 +278,13 @@ def analytic_dsdt(t):

# If the animation is not returned only a static image will be shown in the doc

glob_str = os.path.join("_to_trash", "dump_dustydiffuse_tvi_*.png")
glob_str = os.path.join("_to_trash", "dump_dustydiffuse_tva_*.png")
ani = show_image_sequence(glob_str)

from matplotlib.animation import PillowWriter

writer = PillowWriter(fps=15, metadata=dict(artist="Me"), bitrate=1800)
ani.save("_to_trash/dump_dustydiffuse_tvi.gif", writer=writer)
ani.save("_to_trash/dump_dustydiffuse_tva.gif", writer=writer)

if shamrock.sys.world_rank() == 0:
# Show the animation
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -179,7 +179,7 @@ def setup_model():
alpha_min=0.0, alpha_max=1, sigma_decay=0.1, alpha_u=1, beta_AV=2
)

cfg.set_dust_mode_monofluid_tvi(nvar=ndust)
cfg.set_dust_mode_monofluid_tva(nvar=ndust)
cfg.set_dust_drag_epstein(gamma, mrn_distribution.grain_size, mrn_distribution.rho_grains)
cfg.add_ext_force_vertical_disc_potential(central_mass=1, R0=1)
cfg.add_ext_force_velocity_dissipation(eta=5)
Expand Down Expand Up @@ -781,7 +781,7 @@ def compute_sj_new(patchdata):
from matplotlib.animation import PillowWriter

writer = PillowWriter(fps=15, metadata=dict(artist="Me"), bitrate=1800)
ani.save("_to_trash/dustysettle_vert_slice_tvi.gif", writer=writer)
ani.save("_to_trash/dustysettle_vert_slice_tva.gif", writer=writer)

if shamrock.sys.world_rank() == 0:
# Show the animation
Expand All @@ -794,7 +794,7 @@ def compute_sj_new(patchdata):
from matplotlib.animation import PillowWriter

writer = PillowWriter(fps=15, metadata=dict(artist="Me"), bitrate=1800)
ani.save("_to_trash/dustysettle_vert_slice_s_tvi.gif", writer=writer)
ani.save("_to_trash/dustysettle_vert_slice_s_tva.gif", writer=writer)

if shamrock.sys.world_rank() == 0:
# Show the animation
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -75,7 +75,7 @@ def do_setup(model, cs, delta_v_0):
cfg.set_artif_viscosity_VaryingCD10(
alpha_min=0.0, alpha_max=1, sigma_decay=0.1, alpha_u=1, beta_AV=2
)
cfg.set_dust_mode_monofluid_tvi(nvar=1)
cfg.set_dust_mode_monofluid_tva(nvar=1)
cfg.set_dust_drag_constant([ts])
cfg.set_boundary_periodic()
cfg.set_eos_isothermal(cs)
Expand Down Expand Up @@ -162,7 +162,7 @@ def internal_rho_d(size: int, rho: np.array, eps: np.array) -> np.array:
# Analytics


def dustywave_tvi_matrix(k, cs, ts, eps):
def dustywave_tva_matrix(k, cs, ts, eps):
a = k * cs
b = k * k * ts * cs * cs * eps

Expand All @@ -176,8 +176,8 @@ def dustywave_tvi_matrix(k, cs, ts, eps):
)


def eigensystem_dustywave_tvi(k, cs, ts, eps):
M = dustywave_tvi_matrix(k, cs, ts, eps)
def eigensystem_dustywave_tva(k, cs, ts, eps):
M = dustywave_tva_matrix(k, cs, ts, eps)
vals, vecs = np.linalg.eig(M)
return 1j * vals, vecs

Expand Down Expand Up @@ -295,7 +295,7 @@ def fit_sine_wave(x, y, prev_phi=None):
print(omega_k)

print(f"k={k} cs={cs} ts={ts} epsilon_0={epsilon_0}")
eigval, eigvec = eigensystem_dustywave_tvi(k, cs, ts, epsilon_0)
eigval, eigvec = eigensystem_dustywave_tva(k, cs, ts, epsilon_0)

print(f"eigenval = {eigval}")
print(f"eigenvec = {eigvec}")
Expand Down Expand Up @@ -408,7 +408,7 @@ def fit_sine_wave(x, y, prev_phi=None):
)
plt.legend(loc="upper right")
plt.tight_layout()
plt.savefig(f"_to_trash/dump_dustywave_tvi_{ics:02d}_{i:02d}.png")
plt.savefig(f"_to_trash/dump_dustywave_tva_{ics:02d}_{i:02d}.png")
plt.close()

# if i == 1:
Expand All @@ -434,7 +434,7 @@ def fit_sine_wave(x, y, prev_phi=None):
plt.ylabel("$\delta$ fields [code unit]")
plt.title(f"cs={cs:.6g}")
plt.legend(fontsize=12, loc="upper right")
plt.savefig(f"_to_trash/dustywave_tvi_scan_{ics:04}.png")
plt.savefig(f"_to_trash/dustywave_tva_scan_{ics:04}.png")

# %%
# make gifs
Expand All @@ -445,19 +445,19 @@ def fit_sine_wave(x, y, prev_phi=None):

# %%
# show them the gifs (i have to unroll the loop otherwise the doc does not capture the gifs ...)
ani0 = show_image_sequence(f"_to_trash/dump_dustywave_tvi_{0:02d}_*.png")
ani0 = show_image_sequence(f"_to_trash/dump_dustywave_tva_{0:02d}_*.png")
writer = PillowWriter(fps=15, metadata=dict(artist="Me"), bitrate=1800)
ani0.save(f"_to_trash/dustywave_tvi_scan_{0:04}.gif", writer=writer)
ani0.save(f"_to_trash/dustywave_tva_scan_{0:04}.gif", writer=writer)
plt.show()
# %%
ani1 = show_image_sequence(f"_to_trash/dump_dustywave_tvi_{1:02d}_*.png")
ani1 = show_image_sequence(f"_to_trash/dump_dustywave_tva_{1:02d}_*.png")
writer = PillowWriter(fps=15, metadata=dict(artist="Me"), bitrate=1800)
ani1.save(f"_to_trash/dustywave_tvi_scan_{1:04}.gif", writer=writer)
ani1.save(f"_to_trash/dustywave_tva_scan_{1:04}.gif", writer=writer)
plt.show()
# %%
ani2 = show_image_sequence(f"_to_trash/dump_dustywave_tvi_{2:02d}_*.png")
ani2 = show_image_sequence(f"_to_trash/dump_dustywave_tva_{2:02d}_*.png")
writer = PillowWriter(fps=15, metadata=dict(artist="Me"), bitrate=1800)
ani2.save(f"_to_trash/dustywave_tvi_scan_{2:04}.gif", writer=writer)
ani2.save(f"_to_trash/dustywave_tva_scan_{2:04}.gif", writer=writer)
plt.show()

plt.show()
22 changes: 11 additions & 11 deletions src/shammodels/sph/include/shammodels/sph/SolverConfig.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -115,7 +115,7 @@ namespace shammodels::sph {

struct None {};

struct MonofluidTVI {
struct MonofluidTVA {
u32 ndust;
bool pure_diffusion_mode = false;
};
Expand All @@ -125,28 +125,28 @@ namespace shammodels::sph {
};

/// Variant type to store the EOS configuration
using Variant = std::variant<None, MonofluidTVI, MonofluidComplete>;
using Variant = std::variant<None, MonofluidTVA, MonofluidComplete>;

Variant current_mode = None{};

inline void set_none() { current_mode = None{}; }
inline void set_monofluid_tvi(u32 nvar, bool pure_diffusion_mode = false) {
current_mode = MonofluidTVI{nvar, pure_diffusion_mode};
inline void set_monofluid_tva(u32 nvar, bool pure_diffusion_mode = false) {
current_mode = MonofluidTVA{nvar, pure_diffusion_mode};
}
inline void set_monofluid_complete(u32 nvar) { current_mode = MonofluidComplete{nvar}; }

inline bool is_none() { return std::holds_alternative<None>(current_mode); }
inline bool is_monofluid_tvi() { return bool(std::get_if<MonofluidTVI>(&current_mode)); }
inline bool is_monofluid_tva() { return bool(std::get_if<MonofluidTVA>(&current_mode)); }
inline bool is_monofluid_complete() {
return bool(std::get_if<MonofluidComplete>(&current_mode));
}

inline void mode_to_json(nlohmann::json &j) const {
if (const None *cfg = std::get_if<None>(&current_mode)) {
j = {{"type", "none"}};
} else if (const MonofluidTVI *cfg = std::get_if<MonofluidTVI>(&current_mode)) {
} else if (const MonofluidTVA *cfg = std::get_if<MonofluidTVA>(&current_mode)) {
j
= {{"type", "monofluid_tvi"},
= {{"type", "monofluid_tva"},
{"ndust", cfg->ndust},
{"pure_diffusion_mode", cfg->pure_diffusion_mode}};
} else if (
Expand All @@ -161,8 +161,8 @@ namespace shammodels::sph {
const std::string type = j.at("type").get<std::string>();
if (type == "none") {
set_none();
} else if (type == "monofluid_tvi") {
set_monofluid_tvi(
} else if (type == "monofluid_tva") {
set_monofluid_tva(
j.at("ndust").get<u32>(), j.at("pure_diffusion_mode").get<bool>());
Comment on lines +164 to 166

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

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

medium

To maintain backward compatibility with existing simulation configuration files and checkpoint dumps that contain "type": "monofluid_tvi", we should also accept "monofluid_tvi" and map it to set_monofluid_tva during JSON deserialization.

Suggested change
} else if (type == "monofluid_tva") {
set_monofluid_tva(
j.at("ndust").get<u32>(), j.at("pure_diffusion_mode").get<bool>());
} else if (type == "monofluid_tva" || type == "monofluid_tvi") {
set_monofluid_tva(
j.at("ndust").get<u32>(), j.at("pure_diffusion_mode").get<bool>());

} else if (type == "monofluid_complete") {
set_monofluid_complete(j.at("ndust").get<u32>());
Expand All @@ -172,7 +172,7 @@ namespace shammodels::sph {
}

inline bool has_s_j_field() {
return is_monofluid_tvi(); // S_j = sqrt(\rho \epsilon_j)
return is_monofluid_tva(); // S_j = sqrt(\rho \epsilon_j)
}

inline bool has_epsilon_field() {
Expand All @@ -188,7 +188,7 @@ namespace shammodels::sph {
shambase::throw_with_loc<std::invalid_argument>(
"Querying a dust nvar with no dust as config is ... discutable ...");
return 0;
} else if (MonofluidTVI *cfg = std::get_if<MonofluidTVI>(&current_mode)) {
} else if (MonofluidTVA *cfg = std::get_if<MonofluidTVA>(&current_mode)) {
return cfg->ndust;
} else if (MonofluidComplete *cfg = std::get_if<MonofluidComplete>(&current_mode)) {
return cfg->ndust;
Expand Down
4 changes: 2 additions & 2 deletions src/shammodels/sph/src/pySPHModel.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -266,9 +266,9 @@ void add_instance(py::module &m, std::string name_config, std::string name_model
self.dust_config.set_none();
})
.def(
"set_dust_mode_monofluid_tvi",
"set_dust_mode_monofluid_tva",
[](TConfig &self, u32 nvar, bool pure_diffusion_mode) {
self.dust_config.set_monofluid_tvi(nvar, pure_diffusion_mode);
self.dust_config.set_monofluid_tva(nvar, pure_diffusion_mode);
},
Comment on lines 268 to 272

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

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

medium

To prevent breaking existing Python scripts that call set_dust_mode_monofluid_tvi, we should keep a deprecated alias in the Python bindings that prints a warning and delegates to the new set_dust_mode_monofluid_tva method.

        .def(
            "set_dust_mode_monofluid_tva",
            [](TConfig &self, u32 nvar, bool pure_diffusion_mode) {
                self.dust_config.set_monofluid_tva(nvar, pure_diffusion_mode);
            },
            py::kw_only(),
            py::arg("nvar"),
            py::arg("pure_diffusion_mode") = false)
        .def(
            "set_dust_mode_monofluid_tvi",
            [](TConfig &self, u32 nvar, bool pure_diffusion_mode) {
                ON_RANK_0(shamlog_warn_ln(
                              "SPH",
                              ".set_dust_mode_monofluid_tvi() is deprecated, please use .set_dust_mode_monofluid_tva() instead."));
                self.dust_config.set_monofluid_tva(nvar, pure_diffusion_mode);
            },

py::kw_only(),
py::arg("nvar"),
Expand Down
Loading