From 23af945815449accd3d09ec5ebe5e8f4f9fae3b7 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Timoth=C3=A9e=20David--Cl=C3=A9ris?= Date: Tue, 7 Jul 2026 00:04:47 +0200 Subject: [PATCH 01/24] yolo --- examples/sph/run_dustysettle_tvi.py | 87 +++++++++++++++++++++------ src/shampylib/src/math/pyshammath.cpp | 19 ++++++ 2 files changed, 87 insertions(+), 19 deletions(-) diff --git a/examples/sph/run_dustysettle_tvi.py b/examples/sph/run_dustysettle_tvi.py index 30b55ac704..e161079cfb 100644 --- a/examples/sph/run_dustysettle_tvi.py +++ b/examples/sph/run_dustysettle_tvi.py @@ -136,22 +136,22 @@ def uint_g(r): bmin = (-box / 8, -box / 8, -box) bmax = (box / 8, box / 8, box) -N_target = 1e4 scheduler_split_val = int(2e7) scheduler_merge_val = int(1) -xm, ym, zm = bmin -xM, yM, zM = bmax -vol_b = (xM - xm) * (yM - ym) * (zM - zm) - -part_vol = vol_b / N_target - -# lattice volume -HCP_PACKING_DENSITY = 0.74 -part_vol_lattice = HCP_PACKING_DENSITY * part_vol - -dr = (part_vol_lattice / ((4.0 / 3.0) * np.pi)) ** (1.0 / 3.0) -bmin, bmax = shamrock.math.get_ideal_hcp_box(dr, bmin, bmax) +lx = 12 +ly = 12 +lz = 512 +lmin = (-lx//2,-ly//2,-lz//2) +lmax = (lx//2,ly//2,lz//2) + +# Call with dr = 1 as we will rescale on next call +(xm, ym, zm), (xM, yM, zM) = shamrock.math.get_periodic_hcp_box(1., lmin, lmax) +print(f"base lattice : xM = {xM}, yM = {yM}, zM = {zM}") +dr = 2*box / (zM - zm) +print(f"dr = {dr}") +bmin, bmax = shamrock.math.get_periodic_hcp_box(dr, lmin, lmax) +print(f"new lattice : bmin = {bmin}, bmax = {bmax}") xm, ym, zm = bmin xM, yM, zM = bmax @@ -508,7 +508,7 @@ def setup(self): class ReferenceDustySettleAll: def __init__(self): self.rhomid = get_max_rho() - self.tau = 0.025 + self.tau = 0.1 self.vK = kep_profile(R0) * codeu.to("m") / codeu.to("s") self.OmegaK = omega_k(R0) * codeu.to("s") ** -1 @@ -544,6 +544,32 @@ def evolve_until(self, t): reference_dusty_settle = None +def compute_L2_error(z,field, z_ref, field_ref): + + z_field_sort_args = np.argsort(z) + z_field_sorted = z[z_field_sort_args] + field_field_sorted = field[z_field_sort_args] + + # Interpolate field to the same grid as the reference field + field_interp = np.interp(z_ref, z_field_sorted, field_field_sorted) + + # Compute delta + delta = field_interp - field_ref + + # Compute L2 func + L2_func = delta ** 2 + + # Compute L2 integral # Use the appropriate trapezoidal integration function + if hasattr(np, "trapezoid"): + L2_integral = np.trapezoid(L2_func, z_ref) # NumPy >= 2.0 + else: + L2_integral = np.trapz(L2_func, z_ref) # NumPy < 2.0 + + + return np.sqrt(L2_integral) + +t_l2 = [] +l2_errors = [[] for _ in range(ndust)] def analyse_and_plot(j): @@ -589,6 +615,9 @@ def analyse_and_plot(j): rho_dust_all = np.zeros(len(z)) epsilon_dust_all = np.zeros(len(z)) + if reference_dusty_settle is not None: + t_l2.append(time) + for i in range(ndust): c = dust_colors[i] axs[0].scatter(z, s_j[:, i] ** 2 * to_dens, s=sz, color=c, edgecolors="none") @@ -613,6 +642,10 @@ def analyse_and_plot(j): print(ana_epsilon.max(), ana_epsilon.min()) axs[1].plot(reference_dusty_settle.soluces[i].zbar, ana_epsilon, "--", color="0.0") + L2_error = compute_L2_error(z, s_j[:, i] ** 2 / rho, reference_dusty_settle.soluces[i].zbar, ana_epsilon) + print(f"L2 error for dust {i} = {L2_error}") + l2_errors[i].append(L2_error) + axs[0].scatter(z, rho * to_dens, s=sz, color="0.0", edgecolors="none") axs[0].scatter(z, rho_dust_all, s=sz, color="0.5", edgecolors="none") axs[1].scatter(z, 1 - epsilon_dust_all, s=sz, color="0.0", edgecolors="none") @@ -623,7 +656,7 @@ def analyse_and_plot(j): axs[0].set_xlabel(r"$z$") axs[0].set_yscale("log") axs[0].set_ylim(1e-20, 1e-8) - axs[0].set_xlim(-2.2 * H, 2.2 * H) + axs[0].set_xlim(-box, box) # axs[0].set_ylim(1e-12, 10**2) axs[1].set_ylabel(r"$\epsilon_j$") @@ -631,7 +664,7 @@ def analyse_and_plot(j): axs[1].set_yscale("log") # axs[1].set_ylim(1e-12, 2) # if you want the full range axs[1].set_ylim(1e-4, 1e-1) # if you want to see the dust only - axs[1].set_xlim(-2.2 * H, 2.2 * H) + axs[1].set_xlim(-box, box) gas_handle = Line2D( [0], @@ -785,7 +818,7 @@ def compute_sj_new(patchdata): if shamrock.sys.world_rank() == 0: # Show the animation - plt.show() + plt.close() # %% glob_str = f"{dump_folder}/plots/vert_slice_s_*.png" @@ -798,7 +831,7 @@ def compute_sj_new(patchdata): if shamrock.sys.world_rank() == 0: # Show the animation - plt.show() + plt.close() # %% # Plot the mass history @@ -838,9 +871,25 @@ def compute_sj_new(patchdata): plt.xlabel("t") plt.ylabel("$\|delta M_{dust} / M_{dust,0}$") -plt.yscale("log") +plt.yscale("symlog", linthresh=1e-7) plt.title("Dust mass conservation") plt.legend() plt.tight_layout() plt.savefig(f"{dump_folder}/plots/dust_mass_history.png") plt.show() + +# %% +# Plot the L2 errors + +plt.figure() +for k in range(ndust): + plt.plot(t_l2, l2_errors[k], label=f"dust {k}") + print(f"L2 error for dust {k} = {l2_errors[k]}") +plt.xlabel("t") +plt.ylabel("L2 error") +plt.yscale("log") +plt.title("L2 errors") +plt.legend() +plt.tight_layout() +plt.savefig(f"{dump_folder}/plots/l2_errors.png") +plt.show() \ No newline at end of file diff --git a/src/shampylib/src/math/pyshammath.cpp b/src/shampylib/src/math/pyshammath.cpp index a804de744f..01f272ece8 100644 --- a/src/shampylib/src/math/pyshammath.cpp +++ b/src/shampylib/src/math/pyshammath.cpp @@ -859,4 +859,23 @@ ON_PYTHON_INIT { math_module.def("get_ideal_hcp_box", [](f64 dr, f64_3 box_min, f64_3 box_max) { return shammath::LatticeHCP::get_ideal_hcp_box(dr, {box_min, box_max}); }); + + math_module.def( + "get_periodic_hcp_box", + [](f64 dr, std::array box_min, std::array box_max) { + auto ret = shammath::LatticeHCP::get_periodic_box(dr, box_min, box_max); + return std::tuple{ret.lower, ret.upper}; + }, + py::arg("dr"), + py::arg("box_min"), + py::arg("box_max"), + R"pbdoc( + Get the periodic box corresponding to integer lattice coordinates + this function will throw if the coordinates asked cannot make a periodic lattice + + Args: + dr: the particle spacing in the lattice + box_min: integer triplet for the minimal coordinates on the lattice + box_max: integer triplet for the maximal coordinates on the lattice + )pbdoc"); } From d75539d421726bc6881befa71f922e286703baf5 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?David--Cl=C3=A9ris=20Timoth=C3=A9e?= Date: Tue, 7 Jul 2026 00:27:29 +0200 Subject: [PATCH 02/24] Update examples/sph/run_dustysettle_tvi.py --- examples/sph/run_dustysettle_tvi.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/sph/run_dustysettle_tvi.py b/examples/sph/run_dustysettle_tvi.py index e161079cfb..dcc9ef74cd 100644 --- a/examples/sph/run_dustysettle_tvi.py +++ b/examples/sph/run_dustysettle_tvi.py @@ -508,7 +508,7 @@ def setup(self): class ReferenceDustySettleAll: def __init__(self): self.rhomid = get_max_rho() - self.tau = 0.1 + self.tau = 0.05 self.vK = kep_profile(R0) * codeu.to("m") / codeu.to("s") self.OmegaK = omega_k(R0) * codeu.to("s") ** -1 From c3a674ea20675c1885830b644e7f63b53e70024f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?David--Cl=C3=A9ris=20Timoth=C3=A9e?= Date: Tue, 7 Jul 2026 00:27:38 +0200 Subject: [PATCH 03/24] Update examples/sph/run_dustysettle_tvi.py --- examples/sph/run_dustysettle_tvi.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/sph/run_dustysettle_tvi.py b/examples/sph/run_dustysettle_tvi.py index dcc9ef74cd..41045337c1 100644 --- a/examples/sph/run_dustysettle_tvi.py +++ b/examples/sph/run_dustysettle_tvi.py @@ -141,7 +141,7 @@ def uint_g(r): lx = 12 ly = 12 -lz = 512 +lz = 64 lmin = (-lx//2,-ly//2,-lz//2) lmax = (lx//2,ly//2,lz//2) From 3c73dcf4d5443f1088ff41a8661e58b0b8a6bc63 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?David--Cl=C3=A9ris=20Timoth=C3=A9e?= Date: Tue, 7 Jul 2026 00:27:54 +0200 Subject: [PATCH 04/24] Update examples/sph/run_dustysettle_tvi.py --- examples/sph/run_dustysettle_tvi.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/sph/run_dustysettle_tvi.py b/examples/sph/run_dustysettle_tvi.py index 41045337c1..8e214ad5bd 100644 --- a/examples/sph/run_dustysettle_tvi.py +++ b/examples/sph/run_dustysettle_tvi.py @@ -656,7 +656,7 @@ def analyse_and_plot(j): axs[0].set_xlabel(r"$z$") axs[0].set_yscale("log") axs[0].set_ylim(1e-20, 1e-8) - axs[0].set_xlim(-box, box) + axs[0].set_xlim(-2.2 * H, 2.2 * H) # axs[0].set_ylim(1e-12, 10**2) axs[1].set_ylabel(r"$\epsilon_j$") From 078924839ffc01ad6fe5e77318c68c02dd936261 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?David--Cl=C3=A9ris=20Timoth=C3=A9e?= Date: Tue, 7 Jul 2026 00:28:06 +0200 Subject: [PATCH 05/24] Update examples/sph/run_dustysettle_tvi.py --- examples/sph/run_dustysettle_tvi.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/sph/run_dustysettle_tvi.py b/examples/sph/run_dustysettle_tvi.py index 8e214ad5bd..b0d64f720e 100644 --- a/examples/sph/run_dustysettle_tvi.py +++ b/examples/sph/run_dustysettle_tvi.py @@ -871,7 +871,7 @@ def compute_sj_new(patchdata): plt.xlabel("t") plt.ylabel("$\|delta M_{dust} / M_{dust,0}$") -plt.yscale("symlog", linthresh=1e-7) +plt.yscale("symlog", linthresh=1e-8) plt.title("Dust mass conservation") plt.legend() plt.tight_layout() From f7315914bb0c8b20fbc2a75bb33ca488349e7ab9 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?David--Cl=C3=A9ris=20Timoth=C3=A9e?= Date: Tue, 7 Jul 2026 00:28:13 +0200 Subject: [PATCH 06/24] Update examples/sph/run_dustysettle_tvi.py --- examples/sph/run_dustysettle_tvi.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/sph/run_dustysettle_tvi.py b/examples/sph/run_dustysettle_tvi.py index b0d64f720e..9fa14f10ef 100644 --- a/examples/sph/run_dustysettle_tvi.py +++ b/examples/sph/run_dustysettle_tvi.py @@ -831,7 +831,7 @@ def compute_sj_new(patchdata): if shamrock.sys.world_rank() == 0: # Show the animation - plt.close() + plt.show() # %% # Plot the mass history From 22b65facb7296bf4fa8adc33d8b435af3e2d5551 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?David--Cl=C3=A9ris=20Timoth=C3=A9e?= Date: Tue, 7 Jul 2026 00:28:24 +0200 Subject: [PATCH 07/24] Update examples/sph/run_dustysettle_tvi.py --- examples/sph/run_dustysettle_tvi.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/sph/run_dustysettle_tvi.py b/examples/sph/run_dustysettle_tvi.py index 9fa14f10ef..0fc1946cf0 100644 --- a/examples/sph/run_dustysettle_tvi.py +++ b/examples/sph/run_dustysettle_tvi.py @@ -818,7 +818,7 @@ def compute_sj_new(patchdata): if shamrock.sys.world_rank() == 0: # Show the animation - plt.close() + plt.show() # %% glob_str = f"{dump_folder}/plots/vert_slice_s_*.png" From 96f996f1d0db2d019f112659e9bc35cbc540c565 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?David--Cl=C3=A9ris=20Timoth=C3=A9e?= Date: Tue, 7 Jul 2026 00:28:33 +0200 Subject: [PATCH 08/24] Update examples/sph/run_dustysettle_tvi.py --- examples/sph/run_dustysettle_tvi.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/sph/run_dustysettle_tvi.py b/examples/sph/run_dustysettle_tvi.py index 0fc1946cf0..47f27c209b 100644 --- a/examples/sph/run_dustysettle_tvi.py +++ b/examples/sph/run_dustysettle_tvi.py @@ -664,7 +664,7 @@ def analyse_and_plot(j): axs[1].set_yscale("log") # axs[1].set_ylim(1e-12, 2) # if you want the full range axs[1].set_ylim(1e-4, 1e-1) # if you want to see the dust only - axs[1].set_xlim(-box, box) + axs[1].set_xlim(-2.2 * H, 2.2 * H) gas_handle = Line2D( [0], From b3fafaddc9deb8981a22f15cb3a39547d2e7f00d Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Mon, 6 Jul 2026 22:29:02 +0000 Subject: [PATCH 09/24] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- examples/sph/run_dustysettle_tvi.py | 24 ++++++++++++++---------- 1 file changed, 14 insertions(+), 10 deletions(-) diff --git a/examples/sph/run_dustysettle_tvi.py b/examples/sph/run_dustysettle_tvi.py index 47f27c209b..47c75d5ac9 100644 --- a/examples/sph/run_dustysettle_tvi.py +++ b/examples/sph/run_dustysettle_tvi.py @@ -142,13 +142,13 @@ def uint_g(r): lx = 12 ly = 12 lz = 64 -lmin = (-lx//2,-ly//2,-lz//2) -lmax = (lx//2,ly//2,lz//2) +lmin = (-lx // 2, -ly // 2, -lz // 2) +lmax = (lx // 2, ly // 2, lz // 2) # Call with dr = 1 as we will rescale on next call -(xm, ym, zm), (xM, yM, zM) = shamrock.math.get_periodic_hcp_box(1., lmin, lmax) +(xm, ym, zm), (xM, yM, zM) = shamrock.math.get_periodic_hcp_box(1.0, lmin, lmax) print(f"base lattice : xM = {xM}, yM = {yM}, zM = {zM}") -dr = 2*box / (zM - zm) +dr = 2 * box / (zM - zm) print(f"dr = {dr}") bmin, bmax = shamrock.math.get_periodic_hcp_box(dr, lmin, lmax) print(f"new lattice : bmin = {bmin}, bmax = {bmax}") @@ -544,7 +544,8 @@ def evolve_until(self, t): reference_dusty_settle = None -def compute_L2_error(z,field, z_ref, field_ref): + +def compute_L2_error(z, field, z_ref, field_ref): z_field_sort_args = np.argsort(z) z_field_sorted = z[z_field_sort_args] @@ -557,20 +558,21 @@ def compute_L2_error(z,field, z_ref, field_ref): delta = field_interp - field_ref # Compute L2 func - L2_func = delta ** 2 + L2_func = delta**2 # Compute L2 integral # Use the appropriate trapezoidal integration function if hasattr(np, "trapezoid"): L2_integral = np.trapezoid(L2_func, z_ref) # NumPy >= 2.0 else: - L2_integral = np.trapz(L2_func, z_ref) # NumPy < 2.0 - + L2_integral = np.trapz(L2_func, z_ref) # NumPy < 2.0 return np.sqrt(L2_integral) + t_l2 = [] l2_errors = [[] for _ in range(ndust)] + def analyse_and_plot(j): pmass = model.get_particle_mass() @@ -642,7 +644,9 @@ def analyse_and_plot(j): print(ana_epsilon.max(), ana_epsilon.min()) axs[1].plot(reference_dusty_settle.soluces[i].zbar, ana_epsilon, "--", color="0.0") - L2_error = compute_L2_error(z, s_j[:, i] ** 2 / rho, reference_dusty_settle.soluces[i].zbar, ana_epsilon) + L2_error = compute_L2_error( + z, s_j[:, i] ** 2 / rho, reference_dusty_settle.soluces[i].zbar, ana_epsilon + ) print(f"L2 error for dust {i} = {L2_error}") l2_errors[i].append(L2_error) @@ -892,4 +896,4 @@ def compute_sj_new(patchdata): plt.legend() plt.tight_layout() plt.savefig(f"{dump_folder}/plots/l2_errors.png") -plt.show() \ No newline at end of file +plt.show() From cc140263f4e726e7dac468d7708092c0b818ab06 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Timoth=C3=A9e=20David--Cl=C3=A9ris?= Date: Tue, 7 Jul 2026 22:42:52 +0200 Subject: [PATCH 10/24] reproduce --- examples/sph/run_dustysettle_tvi.py | 21 ++++++++------- .../modules/NodeUpdateDerivsMonofluidTVI.cpp | 27 +++++++++++++++++++ 2 files changed, 38 insertions(+), 10 deletions(-) diff --git a/examples/sph/run_dustysettle_tvi.py b/examples/sph/run_dustysettle_tvi.py index 47c75d5ac9..434ad8d535 100644 --- a/examples/sph/run_dustysettle_tvi.py +++ b/examples/sph/run_dustysettle_tvi.py @@ -59,9 +59,9 @@ epsilon_base = 0.01 -tlist = [0.1 * i for i in range(20)] + [i * 0.1 + 2 for i in range(3000)] -tinject = 2 -iinject = 20 +tlist = [0.01 * i for i in range(1000)] +tinject = 0 +iinject = 0 t_end = 5.0 @@ -141,7 +141,7 @@ def uint_g(r): lx = 12 ly = 12 -lz = 64 +lz = 512 lmin = (-lx // 2, -ly // 2, -lz // 2) lmax = (lx // 2, ly // 2, lz // 2) @@ -215,7 +215,7 @@ def f_remap(r): rn = max(abs(zM), abs(zm)) # print(y, H, H * erfinv(y / rn)) - z = H * erfinv(z / rn) # * (2**0.5) + z = H * erfinv(z / rn) * (2**0.5) z = min(z, zM) z = max(z, zm) @@ -660,7 +660,7 @@ def analyse_and_plot(j): axs[0].set_xlabel(r"$z$") axs[0].set_yscale("log") axs[0].set_ylim(1e-20, 1e-8) - axs[0].set_xlim(-2.2 * H, 2.2 * H) + axs[0].set_xlim(-box, box) # axs[0].set_ylim(1e-12, 10**2) axs[1].set_ylabel(r"$\epsilon_j$") @@ -668,7 +668,7 @@ def analyse_and_plot(j): axs[1].set_yscale("log") # axs[1].set_ylim(1e-12, 2) # if you want the full range axs[1].set_ylim(1e-4, 1e-1) # if you want to see the dust only - axs[1].set_xlim(-2.2 * H, 2.2 * H) + axs[1].set_xlim(-box, box) gas_handle = Line2D( [0], @@ -730,12 +730,12 @@ def analyse_and_plot(j): axs[0].set_ylabel(r"$s_j$") axs[0].set_xlabel(r"$z$") - axs[0].set_xlim(-4 * H, 4 * H) + axs[0].set_xlim(-box, box) axs[0].set_yscale("symlog", linthresh=1e-10) axs[1].set_ylabel(r"$\dot{s}_j$") axs[1].set_xlabel(r"$z$") - axs[1].set_xlim(-4 * H, 4 * H) + axs[1].set_xlim(-box, box) axs[1].set_yscale("symlog", linthresh=1e-10) dust_sm = cm.ScalarMappable(cmap=dust_cmap, norm=dust_norm) @@ -770,6 +770,7 @@ def dust_mass_analysis(): for j in range(1000): if j == iinject: reference_dusty_settle = ReferenceDustySettleAll() + reference_dusty_settle = None if tlist[j] >= t_start: if j > 0: @@ -800,7 +801,7 @@ def compute_sj_new(patchdata): analyse_and_plot(j) - dump_helper.write_dump(j, purge_old_dumps=True, keep_first=1, keep_last=3) + dump_helper.write_dump(j, purge_old_dumps=True, keep_first=1, keep_last=5) if tlist[j] >= t_end: break diff --git a/src/shammodels/sph/src/modules/NodeUpdateDerivsMonofluidTVI.cpp b/src/shammodels/sph/src/modules/NodeUpdateDerivsMonofluidTVI.cpp index b50b8a0259..282bd2f786 100644 --- a/src/shammodels/sph/src/modules/NodeUpdateDerivsMonofluidTVI.cpp +++ b/src/shammodels/sph/src/modules/NodeUpdateDerivsMonofluidTVI.cpp @@ -98,6 +98,18 @@ struct KernelUpdateDerivsMonofluidTVI { Tscal delta_P = P_a - P_b; Tscal Ts_weighted = (Ttilde_sj_a / rho_a + Ttilde_sj_b / rho_b); + + if(id_a == 34559 && jdust == 4){ + printf("--------------------------------\n"); + printf("id_b: %d, s_j_b: %e, Ttilde_sj_b: %e\n", id_b, s_j_b, Ttilde_sj_b); + printf("rho_b: %e, P_b: %e, omega_b: %f\n", rho_b, P_b, omega_b); + printf("Fab_a: %e, Fab_b: %e\n", Fab_a, Fab_b); + printf("F_ab_bar: %e, delta_P: %e, Ts_weighted: %e\n", F_ab_bar, delta_P, Ts_weighted); + printf("rab: %e, rab_inv_sat: %e\n", rab, rab_inv_sat); + printf("v_ab: %f, %f, %f\n", v_ab[0], v_ab[1], v_ab[2]); + printf("r_ab_unit: %f, %f, %f\n", r_ab_unit[0], r_ab_unit[1], r_ab_unit[2]); + } + term1 += (pmass * s_j_b / rho_b) * Ts_weighted * delta_P * F_ab_bar * rab_inv_sat; term2 += pmass * sham::dot(v_ab, r_ab_unit * Fab_a); }); @@ -111,6 +123,21 @@ struct KernelUpdateDerivsMonofluidTVI { // restore it slowly to 0 ds_j_dt_a += (s_j_a < 0) ? -s_j_a / (10 * Ttilde_sj_a) : 0; + if(id_a == 34559){ + printf("id_a: %d, s_j_a: %e, ds_j_dt_a: %e\n", id_a, s_j_a, ds_j_dt_a); + printf("h_a: %e, rho_a: %e, omega_a: %f, Ttilde_sj_a: %f\n", h_a, rho_a, omega_a, Ttilde_sj_a); + printf("term1: %e, term2: %e\n", Tscal{-0.5} * term1, (s_j_a / (2 * rho_a * omega_a)) * term2); + printf("xyz_a: %f, %f, %f\n", xyz_a[0], xyz_a[1], xyz_a[2]); + } + + if(ds_j_dt_a > 1e-1){ + printf("id_a: %d, s_j_a: %e, ds_j_dt_a: %e\n", id_a, s_j_a, ds_j_dt_a); + printf("h_a: %e, rho_a: %e, omega_a: %f, Ttilde_sj_a: %f\n", h_a, rho_a, omega_a, Ttilde_sj_a); + printf("term1: %e, term2: %e\n", Tscal{-0.5} * term1, (s_j_a / (2 * rho_a * omega_a)) * term2); + printf("xyz_a: %f, %f, %f\n", xyz_a[0], xyz_a[1], xyz_a[2]); + throw std::runtime_error("ds_j_dt_a is too large"); + } + ds_j_dt[thread_id] = ds_j_dt_a; } }; From 3a4bc2ce1a0064eb64bd2db2ac697743e6bc1196 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Timoth=C3=A9e=20David--Cl=C3=A9ris?= Date: Wed, 8 Jul 2026 20:15:32 +0200 Subject: [PATCH 11/24] fix --- cmake/ShamConfigureCppTrace.cmake | 1 + 1 file changed, 1 insertion(+) diff --git a/cmake/ShamConfigureCppTrace.cmake b/cmake/ShamConfigureCppTrace.cmake index 7194e737b2..dff5ac4d4e 100644 --- a/cmake/ShamConfigureCppTrace.cmake +++ b/cmake/ShamConfigureCppTrace.cmake @@ -26,5 +26,6 @@ if(SHAMROCK_USE_CPPTRACE) ) set(CPPTRACE_INHERIT_HOST_STANDARD On) + set(CPPTRACE_DISABLE_CXX_20_MODULES On) add_subdirectory(external/cpptrace) endif() From bd35bfd3125a8fbfa9de185e8d7838cff48c8678 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Timoth=C3=A9e=20David--Cl=C3=A9ris?= Date: Wed, 8 Jul 2026 20:16:00 +0200 Subject: [PATCH 12/24] change debug part --- .../src/modules/NodeUpdateDerivsMonofluidTVI.cpp | 15 ++++++++++----- 1 file changed, 10 insertions(+), 5 deletions(-) diff --git a/src/shammodels/sph/src/modules/NodeUpdateDerivsMonofluidTVI.cpp b/src/shammodels/sph/src/modules/NodeUpdateDerivsMonofluidTVI.cpp index 282bd2f786..78de644a9c 100644 --- a/src/shammodels/sph/src/modules/NodeUpdateDerivsMonofluidTVI.cpp +++ b/src/shammodels/sph/src/modules/NodeUpdateDerivsMonofluidTVI.cpp @@ -99,7 +99,7 @@ struct KernelUpdateDerivsMonofluidTVI { Tscal Ts_weighted = (Ttilde_sj_a / rho_a + Ttilde_sj_b / rho_b); - if(id_a == 34559 && jdust == 4){ + if(id_a == 2320 && jdust == 4 && false){ printf("--------------------------------\n"); printf("id_b: %d, s_j_b: %e, Ttilde_sj_b: %e\n", id_b, s_j_b, Ttilde_sj_b); printf("rho_b: %e, P_b: %e, omega_b: %f\n", rho_b, P_b, omega_b); @@ -123,10 +123,15 @@ struct KernelUpdateDerivsMonofluidTVI { // restore it slowly to 0 ds_j_dt_a += (s_j_a < 0) ? -s_j_a / (10 * Ttilde_sj_a) : 0; - if(id_a == 34559){ - printf("id_a: %d, s_j_a: %e, ds_j_dt_a: %e\n", id_a, s_j_a, ds_j_dt_a); - printf("h_a: %e, rho_a: %e, omega_a: %f, Ttilde_sj_a: %f\n", h_a, rho_a, omega_a, Ttilde_sj_a); - printf("term1: %e, term2: %e\n", Tscal{-0.5} * term1, (s_j_a / (2 * rho_a * omega_a)) * term2); + if(id_a == 2320 && jdust == 4){ + printf("s_j_a: %e\n", s_j_a); + printf("ds_j_dt_a: %e\n", ds_j_dt_a); + printf("h_a: %e\n", h_a); + printf("rho_a: %e\n", rho_a); + printf("omega_a: %f\n", omega_a); + printf("Ttilde_sj_a: %f\n", Ttilde_sj_a); + printf("term1: %e\n", Tscal{-0.5} * term1); + printf("term2: %e\n", (s_j_a / (2 * rho_a * omega_a)) * term2); printf("xyz_a: %f, %f, %f\n", xyz_a[0], xyz_a[1], xyz_a[2]); } From 09273b86b9701e4b26166470a4d78fbf7fc5771e Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Timoth=C3=A9e=20David--Cl=C3=A9ris?= Date: Wed, 8 Jul 2026 20:25:11 +0200 Subject: [PATCH 13/24] fix --- .../sph/modules/ComputeCFLDust1Fluid.hpp | 125 ++++++++++++++++++ src/shammodels/sph/src/Solver.cpp | 106 +++++++++++++++ .../modules/NodeUpdateDerivsMonofluidTVI.cpp | 4 +- .../sph/src/modules/UpdateDerivs.cpp | 56 +------- 4 files changed, 235 insertions(+), 56 deletions(-) create mode 100644 src/shammodels/sph/include/shammodels/sph/modules/ComputeCFLDust1Fluid.hpp diff --git a/src/shammodels/sph/include/shammodels/sph/modules/ComputeCFLDust1Fluid.hpp b/src/shammodels/sph/include/shammodels/sph/modules/ComputeCFLDust1Fluid.hpp new file mode 100644 index 0000000000..5686a0c59f --- /dev/null +++ b/src/shammodels/sph/include/shammodels/sph/modules/ComputeCFLDust1Fluid.hpp @@ -0,0 +1,125 @@ +// -------------------------------------------------------// +// +// SHAMROCK code for hydrodynamics +// Copyright (c) 2021-2026 Timothée David--Cléris +// SPDX-License-Identifier: CeCILL Free Software License Agreement v2.1 +// Shamrock is licensed under the CeCILL 2.1 License, see LICENSE for more information +// +// -------------------------------------------------------// + +#pragma once + +/** + * @file ComputeCFLForce.hpp + * @author Timothée David--Cléris (tim.shamrock@proton.me) + * @brief + * + */ + +#include "shambackends/kernel_call_distrib.hpp" +#include "shammodels/sph/math/density.hpp" +#include "shamrock/solvergraph/IFieldSpan.hpp" +#include "shamrock/solvergraph/INode.hpp" +#include "shamrock/solvergraph/Indexes.hpp" +#include "shamrock/solvergraph/ScalarEdge.hpp" +#include "shamsys/NodeInstance.hpp" + +#define NODE_EDGES(X_RO, X_RW) \ + X_RO(shamrock::solvergraph::Indexes, part_counts) \ + X_RO(shamrock::solvergraph::ScalarEdge, C_1_fluid) \ + X_RO(shamrock::solvergraph::ScalarEdge, pmass) \ + X_RO(shamrock::solvergraph::ScalarEdge, hfactd) \ + X_RO(shamrock::solvergraph::IFieldSpan, hpart) \ + X_RO(shamrock::solvergraph::IFieldSpan, soundspeed) \ + X_RO(shamrock::solvergraph::IFieldSpan, s_j) \ + X_RO(shamrock::solvergraph::IFieldSpan, Ts_j) \ + X_RW(shamrock::solvergraph::IFieldSpan, cfl_dt) + +template +class ComputeCFLDust1Fluid : public shamrock::solvergraph::INode { + + using Tscal = shambase::VecComponent; + + u32 nbins; + + public: + ComputeCFLDust1Fluid(u32 nbins) : nbins(nbins) {} + + EXPAND_NODE_EDGES(NODE_EDGES) + + inline void _impl_evaluate_internal() { + auto edges = get_edges(); + + auto dev_sched = shamsys::instance::get_compute_scheduler_ptr(); + + Tscal C_1_fluid = edges.C_1_fluid.value; + Tscal pmass = edges.pmass.value; + Tscal hfactd = edges.hfactd.value; + + sham::distributed_data_kernel_call( + dev_sched, + sham::DDMultiRef{ + edges.hpart.get_spans(), + edges.soundspeed.get_spans(), + edges.s_j.get_spans(), + edges.Ts_j.get_spans()}, + sham::DDMultiRef{edges.cfl_dt.get_spans()}, + edges.part_counts.indexes, + [C_1_fluid, pmass, hfactd, nbins = this->nbins]( + u32 id_a, + const Tscal *hpart, + const Tscal *soundspeed, + const Tscal *s_j, + const Tscal *Ts_j, + Tscal *cfl_dt) { + u32 id_a_d = id_a * nbins; + + Tscal h_a = hpart[id_a]; + Tscal rho_a = shamrock::sph::rho_h(pmass, h_a, hfactd); + + Tscal cs_a = soundspeed[id_a]; + Tscal cs2_a = cs_a * cs_a; + + auto rho_dust = [&](int j) { + auto tmp = s_j[id_a_d + j]; + return tmp * tmp; + }; + + auto epsilon_j = [&](int j) { + return rho_dust(j) / rho_a; + }; + + Tscal sum_eps = 0; + for (int j = 0; j < nbins; j++) { + sum_eps += epsilon_j(j); + } + + Tscal cs_tilde_2_a = cs2_a * (1 - sum_eps); + + Tscal cs4_over_h2 = cs2_a * cs2_a / (h_a * h_a); + + Tscal cfl_tmp = std::numeric_limits::infinity(); + + for (int j = 0; j < nbins; j++) { + Tscal eps_j = epsilon_j(j); + Tscal eps2_j = eps_j * eps_j; + + Tscal Ts_j = Ts_j[id_a_d + j]; + Tscal Ts2_j = Ts_j * Ts_j; + + Tscal dt_j = h_a / sycl::sqrt(cs_tilde_2_a + Ts2_j * eps2_j * cs4_over_h2); + cfl_tmp = sycl::min(cfl_tmp, dt_j); + } + + cfl_tmp *= C_1_fluid; + + cfl_dt[id_a] = sycl::min(cfl_dt[id_a], cfl_tmp); + }); + } + + inline virtual std::string _impl_get_label() const { return "ComputeCFLDust1Fluid"; }; + + inline virtual std::string _impl_get_tex() const { return "C_{1,fluid}"; }; +}; + +#undef NODE_EDGES diff --git a/src/shammodels/sph/src/Solver.cpp b/src/shammodels/sph/src/Solver.cpp index 1313b26aa4..4786907b66 100644 --- a/src/shammodels/sph/src/Solver.cpp +++ b/src/shammodels/sph/src/Solver.cpp @@ -47,6 +47,7 @@ #include "shammodels/sph/modules/BuildTrees.hpp" #include "shammodels/sph/modules/ComputeCFLCourant.hpp" #include "shammodels/sph/modules/ComputeCFLDivBCleaning.hpp" +#include "shammodels/sph/modules/ComputeCFLDust1Fluid.hpp" #include "shammodels/sph/modules/ComputeCFLForce.hpp" #include "shammodels/sph/modules/ComputeEos.hpp" #include "shammodels/sph/modules/ComputeLoadBalanceValue.hpp" @@ -64,6 +65,8 @@ #include "shammodels/sph/modules/LoopSmoothingLengthIter.hpp" #include "shammodels/sph/modules/NeighbourCache.hpp" #include "shammodels/sph/modules/ParticleReordering.hpp" +#include "shammodels/sph/modules/SetDustStoppingTimeConstant.hpp" +#include "shammodels/sph/modules/SetDustStoppingTimeEpstein.hpp" #include "shammodels/sph/modules/SinkParticlesUpdate.hpp" #include "shammodels/sph/modules/UpdateDerivs.hpp" #include "shammodels/sph/modules/UpdateViscosity.hpp" @@ -156,6 +159,9 @@ void shammodels::sph::Solver::init_solver_graph() { if (has_s_j_field) { solver_graph.register_edge("s_j", FieldRefs("s_j", "S_j")); solver_graph.register_edge("ds_j_dt", FieldRefs("ds_j_dt", "dS_j/dt")); + + u32 ndust = solver_config.dust_config.get_dust_nvar(); + solver_graph.register_edge("Ts_j", Field(ndust, "Ts_j", "Ts_j")); } { @@ -1549,9 +1555,98 @@ void shammodels::sph::Solver::prepare_corrector() { template class Kern> void shammodels::sph::Solver::update_derivs(Tscal dt_hydro) { + // if one fluid is enabled time to compute the stopping times + if (solver_config.dust_config.has_s_j_field() ) { + + auto &cfg = solver_config.dust_config; + u32 ndust = cfg.get_dust_nvar(); + + using DustConfig = typename Config::DustConfig; + + using None = typename DustConfig::None; + using ConstantStoppingTimes = typename DustConfig::ConstantStoppingTimes; + using EpsteinDrag = typename DustConfig::EpsteinDrag; + + shamrock::patch::PatchDataLayerLayout &ghost_layout + = shambase::get_check_ref(storage.ghost_layout.get()); + u32 ihpart_interf = ghost_layout.get_field_idx("hpart"); + + auto &part_counts_with_ghost = storage.part_counts_with_ghost; + shambase::DistributedData &mpdats + = storage.merged_patchdata_ghost.get(); + + std::shared_ptr> hpart_refs + = std::make_shared>("hpart", "h"); + { // if was just reset before this call + shambase::get_check_ref(hpart_refs) + .set_refs(mpdats.map>>( + [&](u64 id, shamrock::patch::PatchDataLayer &mpdat) { + return std::ref(mpdat.get_field(ihpart_interf)); + })); + } + + auto gpart_mass + = storage.solver_graph.template get_edge_ptr>( + "gpart_mass"); + auto t_j_field + = storage.solver_graph.template get_edge_ptr>( + "Ts_j"); + + if (std::holds_alternative(cfg.dust_drag_mode)) { + + throw "bro WTF"; + + } else if ( + ConstantStoppingTimes *cfg_drag + = std::get_if(&cfg.dust_drag_mode)) { + + std::shared_ptr>> input_t_j + = std::make_shared>>("", ""); + input_t_j->value = cfg_drag->stopping_times; + + std::shared_ptr> node_set_tj + = std::make_shared>(ndust); + { + node_set_tj->set_edges(input_t_j, part_counts_with_ghost, t_j_field); + } + node_set_tj->evaluate(); + + } else if (EpsteinDrag *cfg_drag = std::get_if(&cfg.dust_drag_mode)) { + + std::shared_ptr> input_gamma + = std::make_shared>("", ""); + input_gamma->value = cfg_drag->gamma; + + std::shared_ptr>> input_sgrain_j + = std::make_shared>>("", ""); + input_sgrain_j->value = cfg_drag->grains_sizes; + + std::shared_ptr>> input_rho_grain_j + = std::make_shared>>("", ""); + input_rho_grain_j->value = cfg_drag->grains_densities; + + std::shared_ptr> node_set_tj + = std::make_shared>(ndust); + { + node_set_tj->set_edges( + gpart_mass, + input_gamma, + input_sgrain_j, + input_rho_grain_j, + part_counts_with_ghost, + hpart_refs, + storage.soundspeed, + t_j_field); + } + node_set_tj->evaluate(); + } + } + modules::UpdateDerivs derivs(context, solver_config, storage); derivs.update_derivs(dt_hydro); + logger::raw_ln("llfaofzefiuzehfiuzehfuizehfiuzehfiuzehfiuzhefuizehfz"); + modules::ExternalForces ext_forces(context, solver_config, storage); ext_forces.add_ext_forces(); } @@ -2607,6 +2702,17 @@ shammodels::sph::TimestepLog shammodels::sph::Solver::evolve_once() storage.part_counts, C_cour_edge, hpart_refs, vclean_dt, cfl_dt); } + // std::shared_ptr> compute_cfl_dust1_fluid; + // if(solver_config.dust_config.has_s_j_field()) { + // u32 ndust = solver_config.dust_config.get_dust_nvar(); + // + // compute_cfl_dust1_fluid = std::make_shared>(ndust); + // + // compute_cfl_dust1_fluid->set_edges( + // storage.part_counts, C_cour_edge, hpart_refs, storage.soundspeed, s_j_refs, + // Ts_j_refs, cfl_dt); + //} + bool show_cfl_detail = solver_config.show_cfl_detail; std::vector> cfl_detail; diff --git a/src/shammodels/sph/src/modules/NodeUpdateDerivsMonofluidTVI.cpp b/src/shammodels/sph/src/modules/NodeUpdateDerivsMonofluidTVI.cpp index 78de644a9c..f7f3e3db39 100644 --- a/src/shammodels/sph/src/modules/NodeUpdateDerivsMonofluidTVI.cpp +++ b/src/shammodels/sph/src/modules/NodeUpdateDerivsMonofluidTVI.cpp @@ -123,7 +123,7 @@ struct KernelUpdateDerivsMonofluidTVI { // restore it slowly to 0 ds_j_dt_a += (s_j_a < 0) ? -s_j_a / (10 * Ttilde_sj_a) : 0; - if(id_a == 2320 && jdust == 4){ + if(id_a == 2320 && jdust == 4&& false){ printf("s_j_a: %e\n", s_j_a); printf("ds_j_dt_a: %e\n", ds_j_dt_a); printf("h_a: %e\n", h_a); @@ -135,7 +135,7 @@ struct KernelUpdateDerivsMonofluidTVI { printf("xyz_a: %f, %f, %f\n", xyz_a[0], xyz_a[1], xyz_a[2]); } - if(ds_j_dt_a > 1e-1){ + if(ds_j_dt_a > 1e-1&& false){ printf("id_a: %d, s_j_a: %e, ds_j_dt_a: %e\n", id_a, s_j_a, ds_j_dt_a); printf("h_a: %e, rho_a: %e, omega_a: %f, Ttilde_sj_a: %f\n", h_a, rho_a, omega_a, Ttilde_sj_a); printf("term1: %e, term2: %e\n", Tscal{-0.5} * term1, (s_j_a / (2 * rho_a * omega_a)) * term2); diff --git a/src/shammodels/sph/src/modules/UpdateDerivs.cpp b/src/shammodels/sph/src/modules/UpdateDerivs.cpp index 31487ba7f0..f00e2d0343 100644 --- a/src/shammodels/sph/src/modules/UpdateDerivs.cpp +++ b/src/shammodels/sph/src/modules/UpdateDerivs.cpp @@ -1151,60 +1151,8 @@ void shammodels::sph::modules::UpdateDerivs::update_derivs_dust })); } - std::shared_ptr> t_j_field - = std::make_shared>(ndust, "t_j", "t_j"); - - using None = typename DustConfig::None; - using ConstantStoppingTimes = typename DustConfig::ConstantStoppingTimes; - using EpsteinDrag = typename DustConfig::EpsteinDrag; - - if (std::holds_alternative(cfg.dust_drag_mode)) { - - throw "bro WTF"; - - } else if ( - ConstantStoppingTimes *cfg_drag = std::get_if(&cfg.dust_drag_mode)) { - - std::shared_ptr>> input_t_j - = std::make_shared>>("", ""); - input_t_j->value = cfg_drag->stopping_times; - - std::shared_ptr> node_set_tj - = std::make_shared>(ndust); - { - node_set_tj->set_edges(input_t_j, part_counts_with_ghost, t_j_field); - } - node_set_tj->evaluate(); - - } else if (EpsteinDrag *cfg_drag = std::get_if(&cfg.dust_drag_mode)) { - - std::shared_ptr> input_gamma - = std::make_shared>("", ""); - input_gamma->value = cfg_drag->gamma; - - std::shared_ptr>> input_sgrain_j - = std::make_shared>>("", ""); - input_sgrain_j->value = cfg_drag->grains_sizes; - - std::shared_ptr>> input_rho_grain_j - = std::make_shared>>("", ""); - input_rho_grain_j->value = cfg_drag->grains_densities; - - std::shared_ptr> node_set_tj - = std::make_shared>(ndust); - { - node_set_tj->set_edges( - gpart_mass, - input_gamma, - input_sgrain_j, - input_rho_grain_j, - part_counts_with_ghost, - hpart_refs, - storage.soundspeed, - t_j_field); - } - node_set_tj->evaluate(); - } + auto t_j_field + = storage.solver_graph.template get_edge_ptr>("Ts_j"); std::shared_ptr> Ttilde_sj_field = std::make_shared>(ndust, "Ttilde_sj", "Ttilde_sj"); From 71122150121f72252e68d9baad58b80f7f403775 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Timoth=C3=A9e=20David--Cl=C3=A9ris?= Date: Wed, 8 Jul 2026 22:26:51 +0200 Subject: [PATCH 14/24] fix --- examples/sph/run_dustysettle_tvi.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/examples/sph/run_dustysettle_tvi.py b/examples/sph/run_dustysettle_tvi.py index 434ad8d535..6ebbd8edd3 100644 --- a/examples/sph/run_dustysettle_tvi.py +++ b/examples/sph/run_dustysettle_tvi.py @@ -186,6 +186,7 @@ def setup_model(): cfg.set_boundary_periodic() cfg.set_units(codeu) cfg.set_eos_isothermal(cs) + cfg.set_show_cfl_detail(True) cfg.print_status() model.set_solver_config(cfg) @@ -224,8 +225,8 @@ def f_remap(r): model.remap_positions(f_remap) model.set_field_value_lambda_f64("uint", uint_g) - model.set_cfl_cour(0.25) - model.set_cfl_force(0.25) + model.set_cfl_cour(0.1) + model.set_cfl_force(0.1) model.timestep() From 4c7fb8d7f5afbd5ce4342d23902fa33013237aa0 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Timoth=C3=A9e=20David--Cl=C3=A9ris?= Date: Wed, 8 Jul 2026 22:32:37 +0200 Subject: [PATCH 15/24] fix --- src/shammodels/sph/src/Solver.cpp | 2 -- 1 file changed, 2 deletions(-) diff --git a/src/shammodels/sph/src/Solver.cpp b/src/shammodels/sph/src/Solver.cpp index 4786907b66..7b9f0295bc 100644 --- a/src/shammodels/sph/src/Solver.cpp +++ b/src/shammodels/sph/src/Solver.cpp @@ -1645,8 +1645,6 @@ void shammodels::sph::Solver::update_derivs(Tscal dt_hydro) { modules::UpdateDerivs derivs(context, solver_config, storage); derivs.update_derivs(dt_hydro); - logger::raw_ln("llfaofzefiuzehfiuzehfuizehfiuzehfiuzehfiuzhefuizehfz"); - modules::ExternalForces ext_forces(context, solver_config, storage); ext_forces.add_ext_forces(); } From b1430680bc115b0e9111d7a5de94cd4d8eaa1f2e Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Timoth=C3=A9e=20David--Cl=C3=A9ris?= Date: Thu, 9 Jul 2026 00:28:38 +0200 Subject: [PATCH 16/24] more tweaks --- examples/sph/run_dustysettle_tvi.py | 6 +- .../sph/modules/ComputeCFLDust1Fluid.hpp | 12 ++-- src/shammodels/sph/src/Solver.cpp | 71 ++++++++++++++++--- .../modules/NodeUpdateDerivsMonofluidTVI.cpp | 38 ++++++---- 4 files changed, 93 insertions(+), 34 deletions(-) diff --git a/examples/sph/run_dustysettle_tvi.py b/examples/sph/run_dustysettle_tvi.py index 6ebbd8edd3..73b21892e1 100644 --- a/examples/sph/run_dustysettle_tvi.py +++ b/examples/sph/run_dustysettle_tvi.py @@ -59,7 +59,7 @@ epsilon_base = 0.01 -tlist = [0.01 * i for i in range(1000)] +tlist = [0.01 * i for i in range(1000)] tinject = 0 iinject = 0 t_end = 5.0 @@ -127,7 +127,7 @@ def uint_g(r): rho_grains_si_edges = np.array([2.3 * 1000 for i in range(ndust + 1)]) -grain_size_si_edges = np.logspace(-5, -3, ndust + 1) +grain_size_si_edges = np.logspace(-5, -2, ndust + 1) mrn_distribution = DustMRNDistribution( codeu, mrn_pow, mrn_cutoff_si, grain_size_si_edges, rho_grains_si_edges @@ -216,7 +216,7 @@ def f_remap(r): rn = max(abs(zM), abs(zm)) # print(y, H, H * erfinv(y / rn)) - z = H * erfinv(z / rn) * (2**0.5) + z = H * erfinv(z / rn) * (2**0.5) z = min(z, zM) z = max(z, zm) diff --git a/src/shammodels/sph/include/shammodels/sph/modules/ComputeCFLDust1Fluid.hpp b/src/shammodels/sph/include/shammodels/sph/modules/ComputeCFLDust1Fluid.hpp index 5686a0c59f..a81687ed4d 100644 --- a/src/shammodels/sph/include/shammodels/sph/modules/ComputeCFLDust1Fluid.hpp +++ b/src/shammodels/sph/include/shammodels/sph/modules/ComputeCFLDust1Fluid.hpp @@ -10,7 +10,7 @@ #pragma once /** - * @file ComputeCFLForce.hpp + * @file ComputeCFLDust1Fluid.hpp * @author Timothée David--Cléris (tim.shamrock@proton.me) * @brief * @@ -101,13 +101,13 @@ class ComputeCFLDust1Fluid : public shamrock::solvergraph::INode { Tscal cfl_tmp = std::numeric_limits::infinity(); for (int j = 0; j < nbins; j++) { - Tscal eps_j = epsilon_j(j); - Tscal eps2_j = eps_j * eps_j; + Tscal eps_j_a = epsilon_j(j); + Tscal eps2_j_a = eps_j_a * eps_j_a; - Tscal Ts_j = Ts_j[id_a_d + j]; - Tscal Ts2_j = Ts_j * Ts_j; + Tscal Ts_j_a = Ts_j[id_a_d + j]; + Tscal Ts2_j_a = Ts_j_a * Ts_j_a; - Tscal dt_j = h_a / sycl::sqrt(cs_tilde_2_a + Ts2_j * eps2_j * cs4_over_h2); + Tscal dt_j = h_a / sycl::sqrt(cs_tilde_2_a + Ts2_j_a * eps2_j_a * cs4_over_h2); cfl_tmp = sycl::min(cfl_tmp, dt_j); } diff --git a/src/shammodels/sph/src/Solver.cpp b/src/shammodels/sph/src/Solver.cpp index ac3b3b18bc..6c68a21e50 100644 --- a/src/shammodels/sph/src/Solver.cpp +++ b/src/shammodels/sph/src/Solver.cpp @@ -2700,16 +2700,43 @@ shammodels::sph::TimestepLog shammodels::sph::Solver::evolve_once() storage.part_counts, C_cour_edge, hpart_refs, vclean_dt, cfl_dt); } - // std::shared_ptr> compute_cfl_dust1_fluid; - // if(solver_config.dust_config.has_s_j_field()) { - // u32 ndust = solver_config.dust_config.get_dust_nvar(); - // - // compute_cfl_dust1_fluid = std::make_shared>(ndust); - // - // compute_cfl_dust1_fluid->set_edges( - // storage.part_counts, C_cour_edge, hpart_refs, storage.soundspeed, s_j_refs, - // Ts_j_refs, cfl_dt); - //} + std::shared_ptr> compute_cfl_dust1_fluid; + std::shared_ptr> s_j_refs; + std::shared_ptr> hfactd_edge; + + if (solver_config.dust_config.has_s_j_field()) { + u32 ndust = solver_config.dust_config.get_dust_nvar(); + + compute_cfl_dust1_fluid = std::make_shared>(ndust); + + auto t_j_field + = storage.solver_graph + .template get_edge_ptr>("Ts_j"); + + auto pmass_edge + = storage.solver_graph + .template get_edge_ptr>( + "gpart_mass"); + + s_j_refs = std::make_shared>("s_j", "s_j"); + + hfactd_edge = std::make_shared>( + "hfactd", "hfactd"); + hfactd_edge->value = Kernel::hfactd; + + map_field_refs(scheduler(), is_j, *s_j_refs); + + compute_cfl_dust1_fluid->set_edges( + storage.part_counts, + C_cour_edge, + pmass_edge, + hfactd_edge, + hpart_refs, + storage.soundspeed, + s_j_refs, + t_j_field, + cfl_dt); + } bool show_cfl_detail = solver_config.show_cfl_detail; std::vector> cfl_detail; @@ -2737,6 +2764,11 @@ shammodels::sph::TimestepLog shammodels::sph::Solver::evolve_once() save_cfl_detail("divB_cleaning"); } + if (solver_config.dust_config.has_s_j_field()) { + compute_cfl_dust1_fluid->evaluate(); + save_cfl_detail("dust1_fluid"); + } + if (!show_cfl_detail) { save_dt_min_to_dt_part(); cfl_detail.push_back({"all SPH", cfl_dt->get_native().compute_rank_min()}); @@ -2871,6 +2903,25 @@ shammodels::sph::TimestepLog shammodels::sph::Solver::evolve_once() reset_merge_ghosts_fields(); reset_eos_fields(); + // restore positivity + if (solver_config.dust_config.has_s_j_field()) { + auto ndust = solver_config.dust_config.get_dust_nvar(); + auto &q = shamsys::instance::get_compute_scheduler().get_queue(); + scheduler().for_each_patchdata_nonempty([&](Patch cur_p, PatchDataLayer &pdat) { + sham::DeviceBuffer &buf_s_j = pdat.get_field_buf_ref(is_j); + sham::kernel_call( + q, + sham::MultiRef{}, + sham::MultiRef{buf_s_j}, + pdat.get_obj_cnt()*ndust, + [](u32 tid, Tscal *s_j) { + if (s_j[tid] < 0) { + s_j[tid] = 0; + } + }); + }); + } + // if delta too big jump to compute force tstep.stop(); diff --git a/src/shammodels/sph/src/modules/NodeUpdateDerivsMonofluidTVI.cpp b/src/shammodels/sph/src/modules/NodeUpdateDerivsMonofluidTVI.cpp index f7f3e3db39..5c5ed7fd7d 100644 --- a/src/shammodels/sph/src/modules/NodeUpdateDerivsMonofluidTVI.cpp +++ b/src/shammodels/sph/src/modules/NodeUpdateDerivsMonofluidTVI.cpp @@ -98,17 +98,17 @@ struct KernelUpdateDerivsMonofluidTVI { Tscal delta_P = P_a - P_b; Tscal Ts_weighted = (Ttilde_sj_a / rho_a + Ttilde_sj_b / rho_b); - - if(id_a == 2320 && jdust == 4 && false){ - printf("--------------------------------\n"); - printf("id_b: %d, s_j_b: %e, Ttilde_sj_b: %e\n", id_b, s_j_b, Ttilde_sj_b); - printf("rho_b: %e, P_b: %e, omega_b: %f\n", rho_b, P_b, omega_b); - printf("Fab_a: %e, Fab_b: %e\n", Fab_a, Fab_b); - printf("F_ab_bar: %e, delta_P: %e, Ts_weighted: %e\n", F_ab_bar, delta_P, Ts_weighted); - printf("rab: %e, rab_inv_sat: %e\n", rab, rab_inv_sat); - printf("v_ab: %f, %f, %f\n", v_ab[0], v_ab[1], v_ab[2]); - printf("r_ab_unit: %f, %f, %f\n", r_ab_unit[0], r_ab_unit[1], r_ab_unit[2]); - } + if (id_a == 2320 && jdust == 4 && false) { + printf("--------------------------------\n"); + printf("id_b: %d, s_j_b: %e, Ttilde_sj_b: %e\n", id_b, s_j_b, Ttilde_sj_b); + printf("rho_b: %e, P_b: %e, omega_b: %f\n", rho_b, P_b, omega_b); + printf("Fab_a: %e, Fab_b: %e\n", Fab_a, Fab_b); + printf( + "F_ab_bar: %e, delta_P: %e, Ts_weighted: %e\n", F_ab_bar, delta_P, Ts_weighted); + printf("rab: %e, rab_inv_sat: %e\n", rab, rab_inv_sat); + printf("v_ab: %f, %f, %f\n", v_ab[0], v_ab[1], v_ab[2]); + printf("r_ab_unit: %f, %f, %f\n", r_ab_unit[0], r_ab_unit[1], r_ab_unit[2]); + } term1 += (pmass * s_j_b / rho_b) * Ts_weighted * delta_P * F_ab_bar * rab_inv_sat; term2 += pmass * sham::dot(v_ab, r_ab_unit * Fab_a); @@ -123,7 +123,7 @@ struct KernelUpdateDerivsMonofluidTVI { // restore it slowly to 0 ds_j_dt_a += (s_j_a < 0) ? -s_j_a / (10 * Ttilde_sj_a) : 0; - if(id_a == 2320 && jdust == 4&& false){ + if (id_a == 2320 && jdust == 4 && false) { printf("s_j_a: %e\n", s_j_a); printf("ds_j_dt_a: %e\n", ds_j_dt_a); printf("h_a: %e\n", h_a); @@ -135,10 +135,18 @@ struct KernelUpdateDerivsMonofluidTVI { printf("xyz_a: %f, %f, %f\n", xyz_a[0], xyz_a[1], xyz_a[2]); } - if(ds_j_dt_a > 1e-1&& false){ + if (ds_j_dt_a > 1e-1 && false) { printf("id_a: %d, s_j_a: %e, ds_j_dt_a: %e\n", id_a, s_j_a, ds_j_dt_a); - printf("h_a: %e, rho_a: %e, omega_a: %f, Ttilde_sj_a: %f\n", h_a, rho_a, omega_a, Ttilde_sj_a); - printf("term1: %e, term2: %e\n", Tscal{-0.5} * term1, (s_j_a / (2 * rho_a * omega_a)) * term2); + printf( + "h_a: %e, rho_a: %e, omega_a: %f, Ttilde_sj_a: %f\n", + h_a, + rho_a, + omega_a, + Ttilde_sj_a); + printf( + "term1: %e, term2: %e\n", + Tscal{-0.5} * term1, + (s_j_a / (2 * rho_a * omega_a)) * term2); printf("xyz_a: %f, %f, %f\n", xyz_a[0], xyz_a[1], xyz_a[2]); throw std::runtime_error("ds_j_dt_a is too large"); } From 3c78519465e74717aaa464bd9b75f74cd1a7c0d8 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Timoth=C3=A9e=20David--Cl=C3=A9ris?= Date: Fri, 10 Jul 2026 00:16:51 +0200 Subject: [PATCH 17/24] welcome to my nightmare --- examples/physics/dustywave_sympy.py | 239 ++++++++++++++++++ examples/sph/run_dustysettle_tvi.py | 2 +- .../common/modules/ForwardEulerPositive.hpp | 96 +++++++ src/shammodels/sph/src/Solver.cpp | 100 +++++++- .../modules/NodeUpdateDerivsMonofluidTVI.cpp | 30 +-- .../include/shamrock/math/integrators.hpp | 10 + .../shamrock/scheduler/SchedulerUtility.hpp | 37 +++ src/shamrock/src/math/integrators.cpp | 45 ++++ 8 files changed, 535 insertions(+), 24 deletions(-) create mode 100644 examples/physics/dustywave_sympy.py create mode 100644 src/shammodels/common/include/shammodels/common/modules/ForwardEulerPositive.hpp diff --git a/examples/physics/dustywave_sympy.py b/examples/physics/dustywave_sympy.py new file mode 100644 index 0000000000..0da1947762 --- /dev/null +++ b/examples/physics/dustywave_sympy.py @@ -0,0 +1,239 @@ +""" +Dustywave TVA dispersion relation +======================================= + +This example shows how to derive the dustywave TVA dispersion relation using SymPy +""" + +import sympy as sp +import numpy as np +import matplotlib.pyplot as plt + +# %% +# Usefull symbols +omega = sp.symbols(r'\omega', complex=True) + +k, cs, ts, eps = sp.symbols( + r'k c_s t_s \epsilon', + positive=True, + real=True, +) + +i = sp.I # imaginary unit + +a = k*cs # s^-1 +b = a**2 * ts*eps # s^-1 + + +# %% +# The perturbation matrix is M = K(omega = 0) +# The eigenfrequencies are the roots of det(M + i omega I) = 0 +K = sp.Matrix([ + [i*omega, 0, -i*a], + [b*(1-eps), i*omega-b, 0], + [-i*a*(1-eps), i*a, i*omega] +]) # s^-1 + + +# %% +# Compute the determinant +det = sp.factor(K.det()) +print(sp.latex(det)) + +# %% +# Let's remove the leading omega mode (omega = 0) +det /= (i * omega) +det = sp.simplify(det) +det = sp.collect(det, omega) +print(sp.latex(det)) + +# %% +# Find the roots of the dispersion relation +r1,r2 = sp.solve(sp.Eq(det, 0), omega) +print(sp.latex(r1)) +print(sp.latex(r2)) + +print(r1) +print(r2) + +# %% +# Function to plot the roots +def get_roots(_r1,_r2,k_list, eps_value, cs_value, ts_value): + + # Substitute all parameters + r1_num = _r1.subs({ + cs: cs_value, + ts: ts_value, + eps: eps_value + }) + + r2_num = _r2.subs({ + cs: cs_value, + ts: ts_value, + eps: eps_value + }) + + r1_num_re = sp.re(r1_num) + r1_num_im = sp.im(r1_num) + r2_num_re = sp.re(r2_num) + r2_num_im = sp.im(r2_num) + + # Lambdify only k remains + r1_re_func = sp.lambdify(k, r1_num_re, modules="numpy") + r1_im_func = sp.lambdify(k, r1_num_im, modules="numpy") + r2_re_func = sp.lambdify(k, r2_num_re, modules="numpy") + r2_im_func = sp.lambdify(k, r2_num_im, modules="numpy") + + # Evaluate + r1_vals_re = r1_re_func(k_list) + r1_vals_im = r1_im_func(k_list) + r2_vals_re = r2_re_func(k_list) + r2_vals_im = r2_im_func(k_list) + + def restore(lst): + # if it is not a numpy array return a np.zeros_like(k_list) + if not isinstance(lst, np.ndarray): + return np.zeros_like(k_list) + return lst + + r1_vals_re = restore(r1_vals_re) + r1_vals_im = restore(r1_vals_im) + r2_vals_re = restore(r2_vals_re) + r2_vals_im = restore(r2_vals_im) + + return r1_vals_re, r1_vals_im, r2_vals_re, r2_vals_im + +def get_roots_LP14(k_list, eps_value, cs_value, ts_value): + _r1 = + cs * sp.sqrt(1-eps) * k - i * ts* k **2 * cs **2 * eps /2 + _r2 = - cs * sp.sqrt(1-eps) * k - i * ts* k **2 * cs **2 * eps /2 + + return get_roots(_r1, _r2, k_list, eps_value, cs_value, ts_value) + + + +def get_overroots_DCL26_simple(k_list, eps_value, cs_value, ts_value): + _r1 = + cs * sp.sqrt(1-eps) * k + i*k**2*cs**2 * eps*ts * (-1 + 1)/2 + _r2 = - cs * sp.sqrt(1-eps) * k + i*k**2*cs**2 * eps*ts *(-1 - 1)/2 + + return get_roots(_r1, _r2, k_list, eps_value, cs_value, ts_value) + +def get_overroots_DCL26(k_list, eps_value, cs_value, ts_value): + D = 4*(1-eps) - eps**2 * cs**2 * ts**2 * k**2 + + sqrtD_real = sp.sqrt(sp.Max(D, 0)) + sqrtD_imag = sp.sqrt(sp.Max(-D, 0)) + + _r1 = cs*k/2 * ( + +sqrtD_real + + i*(sqrtD_imag - eps*cs*k*ts) + ) + + _r2 = cs*k/2 * ( + -sqrtD_real + + i*(-sqrtD_imag - eps*cs*k*ts) + ) + + print(sp.latex(sp.Abs(_r1))) + print(sp.latex(sp.Abs(_r2))) + + return get_roots(_r1, _r2, k_list, eps_value, cs_value, ts_value) + +def plot_case(k_plot, eps_value, cs_value, ts_value): + + r1_vals_re, r1_vals_im, r2_vals_re, r2_vals_im = get_roots(r1,r2,k_plot, eps_value, cs_value, ts_value) + + r1_vals_re_LP14, r1_vals_im_LP14, r2_vals_re_LP14, r2_vals_im_LP14 = get_roots_LP14(k_plot, eps_value, cs_value, ts_value) + + r1_vals_re_DCL26_simple, r1_vals_im_DCL26_simple, r2_vals_re_DCL26_simple, r2_vals_im_DCL26_simple = get_overroots_DCL26_simple(k_plot, eps_value, cs_value, ts_value) + r1_vals_re_DCL26, r1_vals_im_DCL26, r2_vals_re_DCL26, r2_vals_im_DCL26 = get_overroots_DCL26(k_plot, eps_value, cs_value, ts_value) + + + # Create figure + fig, axs = plt.subplots(2, 2, figsize=(8, 8), sharex=True) + + # Real parts + axs[0,0].plot(k_plot, r1_vals_re,color="0",linewidth=2, label="Re($\omega_+$)") + axs[0,0].plot(k_plot, r2_vals_re,color="0",linewidth=2, label="Re($\omega_-$)") + axs[0,0].plot(k_plot, r1_vals_re_LP14,"--", label="Re($\omega_{+,LP14}$)") + axs[0,0].plot(k_plot, r2_vals_re_LP14,"--", label="Re($\omega_{-,LP14}$)") + axs[0,0].plot(k_plot, r1_vals_re_DCL26_simple,linestyle="dotted", label="Re($\omega_{+,approx}$)") + axs[0,0].plot(k_plot, r2_vals_re_DCL26_simple,linestyle="dotted", label="Re($\omega_{-,approx}$)") + #axs[0,0].plot(k_plot, r1_vals_re_DCL26,"--", label="Re($r_1$) DCL26") + #axs[0,0].plot(k_plot, r2_vals_re_DCL26,"--", label="Re($r_2$) DCL26") + axs[0,0].set_ylabel("Real part") + axs[0,0].grid(True) + axs[0,0].legend() + + # Imaginary parts + axs[0,1].plot(k_plot, r1_vals_im,color="0",linewidth=2, label="Im($\omega_+$)") + axs[0,1].plot(k_plot, r2_vals_im,color="0",linewidth=2, label="Im($\omega_-$)") + axs[0,1].plot(k_plot, r1_vals_im_LP14,"--", label="Im($\omega_{+,LP14}$)") + axs[0,1].plot(k_plot, r2_vals_im_LP14,"--", label="Im($\omega_{-,LP14}$)") + axs[0,1].plot(k_plot, r1_vals_im_DCL26_simple,linestyle="dotted", label="Im($\omega_{+,approx}$)") + axs[0,1].plot(k_plot, r2_vals_im_DCL26_simple,linestyle="dotted", label="Im($\omega_{-,approx}$)") + #axs[0,1].plot(k_plot, r1_vals_im_DCL26,"--", label="Im($r_1$) DCL26") + #axs[0,1].plot(k_plot, r2_vals_im_DCL26,"--", label="Im($r_2$) DCL26") + axs[0,1].set_xlabel("$k$") + axs[0,1].set_ylabel("Imaginary part") + axs[0,1].grid(True) + axs[0,1].legend() + + # Abs + r1_vals_abs = np.sqrt(r1_vals_re**2 + r1_vals_im**2) + r2_vals_abs = np.sqrt(r2_vals_re**2 + r2_vals_im**2) + r1_vals_abs_LP14 = np.sqrt(r1_vals_re_LP14**2 + r1_vals_im_LP14**2) + r2_vals_abs_LP14 = np.sqrt(r2_vals_re_LP14**2 + r2_vals_im_LP14**2) + r1_vals_abs_DCL26_simple = np.sqrt(r1_vals_re_DCL26_simple**2 + r1_vals_im_DCL26_simple**2) + r2_vals_abs_DCL26_simple = np.sqrt(r2_vals_re_DCL26_simple**2 + r2_vals_im_DCL26_simple**2) + r1_vals_abs_DCL26 = np.sqrt(r1_vals_re_DCL26**2 + r1_vals_im_DCL26**2) + r2_vals_abs_DCL26 = np.sqrt(r2_vals_re_DCL26**2 + r2_vals_im_DCL26**2) + axs[1,0].plot(k_plot, r1_vals_abs,color="0",linewidth=2, label="Abs($\omega_+$)") + axs[1,0].plot(k_plot, r2_vals_abs,color="0",linewidth=2, label="Abs($\omega_-$)") + axs[1,0].plot(k_plot, r1_vals_abs_LP14,"--", label="Abs($\omega_{+,LP14}$)") + axs[1,0].plot(k_plot, r2_vals_abs_LP14,"--", label="Abs($\omega_{-,LP14}$)") + axs[1,0].plot(k_plot, r1_vals_abs_DCL26_simple,linestyle="dotted", label="Abs($\omega_{+,approx}$)") + axs[1,0].plot(k_plot, r2_vals_abs_DCL26_simple,linestyle="dotted", label="Abs($\omega_{-,approx}$)") + + def approx(_k,_cs,_ts,_eps): + print(type(_k), type(_cs), type(_ts), type(_eps)) + return _cs * _k * np.sqrt((1-_eps) + (_k *_cs *_ts* _eps )**2) + + axs[1,0].plot(k_plot, approx(k_plot, cs_value, ts_value, eps_value), "--", label=r"$max(\vert \omega_{\pm,approx} \vert)$") + + #axs[1,0].plot(k_plot, r1_vals_abs_DCL26,"--", label="Abs($r_1$) DCL26") + #axs[1,0].plot(k_plot, r2_vals_abs_DCL26,"--", label="Abs($r_2$) DCL26") + axs[1,0].set_xlabel("$k$") + axs[1,0].set_ylabel("Abs part") + axs[1,0].grid(True) + axs[1,0].legend() + + # delta with max + r_max = np.maximum(r1_vals_abs,r2_vals_abs) + r_max_LP14 = np.maximum(r1_vals_abs_LP14,r2_vals_abs_LP14) + r_max_DCL26_simple = np.maximum(r1_vals_abs_DCL26_simple,r2_vals_abs_DCL26_simple) + r_max_DCL26 = np.maximum(r1_vals_abs_DCL26,r2_vals_abs_DCL26) + axs[1,1].plot(k_plot, (r_max_LP14 - r_max) / r_max, label="Ana - LP14") + axs[1,1].plot(k_plot, (r_max_DCL26_simple - r_max) / r_max, label="Ana - DCL26 simple") + #axs[1,1].plot(k_plot, (r_max_DCL26 - r_max) / r_max, label="Ana - DCL26") + axs[1,1].set_xlabel("$k$") + axs[1,1].set_ylabel("Abs(Ana) - Abs(Max model) / Abs(Ana)") + axs[1,1].grid(True) + axs[1,1].legend() + + plt.suptitle(f"eps = {eps_value}, cs = {cs_value}, ts = {ts_value}") + + plt.tight_layout() + + +# %% +# Plot the case eps = 0.5, cs = 1.0, ts = 1.0 +k_plot = np.linspace(0, 5, 1000) +plot_case(k_plot, 0.5, 1.0, 1.0) +plt.show() + + +# %% +# Plot the case eps = 0.1, cs = 1.0, ts = 1.0 +k_plot = np.linspace(0, 40, 1000) +plot_case(k_plot, 0.1, 1.0, 1.0) +plt.show() \ No newline at end of file diff --git a/examples/sph/run_dustysettle_tvi.py b/examples/sph/run_dustysettle_tvi.py index 73b21892e1..7385b833cf 100644 --- a/examples/sph/run_dustysettle_tvi.py +++ b/examples/sph/run_dustysettle_tvi.py @@ -59,7 +59,7 @@ epsilon_base = 0.01 -tlist = [0.01 * i for i in range(1000)] +tlist = [0.01 * i for i in range(50)] tinject = 0 iinject = 0 t_end = 5.0 diff --git a/src/shammodels/common/include/shammodels/common/modules/ForwardEulerPositive.hpp b/src/shammodels/common/include/shammodels/common/modules/ForwardEulerPositive.hpp new file mode 100644 index 0000000000..874fff1b80 --- /dev/null +++ b/src/shammodels/common/include/shammodels/common/modules/ForwardEulerPositive.hpp @@ -0,0 +1,96 @@ +// -------------------------------------------------------// +// +// SHAMROCK code for hydrodynamics +// Copyright (c) 2021-2026 Timothée David--Cléris +// SPDX-License-Identifier: CeCILL Free Software License Agreement v2.1 +// Shamrock is licensed under the CeCILL 2.1 License, see LICENSE for more information +// +// -------------------------------------------------------// + +#pragma once + +/** + * @file ForwardEuler.hpp + * @author Timothée David--Cléris (tim.shamrock@proton.me) + * @brief Implements a forward Euler integration step as a solver graph node. + * + */ + +#include "shambase/SourceLocation.hpp" +#include "shambackends/kernel_call_distrib.hpp" +#include "shamcomm/logs.hpp" +#include "shamrock/solvergraph/IDataEdge.hpp" +#include "shamrock/solvergraph/IFieldSpan.hpp" +#include "shamrock/solvergraph/INode.hpp" +#include "shamrock/solvergraph/Indexes.hpp" +#include "shamsys/NodeInstance.hpp" + +#define NODE_EDGES(X_RO, X_RW) \ + /* ------------------- inputs ------------------- */ \ + X_RO(shamrock::solvergraph::IDataEdge, dt) \ + X_RO(shamrock::solvergraph::IFieldSpan, time_derivative) \ + X_RO(shamrock::solvergraph::Indexes, sizes) \ + \ + /* ------------------- outputs ------------------- */ \ + X_RW(shamrock::solvergraph::IFieldSpan, field) + +namespace shammodels::common::modules { + template + class ForwardEulerPositive : public shamrock::solvergraph::INode { + + using Tscal = shambase::VecComponent; + + u32 nvar; + + public: + ForwardEulerPositive(u32 nvar = 1) : nvar(nvar) {} + + EXPAND_NODE_EDGES(NODE_EDGES) + + inline void _impl_evaluate_internal() { + + __shamrock_stack_entry(); + + auto edges = get_edges(); + + edges.field.ensure_sizes(edges.sizes.indexes); + + Tscal dt = edges.dt.data; + + if (nvar == 1) { + + sham::distributed_data_kernel_call( + shamsys::instance::get_compute_scheduler_ptr(), + sham::DDMultiRef{edges.time_derivative.get_spans()}, + sham::DDMultiRef{edges.field.get_spans()}, + edges.sizes.indexes, + [dt](u32 gid, const T *time_derivative, T *field) { + auto tmp = field[gid] + dt * time_derivative[gid]; + field[gid] = std::max(tmp, 0.0); + }); + + } else { + + auto var_count = edges.sizes.indexes.template map([&](u64 id, u32 count) { + return count * nvar; + }); + + sham::distributed_data_kernel_call( + shamsys::instance::get_compute_scheduler_ptr(), + sham::DDMultiRef{edges.time_derivative.get_spans()}, + sham::DDMultiRef{edges.field.get_spans()}, + var_count, + [dt](u32 gid, const T *time_derivative, T *field) { + auto tmp = field[gid] + dt * time_derivative[gid]; + field[gid] = std::max(tmp, 0.0); + }); + } + } + + inline virtual std::string _impl_get_label() const { return "ForwardEulerPositive"; } + + inline virtual std::string _impl_get_tex() const { return "TODO"; } + }; +} // namespace shammodels::common::modules + +#undef NODE_EDGES diff --git a/src/shammodels/sph/src/Solver.cpp b/src/shammodels/sph/src/Solver.cpp index 6c68a21e50..d33f00ef98 100644 --- a/src/shammodels/sph/src/Solver.cpp +++ b/src/shammodels/sph/src/Solver.cpp @@ -35,6 +35,7 @@ #include "shamcomm/wrapper.hpp" #include "shammath/sphkernels.hpp" #include "shammodels/common/modules/ForwardEuler.hpp" +#include "shammodels/common/modules/ForwardEulerPositive.hpp" #include "shammodels/common/timestep_report.hpp" #include "shammodels/sph/BasicSPHGhosts.hpp" #include "shammodels/sph/SPHUtilities.hpp" @@ -453,7 +454,7 @@ void shammodels::sph::Solver::init_solver_graph() { if (has_s_j_field) { u32 ndust = solver_config.dust_config.get_dust_nvar(); auto half_step_s_j = solver_graph.register_node( - prefix + "_s_j", shammodels::common::modules::ForwardEuler(ndust)); + prefix + "_s_j", shammodels::common::modules::ForwardEulerPositive(ndust)); shambase::get_check_ref(half_step_s_j) .set_edges( solver_graph.get_edge_ptr>("dt_half"), @@ -1835,8 +1836,85 @@ shammodels::sph::TimestepLog shammodels::sph::Solver::evolve_once() solver_graph.get_node_ref_base("attach fields to scheduler").evaluate(); } + + // check positivity + if (solver_config.dust_config.has_s_j_field()) { + logger::raw_ln(SourceLocation{}.format_one_line()); + PatchDataLayerLayout &pdl = scheduler().pdl_old(); + const u32 is_j = pdl.get_field_idx("s_j") ; + auto ndust = solver_config.dust_config.get_dust_nvar(); + auto &q = shamsys::instance::get_compute_scheduler().get_queue(); + scheduler().for_each_patchdata_nonempty([&](Patch cur_p, PatchDataLayer &pdat) { + sham::DeviceBuffer &buf_s_j = pdat.get_field_buf_ref(is_j); + + sham::kernel_call( + q, + sham::MultiRef{}, + sham::MultiRef{buf_s_j}, + pdat.get_obj_cnt() * ndust, + [](u32 tid, Tscal *s_j) { + Tscal s_j_tid = s_j[tid]; + if (s_j_tid < 0) { + throw std::runtime_error("s_j is negative"); + } + }); + }); + q.q.wait_and_throw(); + } + sph_prestep(t_current, dt); + + + // check positivity + if (solver_config.dust_config.has_s_j_field()) { + logger::raw_ln(SourceLocation{}.format_one_line()); + PatchDataLayerLayout &pdl = scheduler().pdl_old(); + const u32 is_j = pdl.get_field_idx("s_j") ; + auto ndust = solver_config.dust_config.get_dust_nvar(); + auto &q = shamsys::instance::get_compute_scheduler().get_queue(); + scheduler().for_each_patchdata_nonempty([&](Patch cur_p, PatchDataLayer &pdat) { + sham::DeviceBuffer &buf_s_j = pdat.get_field_buf_ref(is_j); + + sham::kernel_call( + q, + sham::MultiRef{}, + sham::MultiRef{buf_s_j}, + pdat.get_obj_cnt() * ndust, + [](u32 tid, Tscal *s_j) { + Tscal s_j_tid = s_j[tid]; + if (s_j_tid < 0) { + throw std::runtime_error("s_j is negative"); + } + }); + }); + q.q.wait_and_throw(); + } + + + // check positivity + if (solver_config.dust_config.has_s_j_field()) { + PatchDataLayerLayout &pdl = scheduler().pdl_old(); + const u32 is_j = pdl.get_field_idx("s_j") ; + auto ndust = solver_config.dust_config.get_dust_nvar(); + auto &q = shamsys::instance::get_compute_scheduler().get_queue(); + scheduler().for_each_patchdata_nonempty([&](Patch cur_p, PatchDataLayer &pdat) { + sham::DeviceBuffer &buf_s_j = pdat.get_field_buf_ref(is_j); + + sham::kernel_call( + q, + sham::MultiRef{}, + sham::MultiRef{buf_s_j}, + pdat.get_obj_cnt() * ndust, + [](u32 tid, Tscal *s_j) { + Tscal s_j_tid = s_j[tid]; + if (s_j_tid < 0) { + throw std::runtime_error("s_j is negative"); + } + }); + }); + } + using RTree = shamtree::CompressedLeafBVH; // Here we will add self grav to the external forces indep of vel (this will be moved into a @@ -2325,7 +2403,7 @@ shammodels::sph::TimestepLog shammodels::sph::Solver::evolve_once() if (solver_config.dust_config.has_s_j_field()) { ComputeField s_j_s_j_sq = utility.make_compute_field( "s_j s_j^2", solver_config.dust_config.get_dust_nvar()); - utility.fields_leapfrog_corrector( + utility.fields_leapfrog_corrector_positive_only( is_j, ids_j_dt, storage.old_ds_j_dt.get(), s_j_s_j_sq, dt / 2); } @@ -2903,20 +2981,26 @@ shammodels::sph::TimestepLog shammodels::sph::Solver::evolve_once() reset_merge_ghosts_fields(); reset_eos_fields(); - // restore positivity + + + // check positivity if (solver_config.dust_config.has_s_j_field()) { + PatchDataLayerLayout &pdl = scheduler().pdl_old(); + const u32 is_j = pdl.get_field_idx("s_j") ; auto ndust = solver_config.dust_config.get_dust_nvar(); - auto &q = shamsys::instance::get_compute_scheduler().get_queue(); + auto &q = shamsys::instance::get_compute_scheduler().get_queue(); scheduler().for_each_patchdata_nonempty([&](Patch cur_p, PatchDataLayer &pdat) { - sham::DeviceBuffer &buf_s_j = pdat.get_field_buf_ref(is_j); + sham::DeviceBuffer &buf_s_j = pdat.get_field_buf_ref(is_j); + sham::kernel_call( q, sham::MultiRef{}, sham::MultiRef{buf_s_j}, - pdat.get_obj_cnt()*ndust, + pdat.get_obj_cnt() * ndust, [](u32 tid, Tscal *s_j) { - if (s_j[tid] < 0) { - s_j[tid] = 0; + Tscal s_j_tid = s_j[tid]; + if (s_j_tid < 0) { + throw std::runtime_error("s_j is negative"); } }); }); diff --git a/src/shammodels/sph/src/modules/NodeUpdateDerivsMonofluidTVI.cpp b/src/shammodels/sph/src/modules/NodeUpdateDerivsMonofluidTVI.cpp index 5c5ed7fd7d..53379af2b7 100644 --- a/src/shammodels/sph/src/modules/NodeUpdateDerivsMonofluidTVI.cpp +++ b/src/shammodels/sph/src/modules/NodeUpdateDerivsMonofluidTVI.cpp @@ -98,7 +98,7 @@ struct KernelUpdateDerivsMonofluidTVI { Tscal delta_P = P_a - P_b; Tscal Ts_weighted = (Ttilde_sj_a / rho_a + Ttilde_sj_b / rho_b); - if (id_a == 2320 && jdust == 4 && false) { + if (id_a == 29951 && jdust == 4 && false) { printf("--------------------------------\n"); printf("id_b: %d, s_j_b: %e, Ttilde_sj_b: %e\n", id_b, s_j_b, Ttilde_sj_b); printf("rho_b: %e, P_b: %e, omega_b: %f\n", rho_b, P_b, omega_b); @@ -123,7 +123,7 @@ struct KernelUpdateDerivsMonofluidTVI { // restore it slowly to 0 ds_j_dt_a += (s_j_a < 0) ? -s_j_a / (10 * Ttilde_sj_a) : 0; - if (id_a == 2320 && jdust == 4 && false) { + if (id_a == 29951 && jdust == 4 && false) { printf("s_j_a: %e\n", s_j_a); printf("ds_j_dt_a: %e\n", ds_j_dt_a); printf("h_a: %e\n", h_a); @@ -135,20 +135,20 @@ struct KernelUpdateDerivsMonofluidTVI { printf("xyz_a: %f, %f, %f\n", xyz_a[0], xyz_a[1], xyz_a[2]); } - if (ds_j_dt_a > 1e-1 && false) { + if (ds_j_dt_a > 1e-3 && s_j_a < 1e-6 ) { printf("id_a: %d, s_j_a: %e, ds_j_dt_a: %e\n", id_a, s_j_a, ds_j_dt_a); - printf( - "h_a: %e, rho_a: %e, omega_a: %f, Ttilde_sj_a: %f\n", - h_a, - rho_a, - omega_a, - Ttilde_sj_a); - printf( - "term1: %e, term2: %e\n", - Tscal{-0.5} * term1, - (s_j_a / (2 * rho_a * omega_a)) * term2); - printf("xyz_a: %f, %f, %f\n", xyz_a[0], xyz_a[1], xyz_a[2]); - throw std::runtime_error("ds_j_dt_a is too large"); + //printf( + // "h_a: %e, rho_a: %e, omega_a: %f, Ttilde_sj_a: %f\n", + // h_a, + // rho_a, + // omega_a, + // Ttilde_sj_a); + //printf( + // "term1: %e, term2: %e\n", + // Tscal{-0.5} * term1, + // (s_j_a / (2 * rho_a * omega_a)) * term2); + //printf("xyz_a: %f, %f, %f\n", xyz_a[0], xyz_a[1], xyz_a[2]); + //throw std::runtime_error("ds_j_dt_a is too large"); } ds_j_dt[thread_id] = ds_j_dt_a; diff --git a/src/shamrock/include/shamrock/math/integrators.hpp b/src/shamrock/include/shamrock/math/integrators.hpp index 19a0562f05..6888571014 100644 --- a/src/shamrock/include/shamrock/math/integrators.hpp +++ b/src/shamrock/include/shamrock/math/integrators.hpp @@ -65,6 +65,16 @@ namespace shamrock { sham::DeviceBuffer &buf_eps_sq, sycl::range<1> elem_range, flt hdt); + + template + void leapfrog_corrector_positive_only( + sham::DeviceQueue &queue, + sham::DeviceBuffer &buf_val, + sham::DeviceBuffer &buf_der, + sham::DeviceBuffer &buf_der_old, + sham::DeviceBuffer &buf_eps_sq, + sycl::range<1> elem_range, + flt hdt); } // namespace integrators namespace utilities { diff --git a/src/shamrock/include/shamrock/scheduler/SchedulerUtility.hpp b/src/shamrock/include/shamrock/scheduler/SchedulerUtility.hpp index 873e288ca9..4807327222 100644 --- a/src/shamrock/include/shamrock/scheduler/SchedulerUtility.hpp +++ b/src/shamrock/include/shamrock/scheduler/SchedulerUtility.hpp @@ -61,6 +61,22 @@ namespace shamrock { }); } + template + inline void fields_leapfrog_corrector_positive_only( + u32 field_idx, u32 derfield_idx, u32 derfield_old_idx, flt hdt) { + StackEntry stack_loc{}; + using namespace shamrock::patch; + sched.for_each_patchdata_nonempty([&](Patch cur_p, PatchDataLayer &pdat) { + integrators::leapfrog_corrector_positive_only( + shamsys::instance::get_compute_scheduler().get_queue(), + pdat.get_field(field_idx).get_buf(), + pdat.get_field(derfield_idx).get_buf(), + pdat.get_field(derfield_old_idx).get_buf(), + pdat.get_obj_cnt(), + hdt); + }); + } + template inline void fields_leapfrog_corrector( u32 field_idx, @@ -82,6 +98,27 @@ namespace shamrock { }); } + template + inline void fields_leapfrog_corrector_positive_only( + u32 field_idx, + u32 derfield_idx, + ComputeField &derfield_old, + ComputeField &field_epsilon, + flt hdt) { + StackEntry stack_loc{}; + using namespace shamrock::patch; + sched.for_each_patchdata_nonempty([&](Patch cur_p, PatchDataLayer &pdat) { + integrators::leapfrog_corrector_positive_only( + shamsys::instance::get_compute_scheduler().get_queue(), + pdat.get_field(field_idx).get_buf(), + pdat.get_field(derfield_idx).get_buf(), + derfield_old.get_field(cur_p.id_patch).get_buf(), + field_epsilon.get_field(cur_p.id_patch).get_buf(), + pdat.get_obj_cnt() * derfield_old.get_nvar(), + hdt); + }); + } + template inline void fields_apply_periodicity(u32 field_idx, std::pair box) { StackEntry stack_loc{}; diff --git a/src/shamrock/src/math/integrators.cpp b/src/shamrock/src/math/integrators.cpp index f5630484b6..baef01fd01 100644 --- a/src/shamrock/src/math/integrators.cpp +++ b/src/shamrock/src/math/integrators.cpp @@ -18,6 +18,7 @@ #include "shambackends/sycl_utils.hpp" #include "shamrock/math/integrators.hpp" #include +#include "shambackends/math.hpp" namespace integ = shamrock::integrators; namespace util = shamrock::utilities; @@ -117,6 +118,40 @@ void integ::leapfrog_corrector( buf_eps_sq.complete_event_state(e); } +template +void integ::leapfrog_corrector_positive_only( + sham::DeviceQueue &queue, + sham::DeviceBuffer &buf_val, + sham::DeviceBuffer &buf_der, + sham::DeviceBuffer &buf_der_old, + sham::DeviceBuffer &buf_eps_sq, + sycl::range<1> elem_range, + flt hdt) { + + sham::EventList depends_list; + + auto acc_u = buf_val.get_write_access(depends_list); + auto acc_du = buf_der.get_read_access(depends_list); + auto acc_du_old = buf_der_old.get_read_access(depends_list); + auto acc_epsilon_sq = buf_eps_sq.get_write_access(depends_list); + + auto e = queue.submit(depends_list, [&](sycl::handler &cgh) { + cgh.parallel_for(elem_range, [=](sycl::item<1> item) { + u32 gid = (u32) item.get_id(); + + T incr = (hdt) * (acc_du[item] - acc_du_old[item]); + + acc_u[item] = sham::max(acc_u[item] + incr, 0.0); + acc_epsilon_sq[item] = sycl::dot(incr, incr); + }); + }); + + buf_val.complete_event_state(e); + buf_der.complete_event_state(e); + buf_der_old.complete_event_state(e); + buf_eps_sq.complete_event_state(e); +} + #ifndef DOXYGEN template void integ::leapfrog_corrector( sham::DeviceQueue &queue, @@ -153,6 +188,16 @@ template void integ::leapfrog_corrector( sham::DeviceBuffer &buf_eps_sq, sycl::range<1> elem_range, f64 hdt); + +template void integ::leapfrog_corrector_positive_only( + sham::DeviceQueue &queue, + sham::DeviceBuffer &buf_val, + sham::DeviceBuffer &buf_der, + sham::DeviceBuffer &buf_der_old, + sham::DeviceBuffer &buf_eps_sq, + sycl::range<1> elem_range, + f64 hdt); + #endif //////////////////////////////////////////////////////////////////////////////////////////////////// //////////////////////////////////////////////////////////////////////////////////////////////////// From 512c7f7df018923bdfa6b7da856a631446db6078 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Timoth=C3=A9e=20David--Cl=C3=A9ris?= Date: Fri, 10 Jul 2026 20:17:43 +0200 Subject: [PATCH 18/24] fixes --- examples/physics/dustywave_sympy.py | 212 ++++++------ examples/sph/run_dustysettle_tvi.py | 24 +- .../common/modules/ForwardEulerPositive.hpp | 6 +- .../include/shammodels/sph/SolverConfig.hpp | 24 +- .../sph/modules/ComputeCFLDustDeltav.hpp | 86 +++++ src/shammodels/sph/src/Solver.cpp | 305 +++++++++--------- src/shammodels/sph/src/SolverConfig.cpp | 1 + .../modules/NodeUpdateDerivsMonofluidTVI.cpp | 40 --- src/shammodels/sph/src/modules/io/VTKDump.cpp | 108 +++++++ src/shammodels/sph/src/pySPHModel.cpp | 12 +- src/shamrock/src/math/integrators.cpp | 2 +- 11 files changed, 510 insertions(+), 310 deletions(-) create mode 100644 src/shammodels/sph/include/shammodels/sph/modules/ComputeCFLDustDeltav.hpp diff --git a/examples/physics/dustywave_sympy.py b/examples/physics/dustywave_sympy.py index 0da1947762..776994b972 100644 --- a/examples/physics/dustywave_sympy.py +++ b/examples/physics/dustywave_sympy.py @@ -5,34 +5,36 @@ This example shows how to derive the dustywave TVA dispersion relation using SymPy """ -import sympy as sp -import numpy as np import matplotlib.pyplot as plt +import numpy as np +import sympy as sp # %% # Usefull symbols -omega = sp.symbols(r'\omega', complex=True) +omega = sp.symbols(r"\omega", complex=True) k, cs, ts, eps = sp.symbols( - r'k c_s t_s \epsilon', + r"k c_s t_s \epsilon", positive=True, real=True, ) -i = sp.I # imaginary unit +i = sp.I # imaginary unit -a = k*cs # s^-1 -b = a**2 * ts*eps # s^-1 +a = k * cs # s^-1 +b = a**2 * ts * eps # s^-1 # %% # The perturbation matrix is M = K(omega = 0) # The eigenfrequencies are the roots of det(M + i omega I) = 0 -K = sp.Matrix([ - [i*omega, 0, -i*a], - [b*(1-eps), i*omega-b, 0], - [-i*a*(1-eps), i*a, i*omega] -]) # s^-1 +K = sp.Matrix( + [ + [i * omega, 0, -i * a], + [b * (1 - eps), i * omega - b, 0], + [-i * a * (1 - eps), i * a, i * omega], + ] +) # s^-1 # %% @@ -42,36 +44,29 @@ # %% # Let's remove the leading omega mode (omega = 0) -det /= (i * omega) +det /= i * omega det = sp.simplify(det) det = sp.collect(det, omega) print(sp.latex(det)) # %% # Find the roots of the dispersion relation -r1,r2 = sp.solve(sp.Eq(det, 0), omega) +r1, r2 = sp.solve(sp.Eq(det, 0), omega) print(sp.latex(r1)) print(sp.latex(r2)) print(r1) print(r2) + # %% # Function to plot the roots -def get_roots(_r1,_r2,k_list, eps_value, cs_value, ts_value): +def get_roots(_r1, _r2, k_list, eps_value, cs_value, ts_value): # Substitute all parameters - r1_num = _r1.subs({ - cs: cs_value, - ts: ts_value, - eps: eps_value - }) - - r2_num = _r2.subs({ - cs: cs_value, - ts: ts_value, - eps: eps_value - }) + r1_num = _r1.subs({cs: cs_value, ts: ts_value, eps: eps_value}) + + r2_num = _r2.subs({cs: cs_value, ts: ts_value, eps: eps_value}) r1_num_re = sp.re(r1_num) r1_num_im = sp.im(r1_num) @@ -100,83 +95,97 @@ def restore(lst): r1_vals_im = restore(r1_vals_im) r2_vals_re = restore(r2_vals_re) r2_vals_im = restore(r2_vals_im) - + return r1_vals_re, r1_vals_im, r2_vals_re, r2_vals_im + def get_roots_LP14(k_list, eps_value, cs_value, ts_value): - _r1 = + cs * sp.sqrt(1-eps) * k - i * ts* k **2 * cs **2 * eps /2 - _r2 = - cs * sp.sqrt(1-eps) * k - i * ts* k **2 * cs **2 * eps /2 + _r1 = +cs * sp.sqrt(1 - eps) * k - i * ts * k**2 * cs**2 * eps / 2 + _r2 = -cs * sp.sqrt(1 - eps) * k - i * ts * k**2 * cs**2 * eps / 2 return get_roots(_r1, _r2, k_list, eps_value, cs_value, ts_value) - def get_overroots_DCL26_simple(k_list, eps_value, cs_value, ts_value): - _r1 = + cs * sp.sqrt(1-eps) * k + i*k**2*cs**2 * eps*ts * (-1 + 1)/2 - _r2 = - cs * sp.sqrt(1-eps) * k + i*k**2*cs**2 * eps*ts *(-1 - 1)/2 + _r1 = +cs * sp.sqrt(1 - eps) * k + i * k**2 * cs**2 * eps * ts * (-1 + 1) / 2 + _r2 = -cs * sp.sqrt(1 - eps) * k + i * k**2 * cs**2 * eps * ts * (-1 - 1) / 2 return get_roots(_r1, _r2, k_list, eps_value, cs_value, ts_value) + def get_overroots_DCL26(k_list, eps_value, cs_value, ts_value): - D = 4*(1-eps) - eps**2 * cs**2 * ts**2 * k**2 + D = 4 * (1 - eps) - eps**2 * cs**2 * ts**2 * k**2 sqrtD_real = sp.sqrt(sp.Max(D, 0)) sqrtD_imag = sp.sqrt(sp.Max(-D, 0)) - _r1 = cs*k/2 * ( - +sqrtD_real - + i*(sqrtD_imag - eps*cs*k*ts) - ) + _r1 = cs * k / 2 * (+sqrtD_real + i * (sqrtD_imag - eps * cs * k * ts)) - _r2 = cs*k/2 * ( - -sqrtD_real - + i*(-sqrtD_imag - eps*cs*k*ts) - ) + _r2 = cs * k / 2 * (-sqrtD_real + i * (-sqrtD_imag - eps * cs * k * ts)) print(sp.latex(sp.Abs(_r1))) print(sp.latex(sp.Abs(_r2))) return get_roots(_r1, _r2, k_list, eps_value, cs_value, ts_value) -def plot_case(k_plot, eps_value, cs_value, ts_value): - r1_vals_re, r1_vals_im, r2_vals_re, r2_vals_im = get_roots(r1,r2,k_plot, eps_value, cs_value, ts_value) +def plot_case(k_plot, eps_value, cs_value, ts_value): - r1_vals_re_LP14, r1_vals_im_LP14, r2_vals_re_LP14, r2_vals_im_LP14 = get_roots_LP14(k_plot, eps_value, cs_value, ts_value) + r1_vals_re, r1_vals_im, r2_vals_re, r2_vals_im = get_roots( + r1, r2, k_plot, eps_value, cs_value, ts_value + ) - r1_vals_re_DCL26_simple, r1_vals_im_DCL26_simple, r2_vals_re_DCL26_simple, r2_vals_im_DCL26_simple = get_overroots_DCL26_simple(k_plot, eps_value, cs_value, ts_value) - r1_vals_re_DCL26, r1_vals_im_DCL26, r2_vals_re_DCL26, r2_vals_im_DCL26 = get_overroots_DCL26(k_plot, eps_value, cs_value, ts_value) + r1_vals_re_LP14, r1_vals_im_LP14, r2_vals_re_LP14, r2_vals_im_LP14 = get_roots_LP14( + k_plot, eps_value, cs_value, ts_value + ) + ( + r1_vals_re_DCL26_simple, + r1_vals_im_DCL26_simple, + r2_vals_re_DCL26_simple, + r2_vals_im_DCL26_simple, + ) = get_overroots_DCL26_simple(k_plot, eps_value, cs_value, ts_value) + r1_vals_re_DCL26, r1_vals_im_DCL26, r2_vals_re_DCL26, r2_vals_im_DCL26 = get_overroots_DCL26( + k_plot, eps_value, cs_value, ts_value + ) # Create figure fig, axs = plt.subplots(2, 2, figsize=(8, 8), sharex=True) # Real parts - axs[0,0].plot(k_plot, r1_vals_re,color="0",linewidth=2, label="Re($\omega_+$)") - axs[0,0].plot(k_plot, r2_vals_re,color="0",linewidth=2, label="Re($\omega_-$)") - axs[0,0].plot(k_plot, r1_vals_re_LP14,"--", label="Re($\omega_{+,LP14}$)") - axs[0,0].plot(k_plot, r2_vals_re_LP14,"--", label="Re($\omega_{-,LP14}$)") - axs[0,0].plot(k_plot, r1_vals_re_DCL26_simple,linestyle="dotted", label="Re($\omega_{+,approx}$)") - axs[0,0].plot(k_plot, r2_vals_re_DCL26_simple,linestyle="dotted", label="Re($\omega_{-,approx}$)") - #axs[0,0].plot(k_plot, r1_vals_re_DCL26,"--", label="Re($r_1$) DCL26") - #axs[0,0].plot(k_plot, r2_vals_re_DCL26,"--", label="Re($r_2$) DCL26") - axs[0,0].set_ylabel("Real part") - axs[0,0].grid(True) - axs[0,0].legend() + axs[0, 0].plot(k_plot, r1_vals_re, color="0", linewidth=2, label="Re($\omega_+$)") + axs[0, 0].plot(k_plot, r2_vals_re, color="0", linewidth=2, label="Re($\omega_-$)") + axs[0, 0].plot(k_plot, r1_vals_re_LP14, "--", label="Re($\omega_{+,LP14}$)") + axs[0, 0].plot(k_plot, r2_vals_re_LP14, "--", label="Re($\omega_{-,LP14}$)") + axs[0, 0].plot( + k_plot, r1_vals_re_DCL26_simple, linestyle="dotted", label="Re($\omega_{+,approx}$)" + ) + axs[0, 0].plot( + k_plot, r2_vals_re_DCL26_simple, linestyle="dotted", label="Re($\omega_{-,approx}$)" + ) + # axs[0,0].plot(k_plot, r1_vals_re_DCL26,"--", label="Re($r_1$) DCL26") + # axs[0,0].plot(k_plot, r2_vals_re_DCL26,"--", label="Re($r_2$) DCL26") + axs[0, 0].set_ylabel("Real part") + axs[0, 0].grid(True) + axs[0, 0].legend() # Imaginary parts - axs[0,1].plot(k_plot, r1_vals_im,color="0",linewidth=2, label="Im($\omega_+$)") - axs[0,1].plot(k_plot, r2_vals_im,color="0",linewidth=2, label="Im($\omega_-$)") - axs[0,1].plot(k_plot, r1_vals_im_LP14,"--", label="Im($\omega_{+,LP14}$)") - axs[0,1].plot(k_plot, r2_vals_im_LP14,"--", label="Im($\omega_{-,LP14}$)") - axs[0,1].plot(k_plot, r1_vals_im_DCL26_simple,linestyle="dotted", label="Im($\omega_{+,approx}$)") - axs[0,1].plot(k_plot, r2_vals_im_DCL26_simple,linestyle="dotted", label="Im($\omega_{-,approx}$)") - #axs[0,1].plot(k_plot, r1_vals_im_DCL26,"--", label="Im($r_1$) DCL26") - #axs[0,1].plot(k_plot, r2_vals_im_DCL26,"--", label="Im($r_2$) DCL26") - axs[0,1].set_xlabel("$k$") - axs[0,1].set_ylabel("Imaginary part") - axs[0,1].grid(True) - axs[0,1].legend() + axs[0, 1].plot(k_plot, r1_vals_im, color="0", linewidth=2, label="Im($\omega_+$)") + axs[0, 1].plot(k_plot, r2_vals_im, color="0", linewidth=2, label="Im($\omega_-$)") + axs[0, 1].plot(k_plot, r1_vals_im_LP14, "--", label="Im($\omega_{+,LP14}$)") + axs[0, 1].plot(k_plot, r2_vals_im_LP14, "--", label="Im($\omega_{-,LP14}$)") + axs[0, 1].plot( + k_plot, r1_vals_im_DCL26_simple, linestyle="dotted", label="Im($\omega_{+,approx}$)" + ) + axs[0, 1].plot( + k_plot, r2_vals_im_DCL26_simple, linestyle="dotted", label="Im($\omega_{-,approx}$)" + ) + # axs[0,1].plot(k_plot, r1_vals_im_DCL26,"--", label="Im($r_1$) DCL26") + # axs[0,1].plot(k_plot, r2_vals_im_DCL26,"--", label="Im($r_2$) DCL26") + axs[0, 1].set_xlabel("$k$") + axs[0, 1].set_ylabel("Imaginary part") + axs[0, 1].grid(True) + axs[0, 1].legend() # Abs r1_vals_abs = np.sqrt(r1_vals_re**2 + r1_vals_im**2) @@ -187,38 +196,47 @@ def plot_case(k_plot, eps_value, cs_value, ts_value): r2_vals_abs_DCL26_simple = np.sqrt(r2_vals_re_DCL26_simple**2 + r2_vals_im_DCL26_simple**2) r1_vals_abs_DCL26 = np.sqrt(r1_vals_re_DCL26**2 + r1_vals_im_DCL26**2) r2_vals_abs_DCL26 = np.sqrt(r2_vals_re_DCL26**2 + r2_vals_im_DCL26**2) - axs[1,0].plot(k_plot, r1_vals_abs,color="0",linewidth=2, label="Abs($\omega_+$)") - axs[1,0].plot(k_plot, r2_vals_abs,color="0",linewidth=2, label="Abs($\omega_-$)") - axs[1,0].plot(k_plot, r1_vals_abs_LP14,"--", label="Abs($\omega_{+,LP14}$)") - axs[1,0].plot(k_plot, r2_vals_abs_LP14,"--", label="Abs($\omega_{-,LP14}$)") - axs[1,0].plot(k_plot, r1_vals_abs_DCL26_simple,linestyle="dotted", label="Abs($\omega_{+,approx}$)") - axs[1,0].plot(k_plot, r2_vals_abs_DCL26_simple,linestyle="dotted", label="Abs($\omega_{-,approx}$)") - - def approx(_k,_cs,_ts,_eps): + axs[1, 0].plot(k_plot, r1_vals_abs, color="0", linewidth=2, label="Abs($\omega_+$)") + axs[1, 0].plot(k_plot, r2_vals_abs, color="0", linewidth=2, label="Abs($\omega_-$)") + axs[1, 0].plot(k_plot, r1_vals_abs_LP14, "--", label="Abs($\omega_{+,LP14}$)") + axs[1, 0].plot(k_plot, r2_vals_abs_LP14, "--", label="Abs($\omega_{-,LP14}$)") + axs[1, 0].plot( + k_plot, r1_vals_abs_DCL26_simple, linestyle="dotted", label="Abs($\omega_{+,approx}$)" + ) + axs[1, 0].plot( + k_plot, r2_vals_abs_DCL26_simple, linestyle="dotted", label="Abs($\omega_{-,approx}$)" + ) + + def approx(_k, _cs, _ts, _eps): print(type(_k), type(_cs), type(_ts), type(_eps)) - return _cs * _k * np.sqrt((1-_eps) + (_k *_cs *_ts* _eps )**2) + return _cs * _k * np.sqrt((1 - _eps) + (_k * _cs * _ts * _eps) ** 2) - axs[1,0].plot(k_plot, approx(k_plot, cs_value, ts_value, eps_value), "--", label=r"$max(\vert \omega_{\pm,approx} \vert)$") + axs[1, 0].plot( + k_plot, + approx(k_plot, cs_value, ts_value, eps_value), + "--", + label=r"$max(\vert \omega_{\pm,approx} \vert)$", + ) - #axs[1,0].plot(k_plot, r1_vals_abs_DCL26,"--", label="Abs($r_1$) DCL26") - #axs[1,0].plot(k_plot, r2_vals_abs_DCL26,"--", label="Abs($r_2$) DCL26") - axs[1,0].set_xlabel("$k$") - axs[1,0].set_ylabel("Abs part") - axs[1,0].grid(True) - axs[1,0].legend() + # axs[1,0].plot(k_plot, r1_vals_abs_DCL26,"--", label="Abs($r_1$) DCL26") + # axs[1,0].plot(k_plot, r2_vals_abs_DCL26,"--", label="Abs($r_2$) DCL26") + axs[1, 0].set_xlabel("$k$") + axs[1, 0].set_ylabel("Abs part") + axs[1, 0].grid(True) + axs[1, 0].legend() # delta with max - r_max = np.maximum(r1_vals_abs,r2_vals_abs) - r_max_LP14 = np.maximum(r1_vals_abs_LP14,r2_vals_abs_LP14) - r_max_DCL26_simple = np.maximum(r1_vals_abs_DCL26_simple,r2_vals_abs_DCL26_simple) - r_max_DCL26 = np.maximum(r1_vals_abs_DCL26,r2_vals_abs_DCL26) - axs[1,1].plot(k_plot, (r_max_LP14 - r_max) / r_max, label="Ana - LP14") - axs[1,1].plot(k_plot, (r_max_DCL26_simple - r_max) / r_max, label="Ana - DCL26 simple") - #axs[1,1].plot(k_plot, (r_max_DCL26 - r_max) / r_max, label="Ana - DCL26") - axs[1,1].set_xlabel("$k$") - axs[1,1].set_ylabel("Abs(Ana) - Abs(Max model) / Abs(Ana)") - axs[1,1].grid(True) - axs[1,1].legend() + r_max = np.maximum(r1_vals_abs, r2_vals_abs) + r_max_LP14 = np.maximum(r1_vals_abs_LP14, r2_vals_abs_LP14) + r_max_DCL26_simple = np.maximum(r1_vals_abs_DCL26_simple, r2_vals_abs_DCL26_simple) + r_max_DCL26 = np.maximum(r1_vals_abs_DCL26, r2_vals_abs_DCL26) + axs[1, 1].plot(k_plot, (r_max_LP14 - r_max) / r_max, label="Ana - LP14") + axs[1, 1].plot(k_plot, (r_max_DCL26_simple - r_max) / r_max, label="Ana - DCL26 simple") + # axs[1,1].plot(k_plot, (r_max_DCL26 - r_max) / r_max, label="Ana - DCL26") + axs[1, 1].set_xlabel("$k$") + axs[1, 1].set_ylabel("Abs(Ana) - Abs(Max model) / Abs(Ana)") + axs[1, 1].grid(True) + axs[1, 1].legend() plt.suptitle(f"eps = {eps_value}, cs = {cs_value}, ts = {ts_value}") @@ -236,4 +254,4 @@ def approx(_k,_cs,_ts,_eps): # Plot the case eps = 0.1, cs = 1.0, ts = 1.0 k_plot = np.linspace(0, 40, 1000) plot_case(k_plot, 0.1, 1.0, 1.0) -plt.show() \ No newline at end of file +plt.show() diff --git a/examples/sph/run_dustysettle_tvi.py b/examples/sph/run_dustysettle_tvi.py index 7385b833cf..0663d36072 100644 --- a/examples/sph/run_dustysettle_tvi.py +++ b/examples/sph/run_dustysettle_tvi.py @@ -59,7 +59,7 @@ epsilon_base = 0.01 -tlist = [0.01 * i for i in range(50)] +tlist = [0.0025 * i for i in range(500)] tinject = 0 iinject = 0 t_end = 5.0 @@ -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_tvi(nvar=ndust, C_delta_v=10.0) 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) @@ -732,7 +732,8 @@ def analyse_and_plot(j): axs[0].set_ylabel(r"$s_j$") axs[0].set_xlabel(r"$z$") axs[0].set_xlim(-box, box) - axs[0].set_yscale("symlog", linthresh=1e-10) + axs[0].set_yscale("log") + axs[0].set_ylim(1e-20, 1e-1) axs[1].set_ylabel(r"$\dot{s}_j$") axs[1].set_xlabel(r"$z$") @@ -767,6 +768,8 @@ def dust_mass_analysis(): idust_analysis += 1 +ivtk = 0 + tnext = 0 for j in range(1000): if j == iinject: @@ -775,8 +778,19 @@ def dust_mass_analysis(): if tlist[j] >= t_start: if j > 0: - model.evolve_until(tlist[j]) - # model.timestep() + while 1: + res = model.evolve_until(tlist[j], niter_max=10000) + # model.timestep() + + if model.get_time() >= 0.0: + model.do_vtk_dump(f"{dump_folder}/vtk_{ivtk:04d}.vtk", True) + ivtk += 1 + + # if ivtk > 100: + # exit() + + if res.reach_target_time: + break if dust_injected: dust_mass_analysis() diff --git a/src/shammodels/common/include/shammodels/common/modules/ForwardEulerPositive.hpp b/src/shammodels/common/include/shammodels/common/modules/ForwardEulerPositive.hpp index 874fff1b80..ac7ab042d1 100644 --- a/src/shammodels/common/include/shammodels/common/modules/ForwardEulerPositive.hpp +++ b/src/shammodels/common/include/shammodels/common/modules/ForwardEulerPositive.hpp @@ -10,7 +10,7 @@ #pragma once /** - * @file ForwardEuler.hpp + * @file ForwardEulerPositive.hpp * @author Timothée David--Cléris (tim.shamrock@proton.me) * @brief Implements a forward Euler integration step as a solver graph node. * @@ -65,7 +65,7 @@ namespace shammodels::common::modules { sham::DDMultiRef{edges.field.get_spans()}, edges.sizes.indexes, [dt](u32 gid, const T *time_derivative, T *field) { - auto tmp = field[gid] + dt * time_derivative[gid]; + auto tmp = field[gid] + dt * time_derivative[gid]; field[gid] = std::max(tmp, 0.0); }); @@ -81,7 +81,7 @@ namespace shammodels::common::modules { sham::DDMultiRef{edges.field.get_spans()}, var_count, [dt](u32 gid, const T *time_derivative, T *field) { - auto tmp = field[gid] + dt * time_derivative[gid]; + auto tmp = field[gid] + dt * time_derivative[gid]; field[gid] = std::max(tmp, 0.0); }); } diff --git a/src/shammodels/sph/include/shammodels/sph/SolverConfig.hpp b/src/shammodels/sph/include/shammodels/sph/SolverConfig.hpp index 4738dae5ff..fbdf71356b 100644 --- a/src/shammodels/sph/include/shammodels/sph/SolverConfig.hpp +++ b/src/shammodels/sph/include/shammodels/sph/SolverConfig.hpp @@ -118,6 +118,9 @@ namespace shammodels::sph { struct MonofluidTVI { u32 ndust; bool pure_diffusion_mode = false; + + Tscal C_1_fluid = 0.1; + Tscal C_delta_v = 1; }; struct MonofluidComplete { @@ -130,8 +133,12 @@ namespace shammodels::sph { 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_tvi( + u32 nvar, + bool pure_diffusion_mode = false, + Tscal C_1_fluid = 0.1, + Tscal C_delta_v = 0.1) { + current_mode = MonofluidTVI{nvar, pure_diffusion_mode, C_1_fluid, C_delta_v}; } inline void set_monofluid_complete(u32 nvar) { current_mode = MonofluidComplete{nvar}; } @@ -141,6 +148,10 @@ namespace shammodels::sph { return bool(std::get_if(¤t_mode)); } + inline MonofluidTVI get_monofluid_tvi() { + return shambase::get_check_ref(std::get_if(¤t_mode)); + } + inline void mode_to_json(nlohmann::json &j) const { if (const None *cfg = std::get_if(¤t_mode)) { j = {{"type", "none"}}; @@ -148,7 +159,9 @@ namespace shammodels::sph { j = {{"type", "monofluid_tvi"}, {"ndust", cfg->ndust}, - {"pure_diffusion_mode", cfg->pure_diffusion_mode}}; + {"pure_diffusion_mode", cfg->pure_diffusion_mode}, + {"C_1_fluid", cfg->C_1_fluid}, + {"C_delta_v", cfg->C_delta_v}}; } else if ( const MonofluidComplete *cfg = std::get_if(¤t_mode)) { j = {{"type", "monofluid_complete"}, {"ndust", cfg->ndust}}; @@ -163,7 +176,10 @@ namespace shammodels::sph { set_none(); } else if (type == "monofluid_tvi") { set_monofluid_tvi( - j.at("ndust").get(), j.at("pure_diffusion_mode").get()); + j.at("ndust").get(), + j.at("pure_diffusion_mode").get(), + j.at("C_1_fluid").get(), + j.at("C_delta_v").get()); } else if (type == "monofluid_complete") { set_monofluid_complete(j.at("ndust").get()); } else { diff --git a/src/shammodels/sph/include/shammodels/sph/modules/ComputeCFLDustDeltav.hpp b/src/shammodels/sph/include/shammodels/sph/modules/ComputeCFLDustDeltav.hpp new file mode 100644 index 0000000000..82bc3369f8 --- /dev/null +++ b/src/shammodels/sph/include/shammodels/sph/modules/ComputeCFLDustDeltav.hpp @@ -0,0 +1,86 @@ +// -------------------------------------------------------// +// +// SHAMROCK code for hydrodynamics +// Copyright (c) 2021-2026 Timothée David--Cléris +// SPDX-License-Identifier: CeCILL Free Software License Agreement v2.1 +// Shamrock is licensed under the CeCILL 2.1 License, see LICENSE for more information +// +// -------------------------------------------------------// + +#pragma once + +/** + * @file ComputeCFLDustDeltav.hpp + * @author Timothée David--Cléris (tim.shamrock@proton.me) + * @brief + * + */ + +#include "shambackends/kernel_call_distrib.hpp" +#include "shammodels/sph/math/density.hpp" +#include "shamrock/solvergraph/IFieldSpan.hpp" +#include "shamrock/solvergraph/INode.hpp" +#include "shamrock/solvergraph/Indexes.hpp" +#include "shamrock/solvergraph/ScalarEdge.hpp" +#include "shamsys/NodeInstance.hpp" + +#define NODE_EDGES(X_RO, X_RW) \ + X_RO(shamrock::solvergraph::Indexes, part_counts) \ + X_RO(shamrock::solvergraph::ScalarEdge, C_delta_v) \ + X_RO(shamrock::solvergraph::IFieldSpan, hpart) \ + X_RO(shamrock::solvergraph::IFieldSpan, delta_v) \ + X_RW(shamrock::solvergraph::IFieldSpan, cfl_dt) + +template +class ComputeCFLDustDeltav : public shamrock::solvergraph::INode { + + using Tscal = shambase::VecComponent; + + u32 nbins; + + public: + ComputeCFLDustDeltav(u32 nbins) : nbins(nbins) {} + + EXPAND_NODE_EDGES(NODE_EDGES) + + inline void _impl_evaluate_internal() { + auto edges = get_edges(); + + auto dev_sched = shamsys::instance::get_compute_scheduler_ptr(); + + Tscal C_delta_v = edges.C_delta_v.value; + + sham::distributed_data_kernel_call( + dev_sched, + sham::DDMultiRef{edges.hpart.get_spans(), edges.delta_v.get_spans()}, + sham::DDMultiRef{edges.cfl_dt.get_spans()}, + edges.part_counts.indexes, + [C_delta_v, nbins = this->nbins]( + u32 id_a, + const Tscal *hpart, + const Tvec *delta_v, + Tscal *cfl_dt) { + u32 id_a_d = id_a * nbins; + + Tscal h_a = hpart[id_a]; + + Tscal cfl_tmp = std::numeric_limits::infinity(); + + for (int j = 0; j < nbins; j++) { + Tvec delta_v_j_a = delta_v[id_a_d + j]; + Tscal delta_v_j_a_norm = sycl::length(delta_v_j_a); + cfl_tmp = sycl::min(cfl_tmp, h_a / delta_v_j_a_norm); + } + + cfl_tmp *= C_delta_v; + + cfl_dt[id_a] = sycl::min(cfl_dt[id_a], cfl_tmp); + }); + } + + inline virtual std::string _impl_get_label() const { return "ComputeCFLDustDeltav"; }; + + inline virtual std::string _impl_get_tex() const { return "C_{\\Delta_v}"; }; +}; + +#undef NODE_EDGES diff --git a/src/shammodels/sph/src/Solver.cpp b/src/shammodels/sph/src/Solver.cpp index d33f00ef98..6a9fd797f3 100644 --- a/src/shammodels/sph/src/Solver.cpp +++ b/src/shammodels/sph/src/Solver.cpp @@ -49,6 +49,7 @@ #include "shammodels/sph/modules/ComputeCFLCourant.hpp" #include "shammodels/sph/modules/ComputeCFLDivBCleaning.hpp" #include "shammodels/sph/modules/ComputeCFLDust1Fluid.hpp" +#include "shammodels/sph/modules/ComputeCFLDustDeltav.hpp" #include "shammodels/sph/modules/ComputeCFLForce.hpp" #include "shammodels/sph/modules/ComputeEos.hpp" #include "shammodels/sph/modules/ComputeLoadBalanceValue.hpp" @@ -64,7 +65,9 @@ #include "shammodels/sph/modules/IterateSmoothingLengthDensityNeighLim.hpp" #include "shammodels/sph/modules/KillParticles.hpp" #include "shammodels/sph/modules/LoopSmoothingLengthIter.hpp" +#include "shammodels/sph/modules/MonoFluidTVIDeltav.hpp" #include "shammodels/sph/modules/NeighbourCache.hpp" +#include "shammodels/sph/modules/NodeComputePressureGrad.hpp" #include "shammodels/sph/modules/ParticleReordering.hpp" #include "shammodels/sph/modules/SetDustStoppingTimeConstant.hpp" #include "shammodels/sph/modules/SetDustStoppingTimeEpstein.hpp" @@ -454,7 +457,8 @@ void shammodels::sph::Solver::init_solver_graph() { if (has_s_j_field) { u32 ndust = solver_config.dust_config.get_dust_nvar(); auto half_step_s_j = solver_graph.register_node( - prefix + "_s_j", shammodels::common::modules::ForwardEulerPositive(ndust)); + prefix + "_s_j", + shammodels::common::modules::ForwardEulerPositive(ndust)); shambase::get_check_ref(half_step_s_j) .set_edges( solver_graph.get_edge_ptr>("dt_half"), @@ -1553,6 +1557,57 @@ void shammodels::sph::Solver::prepare_corrector() { } } +template +void map_field_refs( + PatchScheduler &sched, u32 field_idx, shamrock::solvergraph::FieldRefs &refs) { + + using namespace shamrock::solvergraph; + using namespace shamrock::patch; + + shamrock::solvergraph::DDPatchDataFieldRef field_refs = {}; + sched.for_each_patchdata_nonempty([&](const Patch p, PatchDataLayer &pdat) { + auto &field = pdat.get_field(field_idx); + field_refs.add_obj(p.id_patch, std::ref(field)); + }); + refs.set_refs(field_refs); +} + +template +void map_field_refs_ext( + PatchScheduler &sched, + shambase::DistributedData &mpdats, + u32 field_idx, + shamrock::solvergraph::FieldRefs &refs) { + + using namespace shamrock::solvergraph; + using namespace shamrock::patch; + + shamrock::solvergraph::DDPatchDataFieldRef field_refs = {}; + sched.for_each_patchdata_nonempty([&](const Patch p, PatchDataLayer &pdat) { + PatchDataLayer &mpdat = mpdats.get(p.id_patch); + auto &field = mpdat.get_field(field_idx); + field_refs.add_obj(p.id_patch, std::ref(field)); + }); + refs.set_refs(field_refs); +} + +template +void map_field_refs_ext( + PatchScheduler &sched, + shamrock::ComputeField &field_data, + shamrock::solvergraph::FieldRefs &refs) { + + using namespace shamrock::solvergraph; + using namespace shamrock::patch; + + shamrock::solvergraph::DDPatchDataFieldRef field_refs = {}; + sched.for_each_patchdata_nonempty([&](const Patch p, PatchDataLayer &pdat) { + auto &field = field_data.get_field(p.id_patch); + field_refs.add_obj(p.id_patch, std::ref(field)); + }); + refs.set_refs(field_refs); +} + template class Kern> void shammodels::sph::Solver::update_derivs(Tscal dt_hydro) { @@ -1570,9 +1625,13 @@ void shammodels::sph::Solver::update_derivs(Tscal dt_hydro) { shamrock::patch::PatchDataLayerLayout &ghost_layout = shambase::get_check_ref(storage.ghost_layout.get()); + shamrock::patch::PatchDataLayerLayout &pdl = scheduler().pdl_old(); + u32 ihpart_interf = ghost_layout.get_field_idx("hpart"); auto &part_counts_with_ghost = storage.part_counts_with_ghost; + auto &part_counts = storage.part_counts; + shambase::DistributedData &mpdats = storage.merged_patchdata_ghost.get(); @@ -1641,6 +1700,61 @@ void shammodels::sph::Solver::update_derivs(Tscal dt_hydro) { } node_set_tj->evaluate(); } + + // delta v computation (for CFL or other uses e.g. COALA) + auto &pressure_field = storage.pressure; + auto &xyz_refs = storage.positions_with_ghosts; + + std::shared_ptr> omega_refs + = std::make_shared>("omega", "omega"); + { + u32 iomega_interf = ghost_layout.get_field_idx("omega"); + shambase::get_check_ref(omega_refs) + .set_refs(mpdats.map>>( + [iomega_interf](u64 id, shamrock::patch::PatchDataLayer &mpdat) { + return std::ref(mpdat.get_field(iomega_interf)); + })); + } + + std::shared_ptr> s_j_refs + = std::make_shared>("s_j", "s_j"); + { + u32 is_j_interf = ghost_layout.get_field_idx("s_j"); + shambase::get_check_ref(s_j_refs).set_refs( + mpdats.map>>( + [is_j_interf](u64 id, shamrock::patch::PatchDataLayer &mpdat) { + return std::ref(mpdat.get_field(is_j_interf)); + })); + } + + std::shared_ptr> grad_P_on_rho + = std::make_shared>(1, "grad P/rho", "grad P/rho"); + + u32 idelta_v = pdl.get_field_idx("delta_v"); + + std::shared_ptr> delta_v + = std::make_shared>("Delta v", "Delta v"); + map_field_refs(scheduler(), idelta_v, *delta_v); + + auto press_grad_node = std::make_shared>(); + auto delta_v_node = std::make_shared>(ndust); + + press_grad_node->set_edges( + gpart_mass, + part_counts, + part_counts_with_ghost, + xyz_refs, + hpart_refs, + omega_refs, + pressure_field, + storage.neigh_cache, + grad_P_on_rho); + + delta_v_node->set_edges( + gpart_mass, part_counts, hpart_refs, grad_P_on_rho, s_j_refs, t_j_field, delta_v); + + press_grad_node->evaluate(); + delta_v_node->evaluate(); } modules::UpdateDerivs derivs(context, solver_config, storage); @@ -1662,57 +1776,6 @@ void shammodels::sph::Solver::update_sync_load_values() { scheduler().scheduler_step(false, false); } -template -void map_field_refs( - PatchScheduler &sched, u32 field_idx, shamrock::solvergraph::FieldRefs &refs) { - - using namespace shamrock::solvergraph; - using namespace shamrock::patch; - - shamrock::solvergraph::DDPatchDataFieldRef field_refs = {}; - sched.for_each_patchdata_nonempty([&](const Patch p, PatchDataLayer &pdat) { - auto &field = pdat.get_field(field_idx); - field_refs.add_obj(p.id_patch, std::ref(field)); - }); - refs.set_refs(field_refs); -} - -template -void map_field_refs_ext( - PatchScheduler &sched, - shambase::DistributedData &mpdats, - u32 field_idx, - shamrock::solvergraph::FieldRefs &refs) { - - using namespace shamrock::solvergraph; - using namespace shamrock::patch; - - shamrock::solvergraph::DDPatchDataFieldRef field_refs = {}; - sched.for_each_patchdata_nonempty([&](const Patch p, PatchDataLayer &pdat) { - PatchDataLayer &mpdat = mpdats.get(p.id_patch); - auto &field = mpdat.get_field(field_idx); - field_refs.add_obj(p.id_patch, std::ref(field)); - }); - refs.set_refs(field_refs); -} - -template -void map_field_refs_ext( - PatchScheduler &sched, - shamrock::ComputeField &field_data, - shamrock::solvergraph::FieldRefs &refs) { - - using namespace shamrock::solvergraph; - using namespace shamrock::patch; - - shamrock::solvergraph::DDPatchDataFieldRef field_refs = {}; - sched.for_each_patchdata_nonempty([&](const Patch p, PatchDataLayer &pdat) { - auto &field = field_data.get_field(p.id_patch); - field_refs.add_obj(p.id_patch, std::ref(field)); - }); - refs.set_refs(field_refs); -} - template class Kern> shammodels::sph::TimestepLog shammodels::sph::Solver::evolve_once() { @@ -1836,85 +1899,8 @@ shammodels::sph::TimestepLog shammodels::sph::Solver::evolve_once() solver_graph.get_node_ref_base("attach fields to scheduler").evaluate(); } - - // check positivity - if (solver_config.dust_config.has_s_j_field()) { - logger::raw_ln(SourceLocation{}.format_one_line()); - PatchDataLayerLayout &pdl = scheduler().pdl_old(); - const u32 is_j = pdl.get_field_idx("s_j") ; - auto ndust = solver_config.dust_config.get_dust_nvar(); - auto &q = shamsys::instance::get_compute_scheduler().get_queue(); - scheduler().for_each_patchdata_nonempty([&](Patch cur_p, PatchDataLayer &pdat) { - sham::DeviceBuffer &buf_s_j = pdat.get_field_buf_ref(is_j); - - sham::kernel_call( - q, - sham::MultiRef{}, - sham::MultiRef{buf_s_j}, - pdat.get_obj_cnt() * ndust, - [](u32 tid, Tscal *s_j) { - Tscal s_j_tid = s_j[tid]; - if (s_j_tid < 0) { - throw std::runtime_error("s_j is negative"); - } - }); - }); - q.q.wait_and_throw(); - } - sph_prestep(t_current, dt); - - - // check positivity - if (solver_config.dust_config.has_s_j_field()) { - logger::raw_ln(SourceLocation{}.format_one_line()); - PatchDataLayerLayout &pdl = scheduler().pdl_old(); - const u32 is_j = pdl.get_field_idx("s_j") ; - auto ndust = solver_config.dust_config.get_dust_nvar(); - auto &q = shamsys::instance::get_compute_scheduler().get_queue(); - scheduler().for_each_patchdata_nonempty([&](Patch cur_p, PatchDataLayer &pdat) { - sham::DeviceBuffer &buf_s_j = pdat.get_field_buf_ref(is_j); - - sham::kernel_call( - q, - sham::MultiRef{}, - sham::MultiRef{buf_s_j}, - pdat.get_obj_cnt() * ndust, - [](u32 tid, Tscal *s_j) { - Tscal s_j_tid = s_j[tid]; - if (s_j_tid < 0) { - throw std::runtime_error("s_j is negative"); - } - }); - }); - q.q.wait_and_throw(); - } - - - // check positivity - if (solver_config.dust_config.has_s_j_field()) { - PatchDataLayerLayout &pdl = scheduler().pdl_old(); - const u32 is_j = pdl.get_field_idx("s_j") ; - auto ndust = solver_config.dust_config.get_dust_nvar(); - auto &q = shamsys::instance::get_compute_scheduler().get_queue(); - scheduler().for_each_patchdata_nonempty([&](Patch cur_p, PatchDataLayer &pdat) { - sham::DeviceBuffer &buf_s_j = pdat.get_field_buf_ref(is_j); - - sham::kernel_call( - q, - sham::MultiRef{}, - sham::MultiRef{buf_s_j}, - pdat.get_obj_cnt() * ndust, - [](u32 tid, Tscal *s_j) { - Tscal s_j_tid = s_j[tid]; - if (s_j_tid < 0) { - throw std::runtime_error("s_j is negative"); - } - }); - }); - } - using RTree = shamtree::CompressedLeafBVH; // Here we will add self grav to the external forces indep of vel (this will be moved into a @@ -2778,14 +2764,14 @@ shammodels::sph::TimestepLog shammodels::sph::Solver::evolve_once() storage.part_counts, C_cour_edge, hpart_refs, vclean_dt, cfl_dt); } - std::shared_ptr> compute_cfl_dust1_fluid; + std::shared_ptr> compute_cfl_dust1_fluid; std::shared_ptr> s_j_refs; std::shared_ptr> hfactd_edge; if (solver_config.dust_config.has_s_j_field()) { u32 ndust = solver_config.dust_config.get_dust_nvar(); - compute_cfl_dust1_fluid = std::make_shared>(ndust); + compute_cfl_dust1_fluid = std::make_shared>(ndust); auto t_j_field = storage.solver_graph @@ -2804,9 +2790,14 @@ shammodels::sph::TimestepLog shammodels::sph::Solver::evolve_once() map_field_refs(scheduler(), is_j, *s_j_refs); + std::shared_ptr> C_1fluid_edge + = std::make_shared>( + "C_1fluid", "C_{1fluid}"); + C_1fluid_edge->value = solver_config.dust_config.get_monofluid_tvi().C_1_fluid; + compute_cfl_dust1_fluid->set_edges( storage.part_counts, - C_cour_edge, + C_1fluid_edge, pmass_edge, hfactd_edge, hpart_refs, @@ -2816,6 +2807,28 @@ shammodels::sph::TimestepLog shammodels::sph::Solver::evolve_once() cfl_dt); } + std::shared_ptr> compute_cfl_dust_deltav; + std::shared_ptr> delta_v_refs; + + if (solver_config.dust_config.has_s_j_field()) { + u32 ndust = solver_config.dust_config.get_dust_nvar(); + + compute_cfl_dust_deltav = std::make_shared>(ndust); + + delta_v_refs = std::make_shared>( + "delta_v", "delta_v"); + const u32 idelta_v = pdl.get_field_idx("delta_v"); + map_field_refs(scheduler(), idelta_v, *delta_v_refs); + + std::shared_ptr> C_delta_v_edge + = std::make_shared>( + "C_delta_v", "C_{delta_v}"); + C_delta_v_edge->value = solver_config.dust_config.get_monofluid_tvi().C_delta_v; + + compute_cfl_dust_deltav->set_edges( + storage.part_counts, C_delta_v_edge, hpart_refs, delta_v_refs, cfl_dt); + } + bool show_cfl_detail = solver_config.show_cfl_detail; std::vector> cfl_detail; @@ -2845,6 +2858,9 @@ shammodels::sph::TimestepLog shammodels::sph::Solver::evolve_once() if (solver_config.dust_config.has_s_j_field()) { compute_cfl_dust1_fluid->evaluate(); save_cfl_detail("dust1_fluid"); + + compute_cfl_dust_deltav->evaluate(); + save_cfl_detail("dust_deltav"); } if (!show_cfl_detail) { @@ -2981,31 +2997,6 @@ shammodels::sph::TimestepLog shammodels::sph::Solver::evolve_once() reset_merge_ghosts_fields(); reset_eos_fields(); - - - // check positivity - if (solver_config.dust_config.has_s_j_field()) { - PatchDataLayerLayout &pdl = scheduler().pdl_old(); - const u32 is_j = pdl.get_field_idx("s_j") ; - auto ndust = solver_config.dust_config.get_dust_nvar(); - auto &q = shamsys::instance::get_compute_scheduler().get_queue(); - scheduler().for_each_patchdata_nonempty([&](Patch cur_p, PatchDataLayer &pdat) { - sham::DeviceBuffer &buf_s_j = pdat.get_field_buf_ref(is_j); - - sham::kernel_call( - q, - sham::MultiRef{}, - sham::MultiRef{buf_s_j}, - pdat.get_obj_cnt() * ndust, - [](u32 tid, Tscal *s_j) { - Tscal s_j_tid = s_j[tid]; - if (s_j_tid < 0) { - throw std::runtime_error("s_j is negative"); - } - }); - }); - } - // if delta too big jump to compute force tstep.stop(); diff --git a/src/shammodels/sph/src/SolverConfig.cpp b/src/shammodels/sph/src/SolverConfig.cpp index b5d3c8d957..a1614e3a21 100644 --- a/src/shammodels/sph/src/SolverConfig.cpp +++ b/src/shammodels/sph/src/SolverConfig.cpp @@ -116,6 +116,7 @@ namespace shammodels::sph { u32 ndust = dust_config.get_dust_nvar(); pdl.add_field("s_j", ndust); pdl.add_field("ds_j_dt", ndust); + pdl.add_field("delta_v", ndust); } } diff --git a/src/shammodels/sph/src/modules/NodeUpdateDerivsMonofluidTVI.cpp b/src/shammodels/sph/src/modules/NodeUpdateDerivsMonofluidTVI.cpp index 53379af2b7..b50b8a0259 100644 --- a/src/shammodels/sph/src/modules/NodeUpdateDerivsMonofluidTVI.cpp +++ b/src/shammodels/sph/src/modules/NodeUpdateDerivsMonofluidTVI.cpp @@ -98,18 +98,6 @@ struct KernelUpdateDerivsMonofluidTVI { Tscal delta_P = P_a - P_b; Tscal Ts_weighted = (Ttilde_sj_a / rho_a + Ttilde_sj_b / rho_b); - if (id_a == 29951 && jdust == 4 && false) { - printf("--------------------------------\n"); - printf("id_b: %d, s_j_b: %e, Ttilde_sj_b: %e\n", id_b, s_j_b, Ttilde_sj_b); - printf("rho_b: %e, P_b: %e, omega_b: %f\n", rho_b, P_b, omega_b); - printf("Fab_a: %e, Fab_b: %e\n", Fab_a, Fab_b); - printf( - "F_ab_bar: %e, delta_P: %e, Ts_weighted: %e\n", F_ab_bar, delta_P, Ts_weighted); - printf("rab: %e, rab_inv_sat: %e\n", rab, rab_inv_sat); - printf("v_ab: %f, %f, %f\n", v_ab[0], v_ab[1], v_ab[2]); - printf("r_ab_unit: %f, %f, %f\n", r_ab_unit[0], r_ab_unit[1], r_ab_unit[2]); - } - term1 += (pmass * s_j_b / rho_b) * Ts_weighted * delta_P * F_ab_bar * rab_inv_sat; term2 += pmass * sham::dot(v_ab, r_ab_unit * Fab_a); }); @@ -123,34 +111,6 @@ struct KernelUpdateDerivsMonofluidTVI { // restore it slowly to 0 ds_j_dt_a += (s_j_a < 0) ? -s_j_a / (10 * Ttilde_sj_a) : 0; - if (id_a == 29951 && jdust == 4 && false) { - printf("s_j_a: %e\n", s_j_a); - printf("ds_j_dt_a: %e\n", ds_j_dt_a); - printf("h_a: %e\n", h_a); - printf("rho_a: %e\n", rho_a); - printf("omega_a: %f\n", omega_a); - printf("Ttilde_sj_a: %f\n", Ttilde_sj_a); - printf("term1: %e\n", Tscal{-0.5} * term1); - printf("term2: %e\n", (s_j_a / (2 * rho_a * omega_a)) * term2); - printf("xyz_a: %f, %f, %f\n", xyz_a[0], xyz_a[1], xyz_a[2]); - } - - if (ds_j_dt_a > 1e-3 && s_j_a < 1e-6 ) { - printf("id_a: %d, s_j_a: %e, ds_j_dt_a: %e\n", id_a, s_j_a, ds_j_dt_a); - //printf( - // "h_a: %e, rho_a: %e, omega_a: %f, Ttilde_sj_a: %f\n", - // h_a, - // rho_a, - // omega_a, - // Ttilde_sj_a); - //printf( - // "term1: %e, term2: %e\n", - // Tscal{-0.5} * term1, - // (s_j_a / (2 * rho_a * omega_a)) * term2); - //printf("xyz_a: %f, %f, %f\n", xyz_a[0], xyz_a[1], xyz_a[2]); - //throw std::runtime_error("ds_j_dt_a is too large"); - } - ds_j_dt[thread_id] = ds_j_dt_a; } }; diff --git a/src/shammodels/sph/src/modules/io/VTKDump.cpp b/src/shammodels/sph/src/modules/io/VTKDump.cpp index a43a1cc4ad..5ed9e1da42 100644 --- a/src/shammodels/sph/src/modules/io/VTKDump.cpp +++ b/src/shammodels/sph/src/modules/io/VTKDump.cpp @@ -125,6 +125,11 @@ namespace shammodels::sph::modules { fnum += ndust; } + if (solver_config.dust_config.has_s_j_field()) { + const u32 ndust = solver_config.dust_config.get_dust_nvar(); + fnum += ndust * 3; // s_j, ds_j_dt and delta_v + } + writer.add_field_data_section(fnum); if (add_patch_world_id) { @@ -236,6 +241,109 @@ namespace shammodels::sph::modules { scheduler(), writer, tmp_deltav, "deltav_" + std::to_string(idust)); } } + + if (solver_config.dust_config.has_s_j_field()) { + const u32 is_j = pdl.get_field_idx("s_j"); + const u32 ndust = solver_config.dust_config.get_dust_nvar(); + + for (u32 idust = 0; idust < ndust; idust++) { + ComputeField tmp_s_j = utility.make_compute_field("tmp_s_j", 1); + + scheduler().for_each_patchdata_nonempty([&](const Patch p, PatchDataLayer &pdat) { + shamlog_debug_ln( + "sph::vtk", "compute extract s_j field with idust =", idust, p.id_patch); + + auto &buf_s_j = pdat.get_field(is_j); + PatchDataFieldSpan span_s_j{buf_s_j, 0, pdat.get_obj_cnt()}; + + auto sptr = shamsys::instance::get_compute_scheduler_ptr(); + auto &q = sptr->get_queue(); + + sham::kernel_call( + q, + sham::MultiRef{span_s_j}, + sham::MultiRef{tmp_s_j.get_buf(p.id_patch)}, + pdat.get_obj_cnt(), + [&, idust](u32 i, auto s_j_field, Tscal *acc_s_j) { + acc_s_j[i] = s_j_field(i, idust); + }); + }); + + vtk_dump_add_compute_field( + scheduler(), writer, tmp_s_j, "s_j_" + std::to_string(idust)); + } + } + + if (solver_config.dust_config.has_s_j_field()) { + const u32 ids_j_dt = pdl.get_field_idx("ds_j_dt"); + const u32 ndust = solver_config.dust_config.get_dust_nvar(); + + for (u32 idust = 0; idust < ndust; idust++) { + ComputeField tmp_ds_j_dt + = utility.make_compute_field("tmp_ds_j_dt", 1); + + scheduler().for_each_patchdata_nonempty([&](const Patch p, PatchDataLayer &pdat) { + shamlog_debug_ln( + "sph::vtk", + "compute extract ds_j_dt field with idust =", + idust, + p.id_patch); + + auto &buf_ds_j_dt = pdat.get_field(ids_j_dt); + PatchDataFieldSpan span_ds_j_dt{buf_ds_j_dt, 0, pdat.get_obj_cnt()}; + + auto sptr = shamsys::instance::get_compute_scheduler_ptr(); + auto &q = sptr->get_queue(); + + sham::kernel_call( + q, + sham::MultiRef{span_ds_j_dt}, + sham::MultiRef{tmp_ds_j_dt.get_buf(p.id_patch)}, + pdat.get_obj_cnt(), + [&, idust](u32 i, auto ds_j_dt_field, Tscal *acc_ds_j_dt) { + acc_ds_j_dt[i] = ds_j_dt_field(i, idust); + }); + }); + + vtk_dump_add_compute_field( + scheduler(), writer, tmp_ds_j_dt, "ds_j_dt_" + std::to_string(idust)); + } + } + + if (solver_config.dust_config.has_s_j_field()) { + const u32 idelta_v = pdl.get_field_idx("delta_v"); + const u32 ndust = solver_config.dust_config.get_dust_nvar(); + + for (u32 idust = 0; idust < ndust; idust++) { + ComputeField tmp_delta_v = utility.make_compute_field("tmp_delta_v", 1); + + scheduler().for_each_patchdata_nonempty([&](const Patch p, PatchDataLayer &pdat) { + shamlog_debug_ln( + "sph::vtk", + "compute extract delta_v field with idust =", + idust, + p.id_patch); + + auto &buf_delta_v = pdat.get_field(idelta_v); + PatchDataFieldSpan span_delta_v{buf_delta_v, 0, pdat.get_obj_cnt()}; + + auto sptr = shamsys::instance::get_compute_scheduler_ptr(); + auto &q = sptr->get_queue(); + + sham::kernel_call( + q, + sham::MultiRef{span_delta_v}, + sham::MultiRef{tmp_delta_v.get_buf(p.id_patch)}, + pdat.get_obj_cnt(), + [&, idust](u32 i, auto delta_v_field, Tvec *acc_delta_v) { + acc_delta_v[i] = delta_v_field(i, idust); + }); + }); + + vtk_dump_add_compute_field( + scheduler(), writer, tmp_delta_v, "delta_v_" + std::to_string(idust)); + } + } } } // namespace shammodels::sph::modules diff --git a/src/shammodels/sph/src/pySPHModel.cpp b/src/shammodels/sph/src/pySPHModel.cpp index e2d40a5b7b..08d192df76 100644 --- a/src/shammodels/sph/src/pySPHModel.cpp +++ b/src/shammodels/sph/src/pySPHModel.cpp @@ -267,12 +267,18 @@ void add_instance(py::module &m, std::string name_config, std::string name_model }) .def( "set_dust_mode_monofluid_tvi", - [](TConfig &self, u32 nvar, bool pure_diffusion_mode) { - self.dust_config.set_monofluid_tvi(nvar, pure_diffusion_mode); + [](TConfig &self, + u32 nvar, + bool pure_diffusion_mode, + Tscal C_1_fluid, + Tscal C_delta_v) { + self.dust_config.set_monofluid_tvi(nvar, pure_diffusion_mode, C_1_fluid, C_delta_v); }, py::kw_only(), py::arg("nvar"), - py::arg("pure_diffusion_mode") = false) + py::arg("pure_diffusion_mode") = false, + py::arg("C_1_fluid") = 0.1, + py::arg("C_delta_v") = 1.0) .def( "set_dust_mode_monofluid_complete", [](TConfig &self, u32 ndust) { diff --git a/src/shamrock/src/math/integrators.cpp b/src/shamrock/src/math/integrators.cpp index baef01fd01..99e8d963e4 100644 --- a/src/shamrock/src/math/integrators.cpp +++ b/src/shamrock/src/math/integrators.cpp @@ -15,10 +15,10 @@ */ #include "shambase/exception.hpp" +#include "shambackends/math.hpp" #include "shambackends/sycl_utils.hpp" #include "shamrock/math/integrators.hpp" #include -#include "shambackends/math.hpp" namespace integ = shamrock::integrators; namespace util = shamrock::utilities; From a04cec82cb5b422e7cc71a0fcdf79154cf93890b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?David--Cl=C3=A9ris=20Timoth=C3=A9e?= Date: Mon, 13 Jul 2026 14:18:42 +0200 Subject: [PATCH 19/24] Apply suggestion from @tdavidcl --- examples/sph/run_dustysettle_tvi.py | 1 - 1 file changed, 1 deletion(-) diff --git a/examples/sph/run_dustysettle_tvi.py b/examples/sph/run_dustysettle_tvi.py index 903c18c261..e64fba7336 100644 --- a/examples/sph/run_dustysettle_tvi.py +++ b/examples/sph/run_dustysettle_tvi.py @@ -879,7 +879,6 @@ def dust_mass_analysis(): for j in range(1000): if j == iinject: reference_dusty_settle = ReferenceDustySettleAll() - reference_dusty_settle = None if tlist[j] >= t_start: if j > 0: From 93b88cfbe7a3af9dc7d3a11106aab00c4bca5cee Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?David--Cl=C3=A9ris=20Timoth=C3=A9e?= Date: Tue, 14 Jul 2026 13:02:29 +0200 Subject: [PATCH 20/24] Apply suggestion from @tdavidcl --- examples/sph/run_dustysettle_tvi.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/sph/run_dustysettle_tvi.py b/examples/sph/run_dustysettle_tvi.py index e64fba7336..306f8c72eb 100644 --- a/examples/sph/run_dustysettle_tvi.py +++ b/examples/sph/run_dustysettle_tvi.py @@ -758,7 +758,7 @@ def analyse_and_plot(j): axs[1].set_ylabel(r"$\dot{s}_j$") axs[1].set_xlabel(r"$z$") - axs[1].set_xlim(-box, box) + axs[1].set_xlim(-4 * H, 4 * H) axs[1].set_yscale("symlog", linthresh=1e-10) dust_sm = cm.ScalarMappable(cmap=dust_cmap, norm=dust_norm) From ee2504b63e064d7447ac6c3708d41db77d0615bf Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?David--Cl=C3=A9ris=20Timoth=C3=A9e?= Date: Tue, 14 Jul 2026 13:02:47 +0200 Subject: [PATCH 21/24] Apply suggestion from @tdavidcl --- examples/sph/run_dustysettle_tvi.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/sph/run_dustysettle_tvi.py b/examples/sph/run_dustysettle_tvi.py index 306f8c72eb..1efe312491 100644 --- a/examples/sph/run_dustysettle_tvi.py +++ b/examples/sph/run_dustysettle_tvi.py @@ -914,7 +914,7 @@ def compute_sj_new(patchdata): analyse_and_plot(j) - dump_helper.write_dump(j, purge_old_dumps=True, keep_first=1, keep_last=5) + dump_helper.write_dump(j, purge_old_dumps=True, keep_first=1, keep_last=3) if tlist[j] >= t_end: break From c08a25ff641831abcc4811d1454a60f9c0bbbdad Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?David--Cl=C3=A9ris=20Timoth=C3=A9e?= Date: Tue, 14 Jul 2026 13:03:16 +0200 Subject: [PATCH 22/24] Apply suggestion from @tdavidcl --- examples/sph/run_dustysettle_tvi.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/sph/run_dustysettle_tvi.py b/examples/sph/run_dustysettle_tvi.py index 1efe312491..f41bf3441e 100644 --- a/examples/sph/run_dustysettle_tvi.py +++ b/examples/sph/run_dustysettle_tvi.py @@ -752,7 +752,7 @@ def analyse_and_plot(j): axs[0].set_ylabel(r"$s_j$") axs[0].set_xlabel(r"$z$") - axs[0].set_xlim(-box, box) + axs[0].set_xlim(-4 * H, 4 * H) axs[0].set_yscale("log") axs[0].set_ylim(1e-20, 1e-1) From 91a005b71cd69cbe2a39a30e10cd6776f743ff07 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Timoth=C3=A9e=20David--Cl=C3=A9ris?= Date: Wed, 15 Jul 2026 15:45:21 +0200 Subject: [PATCH 23/24] add option for positivity --- .../include/shammodels/sph/SolverConfig.hpp | 18 +++++-- src/shammodels/sph/src/Solver.cpp | 48 ++++++++++++++----- src/shammodels/sph/src/pySPHModel.cpp | 13 +++-- 3 files changed, 59 insertions(+), 20 deletions(-) diff --git a/src/shammodels/sph/include/shammodels/sph/SolverConfig.hpp b/src/shammodels/sph/include/shammodels/sph/SolverConfig.hpp index b3484a02e2..2030fe2add 100644 --- a/src/shammodels/sph/include/shammodels/sph/SolverConfig.hpp +++ b/src/shammodels/sph/include/shammodels/sph/SolverConfig.hpp @@ -122,6 +122,8 @@ namespace shammodels::sph { Tscal C_1_fluid = 0.1; Tscal C_delta_v = 1.0; Tscal cfl_density_threshold = shambase::get_epsilon(); + + bool ensure_s_j_positivity = true; }; struct MonofluidComplete { @@ -139,9 +141,15 @@ namespace shammodels::sph { bool pure_diffusion_mode = false, Tscal C_1_fluid = 0.1, Tscal C_delta_v = 1.0, - Tscal cfl_density_threshold = shambase::get_epsilon()) { + Tscal cfl_density_threshold = shambase::get_epsilon(), + bool ensure_s_j_positivity = true) { current_mode = MonofluidTVI{ - nvar, pure_diffusion_mode, C_1_fluid, C_delta_v, cfl_density_threshold}; + nvar, + pure_diffusion_mode, + C_1_fluid, + C_delta_v, + cfl_density_threshold, + ensure_s_j_positivity}; } inline void set_monofluid_complete(u32 nvar) { current_mode = MonofluidComplete{nvar}; } @@ -165,7 +173,8 @@ namespace shammodels::sph { {"pure_diffusion_mode", cfg->pure_diffusion_mode}, {"C_1_fluid", cfg->C_1_fluid}, {"C_delta_v", cfg->C_delta_v}, - {"cfl_density_threshold", cfg->cfl_density_threshold}}; + {"cfl_density_threshold", cfg->cfl_density_threshold}, + {"ensure_s_j_positivity", cfg->ensure_s_j_positivity}}; } else if ( const MonofluidComplete *cfg = std::get_if(¤t_mode)) { j = {{"type", "monofluid_complete"}, {"ndust", cfg->ndust}}; @@ -184,7 +193,8 @@ namespace shammodels::sph { j.at("pure_diffusion_mode").get(), j.at("C_1_fluid").get(), j.at("C_delta_v").get(), - j.at("cfl_density_threshold").get()); + j.at("cfl_density_threshold").get(), + j.at("ensure_s_j_positivity").get()); } else if (type == "monofluid_complete") { set_monofluid_complete(j.at("ndust").get()); } else { diff --git a/src/shammodels/sph/src/Solver.cpp b/src/shammodels/sph/src/Solver.cpp index b321f523ba..d4f26c1d67 100644 --- a/src/shammodels/sph/src/Solver.cpp +++ b/src/shammodels/sph/src/Solver.cpp @@ -455,17 +455,32 @@ void shammodels::sph::Solver::init_solver_graph() { } if (has_s_j_field) { - u32 ndust = solver_config.dust_config.get_dust_nvar(); - auto half_step_s_j = solver_graph.register_node( - prefix + "_s_j", - shammodels::common::modules::ForwardEulerPositive(ndust)); - shambase::get_check_ref(half_step_s_j) - .set_edges( - solver_graph.get_edge_ptr>("dt_half"), - solver_graph.get_edge_ptr>("ds_j_dt"), - solver_graph.get_edge_ptr>("part_counts"), - solver_graph.get_edge_ptr>("s_j")); - half_step_sequence.push_back(half_step_s_j); + u32 ndust = solver_config.dust_config.get_dust_nvar(); + + auto &cfg = solver_config.dust_config.get_monofluid_tvi(); + + if (cfg.ensure_s_j_positivity) { + auto half_step_s_j = solver_graph.register_node( + prefix + "_s_j", + shammodels::common::modules::ForwardEulerPositive(ndust)); + shambase::get_check_ref(half_step_s_j) + .set_edges( + solver_graph.get_edge_ptr>("dt_half"), + solver_graph.get_edge_ptr>("ds_j_dt"), + solver_graph.get_edge_ptr>("part_counts"), + solver_graph.get_edge_ptr>("s_j")); + half_step_sequence.push_back(half_step_s_j); + } else { + auto half_step_s_j = solver_graph.register_node( + prefix + "_s_j", shammodels::common::modules::ForwardEuler(ndust)); + shambase::get_check_ref(half_step_s_j) + .set_edges( + solver_graph.get_edge_ptr>("dt_half"), + solver_graph.get_edge_ptr>("ds_j_dt"), + solver_graph.get_edge_ptr>("part_counts"), + solver_graph.get_edge_ptr>("s_j")); + half_step_sequence.push_back(half_step_s_j); + } } return OperationSequence("half step", std::move(half_step_sequence)); @@ -2389,8 +2404,15 @@ shammodels::sph::TimestepLog shammodels::sph::Solver::evolve_once() if (solver_config.dust_config.has_s_j_field()) { ComputeField s_j_s_j_sq = utility.make_compute_field( "s_j s_j^2", solver_config.dust_config.get_dust_nvar()); - utility.fields_leapfrog_corrector_positive_only( - is_j, ids_j_dt, storage.old_ds_j_dt.get(), s_j_s_j_sq, dt / 2); + bool ensure_positivity + = solver_config.dust_config.get_monofluid_tvi().ensure_s_j_positivity; + if (ensure_positivity) { + utility.fields_leapfrog_corrector_positive_only( + is_j, ids_j_dt, storage.old_ds_j_dt.get(), s_j_s_j_sq, dt / 2); + } else { + utility.fields_leapfrog_corrector( + is_j, ids_j_dt, storage.old_ds_j_dt.get(), s_j_s_j_sq, dt / 2); + } } storage.old_axyz.reset(); diff --git a/src/shammodels/sph/src/pySPHModel.cpp b/src/shammodels/sph/src/pySPHModel.cpp index 381d95d0fc..3cf961c18a 100644 --- a/src/shammodels/sph/src/pySPHModel.cpp +++ b/src/shammodels/sph/src/pySPHModel.cpp @@ -272,16 +272,23 @@ void add_instance(py::module &m, std::string name_config, std::string name_model bool pure_diffusion_mode, Tscal C_1_fluid, Tscal C_delta_v, - Tscal cfl_density_threshold) { + Tscal cfl_density_threshold, + bool ensure_s_j_positivity) { self.dust_config.set_monofluid_tvi( - nvar, pure_diffusion_mode, C_1_fluid, C_delta_v, cfl_density_threshold); + nvar, + pure_diffusion_mode, + C_1_fluid, + C_delta_v, + cfl_density_threshold, + ensure_s_j_positivity); }, py::kw_only(), py::arg("nvar"), py::arg("pure_diffusion_mode") = false, py::arg("C_1_fluid") = 0.1, py::arg("C_delta_v") = 1.0, - py::arg("cfl_density_threshold") = shambase::get_epsilon()) + py::arg("cfl_density_threshold") = shambase::get_epsilon(), + py::arg("ensure_s_j_positivity") = true) .def( "set_dust_mode_monofluid_complete", [](TConfig &self, u32 ndust) { From 3be4e99c6883ce0143427d7f6f7885902a8f30b9 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?David--Cl=C3=A9ris=20Timoth=C3=A9e?= Date: Wed, 15 Jul 2026 20:54:24 +0200 Subject: [PATCH 24/24] Add run_dustywave_sympy.py to examples/physics --- examples/physics/{dustywave_sympy.py => run_dustywave_sympy.py} | 0 1 file changed, 0 insertions(+), 0 deletions(-) rename examples/physics/{dustywave_sympy.py => run_dustywave_sympy.py} (100%) diff --git a/examples/physics/dustywave_sympy.py b/examples/physics/run_dustywave_sympy.py similarity index 100% rename from examples/physics/dustywave_sympy.py rename to examples/physics/run_dustywave_sympy.py