Skip to content

Commit ff836c8

Browse files
Guo-astroGuo-astro
committed
Add GSPH solver and VTK I/O
Adds main solver class, VTK output, and smoothing length iteration for GSPH: Solver.hpp/cpp: - Main GSPH solver orchestrating the simulation - Predictor-corrector time integration (leapfrog) - Ghost particle exchange via SPH's BasicSPHGhostHandler - Neighbor search using compressed BVH trees - Timestep calculation based on CFL and signal velocity - Reuses SPH density utilities (rho_h, newtown_iterate_new_h) - Reuses sph_pressure_symetric for force computation with Riemann p* modules/IterateSmoothingLengthDensity.hpp/cpp: - Outer-loop smoothing length iteration for density convergence - Computes density and omega (grad-h correction) after h converges - Integrates with LoopSmoothingLengthIter from SPH module modules/io/VTKDump.hpp/cpp: - VTK Legacy format output for visualization - Outputs positions, velocity, density, pressure, internal energy - Compatible with ParaView SPH code reuse: - Ghost handling: BasicSPHGhostHandler for all BC types - Forces: sph_pressure_symetric() with Riemann solver pressure - Density: rho_h(), h_rho(), newtown_iterate_new_h() - Safe division: inv_sat_zero() for numerical stability Fixes compatibility with main: - Use get_eos_gamma() instead of gamma member variable - Stub checkpoint methods (GSPHCheckpoint to be added separately)
1 parent 9d3f7f6 commit ff836c8

7 files changed

Lines changed: 2407 additions & 1 deletion

File tree

‎src/shammodels/gsph/CMakeLists.txt‎

Lines changed: 4 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -11,10 +11,13 @@ cmake_minimum_required(VERSION 3.9)
1111

1212
project(Shammodels_gsph CXX C)
1313

14-
# Sources: Core infrastructure + Physics modules
14+
# Sources: GSPH solver and VTK I/O
1515
set(Sources
1616
src/SolverConfig.cpp
17+
src/Solver.cpp
1718
src/modules/UpdateDerivs.cpp
19+
src/modules/IterateSmoothingLengthDensity.cpp
20+
src/modules/io/VTKDump.cpp
1821
)
1922

2023
if(SHAMROCK_USE_SHARED_LIB)

‎src/shammodels/gsph/include/shammodels/gsph/Solver.hpp‎

