diff --git a/exemples/ramses/godunov_disc.py b/exemples/ramses/godunov_disc.py new file mode 100644 index 0000000000..195117b750 --- /dev/null +++ b/exemples/ramses/godunov_disc.py @@ -0,0 +1,260 @@ +import os + +import matplotlib +import matplotlib.animation as animation +import matplotlib.pyplot as plt +import numpy as np +from numba import njit + +import shamrock + +### Parameters + +## Grid setup +resolution = 128 + +## Initial conditions +contrast_factor = 1000 +ref_density = 1 +ref_sound_speed = 0.2 +cs_expo = 0.5 +point_mass_Gmass = 1 +gravity_softening = 0.05 +outer_radius = 0.25 +inner_radius_factor = 0.1 +inner_iso_z_flaring = 2 +gamma0 = 1.6 +smallr = 1e-4 +center = 0.5 +h_over_r = 1 + + +## Initial conditions + + +@njit +def get_cs(rc_soft, z): + inner_radius = np.sqrt( + outer_radius + * inner_radius_factor + * (outer_radius * inner_radius_factor + inner_iso_z_flaring * np.abs(z)) + ) + + if rc_soft > inner_radius: + cs = ref_sound_speed * (rc_soft / outer_radius) ** (-cs_expo) + else: + cs = ref_sound_speed * (inner_radius / outer_radius) ** (-cs_expo) + + return cs, inner_radius + + +@njit +def get_omega(rs_soft, rc_soft, cs, inner_radius): + + if cs_expo == 0.5 or rc_soft > inner_radius: + omega2 = point_mass_Gmass / (rc_soft * rc_soft * rs_soft) - (4.0 - cs_expo) * (cs * cs) / ( + rc_soft * rc_soft + ) + omega = np.sqrt(np.maximum(omega2, 0.0)) + else: + omega2 = point_mass_Gmass / (rs_soft**3) - (3.0 - cs_expo) * (cs * cs) / (rc_soft * rc_soft) + omega = np.sqrt(np.maximum(omega2, 0.0)) + + return omega + + +@njit +def rho_map(rmin, rmax): + + rmin = np.array(rmin) - center + rmax = np.array(rmax) - center + x, y, z = (rmin + rmax) / 2 + + rc_soft = np.sqrt(x * x + y * y + gravity_softening**2) + rs_soft = np.sqrt(x * x + y * y + z * z + gravity_softening**2) + + cs, _ = get_cs(rc_soft, z) + + # Radial power law + dens = ref_density * (outer_radius / rc_soft) ** ((5.0 - 2.0 * cs_expo) / 2.0) + + # Vertical hydrostatic equilibrium + dens *= np.exp((point_mass_Gmass / (cs * cs)) * (1.0 / rs_soft - 1.0 / rc_soft)) + + # Outer taper + if rc_soft > outer_radius or np.abs(z) > 0.5 * outer_radius: + dens /= contrast_factor + + return np.maximum(dens, smallr) + + +@njit +def rhovel_map(rmin, rmax): + + rho = rho_map(rmin, rmax) + + rmin = np.array(rmin) - center + rmax = np.array(rmax) - center + x, y, z = (rmin + rmax) / 2 + + rc_soft = np.sqrt(x * x + y * y + gravity_softening**2) + rs = np.sqrt(x * x + y * y + z * z) + rs_soft = np.sqrt(x * x + y * y + z * z + gravity_softening**2) + + cs, inner_radius = get_cs(rc_soft, z) + omega = get_omega(rs_soft, rc_soft, cs, inner_radius) + + x_soft = x * (rs_soft / rs) + y_soft = y * (rs_soft / rs) + + vx = -omega * y_soft + vy = omega * x_soft + vz = 0.0 + + return (rho * vx, rho * vy, rho * vz) + + +@njit +def rhoe_map(rmin, rmax): + + rho = rho_map(rmin, rmax) + rho_vx, rho_vy, rho_vz = rhovel_map(rmin, rmax) + + rmin = np.array(rmin) - center + rmax = np.array(rmax) - center + x, y, z = (rmin + rmax) / 2 + + rc_soft = np.sqrt(x * x + y * y + gravity_softening**2) + cs, _ = get_cs(rc_soft, z) + + eint = rho * cs * cs / (gamma0 - 1.0) + ekin = 0.5 * (rho_vx * rho_vx + rho_vy * rho_vy + rho_vz * rho_vz) / np.maximum(rho, smallr) + + return eint + ekin + + +def make_cartesian_coords(nx, ny, z_val, min_x, max_x, min_y, max_y): + # Create the cylindrical coordinate grid + x_vals = np.linspace(min_x, max_x, nx) + y_vals = np.linspace(min_y, max_y, ny) + + # Create meshgrid + x_grid, y_grid = np.meshgrid(x_vals, y_vals) + + # Convert to Cartesian coordinates (z = 0 for a disc in the xy-plane) + z_grid = z_val * np.ones_like(x_grid) + + # Flatten and stack to create list of positions + positions = np.column_stack([x_grid.ravel(), y_grid.ravel(), z_grid.ravel()]) + + return [tuple(pos) for pos in positions] + + +nx, ny = 128, 128 + + +positions = make_cartesian_coords(nx, ny, 0.5, 0, 1 - 1e-6, 0, 1 - 1e-6) + + +def plot_rho_slice_cartesian(metadata, arr_rho_pos, iplot, case_name, dpi=200): + ext = metadata["extent"] + + my_cmap = matplotlib.colormaps["gist_heat"].copy() # copy the default cmap + my_cmap.set_bad(color="yellow") + + arr_rho_pos = np.array(arr_rho_pos).reshape(nx, ny) + + plt.figure(dpi=dpi) + res = plt.imshow( + arr_rho_pos, + cmap=my_cmap, + origin="lower", + extent=ext, + aspect="auto", + norm=matplotlib.colors.LogNorm(vmin=1e-3, vmax=6), + ) + plt.xlabel("x") + plt.ylabel("y") + plt.title(f"t = {metadata['time']:0.5f} [seconds]") + cbar = plt.colorbar(res, extend="both") + plt.gca().set_aspect("equal") + cbar.set_label(r"$\rho$ [code unit]") + plt.savefig(os.path.join(".", f"rho_{case_name}_{iplot:04d}.png")) + plt.close() + + +def plot(t, iplot): + metadata = {"extent": [0, 1, 0, 1], "time": t} + arr_rho_pos = model.render_slice("rho", "f64", positions) + plot_rho_slice_cartesian(metadata, arr_rho_pos, iplot, "disk") + pos_XZ = [(x, z, y) for (x, y, z) in positions] + arr_rho_pos = model.render_slice("rho", "f64", pos_XZ) + plot_rho_slice_cartesian(metadata, arr_rho_pos, iplot, "disk_xz") + pos_YZ = [(z, y, x) for (x, y, z) in positions] + arr_rho_pos = model.render_slice("rho", "f64", pos_YZ) + plot_rho_slice_cartesian(metadata, arr_rho_pos, iplot, "disk_yz") + + +def plot_force(t, iplot): + metadata = {"extent": [0, 1, 0, 1], "time": t} + force = model.render_slice("gravitational_force", "f64_3", positions) + force = np.array(force).reshape(nx, ny, 3) + force = np.linalg.norm(force, axis=2) + plot_rho_slice_cartesian(metadata, force, iplot, "force") + + +### Initialisation + +print("initialising model...") + +ctx = shamrock.Context() +ctx.pdata_layout_new() + +model = shamrock.get_Model_Ramses(context=ctx, vector_type="f64_3", grid_repr="i64_3") + + +cell_size = 2 +base = resolution + +cfg = model.gen_default_config() + +scale_fact = 1 / (cell_size * base) +cfg.set_scale_factor(scale_fact) +cfg.set_eos_gamma(1.00001) + +cfg.set_riemann_solver_hll() +cfg.set_slope_lim_minmod() +cfg.set_face_time_interpolation(True) + +model.set_solver_config(cfg) + +model.init_scheduler(int(1e9), 1) # (crit_split - patches, crit_merge - patches) +model.make_base_grid((0, 0, 0), (cell_size, cell_size, cell_size), (base, base, base)) +print("initialising fields...", flush=True) + +print("initial density...", flush=True) +model.set_field_value_lambda_f64("rho", rho_map) +print("initial momentum...", flush=True) +model.set_field_value_lambda_f64_3("rhovel", rhovel_map) +print("initial total energy...", flush=True) +model.set_field_value_lambda_f64("rhoetot", rhoe_map) +print("initialisation done, plotting initial conditions...", flush=True) + +dt = 0.01 +for i in range(0, 30): + t = dt * i + model.evolve_until(dt * i) + d = ctx.collect_data() + print( + f"Step {i}, t = {t:0.5f} seconds, min rho = {d['rho'].min():0.5e}, max rho = {d['rho'].max():0.5e}" + ) + print( + f"Step {i}, t = {t:0.5f} seconds, min rhoe = {d['rhoetot'].min():0.5e}, max rhov = {d['rhoetot'].max():0.5e}" + ) + print( + f"Step {i}, t = {t:0.5f} seconds, min rhov = {d['rhovel'].min():0.5e}, max rhov = {d['rhovel'].max():0.5e}" + ) + plot(t, i) + + +plot(dt * (i + 1), i) diff --git a/src/shammodels/ramses/CMakeLists.txt b/src/shammodels/ramses/CMakeLists.txt index d9e291cec5..1fc604d97d 100644 --- a/src/shammodels/ramses/CMakeLists.txt +++ b/src/shammodels/ramses/CMakeLists.txt @@ -45,6 +45,7 @@ set(Sources src/modules/FuseGhostLayer.cpp src/modules/render/GridRender.cpp src/modules/ComputeCoordinates.cpp + src/modules/ComputeAnalyticalGravity.cpp src/SolverConfig.cpp ) diff --git a/src/shammodels/ramses/include/shammodels/ramses/SolverConfig.hpp b/src/shammodels/ramses/include/shammodels/ramses/SolverConfig.hpp index 9e2f1ab285..c2d84a17dd 100644 --- a/src/shammodels/ramses/include/shammodels/ramses/SolverConfig.hpp +++ b/src/shammodels/ramses/include/shammodels/ramses/SolverConfig.hpp @@ -80,12 +80,26 @@ namespace shammodels::basegodunov { struct GravityConfig { using Tscal = shambase::VecComponent; GravityMode gravity_mode = NoGravity; - bool analytical_gravity = false; // whether to use an external analytical gravity Tscal tol = 1e-6; inline Tscal get_tolerance() { return tol; } inline bool is_gravity_on() { return gravity_mode != NoGravity; } }; + template + struct AnalyticalGravityConfig { + using Tscal = shambase::VecComponent; + + AnalyticalGravityMode analytical_gravity_mode = POINTMASS; + // parameters for analytical gravity can be added here + Tscal point_mass_GM = 1; + Tvec point_mass_pos = Tvec(0.5, 0.5, 0.5); + Tscal epsilon_softening = 0.05; + Tscal min_rho_for_analytical_gravity = 1e-5; // minimum density to apply analytical gravity + inline bool is_analytical_gravity_on() { + return analytical_gravity_mode != NoAnalyticalGravity; + } + }; + template struct SolverStatusVar; @@ -208,9 +222,15 @@ struct shammodels::basegodunov::SolverConfig { auto scal_G = get_constant_G(); return 4 * M_PI * scal_G; } + AnalyticalGravityConfig analytical_gravity_config{}; inline Tscal get_grav_tol() { return gravity_config.get_tolerance(); } inline bool is_gravity_on() { return gravity_config.is_gravity_on(); } - inline bool is_coordinate_field_required() { return gravity_config.analytical_gravity; } + inline bool is_coordinate_field_required() { + return analytical_gravity_config.is_analytical_gravity_on(); + } + inline bool is_analytical_gravity_on() { + return analytical_gravity_config.is_analytical_gravity_on(); + } ////////////////////////////////////////////////////////////////////////////////////////////// // Gravity config (END) diff --git a/src/shammodels/ramses/include/shammodels/ramses/config/enum_GravityMode.hpp b/src/shammodels/ramses/include/shammodels/ramses/config/enum_GravityMode.hpp index 463714c93b..535a644e39 100644 --- a/src/shammodels/ramses/include/shammodels/ramses/config/enum_GravityMode.hpp +++ b/src/shammodels/ramses/include/shammodels/ramses/config/enum_GravityMode.hpp @@ -38,4 +38,14 @@ namespace shammodels::basegodunov { {GravityMode::BICGSTAB, "bicgstab"}, {GravityMode::MULTIGRID, "multigrid"}}); + enum AnalyticalGravityMode { + NoAnalyticalGravity = 0, + POINTMASS = 1, // point mass analytical gravity + }; + + SHAMROCK_JSON_SERIALIZE_ENUM( + AnalyticalGravityMode, + {{AnalyticalGravityMode::NoAnalyticalGravity, "no_analytical_gravity"}, + {AnalyticalGravityMode::POINTMASS, "point_mass"}}); + } // namespace shammodels::basegodunov diff --git a/src/shammodels/ramses/include/shammodels/ramses/modules/ComputeAnalyticalGravity.hpp b/src/shammodels/ramses/include/shammodels/ramses/modules/ComputeAnalyticalGravity.hpp new file mode 100644 index 0000000000..ec5b1ad868 --- /dev/null +++ b/src/shammodels/ramses/include/shammodels/ramses/modules/ComputeAnalyticalGravity.hpp @@ -0,0 +1,65 @@ +// -------------------------------------------------------// +// +// 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 ComputeAnalyticalGravity.hpp + * @author Noé Brucy (noe.brucy@ens-lyon.fr) + * @brief Compute a gravitational acceleration at the center of each cell, + using an analytical formula. + */ + +#include "shambackends/vec.hpp" +#include "shammodels/ramses/SolverConfig.hpp" +#include "shamrock/solvergraph/Field.hpp" +#include "shamrock/solvergraph/INode.hpp" +#include "shamrock/solvergraph/Indexes.hpp" +#include "shamrock/solvergraph/ScalarsEdge.hpp" + +namespace shammodels::basegodunov::modules { + + template + class NodeComputeAnalyticalGravity : public shamrock::solvergraph::INode { + + using AGConfig = typename shammodels::basegodunov::AnalyticalGravityConfig; + using Tscal = shambase::VecComponent; + + u32 block_size; + AGConfig config; + + public: + NodeComputeAnalyticalGravity(u32 block_size, AGConfig config) + : block_size(block_size), config(config) {} + +#define NODE_COMPUTE_ANALYTICAL_GRAVITY(X_RO, X_RW) \ + /* inputs */ \ + X_RO( \ + shamrock::solvergraph::Indexes, \ + sizes) /* number of blocks per patch for all patches on the current MPI process*/ \ + X_RO( \ + shamrock::solvergraph::IFieldSpan, \ + spans_coordinates) /* center coordinates of each cell */ \ + /* outputs */ \ + X_RW(shamrock::solvergraph::IFieldSpan, spans_gravitational_force) + + EXPAND_NODE_EDGES(NODE_COMPUTE_ANALYTICAL_GRAVITY) + +#undef NODE_COMPUTE_ANALYTICAL_GRAVITY + + void _impl_evaluate_internal(); + + void _impl_reset_internal() {}; + + inline virtual std::string _impl_get_label() const { return "ComputeAnalyticalGravity"; }; + + virtual std::string _impl_get_tex() const { return "ComputeAnalyticalGravity"; }; + }; + +} // namespace shammodels::basegodunov::modules diff --git a/src/shammodels/ramses/include/shammodels/ramses/modules/SolverStorage.hpp b/src/shammodels/ramses/include/shammodels/ramses/modules/SolverStorage.hpp index fd7ab07fcf..0ccabc4b47 100644 --- a/src/shammodels/ramses/include/shammodels/ramses/modules/SolverStorage.hpp +++ b/src/shammodels/ramses/include/shammodels/ramses/modules/SolverStorage.hpp @@ -94,6 +94,7 @@ namespace shammodels::basegodunov { std::shared_ptr> block_cell_sizes; std::shared_ptr> cell0block_aabb_lower; std::shared_ptr> coordinates; + std::shared_ptr> gravitational_force; std::shared_ptr> grad_rho; std::shared_ptr> dx_v; diff --git a/src/shammodels/ramses/include/shammodels/ramses/modules/TimeIntegrator.hpp b/src/shammodels/ramses/include/shammodels/ramses/modules/TimeIntegrator.hpp index 538e960993..c8649b1962 100644 --- a/src/shammodels/ramses/include/shammodels/ramses/modules/TimeIntegrator.hpp +++ b/src/shammodels/ramses/include/shammodels/ramses/modules/TimeIntegrator.hpp @@ -46,7 +46,7 @@ namespace shammodels::basegodunov::modules { TimeIntegrator(ShamrockCtx &context, Config &solver_config, Storage &storage) : context(context), solver_config(solver_config), storage(storage) {} - void forward_euler(Tscal dt); + void forward_euler(Tscal dt, bool precall = false); private: inline PatchScheduler &scheduler() { return shambase::get_check_ref(context.sched); } diff --git a/src/shammodels/ramses/src/Solver.cpp b/src/shammodels/ramses/src/Solver.cpp index 899604ce67..f3417ce1df 100644 --- a/src/shammodels/ramses/src/Solver.cpp +++ b/src/shammodels/ramses/src/Solver.cpp @@ -25,6 +25,7 @@ #include "shammodels/ramses/SolverConfig.hpp" #include "shammodels/ramses/modules/AMRGridRefinementHandler.hpp" #include "shammodels/ramses/modules/BlockNeighToCellNeigh.hpp" +#include "shammodels/ramses/modules/ComputeAnalyticalGravity.hpp" #include "shammodels/ramses/modules/ComputeCFL.hpp" #include "shammodels/ramses/modules/ComputeCellAABB.hpp" #include "shammodels/ramses/modules/ComputeCoordinates.hpp" @@ -423,6 +424,12 @@ void shammodels::basegodunov::Solver::init_solver_graph() { "level0_amr", "level0_amr"); } + if (solver_config.is_analytical_gravity_on()) { + // will be filled by NodeComputeAnalyticalGravity + storage.gravitational_force = std::make_shared>( + AMRBlock::block_size, "gravitational_force", "\\mathbf{f}_a"); + } + storage.grad_rho = std::make_shared>( AMRBlock::block_size, "grad_rho", "\\nabla \\rho"); storage.dx_v = std::make_shared>( @@ -707,7 +714,8 @@ void shammodels::basegodunov::Solver::init_solver_graph() { cfg_bc_to_geom(solver_config.bc_config.get_y()), cfg_bc_to_geom(solver_config.bc_config.get_z())}; - // if outflow we want zero gradient so we skip the vector transformation in TransformGhostLayer + // if outflow we want zero gradient so we skip the vector transformation in + // TransformGhostLayer bool transform_vec_x = solver_config.bc_config.get_x() != BCConfig::GhostType::Outflow; bool transform_vec_y = solver_config.bc_config.get_y() != BCConfig::GhostType::Outflow; bool transform_vec_z = solver_config.bc_config.get_z() != BCConfig::GhostType::Outflow; @@ -792,7 +800,8 @@ void shammodels::basegodunov::Solver::init_solver_graph() { // enable this to debug GZ //{ - // auto filename_edge = std::make_shared>( + // auto filename_edge = + // std::make_shared>( // "debug_fuse.vtk", "debug_fuse.vtk"); // filename_edge->data = "debug_fuse.vtk"; // @@ -954,8 +963,8 @@ void shammodels::basegodunov::Solver::init_solver_graph() { std::make_shared(std::move(node_coordinates))); } - if (solver_config.amr_mode.need_level_zero_compute()) { // compute level0 sizes in patch (to be - // enabled later when needed) + if (solver_config.amr_mode.need_level_zero_compute()) { // compute level0 sizes in patch (to + // be enabled later when needed) modules::ComputeLevel0CellSize node_level0_sizes{}; node_level0_sizes.set_edges( graph.get_edge_ptr>>("global_patch_boxes"), @@ -965,6 +974,16 @@ void shammodels::basegodunov::Solver::init_solver_graph() { std::make_shared(std::move(node_level0_sizes))); } + if (solver_config.is_analytical_gravity_on()) { + modules::NodeComputeAnalyticalGravity node_analytical_gravity{ + AMRBlock::block_size, solver_config.analytical_gravity_config}; + node_analytical_gravity.set_edges( + storage.block_counts, storage.coordinates, storage.gravitational_force); + solver_sequence.push_back( + std::make_shared( + std::move(node_analytical_gravity))); + } + if (solver_config.should_compute_rho_mean()) { modules::NodeComputeMass node{AMRBlock::block_size}; node.set_edges( @@ -1528,6 +1547,12 @@ void shammodels::basegodunov::Solver::evolve_once() { shambase::get_check_ref(storage.solver_sequence).evaluate(); } + // apply gravity source term + { + modules::TimeIntegrator dt_integ(context, solver_config, storage); + dt_integ.forward_euler(dt_input, true); + } + // compute dt fields modules::ComputeTimeDerivative dt_compute(context, solver_config, storage); dt_compute.compute_dt_fields(); @@ -1707,8 +1732,9 @@ void shammodels::basegodunov::Solver::do_debug_vtk_dump(std::str */ /* - std::unique_ptr> dtrho = storage.dtrho.get().rankgather_computefield(sched); - writer.write_field("dtrho", dtrho, num_obj * block_size); + std::unique_ptr> dtrho = + storage.dtrho.get().rankgather_computefield(sched); writer.write_field("dtrho", dtrho, + num_obj * block_size); std::unique_ptr> dtrhov = storage.dtrhov.get().rankgather_computefield(sched); diff --git a/src/shammodels/ramses/src/SolverConfig.cpp b/src/shammodels/ramses/src/SolverConfig.cpp index 3d9f0fc068..2418d758d9 100644 --- a/src/shammodels/ramses/src/SolverConfig.cpp +++ b/src/shammodels/ramses/src/SolverConfig.cpp @@ -35,6 +35,7 @@ namespace shammodels::basegodunov { pdl.add_field("rho", AMRBlock::block_size); pdl.add_field("rhovel", AMRBlock::block_size); pdl.add_field("rhoetot", AMRBlock::block_size); + pdl.add_field("gravitational_force", AMRBlock::block_size); if (is_dust_on()) { u32 ndust = dust_config.ndust; diff --git a/src/shammodels/ramses/src/modules/ComputeAnalyticalGravity.cpp b/src/shammodels/ramses/src/modules/ComputeAnalyticalGravity.cpp new file mode 100644 index 0000000000..3b6b26d787 --- /dev/null +++ b/src/shammodels/ramses/src/modules/ComputeAnalyticalGravity.cpp @@ -0,0 +1,75 @@ +// -------------------------------------------------------// +// +// 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 ComputeAnalyticalGravity.cpp + * @author Noé Brucy (noe.brucy@ens-lyon.fr) + * @brief Compute a gravitational acceleration at the center of each cell, + using an analytical formula. + */ + +#include "shammodels/ramses/modules/ComputeAnalyticalGravity.hpp" +#include "shambackends/kernel_call_distrib.hpp" +#include "shamsys/NodeInstance.hpp" + +namespace shammodels::basegodunov::modules { + + template + void NodeComputeAnalyticalGravity::_impl_evaluate_internal() { + using Tscal = shambase::VecComponent; + + auto edges = get_edges(); + + edges.spans_coordinates.check_sizes(edges.sizes.indexes); + edges.spans_gravitational_force.ensure_sizes(edges.sizes.indexes); + + auto &coordinates_spans = edges.spans_coordinates.get_spans(); + auto &gravitational_force_spans = edges.spans_gravitational_force.get_spans(); + + shambase::DistributedData cell_counts + = edges.sizes.indexes.template map([&](u64 id, u32 block_count) { + u32 cell_count = block_count * block_size; + return cell_count; + }); + + static constexpr u32 dim = shambase::VectorProperties::dimension; + + if (config.analytical_gravity_mode == POINTMASS) { + + sham::distributed_data_kernel_call( + shamsys::instance::get_compute_scheduler_ptr(), + sham::DDMultiRef{coordinates_spans}, + sham::DDMultiRef{gravitational_force_spans}, + cell_counts, + [config = this->config, block_size = this->block_size]( + u32 i, const Tvec *__restrict cell_coor, Tvec *__restrict gravitational_force) { + u32 block_id = i / block_size; + Tvec cell_pos = cell_coor[i]; + + // compute vector from the point mass to the cell + Tvec r_vec = cell_pos - config.point_mass_pos; + + // compute distance to the point mass + Tscal r_squared = 0; + for (u32 d = 0; d < dim; d++) { + r_squared += r_vec[d] * r_vec[d]; + } + + Tscal one_over_r_cube_soft = pow( + r_squared + config.epsilon_softening * config.epsilon_softening, -1.5f); + + // compute gravitational acceleration using softened point mass potential + gravitational_force[i] = -config.point_mass_GM * r_vec * one_over_r_cube_soft; + }); + } + } + +} // namespace shammodels::basegodunov::modules + +template class shammodels::basegodunov::modules::NodeComputeAnalyticalGravity; diff --git a/src/shammodels/ramses/src/modules/TimeIntegrator.cpp b/src/shammodels/ramses/src/modules/TimeIntegrator.cpp index 94ac21b428..b0e872b3da 100644 --- a/src/shammodels/ramses/src/modules/TimeIntegrator.cpp +++ b/src/shammodels/ramses/src/modules/TimeIntegrator.cpp @@ -16,7 +16,8 @@ #include "shammodels/ramses/modules/TimeIntegrator.hpp" template -void shammodels::basegodunov::modules::TimeIntegrator::forward_euler(Tscal dt) { +void shammodels::basegodunov::modules::TimeIntegrator::forward_euler( + Tscal dt, bool precall) { StackEntry stack_loc{}; @@ -27,41 +28,98 @@ void shammodels::basegodunov::modules::TimeIntegrator::forward_e shamrock::solvergraph::Field &cfield_dtrho = shambase::get_check_ref(storage.dtrho); shamrock::solvergraph::Field &cfield_dtrhov = shambase::get_check_ref(storage.dtrhov); shamrock::solvergraph::Field &cfield_dtrhoe = shambase::get_check_ref(storage.dtrhoe); + shamrock::solvergraph::Field &cfield_gravitational_force + = shambase::get_check_ref(storage.gravitational_force); // load layout info PatchDataLayerLayout &pdl = scheduler().pdl_old(); - const u32 icell_min = pdl.get_field_idx("cell_min"); - const u32 icell_max = pdl.get_field_idx("cell_max"); - const u32 irho = pdl.get_field_idx("rho"); - const u32 irhoetot = pdl.get_field_idx("rhoetot"); - const u32 irhovel = pdl.get_field_idx("rhovel"); + const u32 icell_min = pdl.get_field_idx("cell_min"); + const u32 icell_max = pdl.get_field_idx("cell_max"); + const u32 irho = pdl.get_field_idx("rho"); + const u32 irhoetot = pdl.get_field_idx("rhoetot"); + const u32 irhovel = pdl.get_field_idx("rhovel"); + const u32 igravforce = pdl.get_field_idx("gravitational_force"); - scheduler().for_each_patchdata_nonempty( - [&, dt](const shamrock::patch::Patch p, shamrock::patch::PatchDataLayer &pdat) { - shamlog_debug_ln("[AMR Flux]", "forward euler integration patch", p.id_patch); + Tscal min_rho_for_analytical_gravity + = solver_config.analytical_gravity_config.min_rho_for_analytical_gravity; - sham::DeviceQueue &q = shamsys::instance::get_compute_scheduler().get_queue(); - u32 id = p.id_patch; + scheduler().for_each_patchdata_nonempty([&, dt, min_rho_for_analytical_gravity, precall]( + const shamrock::patch::Patch p, + shamrock::patch::PatchDataLayer &pdat) { + shamlog_debug_ln("[AMR Flux]", "forward euler integration patch", p.id_patch); - sham::DeviceBuffer &dt_rho_patch = cfield_dtrho.get_buf(id); - sham::DeviceBuffer &dt_rhov_patch = cfield_dtrhov.get_buf(id); - sham::DeviceBuffer &dt_rhoe_patch = cfield_dtrhoe.get_buf(id); + sham::DeviceQueue &q = shamsys::instance::get_compute_scheduler().get_queue(); + u32 id = p.id_patch; + + u32 cell_count = pdat.get_obj_cnt() * AMRBlock::block_size; + + sham::DeviceBuffer &buf_rho = pdat.get_field_buf_ref(irho); + sham::DeviceBuffer &buf_rhov = pdat.get_field_buf_ref(irhovel); + sham::DeviceBuffer &buf_rhoe = pdat.get_field_buf_ref(irhoetot); + + sham::EventList depends_list; - u32 cell_count = pdat.get_obj_cnt() * AMRBlock::block_size; + auto rho = buf_rho.get_write_access(depends_list); + auto rhov = buf_rhov.get_write_access(depends_list); + auto rhoe = buf_rhoe.get_write_access(depends_list); - sham::DeviceBuffer &buf_rho = pdat.get_field_buf_ref(irho); - sham::DeviceBuffer &buf_rhov = pdat.get_field_buf_ref(irhovel); - sham::DeviceBuffer &buf_rhoe = pdat.get_field_buf_ref(irhoetot); + if (precall && solver_config.is_analytical_gravity_on()) { + sham::DeviceBuffer &gravitational_force_patch + = cfield_gravitational_force.get_buf(id); + sham::DeviceBuffer &buf_gravforce = pdat.get_field_buf_ref(igravforce); + auto acc_gravitational_force_patch + = gravitational_force_patch.get_read_access(depends_list); + auto gravforce = buf_gravforce.get_write_access(depends_list); - sham::EventList depends_list; - auto acc_dt_rho_patch = dt_rho_patch.get_read_access(depends_list); - auto acc_dt_rhov_patch = dt_rhov_patch.get_read_access(depends_list); - auto acc_dt_rhoe_patch = dt_rhoe_patch.get_read_access(depends_list); + auto grav_event = q.submit( + depends_list, [&, dt, min_rho_for_analytical_gravity](sycl::handler &cgh) { + shambase::parallel_for( + cgh, cell_count, "add analytical gravity", [=](u32 id_a) { + Tvec g_vec = acc_gravitational_force_patch[id_a]; - auto rho = buf_rho.get_write_access(depends_list); - auto rhov = buf_rhov.get_write_access(depends_list); - auto rhoe = buf_rhoe.get_write_access(depends_list); + rho[id_a] = std::max(min_rho_for_analytical_gravity, rho[id_a]); + Tvec rhov_old = rhov[id_a]; + + Tscal ekin_old = 0; + for (u32 d = 0; d < shambase::VectorProperties::dimension; d++) { + ekin_old += rhov_old[d] * rhov_old[d]; + } + ekin_old *= 0.5 / rho[id_a]; + + Tscal e_other = rhoe[id_a] - ekin_old; + + // Update momentum + rhov[id_a] += dt * g_vec * rho[id_a]; + + Tscal ekin_new = 0; + for (u32 d = 0; d < shambase::VectorProperties::dimension; d++) { + ekin_new += rhov[id_a][d] * rhov[id_a][d]; + } + ekin_new *= 0.5 / rho[id_a]; + + rhoe[id_a] = e_other + ekin_new; + + gravforce[id_a] = g_vec; + }); + }); + + gravitational_force_patch.complete_event_state(grav_event); + + buf_rho.complete_event_state(grav_event); + buf_rhov.complete_event_state(grav_event); + buf_rhoe.complete_event_state(grav_event); + buf_gravforce.complete_event_state(grav_event); + } + + if (!precall) { + + sham::DeviceBuffer &dt_rho_patch = cfield_dtrho.get_buf(id); + sham::DeviceBuffer &dt_rhov_patch = cfield_dtrhov.get_buf(id); + sham::DeviceBuffer &dt_rhoe_patch = cfield_dtrhoe.get_buf(id); + auto acc_dt_rho_patch = dt_rho_patch.get_read_access(depends_list); + auto acc_dt_rhov_patch = dt_rhov_patch.get_read_access(depends_list); + auto acc_dt_rhoe_patch = dt_rhoe_patch.get_read_access(depends_list); auto e = q.submit(depends_list, [&, dt](sycl::handler &cgh) { shambase::parallel_for(cgh, cell_count, "accumulate fluxes", [=](u32 id_a) { @@ -80,55 +138,60 @@ void shammodels::basegodunov::modules::TimeIntegrator::forward_e buf_rho.complete_event_state(e); buf_rhov.complete_event_state(e); buf_rhoe.complete_event_state(e); - }); + } + }); - if (solver_config.is_dust_on()) { + if (!precall) { + if (solver_config.is_dust_on()) { - shamrock::solvergraph::Field &cfield_dtrho_dust - = shambase::get_check_ref(storage.dtrho_dust); - shamrock::solvergraph::Field &cfield_dtrhov_dust - = shambase::get_check_ref(storage.dtrhov_dust); + shamrock::solvergraph::Field &cfield_dtrho_dust + = shambase::get_check_ref(storage.dtrho_dust); + shamrock::solvergraph::Field &cfield_dtrhov_dust + = shambase::get_check_ref(storage.dtrhov_dust); - const u32 irho_dust = pdl.get_field_idx("rho_dust"); - const u32 irhovel_dust = pdl.get_field_idx("rhovel_dust"); + const u32 irho_dust = pdl.get_field_idx("rho_dust"); + const u32 irhovel_dust = pdl.get_field_idx("rhovel_dust"); - scheduler().for_each_patchdata_nonempty([&, dt]( - const shamrock::patch::Patch p, - shamrock::patch::PatchDataLayer &pdat) { - shamlog_debug_ln( - "[AMR Flux]", "forward euler integration patch for dust fields", p.id_patch); + scheduler().for_each_patchdata_nonempty([&, dt]( + const shamrock::patch::Patch p, + shamrock::patch::PatchDataLayer &pdat) { + shamlog_debug_ln( + "[AMR Flux]", "forward euler integration patch for dust fields", p.id_patch); - sham::DeviceQueue &q = shamsys::instance::get_compute_scheduler().get_queue(); - u32 id = p.id_patch; + sham::DeviceQueue &q = shamsys::instance::get_compute_scheduler().get_queue(); + u32 id = p.id_patch; - sham::DeviceBuffer &dt_rho_dust_patch = cfield_dtrho_dust.get_buf(id); - sham::DeviceBuffer &dt_rhov_dust_patch = cfield_dtrhov_dust.get_buf(id); + sham::DeviceBuffer &dt_rho_dust_patch = cfield_dtrho_dust.get_buf(id); + sham::DeviceBuffer &dt_rhov_dust_patch = cfield_dtrhov_dust.get_buf(id); - u32 cell_count = pdat.get_obj_cnt() * AMRBlock::block_size; - u32 ndust = solver_config.dust_config.ndust; + u32 cell_count = pdat.get_obj_cnt() * AMRBlock::block_size; + u32 ndust = solver_config.dust_config.ndust; - sham::DeviceBuffer &buf_rho_dust = pdat.get_field_buf_ref(irho_dust); - sham::DeviceBuffer &buf_rhov_dust = pdat.get_field_buf_ref(irhovel_dust); + sham::DeviceBuffer &buf_rho_dust = pdat.get_field_buf_ref(irho_dust); + sham::DeviceBuffer &buf_rhov_dust + = pdat.get_field_buf_ref(irhovel_dust); - sham::EventList depends_list; - auto acc_dt_rho_dust_patch = dt_rho_dust_patch.get_read_access(depends_list); - auto acc_dt_rhov_dust_patch = dt_rhov_dust_patch.get_read_access(depends_list); + sham::EventList depends_list; + auto acc_dt_rho_dust_patch = dt_rho_dust_patch.get_read_access(depends_list); + auto acc_dt_rhov_dust_patch = dt_rhov_dust_patch.get_read_access(depends_list); - auto rho_dust = buf_rho_dust.get_write_access(depends_list); - auto rhov_dust = buf_rhov_dust.get_write_access(depends_list); + auto rho_dust = buf_rho_dust.get_write_access(depends_list); + auto rhov_dust = buf_rhov_dust.get_write_access(depends_list); - auto e = q.submit(depends_list, [&, dt](sycl::handler &cgh) { - shambase::parallel_for(cgh, ndust * cell_count, "accumulate fluxes", [=](u32 id_a) { - rho_dust[id_a] += dt * acc_dt_rho_dust_patch[id_a]; - rhov_dust[id_a] += dt * acc_dt_rhov_dust_patch[id_a]; + auto e = q.submit(depends_list, [&, dt](sycl::handler &cgh) { + shambase::parallel_for( + cgh, ndust * cell_count, "accumulate fluxes", [=](u32 id_a) { + rho_dust[id_a] += dt * acc_dt_rho_dust_patch[id_a]; + rhov_dust[id_a] += dt * acc_dt_rhov_dust_patch[id_a]; + }); }); - }); - dt_rho_dust_patch.complete_event_state(e); - dt_rhov_dust_patch.complete_event_state(e); - buf_rho_dust.complete_event_state(e); - buf_rhov_dust.complete_event_state(e); - }); + dt_rho_dust_patch.complete_event_state(e); + dt_rhov_dust_patch.complete_event_state(e); + buf_rho_dust.complete_event_state(e); + buf_rhov_dust.complete_event_state(e); + }); + } } }