diff --git a/examples/physics/run_dustywave_sympy.py b/examples/physics/run_dustywave_sympy.py new file mode 100644 index 0000000000..776994b972 --- /dev/null +++ b/examples/physics/run_dustywave_sympy.py @@ -0,0 +1,257 @@ +""" +Dustywave TVA dispersion relation +======================================= + +This example shows how to derive the dustywave TVA dispersion relation using SymPy +""" + +import matplotlib.pyplot as plt +import numpy as np +import sympy as sp + +# %% +# 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() diff --git a/examples/sph/run_dustysettle_tvi.py b/examples/sph/run_dustysettle_tvi.py index 86f7dcc004..f41bf3441e 100644 --- a/examples/sph/run_dustysettle_tvi.py +++ b/examples/sph/run_dustysettle_tvi.py @@ -753,7 +753,8 @@ 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_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$") 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..ac7ab042d1 --- /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 ForwardEulerPositive.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/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 0e30bc7476..d4f26c1d67 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" @@ -454,16 +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::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); + 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)); @@ -2387,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( - 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) { 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..99e8d963e4 100644 --- a/src/shamrock/src/math/integrators.cpp +++ b/src/shamrock/src/math/integrators.cpp @@ -15,6 +15,7 @@ */ #include "shambase/exception.hpp" +#include "shambackends/math.hpp" #include "shambackends/sycl_utils.hpp" #include "shamrock/math/integrators.hpp" #include @@ -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 //////////////////////////////////////////////////////////////////////////////////////////////////// ////////////////////////////////////////////////////////////////////////////////////////////////////