Lines changed: 242 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,242 @@
1+
// -------------------------------------------------------//
2+
//
3+
// SHAMROCK code for hydrodynamics
4+
// Copyright (c) 2021-2025 Timothée David--Cléris <tim.shamrock@proton.me>
5+
// SPDX-License-Identifier: CeCILL Free Software License Agreement v2.1
6+
// Shamrock is licensed under the CeCILL 2.1 License, see LICENSE for more information
7+
//
8+
// -------------------------------------------------------//
9+
10+
#pragma once
11+
12+
/**
13+
* @file Solver.hpp
14+
* @author Guo Yansong (guo.yansong.ngy@gmail.com)
15+
* @author Yona Lapeyre (yona.lapeyre@ens-lyon.fr) --no git blame--
16+
* @brief GSPH Solver class
17+
*
18+
* The GSPH method originated from:
19+
* - Inutsuka, S. (2002) "Reformulation of Smoothed Particle Hydrodynamics
20+
* with Riemann Solver"
21+
*
22+
* This implementation follows:
23+
* - Cha, S.-H. & Whitworth, A.P. (2003) "Implementations and tests of
24+
* Godunov-type particle hydrodynamics"
25+
*/
26+
27+
#include "shambase/exception.hpp"
28+
#include "SolverConfig.hpp"
29+
#include "shambackends/vec.hpp"
30+
#include "shammodels/gsph/modules/SolverStorage.hpp"
31+
#include "shammodels/sph/BasicSPHGhosts.hpp"
32+
#include "shammodels/sph/SPHUtilities.hpp"
33+
#include "shammodels/sph/SolverLog.hpp"
34+
#include "shamrock/patch/PatchDataLayerLayout.hpp"
35+
#include "shamrock/scheduler/ComputeField.hpp"
36+
#include "shamrock/scheduler/InterfacesUtility.hpp"
37+
#include "shamrock/scheduler/SerialPatchTree.hpp"
38+
#include "shamrock/scheduler/ShamrockCtx.hpp"
39+
#include "shamsys/legacy/log.hpp"
40+
#include "shamtree/TreeTraversalCache.hpp"
41+
#include <memory>
42+
#include <stdexcept>
43+
#include <variant>
44+
45+
namespace shammodels::gsph {
46+
47+
struct TimestepLog {
48+
i32 rank;
49+
f64 rate;
50+
u64 npart;
51+
f64 tcompute;
52+
53+
inline f64 rate_sum() { return shamalgs::collective::allreduce_sum(rate); }
54+
inline u64 npart_sum() { return shamalgs::collective::allreduce_sum(npart); }
55+
inline f64 tcompute_max() { return shamalgs::collective::allreduce_max(tcompute); }
56+
};
57+
58+
/**
59+
* @brief The GSPH Solver class
60+
*
61+
* Implements the Godunov SPH method using Riemann solvers at particle
62+
* interfaces instead of artificial viscosity.
63+
*
64+
* @tparam Tvec Vector type (e.g., f64_3)
65+
* @tparam SPHKernel Kernel type (e.g., M4, M6, C2, C4, C6)
66+
*/
67+
template<class Tvec, template<class> class SPHKernel>
68+
class Solver {
69+
public:
70+
using Tscal = shambase::VecComponent<Tvec>;
71+
static constexpr u32 dim = shambase::VectorProperties<Tvec>::dimension;
72+
using Kernel = SPHKernel<Tscal>;
73+
74+
using Config = SolverConfig<Tvec, SPHKernel>;
75+
76+
using u_morton = u32;
77+
78+
static constexpr Tscal Rkern = Kernel::Rkern;
79+
80+
ShamrockCtx &context;
81+
inline PatchScheduler &scheduler() { return shambase::get_check_ref(context.sched); }
82+
83+
SolverStorage<Tvec, u_morton> storage{};
84+
85+
Config solver_config;
86+
sph::SolverLog solve_logs;
87+
88+
inline void init_required_fields() { solver_config.set_layout(context.get_pdl_write()); }
89+
90+
// Serial patch tree control
91+
void gen_serial_patch_tree();
92+
inline void reset_serial_patch_tree() { storage.serial_patch_tree.reset(); }
93+
94+
// Ghost handling - reuse SPH ghost handler
95+
using GhostHandle = sph::BasicSPHGhostHandler<Tvec>;
96+
using GhostHandleCache = typename GhostHandle::CacheMap;
97+
98+
void gen_ghost_handler(Tscal time_val);
99+
inline void reset_ghost_handler() { storage.ghost_handler.reset(); }
100+
101+
void build_ghost_cache();
102+
void clear_ghost_cache();
103+
104+
void merge_position_ghost();
105+
106+
// Tree operations
107+
using RTree = typename Config::RTree;
108+
void build_merged_pos_trees();
109+
void clear_merged_pos_trees();
110+
111+
void compute_presteps_rint();
112+
void reset_presteps_rint();
113+
114+
void start_neighbors_cache();
115+
void reset_neighbors_cache();
116+
117+
void gsph_prestep(Tscal time_val, Tscal dt);
118+
119+
void apply_position_boundary(Tscal time_val);
120+
121+
void do_predictor_leapfrog(Tscal dt);
122+
123+
void init_ghost_layout();
124+
125+
void communicate_merge_ghosts_fields();
126+
void reset_merge_ghosts_fields();
127+
128+
void compute_omega();
129+
void compute_eos_fields();
130+
void reset_eos_fields();
131+
132+
void prepare_corrector();
133+
134+
/**
135+
* @brief Update derivatives using GSPH Riemann solver
136+
*
137+
* This is the core GSPH step: for each particle pair, solve
138+
* the 1D Riemann problem and compute forces from the interface
139+
* pressure p*.
140+
*/
141+
void update_derivs();
142+
143+
/**
144+
* @brief Compute CFL timestep constraint
145+
*
146+
* Computes timestep from:
147+
* - Courant condition: dt_cour = C_cour * h / vsig
148+
* - Force condition: dt_force = C_force * sqrt(h / |a|)
149+
*
150+
* @return Minimum CFL timestep across all particles
151+
*/
152+
Tscal compute_dt_cfl();
153+
154+
bool apply_corrector(Tscal dt, u64 Npart_all);
155+
156+
void update_sync_load_values();
157+
158+
Solver(ShamrockCtx &context) : context(context) {}
159+
160+
void init_solver_graph();
161+
162+
void vtk_do_dump(std::string filename, bool add_patch_world_id);
163+
164+
/**
165+
* @brief Write a checkpoint to disk
166+
*
167+
* Saves the current simulation state to checkpoint files.
168+
*
169+
* @param basename Base filename (without extension)
170+
*/
171+
void write_checkpoint(const std::string &basename);
172+
173+
/**
174+
* @brief Read a checkpoint from disk
175+
*
176+
* Loads simulation state from checkpoint files.
177+
*
178+
* @param basename Base filename (without extension)
179+
*/
180+
void read_checkpoint(const std::string &basename);
181+
182+
/**
183+
* @brief Check if a checkpoint file exists
184+
*
185+
* @param basename Base filename (without extension)
186+
* @return true if checkpoint exists
187+
*/
188+
static bool checkpoint_exists(const std::string &basename);
189+
190+
inline void print_timestep_logs() {
191+
if (shamcomm::world_rank() == 0) {
192+
logger::info_ln(
193+
"GSPH", "iteration since start :", solve_logs.get_iteration_count());
194+
logger::info_ln(
195+
"GSPH", "time since start :", shambase::details::get_wtime(), "(s)");
196+
}
197+
}
198+
199+
TimestepLog evolve_once();
200+
201+
Tscal evolve_once_time_expl(Tscal t_current, Tscal dt_input) {
202+
solver_config.set_time(t_current);
203+
solver_config.set_next_dt(dt_input);
204+
evolve_once();
205+
return solver_config.get_dt();
206+
}
207+
208+
inline bool evolve_until(Tscal target_time, i32 niter_max = -1) {
209+
auto step = [&]() {
210+
Tscal dt = solver_config.get_dt();
211+
Tscal t = solver_config.get_time();
212+
213+
if (t > target_time) {
214+
throw shambase::make_except_with_loc<std::invalid_argument>(
215+
"the target time is higher than the current time");
216+
}
217+
218+
if (t + dt > target_time) {
219+
solver_config.set_next_dt(target_time - t);
220+
}
221+
evolve_once();
222+
};
223+
224+
i32 iter_count = 0;
225+
226+
while (solver_config.get_time() < target_time) {
227+
step();
228+
iter_count++;
229+
230+
if ((iter_count >= niter_max) && (niter_max != -1)) {
231+
logger::info_ln("GSPH", "stopping evolve until because of niter =", iter_count);
232+
return false;
233+
}
234+
}
235+
236+
print_timestep_logs();
237+
238+
return true;
239+
}
240+
};
241+
242+
} // namespace shammodels::gsph

