diff --git a/src/shammodels/ramses/CMakeLists.txt b/src/shammodels/ramses/CMakeLists.txt index eca6794e83..5060c16bc4 100644 --- a/src/shammodels/ramses/CMakeLists.txt +++ b/src/shammodels/ramses/CMakeLists.txt @@ -23,6 +23,8 @@ set(Sources src/modules/ComputeCFL.cpp src/modules/ConsToPrimGas.cpp src/modules/ConsToPrimDust.cpp + src/modules/EulerTimeDerivativeGas.cpp + src/modules/EulerTimeDerivativeDust.cpp src/modules/ComputeCellAABB.cpp src/modules/NodeBuildTrees.cpp src/modules/FindBlockNeigh.cpp diff --git a/src/shammodels/ramses/include/shammodels/ramses/modules/EulerTimeDerivativeDust.hpp b/src/shammodels/ramses/include/shammodels/ramses/modules/EulerTimeDerivativeDust.hpp new file mode 100644 index 0000000000..34a85cfc6b --- /dev/null +++ b/src/shammodels/ramses/include/shammodels/ramses/modules/EulerTimeDerivativeDust.hpp @@ -0,0 +1,69 @@ +// -------------------------------------------------------// +// +// 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 EulerTimeDerivativeDust.hpp + * @author Timothée David--Cléris (tim.shamrock@proton.me) + * @brief Per cell Euler time derivatives of the dust primitive state + * + */ + +#include "shambackends/vec.hpp" +#include "shamrock/solvergraph/IFieldSpan.hpp" +#include "shamrock/solvergraph/Indexes.hpp" +#include "shamsolvergraph/node/INode.hpp" + +#define NODE_EDGES(X_RO, X_RW) \ + /* ------------------- inputs ------------------- */ \ + X_RO(shamrock::solvergraph::Indexes, sizes) \ + X_RO(shamrock::solvergraph::IFieldSpan, spans_rho_dust) \ + X_RO(shamrock::solvergraph::IFieldSpan, spans_vel_dust) \ + X_RO(shamrock::solvergraph::IFieldSpan, spans_grad_rho_dust) \ + X_RO(shamrock::solvergraph::IFieldSpan, spans_dx_v_dust) \ + X_RO(shamrock::solvergraph::IFieldSpan, spans_dy_v_dust) \ + X_RO(shamrock::solvergraph::IFieldSpan, spans_dz_v_dust) \ + \ + /* ------------------- outputs ------------------- */ \ + X_RW(shamrock::solvergraph::IFieldSpan, spans_dt_rho_dust) \ + X_RW(shamrock::solvergraph::IFieldSpan, spans_dt_vel_dust) + +namespace shammodels::basegodunov::modules { + + /** + * @brief Compute the Euler time derivatives of the dust primitive state per cell + * + * Dust is pressureless, so the velocity derivative carries only the + * advection term. Same rationale as the gas counterpart: evaluate once per + * cell instead of once per face side. + */ + template + class NodeEulerTimeDerivativeDust : public shamrock::solvergraph::INode { + using Tscal = shambase::VecComponent; + + u32 block_size; + u32 ndust; + + public: + NodeEulerTimeDerivativeDust(u32 block_size, u32 ndust) + : block_size(block_size), ndust(ndust) {} + + EXPAND_NODE_EDGES(NODE_EDGES) + + void _impl_evaluate_internal(); + + inline virtual std::string _impl_get_label() const { return "EulerTimeDerivativeDust"; }; + + virtual std::string _impl_get_tex() const; + }; + +} // namespace shammodels::basegodunov::modules + +#undef NODE_EDGES diff --git a/src/shammodels/ramses/include/shammodels/ramses/modules/EulerTimeDerivativeGas.hpp b/src/shammodels/ramses/include/shammodels/ramses/modules/EulerTimeDerivativeGas.hpp new file mode 100644 index 0000000000..3cc3273594 --- /dev/null +++ b/src/shammodels/ramses/include/shammodels/ramses/modules/EulerTimeDerivativeGas.hpp @@ -0,0 +1,73 @@ +// -------------------------------------------------------// +// +// 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 EulerTimeDerivativeGas.hpp + * @author Timothée David--Cléris (tim.shamrock@proton.me) + * @brief Per cell Euler time derivatives of the gas primitive state + * + */ + +#include "shambackends/vec.hpp" +#include "shamrock/solvergraph/IFieldSpan.hpp" +#include "shamrock/solvergraph/Indexes.hpp" +#include "shamsolvergraph/node/INode.hpp" + +#define NODE_EDGES(X_RO, X_RW) \ + /* ------------------- inputs ------------------- */ \ + X_RO(shamrock::solvergraph::Indexes, sizes) \ + X_RO(shamrock::solvergraph::IFieldSpan, spans_rho) \ + X_RO(shamrock::solvergraph::IFieldSpan, spans_vel) \ + X_RO(shamrock::solvergraph::IFieldSpan, spans_press) \ + X_RO(shamrock::solvergraph::IFieldSpan, spans_grad_rho) \ + X_RO(shamrock::solvergraph::IFieldSpan, spans_dx_v) \ + X_RO(shamrock::solvergraph::IFieldSpan, spans_dy_v) \ + X_RO(shamrock::solvergraph::IFieldSpan, spans_dz_v) \ + X_RO(shamrock::solvergraph::IFieldSpan, spans_grad_P) \ + \ + /* ------------------- outputs ------------------- */ \ + X_RW(shamrock::solvergraph::IFieldSpan, spans_dt_rho) \ + X_RW(shamrock::solvergraph::IFieldSpan, spans_dt_vel) \ + X_RW(shamrock::solvergraph::IFieldSpan, spans_dt_press) + +namespace shammodels::basegodunov::modules { + + /** + * @brief Compute the Euler time derivatives of the gas primitive state per cell + * + * Those derivatives are the predictor term of the MUSCL-Hancock face + * reconstruction. Evaluating them once per cell here rather than once per + * face side in the interpolation nodes avoids re-fetching the velocity + * gradients for every link. + */ + template + class NodeEulerTimeDerivativeGas : public shamrock::solvergraph::INode { + using Tscal = shambase::VecComponent; + + u32 block_size; + Tscal gamma; + + public: + NodeEulerTimeDerivativeGas(u32 block_size, Tscal gamma) + : block_size(block_size), gamma(gamma) {} + + EXPAND_NODE_EDGES(NODE_EDGES) + + void _impl_evaluate_internal(); + + inline virtual std::string _impl_get_label() const { return "EulerTimeDerivativeGas"; }; + + virtual std::string _impl_get_tex() const; + }; + +} // namespace shammodels::basegodunov::modules + +#undef NODE_EDGES diff --git a/src/shammodels/ramses/include/shammodels/ramses/modules/InterpolateToFace.hpp b/src/shammodels/ramses/include/shammodels/ramses/modules/InterpolateToFace.hpp index e2c1161b56..c99cf87a27 100644 --- a/src/shammodels/ramses/include/shammodels/ramses/modules/InterpolateToFace.hpp +++ b/src/shammodels/ramses/include/shammodels/ramses/modules/InterpolateToFace.hpp @@ -17,13 +17,114 @@ */ #include "shambackends/vec.hpp" -#include "shammodels/ramses/config/enum_SlopeMode.hpp" #include "shammodels/ramses/solvegraph/NeighGraphLinkFieldEdge.hpp" #include "shammodels/ramses/solvegraph/OrientedAMRGraphEdge.hpp" #include "shamrock/solvergraph/IFieldSpan.hpp" #include "shamrock/solvergraph/Indexes.hpp" #include "shamrock/solvergraph/ScalarEdge.hpp" #include "shamsolvergraph/node/INode.hpp" +#include + +// Note on the edge lists below: the Euler time derivatives (dt_rho, dt_vel, +// dt_press and their dust counterparts) used to be recomputed inside these +// kernels for both sides of every link. They are now precomputed per cell by +// NodeEulerTimeDerivativeGas / NodeEulerTimeDerivativeDust and merely loaded +// here, which is why each node only carries the fields its spatial +// reconstruction still needs. + +#define NODE_EDGES_RHO(X_RO, X_RW) \ + /* ------------------- inputs ------------------- */ \ + X_RO(ScalarEdgeScal, dt_interp) \ + X_RO(AMRGraphEdge, cell_neigh_graph) \ + X_RO(shamrock::solvergraph::IFieldSpan, spans_block_cell_sizes) \ + X_RO(shamrock::solvergraph::IFieldSpan, spans_cell0block_aabb_lower) \ + X_RO(shamrock::solvergraph::IFieldSpan, spans_rhos) \ + X_RO(shamrock::solvergraph::IFieldSpan, spans_grad_rho) \ + X_RO(shamrock::solvergraph::IFieldSpan, spans_dt_rho) \ + \ + /* ------------------- outputs ------------------- */ \ + X_RW(LinkFieldScal, rho_face_xp) \ + X_RW(LinkFieldScal, rho_face_xm) \ + X_RW(LinkFieldScal, rho_face_yp) \ + X_RW(LinkFieldScal, rho_face_ym) \ + X_RW(LinkFieldScal, rho_face_zp) \ + X_RW(LinkFieldScal, rho_face_zm) + +#define NODE_EDGES_VEL(X_RO, X_RW) \ + /* ------------------- inputs ------------------- */ \ + X_RO(ScalarEdgeScal, dt_interp) \ + X_RO(AMRGraphEdge, cell_neigh_graph) \ + X_RO(shamrock::solvergraph::IFieldSpan, spans_block_cell_sizes) \ + X_RO(shamrock::solvergraph::IFieldSpan, spans_cell0block_aabb_lower) \ + X_RO(shamrock::solvergraph::IFieldSpan, spans_vel) \ + X_RO(shamrock::solvergraph::IFieldSpan, spans_dx_vel) \ + X_RO(shamrock::solvergraph::IFieldSpan, spans_dy_vel) \ + X_RO(shamrock::solvergraph::IFieldSpan, spans_dz_vel) \ + X_RO(shamrock::solvergraph::IFieldSpan, spans_dt_vel) \ + \ + /* ------------------- outputs ------------------- */ \ + X_RW(LinkFieldVec, vel_face_xp) \ + X_RW(LinkFieldVec, vel_face_xm) \ + X_RW(LinkFieldVec, vel_face_yp) \ + X_RW(LinkFieldVec, vel_face_ym) \ + X_RW(LinkFieldVec, vel_face_zp) \ + X_RW(LinkFieldVec, vel_face_zm) + +#define NODE_EDGES_PRESS(X_RO, X_RW) \ + /* ------------------- inputs ------------------- */ \ + X_RO(ScalarEdgeScal, dt_interp) \ + X_RO(AMRGraphEdge, cell_neigh_graph) \ + X_RO(shamrock::solvergraph::IFieldSpan, spans_block_cell_sizes) \ + X_RO(shamrock::solvergraph::IFieldSpan, spans_cell0block_aabb_lower) \ + X_RO(shamrock::solvergraph::IFieldSpan, spans_press) \ + X_RO(shamrock::solvergraph::IFieldSpan, spans_grad_P) \ + X_RO(shamrock::solvergraph::IFieldSpan, spans_dt_press) \ + \ + /* ------------------- outputs ------------------- */ \ + X_RW(LinkFieldScal, press_face_xp) \ + X_RW(LinkFieldScal, press_face_xm) \ + X_RW(LinkFieldScal, press_face_yp) \ + X_RW(LinkFieldScal, press_face_ym) \ + X_RW(LinkFieldScal, press_face_zp) \ + X_RW(LinkFieldScal, press_face_zm) + +#define NODE_EDGES_RHO_DUST(X_RO, X_RW) \ + /* ------------------- inputs ------------------- */ \ + X_RO(ScalarEdgeScal, dt_interp) \ + X_RO(AMRGraphEdge, cell_neigh_graph) \ + X_RO(shamrock::solvergraph::IFieldSpan, spans_block_cell_sizes) \ + X_RO(shamrock::solvergraph::IFieldSpan, spans_cell0block_aabb_lower) \ + X_RO(shamrock::solvergraph::IFieldSpan, spans_rhos_dust) \ + X_RO(shamrock::solvergraph::IFieldSpan, spans_grad_rho_dust) \ + X_RO(shamrock::solvergraph::IFieldSpan, spans_dt_rho_dust) \ + \ + /* ------------------- outputs ------------------- */ \ + X_RW(LinkFieldScal, rho_dust_face_xp) \ + X_RW(LinkFieldScal, rho_dust_face_xm) \ + X_RW(LinkFieldScal, rho_dust_face_yp) \ + X_RW(LinkFieldScal, rho_dust_face_ym) \ + X_RW(LinkFieldScal, rho_dust_face_zp) \ + X_RW(LinkFieldScal, rho_dust_face_zm) + +#define NODE_EDGES_VEL_DUST(X_RO, X_RW) \ + /* ------------------- inputs ------------------- */ \ + X_RO(ScalarEdgeScal, dt_interp) \ + X_RO(AMRGraphEdge, cell_neigh_graph) \ + X_RO(shamrock::solvergraph::IFieldSpan, spans_block_cell_sizes) \ + X_RO(shamrock::solvergraph::IFieldSpan, spans_cell0block_aabb_lower) \ + X_RO(shamrock::solvergraph::IFieldSpan, spans_vel_dust) \ + X_RO(shamrock::solvergraph::IFieldSpan, spans_dx_vel_dust) \ + X_RO(shamrock::solvergraph::IFieldSpan, spans_dy_vel_dust) \ + X_RO(shamrock::solvergraph::IFieldSpan, spans_dz_vel_dust) \ + X_RO(shamrock::solvergraph::IFieldSpan, spans_dt_vel_dust) \ + \ + /* ------------------- outputs ------------------- */ \ + X_RW(LinkFieldVec, vel_dust_face_xp) \ + X_RW(LinkFieldVec, vel_dust_face_xm) \ + X_RW(LinkFieldVec, vel_dust_face_yp) \ + X_RW(LinkFieldVec, vel_dust_face_ym) \ + X_RW(LinkFieldVec, vel_dust_face_zp) \ + X_RW(LinkFieldVec, vel_dust_face_zm) namespace shammodels::basegodunov::modules { @@ -31,100 +132,17 @@ namespace shammodels::basegodunov::modules { class InterpolateToFaceRho : public shamrock::solvergraph::INode { using Tscal = shambase::VecComponent; - using SlopeMode = shammodels::basegodunov::SlopeMode; + /// Aliases without commas, the edge list macros cannot carry template argument lists + using ScalarEdgeScal = shamrock::solvergraph::ScalarEdge; + using AMRGraphEdge = solvergraph::OrientedAMRGraphEdge; + using LinkFieldScal = solvergraph::NeighGraphLinkFieldEdge>; u32 block_size; public: InterpolateToFaceRho(u32 block_size) : block_size(block_size) {} - struct Edges { - const shamrock::solvergraph::ScalarEdge &dt_interp; - - const solvergraph::OrientedAMRGraphEdge &cell_neigh_graph; - - const shamrock::solvergraph::IFieldSpan &spans_block_cell_sizes; - const shamrock::solvergraph::IFieldSpan &spans_cell0block_aabb_lower; - const shamrock::solvergraph::IFieldSpan &spans_rhos; - const shamrock::solvergraph::IFieldSpan &spans_grad_rho; - const shamrock::solvergraph::IFieldSpan &spans_vel; - const shamrock::solvergraph::IFieldSpan &spans_dx_vel; - const shamrock::solvergraph::IFieldSpan &spans_dy_vel; - const shamrock::solvergraph::IFieldSpan &spans_dz_vel; - - solvergraph::NeighGraphLinkFieldEdge> &rho_face_xp; - solvergraph::NeighGraphLinkFieldEdge> &rho_face_xm; - solvergraph::NeighGraphLinkFieldEdge> &rho_face_yp; - solvergraph::NeighGraphLinkFieldEdge> &rho_face_ym; - solvergraph::NeighGraphLinkFieldEdge> &rho_face_zp; - solvergraph::NeighGraphLinkFieldEdge> &rho_face_zm; - }; - - inline void set_edges( - std::shared_ptr> dt_interp, - std::shared_ptr> cell_neigh_graph, - std::shared_ptr> spans_block_cell_sizes, - std::shared_ptr> spans_cell0block_aabb_lower, - std::shared_ptr> spans_rhos, - std::shared_ptr> spans_grad_rho, - std::shared_ptr> spans_vel, - std::shared_ptr> spans_dx_vel, - std::shared_ptr> spans_dy_vel, - std::shared_ptr> spans_dz_vel, - std::shared_ptr>> - &rho_face_xp, - std::shared_ptr>> - &rho_face_xm, - std::shared_ptr>> - &rho_face_yp, - std::shared_ptr>> - &rho_face_ym, - std::shared_ptr>> - &rho_face_zp, - std::shared_ptr>> - &rho_face_zm) { - __internal_set_ro_edges({ - dt_interp, - cell_neigh_graph, - spans_block_cell_sizes, - spans_cell0block_aabb_lower, - spans_rhos, - spans_grad_rho, - spans_vel, - spans_dx_vel, - spans_dy_vel, - spans_dz_vel, - }); - __internal_set_rw_edges({ - rho_face_xp, - rho_face_xm, - rho_face_yp, - rho_face_ym, - rho_face_zp, - rho_face_zm, - }); - } - - inline Edges get_edges() { - return Edges{ - get_ro_edge>(0), - get_ro_edge>(1), - get_ro_edge>(2), - get_ro_edge>(3), - get_ro_edge>(4), - get_ro_edge>(5), - get_ro_edge>(6), - get_ro_edge>(7), - get_ro_edge>(8), - get_ro_edge>(9), - get_rw_edge>>(0), - get_rw_edge>>(1), - get_rw_edge>>(2), - get_rw_edge>>(3), - get_rw_edge>>(4), - get_rw_edge>>(5), - }; - } + EXPAND_NODE_EDGES(NODE_EDGES_RHO) void _impl_evaluate_internal(); @@ -137,96 +155,17 @@ namespace shammodels::basegodunov::modules { class InterpolateToFaceVel : public shamrock::solvergraph::INode { using Tscal = shambase::VecComponent; - using SlopeMode = shammodels::basegodunov::SlopeMode; + /// Aliases without commas, the edge list macros cannot carry template argument lists + using ScalarEdgeScal = shamrock::solvergraph::ScalarEdge; + using AMRGraphEdge = solvergraph::OrientedAMRGraphEdge; + using LinkFieldVec = solvergraph::NeighGraphLinkFieldEdge>; u32 block_size; public: InterpolateToFaceVel(u32 block_size) : block_size(block_size) {} - struct Edges { - const shamrock::solvergraph::ScalarEdge &dt_interp; - - const solvergraph::OrientedAMRGraphEdge &cell_neigh_graph; - - const shamrock::solvergraph::IFieldSpan &spans_block_cell_sizes; - const shamrock::solvergraph::IFieldSpan &spans_cell0block_aabb_lower; - const shamrock::solvergraph::IFieldSpan &spans_rhos; - const shamrock::solvergraph::IFieldSpan &spans_grad_P; - const shamrock::solvergraph::IFieldSpan &spans_vel; - const shamrock::solvergraph::IFieldSpan &spans_dx_vel; - const shamrock::solvergraph::IFieldSpan &spans_dy_vel; - const shamrock::solvergraph::IFieldSpan &spans_dz_vel; - - solvergraph::NeighGraphLinkFieldEdge> &vel_face_xp; - solvergraph::NeighGraphLinkFieldEdge> &vel_face_xm; - solvergraph::NeighGraphLinkFieldEdge> &vel_face_yp; - solvergraph::NeighGraphLinkFieldEdge> &vel_face_ym; - solvergraph::NeighGraphLinkFieldEdge> &vel_face_zp; - solvergraph::NeighGraphLinkFieldEdge> &vel_face_zm; - }; - - inline void set_edges( - std::shared_ptr> dt_interp, - std::shared_ptr> cell_neigh_graph, - std::shared_ptr> spans_block_cell_sizes, - std::shared_ptr> spans_cell0block_aabb_lower, - std::shared_ptr> spans_rhos, - std::shared_ptr> spans_grad_P, - std::shared_ptr> spans_vel, - std::shared_ptr> spans_dx_vel, - std::shared_ptr> spans_dy_vel, - std::shared_ptr> spans_dz_vel, - - std::shared_ptr>> &vel_face_xp, - std::shared_ptr>> &vel_face_xm, - std::shared_ptr>> &vel_face_yp, - std::shared_ptr>> &vel_face_ym, - std::shared_ptr>> &vel_face_zp, - std::shared_ptr>> - &vel_face_zm) { - __internal_set_ro_edges({ - dt_interp, - cell_neigh_graph, - spans_block_cell_sizes, - spans_cell0block_aabb_lower, - spans_rhos, - spans_grad_P, - spans_vel, - spans_dx_vel, - spans_dy_vel, - spans_dz_vel, - }); - __internal_set_rw_edges({ - vel_face_xp, - vel_face_xm, - vel_face_yp, - vel_face_ym, - vel_face_zp, - vel_face_zm, - }); - } - - inline Edges get_edges() { - return Edges{ - get_ro_edge>(0), - get_ro_edge>(1), - get_ro_edge>(2), - get_ro_edge>(3), - get_ro_edge>(4), - get_ro_edge>(5), - get_ro_edge>(6), - get_ro_edge>(7), - get_ro_edge>(8), - get_ro_edge>(9), - get_rw_edge>>(0), - get_rw_edge>>(1), - get_rw_edge>>(2), - get_rw_edge>>(3), - get_rw_edge>>(4), - get_rw_edge>>(5), - }; - } + EXPAND_NODE_EDGES(NODE_EDGES_VEL) void _impl_evaluate_internal(); @@ -239,103 +178,17 @@ namespace shammodels::basegodunov::modules { class InterpolateToFacePress : public shamrock::solvergraph::INode { using Tscal = shambase::VecComponent; - using SlopeMode = shammodels::basegodunov::SlopeMode; + /// Aliases without commas, the edge list macros cannot carry template argument lists + using ScalarEdgeScal = shamrock::solvergraph::ScalarEdge; + using AMRGraphEdge = solvergraph::OrientedAMRGraphEdge; + using LinkFieldScal = solvergraph::NeighGraphLinkFieldEdge>; u32 block_size; - Tscal gamma; public: - InterpolateToFacePress(u32 block_size, Tscal gamma) - : block_size(block_size), gamma(gamma) {} - - struct Edges { - const shamrock::solvergraph::ScalarEdge &dt_interp; - - const solvergraph::OrientedAMRGraphEdge &cell_neigh_graph; - - const shamrock::solvergraph::IFieldSpan &spans_block_cell_sizes; - const shamrock::solvergraph::IFieldSpan &spans_cell0block_aabb_lower; - const shamrock::solvergraph::IFieldSpan &spans_press; - const shamrock::solvergraph::IFieldSpan &spans_grad_P; - const shamrock::solvergraph::IFieldSpan &spans_vel; - const shamrock::solvergraph::IFieldSpan &spans_dx_vel; - const shamrock::solvergraph::IFieldSpan &spans_dy_vel; - const shamrock::solvergraph::IFieldSpan &spans_dz_vel; - - solvergraph::NeighGraphLinkFieldEdge> &press_face_xp; - solvergraph::NeighGraphLinkFieldEdge> &press_face_xm; - solvergraph::NeighGraphLinkFieldEdge> &press_face_yp; - solvergraph::NeighGraphLinkFieldEdge> &press_face_ym; - solvergraph::NeighGraphLinkFieldEdge> &press_face_zp; - solvergraph::NeighGraphLinkFieldEdge> &press_face_zm; - }; + InterpolateToFacePress(u32 block_size) : block_size(block_size) {} - inline void set_edges( - std::shared_ptr> dt_interp, - std::shared_ptr> cell_neigh_graph, - std::shared_ptr> spans_block_cell_sizes, - std::shared_ptr> spans_cell0block_aabb_lower, - std::shared_ptr> spans_press, - std::shared_ptr> spans_grad_P, - std::shared_ptr> spans_vel, - std::shared_ptr> spans_dx_vel, - std::shared_ptr> spans_dy_vel, - std::shared_ptr> spans_dz_vel, - - std::shared_ptr>> - &press_face_xp, - std::shared_ptr>> - &press_face_xm, - std::shared_ptr>> - &press_face_yp, - std::shared_ptr>> - &press_face_ym, - std::shared_ptr>> - &press_face_zp, - std::shared_ptr>> - &press_face_zm) { - __internal_set_ro_edges({ - dt_interp, - cell_neigh_graph, - spans_block_cell_sizes, - spans_cell0block_aabb_lower, - spans_press, - spans_grad_P, - spans_vel, - spans_dx_vel, - spans_dy_vel, - spans_dz_vel, - }); - __internal_set_rw_edges({ - press_face_xp, - press_face_xm, - press_face_yp, - press_face_ym, - press_face_zp, - press_face_zm, - }); - } - - inline Edges get_edges() { - return Edges{ - get_ro_edge>(0), - get_ro_edge>(1), - get_ro_edge>(2), - get_ro_edge>(3), - get_ro_edge>(4), - get_ro_edge>(5), - get_ro_edge>(6), - get_ro_edge>(7), - get_ro_edge>(8), - get_ro_edge>(9), - get_rw_edge>>(0), - get_rw_edge>>(1), - get_rw_edge>>(2), - get_rw_edge>>(3), - get_rw_edge>>(4), - get_rw_edge>>(5), - }; - } + EXPAND_NODE_EDGES(NODE_EDGES_PRESS) void _impl_evaluate_internal(); @@ -350,7 +203,10 @@ namespace shammodels::basegodunov::modules { class InterpolateToFaceRhoDust : public shamrock::solvergraph::INode { using Tscal = shambase::VecComponent; - using SlopeMode = shammodels::basegodunov::SlopeMode; + /// Aliases without commas, the edge list macros cannot carry template argument lists + using ScalarEdgeScal = shamrock::solvergraph::ScalarEdge; + using AMRGraphEdge = solvergraph::OrientedAMRGraphEdge; + using LinkFieldScal = solvergraph::NeighGraphLinkFieldEdge>; u32 block_size; u32 ndust; @@ -359,97 +215,13 @@ namespace shammodels::basegodunov::modules { InterpolateToFaceRhoDust(u32 block_size, u32 ndust) : block_size(block_size), ndust(ndust) {} - struct Edges { - const shamrock::solvergraph::ScalarEdge &dt_interp; - - const solvergraph::OrientedAMRGraphEdge &cell_neigh_graph; - - const shamrock::solvergraph::IFieldSpan &spans_block_cell_sizes; - const shamrock::solvergraph::IFieldSpan &spans_cell0block_aabb_lower; - const shamrock::solvergraph::IFieldSpan &spans_rhos_dust; - const shamrock::solvergraph::IFieldSpan &spans_grad_rho_dust; - const shamrock::solvergraph::IFieldSpan &spans_vel_dust; - const shamrock::solvergraph::IFieldSpan &spans_dx_vel_dust; - const shamrock::solvergraph::IFieldSpan &spans_dy_vel_dust; - const shamrock::solvergraph::IFieldSpan &spans_dz_vel_dust; - - solvergraph::NeighGraphLinkFieldEdge> &rho_dust_face_xp; - solvergraph::NeighGraphLinkFieldEdge> &rho_dust_face_xm; - solvergraph::NeighGraphLinkFieldEdge> &rho_dust_face_yp; - solvergraph::NeighGraphLinkFieldEdge> &rho_dust_face_ym; - solvergraph::NeighGraphLinkFieldEdge> &rho_dust_face_zp; - solvergraph::NeighGraphLinkFieldEdge> &rho_dust_face_zm; - }; - - inline void set_edges( - std::shared_ptr> dt_interp, - std::shared_ptr> cell_neigh_graph, - std::shared_ptr> spans_block_cell_sizes, - std::shared_ptr> spans_cell0block_aabb_lower, - std::shared_ptr> spans_rhos_dust, - std::shared_ptr> spans_grad_rho_dust, - std::shared_ptr> spans_vel_dust, - std::shared_ptr> spans_dx_vel_dust, - std::shared_ptr> spans_dy_vel_dust, - std::shared_ptr> spans_dz_vel_dust, - std::shared_ptr>> - &rho_dust_face_xp, - std::shared_ptr>> - &rho_dust_face_xm, - std::shared_ptr>> - &rho_dust_face_yp, - std::shared_ptr>> - &rho_dust_face_ym, - std::shared_ptr>> - &rho_dust_face_zp, - std::shared_ptr>> - &rho_dust_face_zm) { - __internal_set_ro_edges({ - dt_interp, - cell_neigh_graph, - spans_block_cell_sizes, - spans_cell0block_aabb_lower, - spans_rhos_dust, - spans_grad_rho_dust, - spans_vel_dust, - spans_dx_vel_dust, - spans_dy_vel_dust, - spans_dz_vel_dust, - }); - __internal_set_rw_edges({ - rho_dust_face_xp, - rho_dust_face_xm, - rho_dust_face_yp, - rho_dust_face_ym, - rho_dust_face_zp, - rho_dust_face_zm, - }); - } - - inline Edges get_edges() { - return Edges{ - get_ro_edge>(0), - get_ro_edge>(1), - get_ro_edge>(2), - get_ro_edge>(3), - get_ro_edge>(4), - get_ro_edge>(5), - get_ro_edge>(6), - get_ro_edge>(7), - get_ro_edge>(8), - get_ro_edge>(9), - get_rw_edge>>(0), - get_rw_edge>>(1), - get_rw_edge>>(2), - get_rw_edge>>(3), - get_rw_edge>>(4), - get_rw_edge>>(5), - }; - } + EXPAND_NODE_EDGES(NODE_EDGES_RHO_DUST) void _impl_evaluate_internal(); - inline virtual std::string _impl_get_label() const { return "InterpolateRhoToFaceRho"; }; + inline virtual std::string _impl_get_label() const { + return "InterpolateRhoDustToFaceRhoDust"; + }; virtual std::string _impl_get_tex() const; }; @@ -458,7 +230,10 @@ namespace shammodels::basegodunov::modules { class InterpolateToFaceVelDust : public shamrock::solvergraph::INode { using Tscal = shambase::VecComponent; - using SlopeMode = shammodels::basegodunov::SlopeMode; + /// Aliases without commas, the edge list macros cannot carry template argument lists + using ScalarEdgeScal = shamrock::solvergraph::ScalarEdge; + using AMRGraphEdge = solvergraph::OrientedAMRGraphEdge; + using LinkFieldVec = solvergraph::NeighGraphLinkFieldEdge>; u32 block_size; u32 ndust; @@ -467,96 +242,21 @@ namespace shammodels::basegodunov::modules { InterpolateToFaceVelDust(u32 block_size, u32 ndust) : block_size(block_size), ndust(ndust) {} - struct Edges { - const shamrock::solvergraph::ScalarEdge &dt_interp; - - const solvergraph::OrientedAMRGraphEdge &cell_neigh_graph; - - const shamrock::solvergraph::IFieldSpan &spans_block_cell_sizes; - const shamrock::solvergraph::IFieldSpan &spans_cell0block_aabb_lower; - const shamrock::solvergraph::IFieldSpan &spans_rhos_dust; - const shamrock::solvergraph::IFieldSpan &spans_vel_dust; - const shamrock::solvergraph::IFieldSpan &spans_dx_vel_dust; - const shamrock::solvergraph::IFieldSpan &spans_dy_vel_dust; - const shamrock::solvergraph::IFieldSpan &spans_dz_vel_dust; - - solvergraph::NeighGraphLinkFieldEdge> &vel_dust_face_xp; - solvergraph::NeighGraphLinkFieldEdge> &vel_dust_face_xm; - solvergraph::NeighGraphLinkFieldEdge> &vel_dust_face_yp; - solvergraph::NeighGraphLinkFieldEdge> &vel_dust_face_ym; - solvergraph::NeighGraphLinkFieldEdge> &vel_dust_face_zp; - solvergraph::NeighGraphLinkFieldEdge> &vel_dust_face_zm; - }; - - inline void set_edges( - std::shared_ptr> dt_interp, - std::shared_ptr> cell_neigh_graph, - std::shared_ptr> spans_block_cell_sizes, - std::shared_ptr> spans_cell0block_aabb_lower, - std::shared_ptr> spans_rhos_dust, - std::shared_ptr> spans_vel_dust, - std::shared_ptr> spans_dx_vel_dust, - std::shared_ptr> spans_dy_vel_dust, - std::shared_ptr> spans_dz_vel_dust, - - std::shared_ptr>> - &vel_dust_face_xp, - std::shared_ptr>> - &vel_dust_face_xm, - std::shared_ptr>> - &vel_dust_face_yp, - std::shared_ptr>> - &vel_dust_face_ym, - std::shared_ptr>> - &vel_dust_face_zp, - std::shared_ptr>> - &vel_dust_face_zm) { - __internal_set_ro_edges({ - dt_interp, - cell_neigh_graph, - spans_block_cell_sizes, - spans_cell0block_aabb_lower, - spans_rhos_dust, - spans_vel_dust, - spans_dx_vel_dust, - spans_dy_vel_dust, - spans_dz_vel_dust, - }); - __internal_set_rw_edges({ - vel_dust_face_xp, - vel_dust_face_xm, - vel_dust_face_yp, - vel_dust_face_ym, - vel_dust_face_zp, - vel_dust_face_zm, - }); - } - - inline Edges get_edges() { - return Edges{ - get_ro_edge>(0), - get_ro_edge>(1), - get_ro_edge>(2), - get_ro_edge>(3), - get_ro_edge>(4), - get_ro_edge>(5), - get_ro_edge>(6), - get_ro_edge>(7), - get_ro_edge>(8), - get_rw_edge>>(0), - get_rw_edge>>(1), - get_rw_edge>>(2), - get_rw_edge>>(3), - get_rw_edge>>(4), - get_rw_edge>>(5), - }; - } + EXPAND_NODE_EDGES(NODE_EDGES_VEL_DUST) void _impl_evaluate_internal(); - inline virtual std::string _impl_get_label() const { return "InterpolateVelToFaceVel"; }; + inline virtual std::string _impl_get_label() const { + return "InterpolateVelDustToFaceVelDust"; + }; virtual std::string _impl_get_tex() const; }; } // namespace shammodels::basegodunov::modules + +#undef NODE_EDGES_RHO +#undef NODE_EDGES_VEL +#undef NODE_EDGES_PRESS +#undef NODE_EDGES_RHO_DUST +#undef NODE_EDGES_VEL_DUST diff --git a/src/shammodels/ramses/include/shammodels/ramses/modules/SolverStorage.hpp b/src/shammodels/ramses/include/shammodels/ramses/modules/SolverStorage.hpp index 96103dccb2..1bfe4cc056 100644 --- a/src/shammodels/ramses/include/shammodels/ramses/modules/SolverStorage.hpp +++ b/src/shammodels/ramses/include/shammodels/ramses/modules/SolverStorage.hpp @@ -118,6 +118,17 @@ namespace shammodels::basegodunov { /// dust fields gradients (d vdust / d z) std::shared_ptr> dz_v_dust; + /// Euler time derivative of the gas primitive density (predictor term) + std::shared_ptr> euler_dt_rho; + /// Euler time derivative of the gas primitive velocity (predictor term) + std::shared_ptr> euler_dt_vel; + /// Euler time derivative of the gas primitive pressure (predictor term) + std::shared_ptr> euler_dt_press; + /// Euler time derivative of the dust primitive density (predictor term) + std::shared_ptr> euler_dt_rho_dust; + /// Euler time derivative of the dust primitive velocity (predictor term) + std::shared_ptr> euler_dt_vel_dust; + std::shared_ptr> rho_mean; std::shared_ptr> simulation_volume; std::shared_ptr> cell_mass; diff --git a/src/shammodels/ramses/src/Solver.cpp b/src/shammodels/ramses/src/Solver.cpp index 44a1445f3c..87b04f7bb9 100644 --- a/src/shammodels/ramses/src/Solver.cpp +++ b/src/shammodels/ramses/src/Solver.cpp @@ -35,6 +35,8 @@ #include "shammodels/ramses/modules/ConsToPrimDust.hpp" #include "shammodels/ramses/modules/ConsToPrimGas.hpp" #include "shammodels/ramses/modules/DragIntegrator.hpp" +#include "shammodels/ramses/modules/EulerTimeDerivativeDust.hpp" +#include "shammodels/ramses/modules/EulerTimeDerivativeGas.hpp" #include "shammodels/ramses/modules/ExtractGhostLayer.hpp" #include "shammodels/ramses/modules/FindBlockNeigh.hpp" #include "shammodels/ramses/modules/FindGhostLayerIndices.hpp" @@ -464,6 +466,14 @@ void shammodels::basegodunov::Solver::init_solver_graph() { storage.grad_P = std::make_shared>( AMRBlock::block_size, "grad_P", "\\nabla P"); + // will be filled by NodeEulerTimeDerivativeGas + storage.euler_dt_rho = std::make_shared>( + AMRBlock::block_size, "euler_dt_rho", "\\partial_t \\rho"); + storage.euler_dt_vel = std::make_shared>( + AMRBlock::block_size, "euler_dt_vel", "\\partial_t \\mathbf{v}"); + storage.euler_dt_press = std::make_shared>( + AMRBlock::block_size, "euler_dt_press", "\\partial_t P"); + if (solver_config.is_dust_on()) { u32 ndust = solver_config.dust_config.ndust; storage.grad_rho_dust = std::make_shared>( @@ -474,6 +484,14 @@ void shammodels::basegodunov::Solver::init_solver_graph() { AMRBlock::block_size * ndust, "dy_v_dust", "\\nabla_y \\mathbf{v}_{\\rm dust}"); storage.dz_v_dust = std::make_shared>( AMRBlock::block_size * ndust, "dz_v_dust", "\\nabla_z \\mathbf{v}_{\\rm dust}"); + + // will be filled by NodeEulerTimeDerivativeDust + storage.euler_dt_rho_dust = std::make_shared>( + AMRBlock::block_size * ndust, "euler_dt_rho_dust", "\\partial_t \\rho_{\\rm dust}"); + storage.euler_dt_vel_dust = std::make_shared>( + AMRBlock::block_size * ndust, + "euler_dt_vel_dust", + "\\partial_t \\mathbf{v}_{\\rm dust}"); } { @@ -1123,6 +1141,52 @@ void shammodels::basegodunov::Solver::init_solver_graph() { solver_sequence.push_back(std::make_shared(std::move(seq))); } + { // Euler time derivatives of the primitive state + // Hoisted out of the face interpolation nodes: they only depend on cell + // local quantities, so computing them once per cell here avoids + // re-fetching the velocity gradients for every face link. + std::vector> dt_prim_sequence; + + { + modules::NodeEulerTimeDerivativeGas node{ + AMRBlock::block_size, solver_config.eos_gamma}; + node.set_edges( + storage.block_counts_with_ghost, + storage.refs_rho, + storage.vel, + storage.press, + storage.grad_rho, + storage.dx_v, + storage.dy_v, + storage.dz_v, + storage.grad_P, + storage.euler_dt_rho, + storage.euler_dt_vel, + storage.euler_dt_press); + dt_prim_sequence.push_back(std::make_shared(std::move(node))); + } + + if (solver_config.is_dust_on()) { + u32 ndust = solver_config.dust_config.ndust; + modules::NodeEulerTimeDerivativeDust node{AMRBlock::block_size, ndust}; + node.set_edges( + storage.block_counts_with_ghost, + storage.refs_rho_dust, + storage.vel_dust, + storage.grad_rho_dust, + storage.dx_v_dust, + storage.dy_v_dust, + storage.dz_v_dust, + storage.euler_dt_rho_dust, + storage.euler_dt_vel_dust); + dt_prim_sequence.push_back(std::make_shared(std::move(node))); + } + + shamrock::solvergraph::OperationSequence seq( + "Euler time derivatives", std::move(dt_prim_sequence)); + solver_sequence.push_back(std::make_shared(std::move(seq))); + } + { // interpolate to face std::vector> interp_sequence; { @@ -1134,10 +1198,7 @@ void shammodels::basegodunov::Solver::init_solver_graph() { storage.cell0block_aabb_lower, storage.refs_rho, storage.grad_rho, - storage.vel, - storage.dx_v, - storage.dy_v, - storage.dz_v, + storage.euler_dt_rho, storage.rho_face_xp, storage.rho_face_xm, storage.rho_face_yp, @@ -1154,12 +1215,11 @@ void shammodels::basegodunov::Solver::init_solver_graph() { storage.cell_graph_edge, storage.block_cell_sizes, storage.cell0block_aabb_lower, - storage.refs_rho, - storage.grad_P, storage.vel, storage.dx_v, storage.dy_v, storage.dz_v, + storage.euler_dt_vel, storage.vel_face_xp, storage.vel_face_xm, storage.vel_face_yp, @@ -1170,8 +1230,7 @@ void shammodels::basegodunov::Solver::init_solver_graph() { } { - modules::InterpolateToFacePress node{ - AMRBlock::block_size, solver_config.eos_gamma}; + modules::InterpolateToFacePress node{AMRBlock::block_size}; node.set_edges( storage.dt_over2, storage.cell_graph_edge, @@ -1179,10 +1238,7 @@ void shammodels::basegodunov::Solver::init_solver_graph() { storage.cell0block_aabb_lower, storage.press, storage.grad_P, - storage.vel, - storage.dx_v, - storage.dy_v, - storage.dz_v, + storage.euler_dt_press, storage.press_face_xp, storage.press_face_xm, storage.press_face_yp, @@ -1202,10 +1258,7 @@ void shammodels::basegodunov::Solver::init_solver_graph() { storage.cell0block_aabb_lower, storage.refs_rho_dust, storage.grad_rho_dust, - storage.vel_dust, - storage.dx_v_dust, - storage.dy_v_dust, - storage.dz_v_dust, + storage.euler_dt_rho_dust, storage.rho_dust_face_xp, storage.rho_dust_face_xm, storage.rho_dust_face_yp, @@ -1223,11 +1276,11 @@ void shammodels::basegodunov::Solver::init_solver_graph() { storage.cell_graph_edge, storage.block_cell_sizes, storage.cell0block_aabb_lower, - storage.refs_rho_dust, storage.vel_dust, storage.dx_v_dust, storage.dy_v_dust, storage.dz_v_dust, + storage.euler_dt_vel_dust, storage.vel_dust_face_xp, storage.vel_dust_face_xm, storage.vel_dust_face_yp, diff --git a/src/shammodels/ramses/src/modules/EulerTimeDerivativeDust.cpp b/src/shammodels/ramses/src/modules/EulerTimeDerivativeDust.cpp new file mode 100644 index 0000000000..78f404ff41 --- /dev/null +++ b/src/shammodels/ramses/src/modules/EulerTimeDerivativeDust.cpp @@ -0,0 +1,173 @@ +// -------------------------------------------------------// +// +// 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 EulerTimeDerivativeDust.cpp + * @author Timothée David--Cléris (tim.shamrock@proton.me) + * @brief + * + */ + +#include "shambase/string.hpp" +#include "shambackends/kernel_call_distrib.hpp" +#include "shambackends/math.hpp" +#include "shammodels/ramses/modules/EulerTimeDerivativeDust.hpp" +#include "shamrock/patch/PatchDataField.hpp" +#include "shamsys/NodeInstance.hpp" + +namespace { + + template + struct KernelEulerTimeDerivativeDust { + using Tscal = shambase::VecComponent; + + inline static void kernel( + const shambase::DistributedData> + &spans_rho_dust, + const shambase::DistributedData> + &spans_vel_dust, + const shambase::DistributedData> + &spans_grad_rho_dust, + const shambase::DistributedData> + &spans_dx_v_dust, + const shambase::DistributedData> + &spans_dy_v_dust, + const shambase::DistributedData> + &spans_dz_v_dust, + + shambase::DistributedData> + &spans_dt_rho_dust, + shambase::DistributedData> &spans_dt_vel_dust, + const shambase::DistributedData &sizes, + u32 block_size, + u32 ndust) { + + shambase::DistributedData cell_counts + = sizes.map([&](u64 id, u32 block_count) { + u32 cell_count = block_count * block_size * ndust; + return cell_count; + }); + + sham::distributed_data_kernel_call( + shamsys::instance::get_compute_scheduler_ptr(), + sham::DDMultiRef{ + spans_rho_dust, + spans_vel_dust, + spans_grad_rho_dust, + spans_dx_v_dust, + spans_dy_v_dust, + spans_dz_v_dust}, + sham::DDMultiRef{spans_dt_rho_dust, spans_dt_vel_dust}, + cell_counts, + [](u32 i, + const Tscal *__restrict rho_dust, + const Tvec *__restrict vel_dust, + const Tvec *__restrict grad_rho_dust, + const Tvec *__restrict dx_v_dust, + const Tvec *__restrict dy_v_dust, + const Tvec *__restrict dz_v_dust, + Tscal *__restrict dt_rho_dust, + Tvec *__restrict dt_vel_dust) { + Tscal rho_dust_i = rho_dust[i]; + Tvec v_dust_i = vel_dust[i]; + Tvec grad_rho_dust_i = grad_rho_dust[i]; + Tvec dx_v_dust_i = dx_v_dust[i]; + Tvec dy_v_dust_i = dy_v_dust[i]; + Tvec dz_v_dust_i = dz_v_dust[i]; + + dt_rho_dust[i] + = -(sham::dot(v_dust_i, grad_rho_dust_i) + + rho_dust_i * (dx_v_dust_i[0] + dy_v_dust_i[1] + dz_v_dust_i[2])); + + // Dust is pressureless, only the advection term remains + dt_vel_dust[i] + = -(v_dust_i[0] * dx_v_dust_i + v_dust_i[1] * dy_v_dust_i + + v_dust_i[2] * dz_v_dust_i); + }); + } + }; + +} // namespace + +namespace shammodels::basegodunov::modules { + + template + void NodeEulerTimeDerivativeDust::_impl_evaluate_internal() { + auto edges = get_edges(); + + edges.spans_rho_dust.check_sizes(edges.sizes.indexes); + edges.spans_vel_dust.check_sizes(edges.sizes.indexes); + edges.spans_grad_rho_dust.check_sizes(edges.sizes.indexes); + edges.spans_dx_v_dust.check_sizes(edges.sizes.indexes); + edges.spans_dy_v_dust.check_sizes(edges.sizes.indexes); + edges.spans_dz_v_dust.check_sizes(edges.sizes.indexes); + + edges.spans_dt_rho_dust.ensure_sizes(edges.sizes.indexes); + edges.spans_dt_vel_dust.ensure_sizes(edges.sizes.indexes); + + KernelEulerTimeDerivativeDust::kernel( + edges.spans_rho_dust.get_spans(), + edges.spans_vel_dust.get_spans(), + edges.spans_grad_rho_dust.get_spans(), + edges.spans_dx_v_dust.get_spans(), + edges.spans_dy_v_dust.get_spans(), + edges.spans_dz_v_dust.get_spans(), + edges.spans_dt_rho_dust.get_spans(), + edges.spans_dt_vel_dust.get_spans(), + edges.sizes.indexes, + block_size, + ndust); + } + + template + std::string NodeEulerTimeDerivativeDust::_impl_get_tex() const { + + auto block_count = get_ro_edge_base(0).get_tex_symbol(); + auto rho_dust = get_ro_edge_base(1).get_tex_symbol(); + auto vel_dust = get_ro_edge_base(2).get_tex_symbol(); + auto grad_rho_dust = get_ro_edge_base(3).get_tex_symbol(); + auto dx_v_dust = get_ro_edge_base(4).get_tex_symbol(); + auto dy_v_dust = get_ro_edge_base(5).get_tex_symbol(); + auto dz_v_dust = get_ro_edge_base(6).get_tex_symbol(); + auto dt_rho_dust = get_rw_edge_base(0).get_tex_symbol(); + auto dt_vel_dust = get_rw_edge_base(1).get_tex_symbol(); + + std::string tex = R"tex( + Euler time derivatives of the dust primitive state (pressureless) + + \begin{align} + {dt_rho_dust}_i &= - \left( {vel_dust}_i \cdot {grad_rho_dust}_i + + {rho_dust}_i \left( {dx_v_dust}_{i,x} + {dy_v_dust}_{i,y} + + {dz_v_dust}_{i,z} \right) \right) \\ + {dt_vel_dust}_i &= - \left( {vel_dust}_{i,x} {dx_v_dust}_i + + {vel_dust}_{i,y} {dy_v_dust}_i + {vel_dust}_{i,z} {dz_v_dust}_i \right) \\ + i &\in [0,{block_count} * N_{\rm cell/block} * N_{\rm dust}) \\ + N_{\rm cell/block} & = {block_size} \\ + N_{\rm dust} & = {ndust} + \end{align} + )tex"; + + shambase::replace_all(tex, "{dt_rho_dust}", dt_rho_dust); + shambase::replace_all(tex, "{dt_vel_dust}", dt_vel_dust); + shambase::replace_all(tex, "{grad_rho_dust}", grad_rho_dust); + shambase::replace_all(tex, "{rho_dust}", rho_dust); + shambase::replace_all(tex, "{vel_dust}", vel_dust); + shambase::replace_all(tex, "{dx_v_dust}", dx_v_dust); + shambase::replace_all(tex, "{dy_v_dust}", dy_v_dust); + shambase::replace_all(tex, "{dz_v_dust}", dz_v_dust); + shambase::replace_all(tex, "{block_count}", block_count); + shambase::replace_all(tex, "{block_size}", sham::format("{}", block_size)); + shambase::replace_all(tex, "{ndust}", sham::format("{}", ndust)); + + return tex; + } + +} // namespace shammodels::basegodunov::modules + +template class shammodels::basegodunov::modules::NodeEulerTimeDerivativeDust; diff --git a/src/shammodels/ramses/src/modules/EulerTimeDerivativeGas.cpp b/src/shammodels/ramses/src/modules/EulerTimeDerivativeGas.cpp new file mode 100644 index 0000000000..6b95ea29c3 --- /dev/null +++ b/src/shammodels/ramses/src/modules/EulerTimeDerivativeGas.cpp @@ -0,0 +1,195 @@ +// -------------------------------------------------------// +// +// 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 EulerTimeDerivativeGas.cpp + * @author Timothée David--Cléris (tim.shamrock@proton.me) + * @brief + * + */ + +#include "shambase/string.hpp" +#include "shambackends/kernel_call_distrib.hpp" +#include "shambackends/math.hpp" +#include "shammodels/ramses/modules/EulerTimeDerivativeGas.hpp" +#include "shamrock/patch/PatchDataField.hpp" +#include "shamsys/NodeInstance.hpp" + +namespace { + + template + struct KernelEulerTimeDerivativeGas { + using Tscal = shambase::VecComponent; + + inline static void kernel( + const shambase::DistributedData> &spans_rho, + const shambase::DistributedData> &spans_vel, + const shambase::DistributedData> + &spans_press, + const shambase::DistributedData> + &spans_grad_rho, + const shambase::DistributedData> &spans_dx_v, + const shambase::DistributedData> &spans_dy_v, + const shambase::DistributedData> &spans_dz_v, + const shambase::DistributedData> + &spans_grad_P, + + shambase::DistributedData> &spans_dt_rho, + shambase::DistributedData> &spans_dt_vel, + shambase::DistributedData> &spans_dt_press, + const shambase::DistributedData &sizes, + u32 block_size, + Tscal gamma) { + + shambase::DistributedData cell_counts + = sizes.map([&](u64 id, u32 block_count) { + u32 cell_count = block_count * block_size; + return cell_count; + }); + + sham::distributed_data_kernel_call( + shamsys::instance::get_compute_scheduler_ptr(), + sham::DDMultiRef{ + spans_rho, + spans_vel, + spans_press, + spans_grad_rho, + spans_dx_v, + spans_dy_v, + spans_dz_v, + spans_grad_P}, + sham::DDMultiRef{spans_dt_rho, spans_dt_vel, spans_dt_press}, + cell_counts, + [gamma]( + u32 i, + const Tscal *__restrict rho, + const Tvec *__restrict vel, + const Tscal *__restrict press, + const Tvec *__restrict grad_rho, + const Tvec *__restrict dx_v, + const Tvec *__restrict dy_v, + const Tvec *__restrict dz_v, + const Tvec *__restrict grad_P, + Tscal *__restrict dt_rho, + Tvec *__restrict dt_vel, + Tscal *__restrict dt_press) { + Tscal rho_i = rho[i]; + Tvec v_i = vel[i]; + Tscal P_i = press[i]; + Tvec grad_rho_i = grad_rho[i]; + Tvec dx_v_i = dx_v[i]; + Tvec dy_v_i = dy_v[i]; + Tvec dz_v_i = dz_v[i]; + Tvec grad_P_i = grad_P[i]; + + dt_rho[i] = -( + sham::dot(v_i, grad_rho_i) + rho_i * (dx_v_i[0] + dy_v_i[1] + dz_v_i[2])); + + dt_vel[i] + = -(v_i[0] * dx_v_i + v_i[1] * dy_v_i + v_i[2] * dz_v_i + grad_P_i / rho_i); + + dt_press[i] + = -(gamma * P_i * (dx_v_i[0] + dy_v_i[1] + dz_v_i[2]) + + sham::dot(v_i, grad_P_i)); + }); + } + }; + +} // namespace + +namespace shammodels::basegodunov::modules { + + template + void NodeEulerTimeDerivativeGas::_impl_evaluate_internal() { + auto edges = get_edges(); + + edges.spans_rho.check_sizes(edges.sizes.indexes); + edges.spans_vel.check_sizes(edges.sizes.indexes); + edges.spans_press.check_sizes(edges.sizes.indexes); + edges.spans_grad_rho.check_sizes(edges.sizes.indexes); + edges.spans_dx_v.check_sizes(edges.sizes.indexes); + edges.spans_dy_v.check_sizes(edges.sizes.indexes); + edges.spans_dz_v.check_sizes(edges.sizes.indexes); + edges.spans_grad_P.check_sizes(edges.sizes.indexes); + + edges.spans_dt_rho.ensure_sizes(edges.sizes.indexes); + edges.spans_dt_vel.ensure_sizes(edges.sizes.indexes); + edges.spans_dt_press.ensure_sizes(edges.sizes.indexes); + + KernelEulerTimeDerivativeGas::kernel( + edges.spans_rho.get_spans(), + edges.spans_vel.get_spans(), + edges.spans_press.get_spans(), + edges.spans_grad_rho.get_spans(), + edges.spans_dx_v.get_spans(), + edges.spans_dy_v.get_spans(), + edges.spans_dz_v.get_spans(), + edges.spans_grad_P.get_spans(), + edges.spans_dt_rho.get_spans(), + edges.spans_dt_vel.get_spans(), + edges.spans_dt_press.get_spans(), + edges.sizes.indexes, + block_size, + gamma); + } + + template + std::string NodeEulerTimeDerivativeGas::_impl_get_tex() const { + + auto block_count = get_ro_edge_base(0).get_tex_symbol(); + auto rho = get_ro_edge_base(1).get_tex_symbol(); + auto vel = get_ro_edge_base(2).get_tex_symbol(); + auto press = get_ro_edge_base(3).get_tex_symbol(); + auto grad_rho = get_ro_edge_base(4).get_tex_symbol(); + auto dx_v = get_ro_edge_base(5).get_tex_symbol(); + auto dy_v = get_ro_edge_base(6).get_tex_symbol(); + auto dz_v = get_ro_edge_base(7).get_tex_symbol(); + auto grad_P = get_ro_edge_base(8).get_tex_symbol(); + auto dt_rho = get_rw_edge_base(0).get_tex_symbol(); + auto dt_vel = get_rw_edge_base(1).get_tex_symbol(); + auto dt_press = get_rw_edge_base(2).get_tex_symbol(); + + std::string tex = R"tex( + Euler time derivatives of the gas primitive state + + \begin{align} + {dt_rho}_i &= - \left( {vel}_i \cdot {grad_rho}_i + + {rho}_i \left( {dx_v}_{i,x} + {dy_v}_{i,y} + {dz_v}_{i,z} \right) \right) \\ + {dt_vel}_i &= - \left( {vel}_{i,x} {dx_v}_i + {vel}_{i,y} {dy_v}_i + + {vel}_{i,z} {dz_v}_i + \frac{ {grad_P}_i }{ {rho}_i } \right) \\ + {dt_press}_i &= - \left( \gamma {press}_i + \left( {dx_v}_{i,x} + {dy_v}_{i,y} + {dz_v}_{i,z} \right) + + {vel}_i \cdot {grad_P}_i \right) \\ + i &\in [0,{block_count} * N_{\rm cell/block}) \\ + \gamma &= {gamma} \\ + N_{\rm cell/block} & = {block_size} + \end{align} + )tex"; + + shambase::replace_all(tex, "{dt_rho}", dt_rho); + shambase::replace_all(tex, "{dt_vel}", dt_vel); + shambase::replace_all(tex, "{dt_press}", dt_press); + shambase::replace_all(tex, "{rho}", rho); + shambase::replace_all(tex, "{vel}", vel); + shambase::replace_all(tex, "{press}", press); + shambase::replace_all(tex, "{grad_rho}", grad_rho); + shambase::replace_all(tex, "{dx_v}", dx_v); + shambase::replace_all(tex, "{dy_v}", dy_v); + shambase::replace_all(tex, "{dz_v}", dz_v); + shambase::replace_all(tex, "{grad_P}", grad_P); + shambase::replace_all(tex, "{block_count}", block_count); + shambase::replace_all(tex, "{gamma}", sham::format("{}", gamma)); + shambase::replace_all(tex, "{block_size}", sham::format("{}", block_size)); + + return tex; + } + +} // namespace shammodels::basegodunov::modules + +template class shammodels::basegodunov::modules::NodeEulerTimeDerivativeGas; diff --git a/src/shammodels/ramses/src/modules/InterpolateToFace.cpp b/src/shammodels/ramses/src/modules/InterpolateToFace.cpp index 51f6c682c1..d5d1b2978d 100644 --- a/src/shammodels/ramses/src/modules/InterpolateToFace.cpp +++ b/src/shammodels/ramses/src/modules/InterpolateToFace.cpp @@ -82,10 +82,7 @@ namespace { shamrock::PatchDataFieldSpanPointer grad_rho_cell; // For time interpolation Tscal dt_interp; - shamrock::PatchDataFieldSpanPointer vel_cell; - shamrock::PatchDataFieldSpanPointer dx_v_cell; - shamrock::PatchDataFieldSpanPointer dy_v_cell; - shamrock::PatchDataFieldSpanPointer dz_v_cell; + shamrock::PatchDataFieldSpanPointer dt_rho_cell; class acc { public: @@ -95,10 +92,7 @@ namespace { const Tvec *acc_grad_rho_cell; // For time interpolation - const Tvec *acc_vel_cell; - const Tvec *acc_dx_v_cell; - const Tvec *acc_dy_v_cell; - const Tvec *acc_dz_v_cell; + const Tscal *acc_dt_rho_cell; Tscal dt_interp; @@ -108,18 +102,10 @@ namespace { const Tvec *grad_rho_cell, // For time interpolation Tscal dt_interp, - const Tvec *vel_cell, - const Tvec *dx_v_cell, - const Tvec *dy_v_cell, - const Tvec *dz_v_cell) + const Tscal *dt_rho_cell) : shift_get(aabb_block_lower, aabb_cell_size), acc_rho_cell{rho_cell}, - acc_grad_rho_cell{grad_rho_cell}, dt_interp(dt_interp), acc_vel_cell{vel_cell}, - acc_dx_v_cell{dx_v_cell}, acc_dy_v_cell{dy_v_cell}, acc_dz_v_cell{dz_v_cell} {} - - Tscal get_dt_rho( - Tscal rho, Tvec v, Tvec grad_rho, Tvec dx_v, Tvec dy_v, Tvec dz_v) const { - return -(sham::dot(v, grad_rho) + rho * (dx_v[0] + dy_v[1] + dz_v[2])); - } + acc_grad_rho_cell{grad_rho_cell}, acc_dt_rho_cell{dt_rho_cell}, + dt_interp(dt_interp) {} std::array get_link_field_val(u32 id_a, u32 id_b) const { @@ -130,24 +116,13 @@ namespace { Tscal rho_b = acc_rho_cell[id_b]; Tvec grad_rho_b = acc_grad_rho_cell[id_b]; - Tvec vel_a = acc_vel_cell[id_a]; - Tvec dx_v_a = acc_dx_v_cell[id_a]; - Tvec dy_v_a = acc_dy_v_cell[id_a]; - Tvec dz_v_a = acc_dz_v_cell[id_a]; - Tvec vel_b = acc_vel_cell[id_b]; - Tvec dx_v_b = acc_dx_v_cell[id_b]; - Tvec dy_v_b = acc_dy_v_cell[id_b]; - Tvec dz_v_b = acc_dz_v_cell[id_b]; - // Spatial interpolate Tscal rho_face_a = rho_a + sycl::dot(grad_rho_a, shift_a); Tscal rho_face_b = rho_b + sycl::dot(grad_rho_b, shift_b); // Interpolate also to half a timestep - rho_face_a - += get_dt_rho(rho_a, vel_a, grad_rho_a, dx_v_a, dy_v_a, dz_v_a) * dt_interp; - rho_face_b - += get_dt_rho(rho_b, vel_b, grad_rho_b, dx_v_b, dy_v_b, dz_v_b) * dt_interp; + rho_face_a += acc_dt_rho_cell[id_a] * dt_interp; + rho_face_b += acc_dt_rho_cell[id_b] * dt_interp; return {rho_face_a, rho_face_b}; } @@ -161,10 +136,7 @@ namespace { grad_rho_cell.get_read_access(deps), // For time interpolation dt_interp, - vel_cell.get_read_access(deps), - dx_v_cell.get_read_access(deps), - dy_v_cell.get_read_access(deps), - dz_v_cell.get_read_access(deps)); + dt_rho_cell.get_read_access(deps)); } inline void complete_event_state(sycl::event e) { @@ -172,10 +144,7 @@ namespace { aabb_cell_size.complete_event_state(e); rho_cell.complete_event_state(e); grad_rho_cell.complete_event_state(e); - vel_cell.complete_event_state(e); - dx_v_cell.complete_event_state(e); - dy_v_cell.complete_event_state(e); - dz_v_cell.complete_event_state(e); + dt_rho_cell.complete_event_state(e); } }; @@ -193,8 +162,7 @@ namespace { shamrock::PatchDataFieldSpanPointer dz_v_cell; // For time interpolation Tscal dt_interp; - shamrock::PatchDataFieldSpanPointer rho_cell; - shamrock::PatchDataFieldSpanPointer grad_P_cell; + shamrock::PatchDataFieldSpanPointer dt_vel_cell; class acc { public: @@ -206,8 +174,7 @@ namespace { const Tvec *acc_dz_v_cell; // For time interpolation - const Tscal *acc_rho_cell; - const Tvec *acc_grad_P_cell; + const Tvec *acc_dt_vel_cell; Tscal dt_interp; @@ -219,15 +186,10 @@ namespace { const Tvec *dz_v_cell, // For time interpolation Tscal dt_interp, - const Tscal *rho_cell, - const Tvec *grad_P_cell) + const Tvec *dt_vel_cell) : shift_get(aabb_block_lower, aabb_cell_size), acc_vel_cell{vel_cell}, acc_dx_v_cell{dx_v_cell}, acc_dy_v_cell{dy_v_cell}, acc_dz_v_cell{dz_v_cell}, - dt_interp(dt_interp), acc_rho_cell{rho_cell}, acc_grad_P_cell{grad_P_cell} {} - - Tvec get_dt_v(Tvec v, Tvec dx_v, Tvec dy_v, Tvec dz_v, Tscal rho, Tvec grad_P) const { - return -(v[0] * dx_v + v[1] * dy_v + v[2] * dz_v + grad_P / rho); - } + acc_dt_vel_cell{dt_vel_cell}, dt_interp(dt_interp) {} std::array get_link_field_val(u32 id_a, u32 id_b) const { @@ -243,18 +205,13 @@ namespace { Tvec dy_vel_b = acc_dy_v_cell[id_b]; Tvec dz_vel_b = acc_dz_v_cell[id_b]; - Tscal rho_a = acc_rho_cell[id_a]; - Tvec grad_P_a = acc_grad_P_cell[id_a]; - Tscal rho_b = acc_rho_cell[id_b]; - Tvec grad_P_b = acc_grad_P_cell[id_b]; - Tvec dx_v_a_dot_shift = shift_a.x() * dx_vel_a + shift_a.y() * dy_vel_a + shift_a.z() * dz_vel_a; Tvec dx_v_b_dot_shift = shift_b.x() * dx_vel_b + shift_b.y() * dy_vel_b + shift_b.z() * dz_vel_b; - Tvec dt_v_a = get_dt_v(v_a, dx_vel_a, dy_vel_a, dz_vel_a, rho_a, grad_P_a); - Tvec dt_v_b = get_dt_v(v_b, dx_vel_b, dy_vel_b, dz_vel_b, rho_b, grad_P_b); + Tvec dt_v_a = acc_dt_vel_cell[id_a]; + Tvec dt_v_b = acc_dt_vel_cell[id_b]; Tvec vel_face_a = v_a + dx_v_a_dot_shift + dt_v_a * dt_interp; Tvec vel_face_b = v_b + dx_v_b_dot_shift + dt_v_b * dt_interp; @@ -273,8 +230,7 @@ namespace { dz_v_cell.get_read_access(deps), // For time interpolation dt_interp, - rho_cell.get_read_access(deps), - grad_P_cell.get_read_access(deps)); + dt_vel_cell.get_read_access(deps)); } inline void complete_event_state(sycl::event e) { @@ -284,8 +240,7 @@ namespace { dx_v_cell.complete_event_state(e); dy_v_cell.complete_event_state(e); dz_v_cell.complete_event_state(e); - rho_cell.complete_event_state(e); - grad_P_cell.complete_event_state(e); + dt_vel_cell.complete_event_state(e); } }; @@ -300,11 +255,7 @@ namespace { shamrock::PatchDataFieldSpanPointer grad_P_cell; // For time interpolation Tscal dt_interp; - Tscal gamma; - shamrock::PatchDataFieldSpanPointer vel_cell; - shamrock::PatchDataFieldSpanPointer dx_v_cell; - shamrock::PatchDataFieldSpanPointer dy_v_cell; - shamrock::PatchDataFieldSpanPointer dz_v_cell; + shamrock::PatchDataFieldSpanPointer dt_P_cell; class acc { public: @@ -314,12 +265,8 @@ namespace { const Tvec *acc_grad_P_cell; // For time interpolation - const Tvec *acc_vel_cell; - const Tvec *acc_dx_v_cell; - const Tvec *acc_dy_v_cell; - const Tvec *acc_dz_v_cell; + const Tscal *acc_dt_P_cell; - Tscal gamma; Tscal dt_interp; acc(const Tvec *aabb_block_lower, @@ -328,20 +275,9 @@ namespace { const Tvec *grad_P_cell, // For time interpolation Tscal dt_interp, - Tscal gamma, - const Tvec *vel_cell, - const Tvec *dx_v_cell, - const Tvec *dy_v_cell, - const Tvec *dz_v_cell) + const Tscal *dt_P_cell) : shift_get(aabb_block_lower, aabb_cell_size), acc_P_cell{P_cell}, - acc_grad_P_cell{grad_P_cell}, dt_interp(dt_interp), gamma(gamma), - acc_vel_cell{vel_cell}, acc_dx_v_cell{dx_v_cell}, acc_dy_v_cell{dy_v_cell}, - acc_dz_v_cell{dz_v_cell} {} - - Tscal get_dt_P( - Tscal P, Tvec grad_P, Tvec v, Tvec dx_v, Tvec dy_v, Tvec dz_v, Tscal gamma) const { - return -(gamma * P * (dx_v[0] + dy_v[1] + dz_v[2]) + sham::dot(v, grad_P)); - } + acc_grad_P_cell{grad_P_cell}, acc_dt_P_cell{dt_P_cell}, dt_interp(dt_interp) {} std::array get_link_field_val(u32 id_a, u32 id_b) const { @@ -352,17 +288,8 @@ namespace { Tscal P_b = acc_P_cell[id_b]; Tvec grad_P_b = acc_grad_P_cell[id_b]; - Tvec v_a = acc_vel_cell[id_a]; - Tvec dx_v_a = acc_dx_v_cell[id_a]; - Tvec dy_v_a = acc_dy_v_cell[id_a]; - Tvec dz_v_a = acc_dz_v_cell[id_a]; - Tvec v_b = acc_vel_cell[id_b]; - Tvec dx_v_b = acc_dx_v_cell[id_b]; - Tvec dy_v_b = acc_dy_v_cell[id_b]; - Tvec dz_v_b = acc_dz_v_cell[id_b]; - - Tscal dtP_cell_a = get_dt_P(P_a, grad_P_a, v_a, dx_v_a, dy_v_a, dz_v_a, gamma); - Tscal dtP_cell_b = get_dt_P(P_b, grad_P_b, v_b, dx_v_b, dy_v_b, dz_v_b, gamma); + Tscal dtP_cell_a = acc_dt_P_cell[id_a]; + Tscal dtP_cell_b = acc_dt_P_cell[id_b]; Tscal P_face_a = P_a + sycl::dot(grad_P_a, shift_a) + dtP_cell_a * dt_interp; Tscal P_face_b = P_b + sycl::dot(grad_P_b, shift_b) + dtP_cell_b * dt_interp; @@ -381,11 +308,7 @@ namespace { P_cell.get_read_access(deps), grad_P_cell.get_read_access(deps), dt_interp, - gamma, - vel_cell.get_read_access(deps), - dx_v_cell.get_read_access(deps), - dy_v_cell.get_read_access(deps), - dz_v_cell.get_read_access(deps)); + dt_P_cell.get_read_access(deps)); } inline void complete_event_state(sycl::event e) { @@ -393,10 +316,7 @@ namespace { aabb_cell_size.complete_event_state(e); P_cell.complete_event_state(e); grad_P_cell.complete_event_state(e); - vel_cell.complete_event_state(e); - dx_v_cell.complete_event_state(e); - dy_v_cell.complete_event_state(e); - dz_v_cell.complete_event_state(e); + dt_P_cell.complete_event_state(e); } }; @@ -412,10 +332,7 @@ namespace { shamrock::PatchDataFieldSpanPointer grad_rho_dust_cell; // For time interpolation Tscal dt_interp; - shamrock::PatchDataFieldSpanPointer vel_dust_cell; - shamrock::PatchDataFieldSpanPointer dx_v_dust_cell; - shamrock::PatchDataFieldSpanPointer dy_v_dust_cell; - shamrock::PatchDataFieldSpanPointer dz_v_dust_cell; + shamrock::PatchDataFieldSpanPointer dt_rho_dust_cell; class acc { public: @@ -426,10 +343,7 @@ namespace { const Tvec *acc_grad_rho_dust_cell; // For time interpolation - const Tvec *acc_vel_dust_cell; - const Tvec *acc_dx_v_dust_cell; - const Tvec *acc_dy_v_dust_cell; - const Tvec *acc_dz_v_dust_cell; + const Tscal *acc_dt_rho_dust_cell; Tscal dt_interp; @@ -440,27 +354,10 @@ namespace { const Tvec *grad_rho_dust_cell, // For time interpolation Tscal dt_interp, - const Tvec *vel_dust_cell, - const Tvec *dx_v_dust_cell, - const Tvec *dy_v_dust_cell, - const Tvec *dz_v_dust_cell) + const Tscal *dt_rho_dust_cell) : shift_get(aabb_block_lower, aabb_cell_size), nvar(nvar), acc_rho_dust_cell{rho_dust_cell}, acc_grad_rho_dust_cell{grad_rho_dust_cell}, - dt_interp(dt_interp), acc_vel_dust_cell{vel_dust_cell}, - acc_dx_v_dust_cell{dx_v_dust_cell}, acc_dy_v_dust_cell{dy_v_dust_cell}, - acc_dz_v_dust_cell{dz_v_dust_cell} {} - - Tscal get_dt_rho_dust( - Tscal rho_dust, - Tvec v_dust, - Tvec grad_rho_dust, - Tvec dx_v_dust, - Tvec dy_v_dust, - Tvec dz_v_dust) const { - return -( - sham::dot(v_dust, grad_rho_dust) - + rho_dust * (dx_v_dust[0] + dy_v_dust[1] + dz_v_dust[2])); - } + acc_dt_rho_dust_cell{dt_rho_dust_cell}, dt_interp(dt_interp) {} std::array get_link_field_val(u32 id_a, u32 id_b) const { const u32 icell_a = id_a / nvar; @@ -473,34 +370,11 @@ namespace { Tscal rho_dust_b = acc_rho_dust_cell[id_b]; Tvec grad_rho_dust_b = acc_grad_rho_dust_cell[id_b]; - Tvec vel_dust_a = acc_vel_dust_cell[id_a]; - Tvec dx_v_dust_a = acc_dx_v_dust_cell[id_a]; - Tvec dy_v_dust_a = acc_dy_v_dust_cell[id_a]; - Tvec dz_v_dust_a = acc_dz_v_dust_cell[id_a]; - Tvec vel_dust_b = acc_vel_dust_cell[id_b]; - Tvec dx_v_dust_b = acc_dx_v_dust_cell[id_b]; - Tvec dy_v_dust_b = acc_dy_v_dust_cell[id_b]; - Tvec dz_v_dust_b = acc_dz_v_dust_cell[id_b]; - Tscal rho_dust_face_a = rho_dust_a + sycl::dot(grad_rho_dust_a, shift_a); Tscal rho_dust_face_b = rho_dust_b + sycl::dot(grad_rho_dust_b, shift_b); - rho_dust_face_a += get_dt_rho_dust( - rho_dust_a, - vel_dust_a, - grad_rho_dust_a, - dx_v_dust_a, - dy_v_dust_a, - dz_v_dust_a) - * dt_interp; - rho_dust_face_b += get_dt_rho_dust( - rho_dust_b, - vel_dust_b, - grad_rho_dust_b, - dx_v_dust_b, - dy_v_dust_b, - dz_v_dust_b) - * dt_interp; + rho_dust_face_a += acc_dt_rho_dust_cell[id_a] * dt_interp; + rho_dust_face_b += acc_dt_rho_dust_cell[id_b] * dt_interp; return {rho_dust_face_a, rho_dust_face_b}; } @@ -515,10 +389,7 @@ namespace { grad_rho_dust_cell.get_read_access(deps), // For time interpolation dt_interp, - vel_dust_cell.get_read_access(deps), - dx_v_dust_cell.get_read_access(deps), - dy_v_dust_cell.get_read_access(deps), - dz_v_dust_cell.get_read_access(deps)); + dt_rho_dust_cell.get_read_access(deps)); } inline void complete_event_state(sycl::event e) { @@ -526,10 +397,7 @@ namespace { aabb_cell_size.complete_event_state(e); rho_dust_cell.complete_event_state(e); grad_rho_dust_cell.complete_event_state(e); - vel_dust_cell.complete_event_state(e); - dx_v_dust_cell.complete_event_state(e); - dy_v_dust_cell.complete_event_state(e); - dz_v_dust_cell.complete_event_state(e); + dt_rho_dust_cell.complete_event_state(e); } }; @@ -547,7 +415,7 @@ namespace { shamrock::PatchDataFieldSpanPointer dz_v_dust_cell; // For time interpolation Tscal dt_interp; - shamrock::PatchDataFieldSpanPointer rho_dust_cell; + shamrock::PatchDataFieldSpanPointer dt_vel_dust_cell; class acc { public: @@ -560,7 +428,7 @@ namespace { const Tvec *acc_dz_v_dust_cell; // For time interpolation - const Tscal *acc_rho_dust_cell; + const Tvec *acc_dt_vel_dust_cell; Tscal dt_interp; @@ -573,15 +441,11 @@ namespace { const Tvec *dz_v_dust_cell, // For time interpolation Tscal dt_interp, - const Tscal *rho_dust_cell) + const Tvec *dt_vel_dust_cell) : shift_get(aabb_block_lower, aabb_cell_size), nvar(nvar), acc_vel_dust_cell{vel_dust_cell}, acc_dx_v_dust_cell{dx_v_dust_cell}, acc_dy_v_dust_cell{dy_v_dust_cell}, acc_dz_v_dust_cell{dz_v_dust_cell}, - dt_interp(dt_interp), acc_rho_dust_cell{rho_dust_cell} {} - - Tvec get_dt_v_dust(Tvec v, Tvec dx_v, Tvec dy_v, Tvec dz_v, Tscal rho) const { - return -(v[0] * dx_v + v[1] * dy_v + v[2] * dz_v); - } + acc_dt_vel_dust_cell{dt_vel_dust_cell}, dt_interp(dt_interp) {} std::array get_link_field_val(u32 id_a, u32 id_b) const { const u32 icell_a = id_a / nvar; @@ -599,9 +463,6 @@ namespace { Tvec dy_vel_dust_b = acc_dy_v_dust_cell[id_b]; Tvec dz_vel_dust_b = acc_dz_v_dust_cell[id_b]; - Tscal rho_dust_a = acc_rho_dust_cell[id_a]; - Tscal rho_dust_b = acc_rho_dust_cell[id_b]; - Tvec dx_v_dust_a_dot_shift = shift_a.x() * dx_vel_dust_a + shift_a.y() * dy_vel_dust_a + shift_a.z() * dz_vel_dust_a; @@ -609,10 +470,8 @@ namespace { + shift_b.y() * dy_vel_dust_b + shift_b.z() * dz_vel_dust_b; - Tvec dt_v_dust_a = get_dt_v_dust( - v_dust_a, dx_vel_dust_a, dy_vel_dust_a, dz_vel_dust_a, rho_dust_a); - Tvec dt_v_dust_b = get_dt_v_dust( - v_dust_b, dx_vel_dust_b, dy_vel_dust_b, dz_vel_dust_b, rho_dust_b); + Tvec dt_v_dust_a = acc_dt_vel_dust_cell[id_a]; + Tvec dt_v_dust_b = acc_dt_vel_dust_cell[id_b]; Tvec vel_dust_face_a = v_dust_a + dx_v_dust_a_dot_shift + dt_v_dust_a * dt_interp; Tvec vel_dust_face_b = v_dust_b + dx_v_dust_b_dot_shift + dt_v_dust_b * dt_interp; @@ -632,7 +491,7 @@ namespace { dz_v_dust_cell.get_read_access(deps), // For time interpolation dt_interp, - rho_dust_cell.get_read_access(deps)); + dt_vel_dust_cell.get_read_access(deps)); } inline void complete_event_state(sycl::event e) { @@ -642,7 +501,7 @@ namespace { dx_v_dust_cell.complete_event_state(e); dy_v_dust_cell.complete_event_state(e); dz_v_dust_cell.complete_event_state(e); - rho_dust_cell.complete_event_state(e); + dt_vel_dust_cell.complete_event_state(e); } }; @@ -682,10 +541,7 @@ void shammodels::basegodunov::modules::InterpolateToFaceRho:: auto spans_cell0block_aabb_lower = edges.spans_cell0block_aabb_lower.get_spans(); auto spans_rhos = edges.spans_rhos.get_spans(); auto spans_grad_rho = edges.spans_grad_rho.get_spans(); - auto spans_vel = edges.spans_vel.get_spans(); - auto spans_dx_vel = edges.spans_dx_vel.get_spans(); - auto spans_dy_vel = edges.spans_dy_vel.get_spans(); - auto spans_dz_vel = edges.spans_dz_vel.get_spans(); + auto spans_dt_rho = edges.spans_dt_rho.get_spans(); using Interp = RhoInterpolate; auto interpolators @@ -696,10 +552,7 @@ void shammodels::basegodunov::modules::InterpolateToFaceRho:: spans_rhos.get(id), spans_grad_rho.get(id), dt_interp, - spans_vel.get(id), - spans_dx_vel.get(id), - spans_dy_vel.get(id), - spans_dz_vel.get(id)}; + spans_dt_rho.get(id)}; }); auto graphs_xp = edges.cell_neigh_graph.get_refs_dir(Direction::xp); @@ -836,12 +689,11 @@ void shammodels::basegodunov::modules::InterpolateToFaceVel:: auto spans_block_cell_sizes = edges.spans_block_cell_sizes.get_spans(); auto spans_cell0block_aabb_lower = edges.spans_cell0block_aabb_lower.get_spans(); - auto spans_rhos = edges.spans_rhos.get_spans(); - auto spans_grad_P = edges.spans_grad_P.get_spans(); auto spans_vel = edges.spans_vel.get_spans(); auto spans_dx_vel = edges.spans_dx_vel.get_spans(); auto spans_dy_vel = edges.spans_dy_vel.get_spans(); auto spans_dz_vel = edges.spans_dz_vel.get_spans(); + auto spans_dt_vel = edges.spans_dt_vel.get_spans(); using Interp = VelInterpolate; auto interpolators @@ -854,8 +706,7 @@ void shammodels::basegodunov::modules::InterpolateToFaceVel:: spans_dy_vel.get(id), spans_dz_vel.get(id), dt_interp, - spans_rhos.get(id), - spans_grad_P.get(id)}; + spans_dt_vel.get(id)}; }); auto graphs_xp = edges.cell_neigh_graph.get_refs_dir(Direction::xp); @@ -994,10 +845,7 @@ void shammodels::basegodunov::modules::InterpolateToFacePress:: auto spans_cell0block_aabb_lower = edges.spans_cell0block_aabb_lower.get_spans(); auto spans_press = edges.spans_press.get_spans(); auto spans_grad_P = edges.spans_grad_P.get_spans(); - auto spans_vel = edges.spans_vel.get_spans(); - auto spans_dx_vel = edges.spans_dx_vel.get_spans(); - auto spans_dy_vel = edges.spans_dy_vel.get_spans(); - auto spans_dz_vel = edges.spans_dz_vel.get_spans(); + auto spans_dt_press = edges.spans_dt_press.get_spans(); using Interp = PressInterpolate; auto interpolators @@ -1008,11 +856,7 @@ void shammodels::basegodunov::modules::InterpolateToFacePress:: spans_press.get(id), spans_grad_P.get(id), dt_interp, - gamma, - spans_vel.get(id), - spans_dx_vel.get(id), - spans_dy_vel.get(id), - spans_dz_vel.get(id)}; + spans_dt_press.get(id)}; }); auto graphs_xp = edges.cell_neigh_graph.get_refs_dir(Direction::xp); @@ -1158,10 +1002,7 @@ void shammodels::basegodunov::modules::InterpolateToFaceRhoDust: auto spans_cell0block_aabb_lower = edges.spans_cell0block_aabb_lower.get_spans(); auto spans_rhos_dust = edges.spans_rhos_dust.get_spans(); auto spans_grad_rho_dust = edges.spans_grad_rho_dust.get_spans(); - auto spans_vel_dust = edges.spans_vel_dust.get_spans(); - auto spans_dx_vel_dust = edges.spans_dx_vel_dust.get_spans(); - auto spans_dy_vel_dust = edges.spans_dy_vel_dust.get_spans(); - auto spans_dz_vel_dust = edges.spans_dz_vel_dust.get_spans(); + auto spans_dt_rho_dust = edges.spans_dt_rho_dust.get_spans(); using Interp = RhoDustInterpolate; auto interpolators @@ -1173,10 +1014,7 @@ void shammodels::basegodunov::modules::InterpolateToFaceRhoDust: spans_rhos_dust.get(id), spans_grad_rho_dust.get(id), dt_interp, - spans_vel_dust.get(id), - spans_dx_vel_dust.get(id), - spans_dy_vel_dust.get(id), - spans_dz_vel_dust.get(id)}; + spans_dt_rho_dust.get(id)}; }); auto graphs_xp = edges.cell_neigh_graph.get_refs_dir(Direction::xp); @@ -1338,11 +1176,11 @@ void shammodels::basegodunov::modules::InterpolateToFaceVelDust: auto spans_block_cell_sizes = edges.spans_block_cell_sizes.get_spans(); auto spans_cell0block_aabb_lower = edges.spans_cell0block_aabb_lower.get_spans(); - auto spans_rhos_dust = edges.spans_rhos_dust.get_spans(); auto spans_vel_dust = edges.spans_vel_dust.get_spans(); auto spans_dx_vel_dust = edges.spans_dx_vel_dust.get_spans(); auto spans_dy_vel_dust = edges.spans_dy_vel_dust.get_spans(); auto spans_dz_vel_dust = edges.spans_dz_vel_dust.get_spans(); + auto spans_dt_vel_dust = edges.spans_dt_vel_dust.get_spans(); using Interp = VelDustInterpolate; auto interpolators @@ -1356,7 +1194,7 @@ void shammodels::basegodunov::modules::InterpolateToFaceVelDust: spans_dy_vel_dust.get(id), spans_dz_vel_dust.get(id), dt_interp, - spans_rhos_dust.get(id)}; + spans_dt_vel_dust.get(id)}; }); auto graphs_xp = edges.cell_neigh_graph.get_refs_dir(Direction::xp);