diff --git a/src/shammodels/sph/CMakeLists.txt b/src/shammodels/sph/CMakeLists.txt index aabe7f345..da1dc5cf8 100644 --- a/src/shammodels/sph/CMakeLists.txt +++ b/src/shammodels/sph/CMakeLists.txt @@ -49,6 +49,7 @@ set(Sources src/modules/NodeUpdateDerivsVaryingAlphaAV.cpp src/modules/NodeComputePressureGrad.cpp src/modules/NodeEvolveDustCOALASourceTerm.cpp + src/modules/NodeEvolveDustCOALASourceTermOrmel.cpp src/modules/setup/GeneratorMCDisc.cpp src/modules/setup/ModifierApplyDiscWarp.cpp src/modules/setup/ModifierApplyCustomWarp.cpp diff --git a/src/shammodels/sph/include/shammodels/sph/modules/NodeEvolveDustCOALASourceTerm.hpp b/src/shammodels/sph/include/shammodels/sph/modules/NodeEvolveDustCOALASourceTerm.hpp index 8f804bdd6..412c9d7fc 100644 --- a/src/shammodels/sph/include/shammodels/sph/modules/NodeEvolveDustCOALASourceTerm.hpp +++ b/src/shammodels/sph/include/shammodels/sph/modules/NodeEvolveDustCOALASourceTerm.hpp @@ -65,3 +65,5 @@ namespace shammodels::sph::modules { inline virtual std::string _impl_get_tex() const; }; } // namespace shammodels::sph::modules + +#undef NODE_EVOLVE_DUST_COALA_SOURCE_TERM_EDGES diff --git a/src/shammodels/sph/include/shammodels/sph/modules/NodeEvolveDustCOALASourceTermOrmel.hpp b/src/shammodels/sph/include/shammodels/sph/modules/NodeEvolveDustCOALASourceTermOrmel.hpp new file mode 100644 index 000000000..e6371d4c0 --- /dev/null +++ b/src/shammodels/sph/include/shammodels/sph/modules/NodeEvolveDustCOALASourceTermOrmel.hpp @@ -0,0 +1,77 @@ +// -------------------------------------------------------// +// +// 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 NodeEvolveDustCOALASourceTermOrmel.hpp + * @author Timothée David--Cléris (tim.shamrock@proton.me) + * @brief + * + */ + +#include "shambackends/vec.hpp" +#include "shamrock/solvergraph/IFieldSpan.hpp" +#include "shamrock/solvergraph/INode.hpp" +#include "shamrock/solvergraph/Indexes.hpp" +#include "shamrock/solvergraph/ScalarEdge.hpp" +#include +#include + +#define NODE_EVOLVE_DUST_COALA_SOURCE_TERM_EDGES(X_RO, X_RW) \ + /* scalars */ \ + X_RO(shamrock::solvergraph::ScalarEdge, rhodust_eps) \ + X_RO(shamrock::solvergraph::ScalarEdge>, massgrid) \ + X_RO(shamrock::solvergraph::ScalarEdge>, tensor_tabflux_coag) \ + X_RO(shamrock::solvergraph::ScalarEdge, alpha_turb) \ + X_RO(shamrock::solvergraph::ScalarEdge, mu) \ + X_RO(shamrock::solvergraph::ScalarEdge, mh) \ + X_RO(shamrock::solvergraph::ScalarEdge, kb) \ + X_RO(shamrock::solvergraph::ScalarEdge, gpart_mass) \ + \ + X_RO(shamrock::solvergraph::IFieldSpan, hpart) \ + X_RO(shamrock::solvergraph::IFieldSpan, pressure) \ + X_RO(shamrock::solvergraph::IFieldSpan, cs) \ + X_RO(shamrock::solvergraph::IFieldSpan, ts) \ + X_RO(shamrock::solvergraph::IFieldSpan, a_ext) \ + \ + /* counts */ \ + X_RO(shamrock::solvergraph::Indexes, part_counts) \ + \ + /* to get rho_dust_j */ \ + X_RO(shamrock::solvergraph::IFieldSpan, s_j) \ + \ + /* outputs */ \ + X_RW(shamrock::solvergraph::IFieldSpan, S_coag) + +namespace shammodels::sph::modules { + + template class SPHKernel> + class NodeEvolveDustCOALASourceTermOrmel : public shamrock::solvergraph::INode { + + using Tscal = shambase::VecComponent; + + u32 nbins; + + public: + NodeEvolveDustCOALASourceTermOrmel(u32 nbins) : nbins(nbins) {} + + EXPAND_NODE_EDGES(NODE_EVOLVE_DUST_COALA_SOURCE_TERM_EDGES) + + void _impl_evaluate_internal(); + + inline virtual std::string _impl_get_label() const { + return "NodeEvolveDustCOALASourceTermOrmel"; + }; + + inline virtual std::string _impl_get_tex() const; + }; +} // namespace shammodels::sph::modules + +#undef NODE_EVOLVE_DUST_COALA_SOURCE_TERM_EDGES diff --git a/src/shammodels/sph/src/modules/NodeEvolveDustCOALASourceTermOrmel.cpp b/src/shammodels/sph/src/modules/NodeEvolveDustCOALASourceTermOrmel.cpp new file mode 100644 index 000000000..edb566caa --- /dev/null +++ b/src/shammodels/sph/src/modules/NodeEvolveDustCOALASourceTermOrmel.cpp @@ -0,0 +1,349 @@ +// -------------------------------------------------------// +// +// 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 +// +// -------------------------------------------------------// + +/** + * @file NodeEvolveDustCOALASourceTermOrmel.cpp + * @author Timothée David--Cléris (tim.shamrock@proton.me) + * @brief + * + */ + +#include "shambase/memory.hpp" +#include "shambase/stacktrace.hpp" +#include "shambase/string.hpp" +#include "shambackends/DeviceBuffer.hpp" +#include "shambackends/kernel_call.hpp" +#include "shambackends/math.hpp" +#include "shambackends/vec.hpp" +#include "shamcomm/logs.hpp" +#include "shammath/sphkernels.hpp" +#include "shammodels/sph/math/density.hpp" +#include "shammodels/sph/modules/NodeEvolveDustCOALASourceTermOrmel.hpp" +#include "shamphys/coala_interface.hpp" +#include "shamrock/patch/PatchDataField.hpp" // IWYU pragma: keep +#include "shamsys/NodeInstance.hpp" +#include +#include + +namespace shammodels::sph::modules { + + template class SPHKernel> + struct KernelGenCoala_k0_Ormel { + + using Tscal = shambase::VecComponent; + + using Kernel = SPHKernel; + static constexpr Tscal hfactd = Kernel::hfactd; + static constexpr Tscal Rkern = Kernel::Rkern; + static constexpr Tscal Rker2 = Rkern * Rkern; + + using mdspan_rank_1 = std::mdspan>; + using mdspan_rank_3 = std::mdspan>; + + using const_mdspan_rank_1 = std::mdspan>; + using const_mdspan_rank_3 = std::mdspan>; + + u32 nbins; + Tscal rho_eps; + u32 corrected_len; + u32 group_size; + u32 true_size; + + Tscal pmass; + Tscal mu; + Tscal mh; + Tscal kb; + Tscal alpha_turb; + + auto operator()( + u32 /**/, + // common to all kernel calls + const Tscal *__restrict massgrid_ptr, + const Tscal *__restrict tensor_tabflux_coag, + // field specific data + const Tscal *__restrict s_j, + const Tscal *__restrict hpart, + const Tscal *__restrict pressure, + const Tscal *__restrict cs, + const Tscal *__restrict ts, + + const Tvec *__restrict a_ext, + Tscal *__restrict S_coag) const { + + auto range = sycl::nd_range<1>{corrected_len, group_size}; + + auto local_acc_sz_nbins = sycl::range<1>{group_size * nbins}; + + auto true_size = this->true_size; + auto rho_eps = this->rho_eps; + + Tscal pmass = this->pmass; + Tscal mu = this->mu; + Tscal mh = this->mh; + Tscal kb = this->kb; + Tscal alpha_turb = this->alpha_turb; + + return [=, nbins = this->nbins](sycl::handler &cgh) { + auto gij_acc = sycl::local_accessor{local_acc_sz_nbins, cgh}; + auto flux_acc = sycl::local_accessor{local_acc_sz_nbins, cgh}; + + cgh.parallel_for(range, [=](sycl::nd_item<1> tid) { + const u64 id_a = tid.get_global_linear_id(); + const u64 lid = tid.get_local_linear_id(); + + if (id_a >= true_size) { + return; + } + + u32 id_a_d = id_a * nbins; + + /* inputs */ + const_mdspan_rank_3 tabflux_coag(tensor_tabflux_coag, nbins, nbins, nbins); + const_mdspan_rank_1 massgrid(massgrid_ptr, nbins + 1); + + /* internal */ + auto gij_loc = &(gij_acc[nbins * lid]); + auto flux_loc = &(flux_acc[nbins * lid]); + + mdspan_rank_1 gij(gij_loc, nbins); + mdspan_rank_1 flux(flux_loc, nbins); + + /* output */ + mdspan_rank_1 S_coag_span(S_coag + id_a_d, nbins); + + /* lambda getters */ + auto rho_dust = [&](int j) { + auto tmp = s_j[id_a_d + j]; + return tmp * tmp; + }; + + // auto dv = [delta_v = delta_v_j + id_a_d](int i, int j) { + // // dv_ij = v_dust_j - v_dust_i = delta_v_j[j] - delta_v_j[i] + // return sycl::length(delta_v[j] - delta_v[i]); + // }; + + Tscal cs_a = cs[id_a]; + Tscal a_ext_a = sycl::length(a_ext[id_a]); + Tscal h_a = hpart[id_a]; + Tscal P_a = pressure[id_a]; + + // TODO use true gas density instead of combined dust + gas + Tscal rho_gas_a = shamrock::sph::rho_h(pmass, h_a, Kernel::hfactd); + + auto temperature = [](Tscal P, Tscal rho, Tscal mu, Tscal mh, Tscal kb) { + return mu * mh * P / (rho * kb); + }; + + Tscal T_a = temperature(P_a, rho_gas_a, mu, mh, kb); + + auto dv = [&](int i, int j) { + // dv_ij = v_dust_j - v_dust_i = delta_v_j[j] - delta_v_j[i] + + Tscal ts_i = ts[i]; + Tscal ts_j = ts[j]; + + // Ormel model for differential velocities + Tscal ts_1 = 0.; + Tscal St_1 = 0.; + Tscal St_2 = 0.; + Tscal x_St = 0.; + Tscal beta_St = 0.; + Tscal res = 0.; + + Tscal eps_diff = 1e-16; + + Tscal t_L = cs_a / a_ext_a; + Tscal nh = rho_gas_a / (mu * mh); + Tscal Re = 62e6 * sycl::sqrt(nh / 1e5) * sycl::sqrt(T_a / 10.); + Tscal t_eta = t_L / sycl::sqrt(Re); + + // to symmetrize dv + if (j > i) { + ts_1 = ts_j; + St_1 = ts_j / t_L; + St_2 = ts_i / t_L; + } else { + ts_1 = ts_i; + St_1 = ts_i / t_L; + St_2 = ts_j / t_L; + } + + x_St = St_1 / St_2; + beta_St + = 3.2 - (1. + x_St) + + 2. / (1. + x_St) * (1. / 2.6 + sycl::pow(x_St, 3) / (1.6 + x_St)); + + if (ts_1 < t_eta) { + + if (sham::abs(St_1 - St_2) < eps_diff) { + res = 0.; + } else { + res = alpha_turb * sycl::pow(cs_a, 2) * (St_1 - St_2) + / (St_1 + St_2) + * (sycl::pow(St_1, 2) / (St_1 + 1. / sycl::sqrt(Re)) + + sycl::pow(St_2, 2) / (St_2 + 1. / sycl::sqrt(Re))); + } + + } else if ((t_eta <= ts_1) && (ts_1 < t_L)) { + + res = alpha_turb * sycl::pow(cs_a, 2) * beta_St * St_1; + + } else { + + res = alpha_turb * sycl::pow(cs_a, 2) + * (1. / (St_1 + 1.) + 1. / (St_2 + 1.)); + } + + return sycl::sqrt(res); + }; + + // should implement the same content as + // src/pylib/shamrock/external/coala/interface_coala_shamrock.py + + shamphys::coala_k0_source_term( + nbins, + dv, + rho_dust, + rho_eps, + massgrid, + tabflux_coag, + gij, + flux, + S_coag_span); + }); + }; + } + }; + + template class SPHKernel> + inline void NodeEvolveDustCOALASourceTermOrmel::_impl_evaluate_internal() { + + __shamrock_stack_entry(); + + auto edges = get_edges(); + + auto s_j_spans = edges.s_j.get_spans(); + auto hpart_spans = edges.hpart.get_spans(); + auto pressure_spans = edges.pressure.get_spans(); + auto cs_spans = edges.cs.get_spans(); + auto ts_spans = edges.ts.get_spans(); + auto a_ext_spans = edges.a_ext.get_spans(); + + auto counts = edges.part_counts.indexes; + + edges.S_coag.ensure_sizes(counts); + auto S_coag_spans = edges.S_coag.get_spans(); + + Tscal rho_eps = edges.rhodust_eps.value; + const std::vector &massgrid = edges.massgrid.value; + const std::vector &tensor_tabflux_coag = edges.tensor_tabflux_coag.value; + Tscal alpha_turb = edges.alpha_turb.value; + Tscal mu = edges.mu.value; + Tscal mh = edges.mh.value; + Tscal kb = edges.kb.value; + Tscal pmass = edges.gpart_mass.value; + + auto dev_sched = shamsys::instance::get_compute_scheduler_ptr(); + auto &q = shambase::get_check_ref(dev_sched).get_queue(); + + sham::DeviceBuffer massgrid_buf(nbins + 1, dev_sched); + massgrid_buf.copy_from_stdvec(massgrid); + + sham::DeviceBuffer tensor_tabflux_coag_buf(nbins * nbins * nbins, dev_sched); + tensor_tabflux_coag_buf.copy_from_stdvec(tensor_tabflux_coag); + + u32 group_size = 64; + + counts.for_each([&](u64 id_patch, u64 count) { + u32 group_cnt = shambase::group_count(count, group_size); + u32 corrected_len = group_cnt * group_size; + + sham::kernel_call_hndl( + q, + sham::MultiRef{ + massgrid_buf, + tensor_tabflux_coag_buf, + s_j_spans.get(id_patch), + hpart_spans.get(id_patch), + pressure_spans.get(id_patch), + cs_spans.get(id_patch), + ts_spans.get(id_patch), + a_ext_spans.get(id_patch)}, + sham::MultiRef{S_coag_spans.get(id_patch)}, + count, + KernelGenCoala_k0_Ormel{ + nbins, + rho_eps, + corrected_len, + group_size, + u32(count), + pmass, + mu, + mh, + kb, + alpha_turb}); + }); + } + + template class SPHKernel> + inline std::string NodeEvolveDustCOALASourceTermOrmel::_impl_get_tex() const { + + auto rhodust_eps = get_ro_edge_base(0).get_tex_symbol(); + auto massgrid = get_ro_edge_base(1).get_tex_symbol(); + auto tensor_tabflux_coag = get_ro_edge_base(2).get_tex_symbol(); + auto part_counts = get_ro_edge_base(3).get_tex_symbol(); + auto s_j = get_ro_edge_base(5).get_tex_symbol(); + auto delta_v_j = get_ro_edge_base(6).get_tex_symbol(); + auto S_coag = get_rw_edge_base(0).get_tex_symbol(); + + std::string tex = R"tex( + COALA dust coagulation source term, DG $k=0$ (Lombart et al., 2021) + + Per gas particle $a$ and mass bin $j$ (monofluid: $\rho_{{\rm d},j,a}} = {s_j}_{j,a}^2$): + + \begin{align} + \rho_{{\rm d},j,a} &= {s_j}_{j,a}^2 \\ + \Delta m_j &= {massgrid}_{j+1} - {massgrid}_j \\ + g_{j,a} &= \begin{cases} + \rho_{{\rm d},j,a} / \Delta m_j & \rho_{{\rm d},j,a} > \rho_{\rm eps} \\ + 0 & \text{otherwise} + \end{cases} \\ + \mathrm{dv}_{l,m,a} &= \left| {delta_v_j}_{m,a} - {delta_v_j}_{l,a} \right| \\ + \mathrm{flux}_{j,a} &= \sum_{l,m} + {tensor_tabflux_coag}_{j,l,m}\, + \mathrm{dv}_{l,m,a}\, g_{l,a}\, g_{m,a} \\ + {S_coag}_{0,a} &= -\mathrm{flux}_{0,a}, \quad + {S_coag}_{j,a} = \mathrm{flux}_{j-1,a} - \mathrm{flux}_{j,a} + \quad (j \ge 1) \\ + a &\in [0, {part_counts}), \quad j,l,m \in [0, N_{\rm bins}) \\ + \rho_{\rm eps} &= {rhodust_eps}, \quad N_{\rm bins} = {nbins} + \end{align} + )tex"; + + shambase::replace_all(tex, "{rhodust_eps}", rhodust_eps); + shambase::replace_all(tex, "{massgrid}", massgrid); + shambase::replace_all(tex, "{tensor_tabflux_coag}", tensor_tabflux_coag); + shambase::replace_all(tex, "{part_counts}", part_counts); + shambase::replace_all(tex, "{s_j}", s_j); + shambase::replace_all(tex, "{delta_v_j}", delta_v_j); + shambase::replace_all(tex, "{S_coag}", S_coag); + shambase::replace_all(tex, "{nbins}", shambase::format("{}", nbins)); + + return tex; + } +} // namespace shammodels::sph::modules + +using namespace shammath; +template class shammodels::sph::modules::NodeEvolveDustCOALASourceTermOrmel; +template class shammodels::sph::modules::NodeEvolveDustCOALASourceTermOrmel; +template class shammodels::sph::modules::NodeEvolveDustCOALASourceTermOrmel; + +template class shammodels::sph::modules::NodeEvolveDustCOALASourceTermOrmel; +template class shammodels::sph::modules::NodeEvolveDustCOALASourceTermOrmel; +template class shammodels::sph::modules::NodeEvolveDustCOALASourceTermOrmel;