‎src/shammodels/gsph/include/shammodels/gsph/modules/IterateSmoothingLengthDensity.hpp‎

Lines changed: 105 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,105 @@
1+
// -------------------------------------------------------//
2+
//
3+
// SHAMROCK code for hydrodynamics
4+
// Copyright (c) 2021-2025 Timothée David--Cléris <tim.shamrock@proton.me>
5+
// SPDX-License-Identifier: CeCILL Free Software License Agreement v2.1
6+
// Shamrock is licensed under the CeCILL 2.1 License, see LICENSE for more information
7+
//
8+
// -------------------------------------------------------//
9+
10+
#pragma once
11+
12+
/**
13+
* @file IterateSmoothingLengthDensity.hpp
14+
* @author Guo Yansong (guo.yansong.ngy@gmail.com)
15+
* @brief Declares the IterateSmoothingLengthDensity module for GSPH.
16+
*
17+
* This module performs ONE Newton-Raphson iteration step for the smoothing length,
18+
* computing density and omega (grad-h correction) in the process.
19+
*
20+
* Unlike the buggy inner-loop approach, this module is designed to be called
21+
* from an outer loop (LoopSmoothingLengthIter) which can rebuild the neighbor
22+
* cache if h changes too much.
23+
*/
24+
25+
#include "shambackends/vec.hpp"
26+
#include "shammodels/sph/solvergraph/NeighCache.hpp"
27+
#include "shamrock/solvergraph/Field.hpp"
28+
#include "shamrock/solvergraph/IFieldSpan.hpp"
29+
#include "shamrock/solvergraph/INode.hpp"
30+
#include "shamrock/solvergraph/Indexes.hpp"
31+
#include <memory>
32+
33+
namespace shammodels::gsph::modules {
34+
35+
/**
36+
* @brief GSPH smoothing length iteration module.
37+
*
38+
* Performs one Newton-Raphson step to find h that satisfies:
39+
* rho_sum = rho_h(m, h) = m * (hfact/h)^3
40+
*
41+
* Also computes:
42+
* - density (SPH summation)
43+
* - omega (grad-h correction factor)
44+
*
45+
* @tparam Tvec Vector type (e.g., f64_3)
46+
* @tparam SPHKernel SPH kernel type (e.g., M6<f64>)
47+
*/
48+
template<class Tvec, class SPHKernel>
49+
class IterateSmoothingLengthDensity : public shamrock::solvergraph::INode {
50+
51+
using Tscal = shambase::VecComponent<Tvec>;
52+
53+
Tscal gpart_mass;
54+
Tscal h_evol_max;
55+
Tscal h_evol_iter_max;
56+
57+
public:
58+
IterateSmoothingLengthDensity(Tscal gpart_mass, Tscal h_evol_max, Tscal h_evol_iter_max)
59+
: gpart_mass(gpart_mass), h_evol_max(h_evol_max), h_evol_iter_max(h_evol_iter_max) {}
60+
61+
struct Edges {
62+
const shamrock::solvergraph::Indexes<u32> &sizes;
63+
const shammodels::sph::solvergraph::NeighCache &neigh_cache;
64+
const shamrock::solvergraph::IFieldSpan<Tvec> &positions;
65+
const shamrock::solvergraph::IFieldSpan<Tscal> &old_h;
66+
shamrock::solvergraph::IFieldSpan<Tscal> &new_h;
67+
shamrock::solvergraph::IFieldSpan<Tscal> &eps_h;
68+
shamrock::solvergraph::Field<Tscal> &density;
69+
shamrock::solvergraph::Field<Tscal> &omega;
70+
};
71+
72+
inline void set_edges(
73+
std::shared_ptr<shamrock::solvergraph::Indexes<u32>> sizes,
74+
std::shared_ptr<shammodels::sph::solvergraph::NeighCache> neigh_cache,
75+
std::shared_ptr<shamrock::solvergraph::IFieldSpan<Tvec>> positions,
76+
std::shared_ptr<shamrock::solvergraph::IFieldSpan<Tscal>> old_h,
77+
std::shared_ptr<shamrock::solvergraph::IFieldSpan<Tscal>> new_h,
78+
std::shared_ptr<shamrock::solvergraph::IFieldSpan<Tscal>> eps_h,
79+
std::shared_ptr<shamrock::solvergraph::Field<Tscal>> density,
80+
std::shared_ptr<shamrock::solvergraph::Field<Tscal>> omega) {
81+
__internal_set_ro_edges({sizes, neigh_cache, positions, old_h});
82+
__internal_set_rw_edges({new_h, eps_h, density, omega});
83+
}
84+
85+
inline Edges get_edges() {
86+
return Edges{
87+
get_ro_edge<shamrock::solvergraph::Indexes<u32>>(0),
88+
get_ro_edge<shammodels::sph::solvergraph::NeighCache>(1),
89+
get_ro_edge<shamrock::solvergraph::IFieldSpan<Tvec>>(2),
90+
get_ro_edge<shamrock::solvergraph::IFieldSpan<Tscal>>(3),
91+
get_rw_edge<shamrock::solvergraph::IFieldSpan<Tscal>>(0),
92+
get_rw_edge<shamrock::solvergraph::IFieldSpan<Tscal>>(1),
93+
get_rw_edge<shamrock::solvergraph::Field<Tscal>>(2),
94+
get_rw_edge<shamrock::solvergraph::Field<Tscal>>(3)};
95+
}
96+
97+
void _impl_evaluate_internal();
98+
99+
inline virtual std::string _impl_get_label() const {
100+
return "GSPHIterateSmoothingLengthDensity";
101+
};
102+
103+
virtual std::string _impl_get_tex() const;
104+
};
105+
} // namespace shammodels::gsph::modules

0 commit comments

Comments
 (0)