Skip to content

Commit 5b2f0bb

Browse files
tdavidcltdavidcl
authored
[SPH] add q_av dust corrections (#1950)
see discussion in danieljprice/phantom#857
1 parent 430c187 commit 5b2f0bb

6 files changed

Lines changed: 379 additions & 28 deletions

File tree

‎src/shammodels/sph/CMakeLists.txt‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -47,6 +47,7 @@ set(Sources
4747
src/modules/LoopSmoothingLengthIter.cpp
4848
src/modules/ComputeNeighStats.cpp
4949
src/modules/NodeUpdateDerivsVaryingAlphaAV.cpp
50+
src/modules/NodeUpdateDerivsVaryingAlphaAVDustTVA.cpp
5051
src/modules/NodeUpdateDerivsMonofluidTVA.cpp
5152
src/modules/NodeComputePressureGrad.cpp
5253
src/modules/NodeEvolveDustCOALASourceTerm.cpp

‎src/shammodels/sph/include/shammodels/sph/SolverConfig.hpp‎

Lines changed: 18 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -116,6 +116,9 @@ namespace shammodels::sph {
116116
Tscal cfl_density_threshold = shambase::get_epsilon<Tscal>();
117117

118118
bool ensure_s_j_positivity = true;
119+
120+
// use the corrected q_AV from Hutchison 2018 & Price Laibe 15
121+
bool dust_corrected_av = false;
119122
};
120123

121124
struct MonofluidComplete {
@@ -134,14 +137,16 @@ namespace shammodels::sph {
134137
Tscal C_1_fluid = 0.1,
135138
Tscal C_drift = 1.0,
136139
Tscal cfl_density_threshold = shambase::get_epsilon<Tscal>(),
137-
bool ensure_s_j_positivity = true) {
140+
bool ensure_s_j_positivity = true,
141+
bool dust_corrected_av = false) {
138142
current_mode = MonofluidTVA{
139143
nvar,
140144
pure_diffusion_mode,
141145
C_1_fluid,
142146
C_drift,
143147
cfl_density_threshold,
144-
ensure_s_j_positivity};
148+
ensure_s_j_positivity,
149+
dust_corrected_av};
145150
}
146151
inline void set_monofluid_complete(u32 nvar) { current_mode = MonofluidComplete{nvar}; }
147152

@@ -166,7 +171,8 @@ namespace shammodels::sph {
166171
{"C_1_fluid", cfg->C_1_fluid},
167172
{"C_drift", cfg->C_drift},
168173
{"cfl_density_threshold", cfg->cfl_density_threshold},
169-
{"ensure_s_j_positivity", cfg->ensure_s_j_positivity}};
174+
{"ensure_s_j_positivity", cfg->ensure_s_j_positivity},
175+
{"dust_corrected_av", cfg->dust_corrected_av}};
170176
} else if (
171177
const MonofluidComplete *cfg = std::get_if<MonofluidComplete>(&current_mode)) {
172178
j = {{"type", "monofluid_complete"}, {"ndust", cfg->ndust}};
@@ -186,7 +192,8 @@ namespace shammodels::sph {
186192
j.at("C_1_fluid").get<Tscal>(),
187193
j.at("C_drift").get<Tscal>(),
188194
j.at("cfl_density_threshold").get<Tscal>(),
189-
j.at("ensure_s_j_positivity").get<bool>());
195+
j.at("ensure_s_j_positivity").get<bool>(),
196+
j.value("dust_corrected_av", false));
190197
} else if (type == "monofluid_complete") {
191198
set_monofluid_complete(j.at("ndust").get<u32>());
192199
} else {
@@ -198,6 +205,13 @@ namespace shammodels::sph {
198205
return is_monofluid_tva(); // S_j = sqrt(\rho \epsilon_j)
199206
}
200207

208+
inline bool should_use_dust_av() {
209+
if (!is_monofluid_tva()) {
210+
return false;
211+
}
212+
return get_monofluid_tva().dust_corrected_av;
213+
}
214+
201215
inline bool has_epsilon_field() {
202216
return bool(std::get_if<MonofluidComplete>(&current_mode));
203217
}

‎src/shammodels/sph/include/shammodels/sph/modules/NodeUpdateDerivsVaryingAlphaAVDustTVA.hpp‎

Lines changed: 80 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,80 @@
1+
// -------------------------------------------------------//
2+
//
3+
// SHAMROCK code for hydrodynamics
4+
// Copyright (c) 2021-2026 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 NodeUpdateDerivsVaryingAlphaAVDustTVA.hpp
14+
* @author Timothée David--Cléris (tim.shamrock@proton.me)
15+
* @brief
16+
*
17+
*/
18+
19+
#include "shambackends/vec.hpp"
20+
#include "shammodels/sph/solvergraph/NeighCache.hpp"
21+
#include "shamrock/solvergraph/IFieldSpan.hpp"
22+
#include "shamrock/solvergraph/INode.hpp"
23+
#include "shamrock/solvergraph/Indexes.hpp"
24+
#include "shamrock/solvergraph/ScalarEdge.hpp"
25+
26+
#define NODE_EDGES(X_RO, X_RW) \
27+
/* scalars */ \
28+
X_RO(shamrock::solvergraph::ScalarEdge<Tscal>, gpart_mass) \
29+
X_RO(shamrock::solvergraph::ScalarEdge<Tscal>, alpha_u) \
30+
X_RO(shamrock::solvergraph::ScalarEdge<Tscal>, beta_AV) \
31+
\
32+
/* counts */ \
33+
X_RO(shamrock::solvergraph::Indexes<u32>, part_counts) \
34+
X_RO(shamrock::solvergraph::Indexes<u32>, part_counts_with_ghost) \
35+
\
36+
/* fields */ \
37+
X_RO(shamrock::solvergraph::IFieldSpan<Tvec>, xyz) \
38+
X_RO(shamrock::solvergraph::IFieldSpan<Tscal>, hpart) \
39+
X_RO(shamrock::solvergraph::IFieldSpan<Tvec>, vxyz) \
40+
X_RO(shamrock::solvergraph::IFieldSpan<Tscal>, uint) \
41+
X_RO(shamrock::solvergraph::IFieldSpan<Tscal>, omega) \
42+
X_RO(shamrock::solvergraph::IFieldSpan<Tscal>, pressure) \
43+
X_RO(shamrock::solvergraph::IFieldSpan<Tscal>, cs) \
44+
X_RO(shamrock::solvergraph::IFieldSpan<Tscal>, alpha_AV) \
45+
X_RO(shamrock::solvergraph::IFieldSpan<Tscal>, s_j) \
46+
\
47+
/* neigh */ \
48+
X_RO(shammodels::sph::solvergraph::NeighCache, neigh_cache) \
49+
\
50+
/* outputs */ \
51+
X_RW(shamrock::solvergraph::IFieldSpan<Tvec>, axyz) \
52+
X_RW(shamrock::solvergraph::IFieldSpan<Tscal>, duint)
53+
54+
namespace shammodels::sph::modules {
55+
56+
template<class Tvec, template<class> class SPHKernel>
57+
class NodeUpdateDerivsVaryingAlphaAVDustTVA : public shamrock::solvergraph::INode {
58+
59+
using Tscal = shambase::VecComponent<Tvec>;
60+
61+
static constexpr Tscal kernel_radius = SPHKernel<Tscal>::Rkern;
62+
63+
u32 ndust;
64+
65+
public:
66+
NodeUpdateDerivsVaryingAlphaAVDustTVA(u32 ndust) : ndust(ndust) {}
67+
68+
EXPAND_NODE_EDGES(NODE_EDGES)
69+
70+
void _impl_evaluate_internal();
71+
72+
inline virtual std::string _impl_get_label() const {
73+
return "UpdateDerivsVaryingAlphaAVDustTVA";
74+
};
75+
76+
inline virtual std::string _impl_get_tex() const { return "TODO"; };
77+
};
78+
} // namespace shammodels::sph::modules
79+
80+
#undef NODE_EDGES

‎src/shammodels/sph/src/modules/NodeUpdateDerivsVaryingAlphaAVDustTVA.cpp‎

Lines changed: 213 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,213 @@
1+
// -------------------------------------------------------//
2+
//
3+
// SHAMROCK code for hydrodynamics
4+
// Copyright (c) 2021-2026 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+
/**
11+
* @file NodeUpdateDerivsVaryingAlphaAVDustTVA.cpp
12+
* @author Timothée David--Cléris (tim.shamrock@proton.me)
13+
* @brief
14+
*
15+
*/
16+
17+
#include "shammodels/sph/modules/NodeUpdateDerivsVaryingAlphaAVDustTVA.hpp"
18+
#include "shambackends/kernel_call_distrib.hpp"
19+
#include "shammath/sphkernels.hpp"
20+
#include "shammodels/sph/math/density.hpp"
21+
#include "shammodels/sph/math/forces.hpp"
22+
#include "shammodels/sph/math/q_ab.hpp"
23+
#include "shamrock/patch/PatchDataField.hpp"
24+
25+
template<class Tvec, template<class> class SPHKernel>
26+
struct KernelUpdateDerivsVaryingAlphaAVDustTVA {
27+
using Tscal = shambase::VecComponent<Tvec>;
28+
using Kernel = SPHKernel<Tscal>;
29+
static constexpr Tscal hfactd = Kernel::hfactd;
30+
static constexpr Tscal Rkern = Kernel::Rkern;
31+
static constexpr Tscal Rker2 = Rkern * Rkern;
32+
33+
Tscal pmass;
34+
Tscal alpha_u;
35+
Tscal beta_AV;
36+
u32 ndust;
37+
38+
inline void operator()(
39+
unsigned int id_a,
40+
const Tvec *__restrict xyz,
41+
const Tscal *__restrict hpart,
42+
const Tvec *__restrict vxyz,
43+
const Tscal *__restrict uint,
44+
const Tscal *__restrict omega,
45+
const Tscal *__restrict pressure,
46+
const Tscal *__restrict cs,
47+
const Tscal *__restrict alpha_AV,
48+
const Tscal *__restrict s_j,
49+
shamrock::tree::ObjectCache::ptrs_read ploop_ptrs,
50+
Tvec *__restrict axyz,
51+
Tscal *__restrict duint) const {
52+
53+
using namespace shamrock::sph;
54+
55+
shamrock::tree::ObjectCacheIterator particle_looper(ploop_ptrs);
56+
57+
Tvec xyz_a = xyz[id_a];
58+
Tscal h_a = hpart[id_a];
59+
Tvec vxyz_a = vxyz[id_a];
60+
Tscal u_a = uint[id_a];
61+
Tscal omega_a = omega[id_a];
62+
Tscal P_a = pressure[id_a];
63+
Tscal cs_a = cs[id_a];
64+
Tscal alpha_a = alpha_AV[id_a];
65+
66+
Tscal rho_a = rho_h(pmass, h_a, hfactd);
67+
Tscal rho_a_sq = rho_a * rho_a;
68+
Tscal rho_a_inv = 1. / rho_a;
69+
70+
Tscal epsilon_sum_a = 0;
71+
for (u32 j = 0; j < ndust; j++) {
72+
Tscal s = s_j[id_a * ndust + j];
73+
epsilon_sum_a += s * s / rho_a;
74+
}
75+
76+
Tscal omega_a_rho_a_inv = 1 / (omega_a * rho_a);
77+
78+
Tvec force_pressure = Tvec{0, 0, 0};
79+
Tscal tmpdU_pressure = Tscal{0};
80+
81+
particle_looper.for_each_object(id_a, [&](u32 id_b) {
82+
Tvec dr = xyz_a - xyz[id_b];
83+
Tscal rab2 = sycl::dot(dr, dr);
84+
Tscal h_b = hpart[id_b];
85+
86+
if (rab2 > h_a * h_a * Rker2 && rab2 > h_b * h_b * Rker2) {
87+
return;
88+
}
89+
90+
Tvec vxyz_b = vxyz[id_b];
91+
const Tscal u_b = uint[id_b];
92+
Tscal P_b = pressure[id_b];
93+
Tscal omega_b = omega[id_b];
94+
const Tscal alpha_b = alpha_AV[id_b];
95+
Tscal cs_b = cs[id_b];
96+
97+
Tscal rab = sycl::sqrt(rab2);
98+
99+
Tscal rho_b = rho_h(pmass, h_b, hfactd);
100+
101+
Tscal epsilon_sum_b = 0;
102+
for (u32 j = 0; j < ndust; j++) {
103+
Tscal s = s_j[id_b * ndust + j];
104+
epsilon_sum_b += s * s / rho_b;
105+
}
106+
107+
Tscal Fab_a = Kernel::dW_3d(rab, h_a);
108+
Tscal Fab_b = Kernel::dW_3d(rab, h_b);
109+
110+
Tvec v_ab = vxyz_a - vxyz_b;
111+
112+
Tvec r_ab_unit = dr * sham::inv_sat_positive(rab);
113+
114+
Tscal v_ab_r_ab = sycl::dot(v_ab, r_ab_unit);
115+
Tscal abs_v_ab_r_ab = sycl::fabs(v_ab_r_ab);
116+
117+
Tscal vsig_a = alpha_a * cs_a + beta_AV * abs_v_ab_r_ab;
118+
Tscal vsig_b = alpha_b * cs_b + beta_AV * abs_v_ab_r_ab;
119+
120+
Tscal vsig_u = shamrock::sph::vsig_u(P_a, P_b, rho_a, rho_b);
121+
122+
Tscal qa_ab = shamrock::sph::q_av(rho_a * (1 - epsilon_sum_a), vsig_a, v_ab_r_ab);
123+
Tscal qb_ab = shamrock::sph::q_av(rho_b * (1 - epsilon_sum_b), vsig_b, v_ab_r_ab);
124+
125+
add_to_derivs_sph_artif_visco_cond(
126+
pmass,
127+
rho_a_sq,
128+
omega_a_rho_a_inv,
129+
rho_a_inv,
130+
rho_b,
131+
omega_a,
132+
omega_b,
133+
Fab_a,
134+
Fab_b,
135+
u_a,
136+
u_b,
137+
P_a,
138+
P_b,
139+
alpha_u,
140+
v_ab,
141+
r_ab_unit,
142+
vsig_u,
143+
qa_ab,
144+
qb_ab,
145+
146+
force_pressure,
147+
tmpdU_pressure);
148+
});
149+
150+
axyz[id_a] = force_pressure;
151+
duint[id_a] = tmpdU_pressure;
152+
}
153+
};
154+
155+
template<class Tvec, template<class> class SPHKernel>
156+
void shammodels::sph::modules::NodeUpdateDerivsVaryingAlphaAVDustTVA<Tvec, SPHKernel>::
157+
_impl_evaluate_internal() {
158+
159+
__shamrock_stack_entry();
160+
161+
auto edges = get_edges();
162+
163+
auto &part_counts_with_ghost = edges.part_counts_with_ghost.indexes;
164+
auto &part_counts = edges.part_counts.indexes;
165+
166+
// check that all input edges have the particles with ghosts zones
167+
edges.xyz.check_sizes(part_counts_with_ghost);
168+
edges.hpart.check_sizes(part_counts_with_ghost);
169+
edges.vxyz.check_sizes(part_counts_with_ghost);
170+
edges.uint.check_sizes(part_counts_with_ghost);
171+
edges.omega.check_sizes(part_counts_with_ghost);
172+
edges.pressure.check_sizes(part_counts_with_ghost);
173+
edges.cs.check_sizes(part_counts_with_ghost);
174+
edges.alpha_AV.check_sizes(part_counts_with_ghost);
175+
edges.s_j.check_sizes(part_counts_with_ghost);
176+
177+
// ensure that the output edges are of size part_counts (output without ghosts zones)
178+
edges.axyz.ensure_sizes(part_counts);
179+
edges.duint.ensure_sizes(part_counts);
180+
181+
const Tscal pmass = edges.gpart_mass.value;
182+
const Tscal alpha_u = edges.alpha_u.value;
183+
const Tscal beta_AV = edges.beta_AV.value;
184+
185+
using ComputeKernel = KernelUpdateDerivsVaryingAlphaAVDustTVA<Tvec, SPHKernel>;
186+
187+
// call the kernel for each patches with part_counts.get(id_patch) threads of patch id_patch
188+
sham::distributed_data_kernel_call(
189+
shamsys::instance::get_compute_scheduler_ptr(),
190+
sham::DDMultiRef{
191+
edges.xyz.get_spans(),
192+
edges.hpart.get_spans(),
193+
edges.vxyz.get_spans(),
194+
edges.uint.get_spans(),
195+
edges.omega.get_spans(),
196+
edges.pressure.get_spans(),
197+
edges.cs.get_spans(),
198+
edges.alpha_AV.get_spans(),
199+
edges.s_j.get_spans(),
200+
edges.neigh_cache},
201+
sham::DDMultiRef{edges.axyz.get_spans(), edges.duint.get_spans()},
202+
part_counts,
203+
ComputeKernel{pmass, alpha_u, beta_AV, ndust});
204+
}
205+
206+
using namespace shammath;
207+
template class shammodels::sph::modules::NodeUpdateDerivsVaryingAlphaAVDustTVA<f64_3, M4>;
208+
template class shammodels::sph::modules::NodeUpdateDerivsVaryingAlphaAVDustTVA<f64_3, M6>;
209+
template class shammodels::sph::modules::NodeUpdateDerivsVaryingAlphaAVDustTVA<f64_3, M8>;
210+
211+
template class shammodels::sph::modules::NodeUpdateDerivsVaryingAlphaAVDustTVA<f64_3, C2>;
212+
template class shammodels::sph::modules::NodeUpdateDerivsVaryingAlphaAVDustTVA<f64_3, C4>;
213+
template class shammodels::sph::modules::NodeUpdateDerivsVaryingAlphaAVDustTVA<f64_3, C6>;

0 commit comments

Comments
 (0)