Peano
Loading...
Searching...
No Matches
main.cpp
Go to the documentation of this file.
1//#define DISABLE_TRACING
2
5#include <spacetree/sort.h>
6
7#include <lang/trace.h>
8#include <io/writer/dma.h>
9#include <vtk/writer/point.h>
10
12
13#include "align.h"
14#include "hydro_part.h"
15#include "hydro_iact.h"
16
17#include "setup.h"
18
19#define dt_therm DT_DRIFT
20#define dt_kick_therm DT_KICK_HYDRO
21
23
25 float h_max = 0;
26 float h_max_active = 0;
28
30 .eta_neighbours = 1.2348,
31 .h_tolerance = 1e-4,
32 .h_max = 4.0 / N_PART,
33 .h_min = 1e-5,
34 .max_smoothing_iterations = 30,
35 .use_mass_weighted_num_ngb = 0,
36 .minimal_internal_energy = 0.,
37 .viscosity = {
39 },
40};
41
42const auto COSMO = cosmology{
43 .a2_inv = 1.0,
44 .a_factor_internal_energy = 1.0,
45 .a_factor_Balsara_eps = 1.0,
46 .H = 0.0,
47};
48
49#include "runner_do_ghost.h"
50
51auto GetDefaultParticle = [](const Geo::Point &pos) {
52 auto v = pos - 0.5;
53 auto v_norm = std::sqrt(v[0] * v[0] + v[1] * v[1]);
54 v[0] /= v_norm;
55 v[1] /= v_norm;
56 auto p = part{
57 ._x = pos,
58 .v = {(float) -v[0], (float) -v[1], 0},
59 .mass = 1.0 / N_PART / N_PART,
60 .h = 4.f / N_PART,
61 .u = 1.0e-6 / (5.0 / 3 - 1),
62 };
63 hydro_first_init_part(&p, &p.xpart);
64 hydro_convert_quantities(&p, &p.xpart, &COSMO, &HYDRO_PROPS, nullptr);
65 return p;
66};
67
68auto createTree() {
69 auto tree = Tree({Geo::Region::fromInclusiveBounds(Geo::Point(0), Geo::Point(1))});
70 tree.split(3);
71
72 tree.onNode(InitialConditionsTask<Tree, GetDefaultParticle>(Geo::Region::fromInclusiveBounds(Geo::Point(0.005), Geo::Point(0.995)), N_PART), 8u);
73
74 return tree;
75}
76
77auto ReferenceDensityKernel = [](Tree *leaf) {
78 for (int i = 0; i < leaf->items.size(); i++) {
79 auto &localPart = leaf->items[i];
80 hydro_init_part(&localPart, nullptr);
81 }
82
83 do_density(leaf);
84
86};
87
88auto ReferenceHydroKernel = [](Tree *supernode, Tree *leaf) {
89 constexpr float a = 1.0f; // SWIFT parameter, ignored by the Minimal SPH scheme
90 constexpr float H = 1.0f; // SWIFT parameter, ignored by the Minimal SPH scheme
91
92 if (leaf->items.empty()) [[unlikely]] return;
93
94 auto activeRegion = Tree::Geo::Region(leaf->range.low - 0.001, leaf->range.high + 0.001);
95 auto *activeRegionParent = leaf->getRegionParent(activeRegion, true);
96 auto activeRegionIter = Tree::InsideRegionLeafsIterator(activeRegionParent, activeRegion);
97
98 auto nonsymIact = [](Tree * __restrict__ leaf, Tree * __restrict__ activeCell) {
99 for (auto localIdx = 0; localIdx < leaf->items.size(); localIdx++) {
100 auto *pi = &leaf->items[localIdx];
101 auto hi = pi->h;
102 auto hig = hi * kernel_gamma;
103 auto hig2 = hig * hig;
104
105 if (!activeCell->range.overlaps({pi->_x - hig, pi->_x + hig})) continue;
106
107 for (auto activeIdx = 0; activeIdx < activeCell->items.size(); activeIdx++) {
108 auto *pj = &activeCell->items[activeIdx];
109 auto r2 = 0.0f;
110 float dx[3] = {0};
111
112 for (int k = 0; k < Geo::D::N; k++) {
113 dx[k] = pi->_x[k] - pj->_x[k];
114 r2 += dx[k] * dx[k];
115 }
116
117 auto hjg = pj->h * kernel_gamma;
118 auto hjg2 = hjg * hjg;
119
120 /* Hit or miss? */
121 if (r2 < hig2 | r2 < hjg2) {
122
123 /* Interact */
124 runner_iact_nonsym_force(r2, dx, hi, pj->h, pi, pj, a, H);
125 }
126 }
127 }
128 };
129
130 auto symIact = [](Tree * __restrict__ leaf, Tree * __restrict__ activeCell) {
131 for (auto localIdx = 0; localIdx < leaf->items.size(); localIdx++) {
132 auto *pi = &leaf->items[localIdx];
133 auto hi = pi->h;
134 auto hig = hi * kernel_gamma;
135 auto hig2 = hig * hig;
136
137 if (!activeCell->range.overlaps({pi->_x - hig, pi->_x + hig})) continue;
138
139 for (auto activeIdx = 0; activeIdx < activeCell->items.size(); activeIdx++) {
140 auto *pj = &activeCell->items[activeIdx];
141 auto r2 = 0.0f;
142 float dx[3] = {0};
143
144 for (int k = 0; k < Geo::D::N; k++) {
145 dx[k] = pi->_x[k] - pj->_x[k];
146 r2 += dx[k] * dx[k];
147 }
148
149 auto hjg = pj->h * kernel_gamma;
150 auto hjg2 = hjg * hjg;
151
152 assert(r2 != 0)
153
154 /* Hit or miss? */
155 if (r2 < hig2 | r2 < hjg2) {
156
157 /* Interact */
158 runner_iact_force(r2, dx, hi, pj->h, pi, pj, a, H);
159 }
160 }
161 }
162 };
163
164 auto selfIact = [](Tree *leaf) {
165 for (auto localIdx = 0; localIdx < leaf->items.size(); localIdx++) {
166 auto *pi = &leaf->items[localIdx];
167 auto hi = pi->h;
168 auto hig = hi * kernel_gamma;
169 auto hig2 = hig * hig;
170
171 for (auto activeIdx = localIdx + 1; activeIdx < leaf->items.size(); activeIdx++) {
172 auto *pj = &leaf->items[activeIdx];
173 auto r2 = 0.0f;
174 float dx[3] = {0};
175
176 for (int k = 0; k < Geo::D::N; k++) {
177 dx[k] = pi->_x[k] - pj->_x[k];
178 r2 += dx[k] * dx[k];
179 }
180
181 assert(r2 != 0)
182
183 auto hjg = pj->h * kernel_gamma;
184 auto hjg2 = hjg * hjg;
185
186 /* Hit or miss? */
187 if (r2 < hig2 | r2 < hjg2) {
188
189 /* Interact */
190 runner_iact_force(r2, dx, hi, pj->h, pi, pj, a, H);
191 }
192 }
193 }
194 };
195
196 auto *activeCell = (Tree *) nullptr;
197 while ((activeCell = activeRegionIter.next().value_or(nullptr))) {
198 auto isActiveCellWithinSupernode = supernode->range.contains(activeCell->range);
199 if (!isActiveCellWithinSupernode) [[unlikely]] {
200 nonsymIact(leaf, activeCell);
201 continue;
202 }
203
204 auto isActiveCellLocalCell = leaf == activeCell;
205 if (isActiveCellLocalCell) [[unlikely]] {
206 selfIact(leaf);
207 continue;
208 }
209
210 auto isLocalLowerAddressThanActive = (u64) leaf < (u64) activeCell;
211 if (isLocalLowerAddressThanActive) symIact(leaf, activeCell);
212 }
213};
214
215auto Kick1Kernel = [](Tree *leaf) {
216 for (int i = 0; i < leaf->items.size(); i++) {
217 auto &localPart = leaf->items[i];
218
219 localPart.xpart.v_full[0] += localPart.a_hydro[0] * DT_KICK_HYDRO;
220 localPart.xpart.v_full[1] += localPart.a_hydro[1] * DT_KICK_HYDRO;
221
222 hydro_kick_extra(&localPart, &localPart.xpart, dt_kick_therm, 0, 0, 0, 0, &COSMO, &HYDRO_PROPS, nullptr);
223 }
224};
225
226auto Kick2Kernel = [](Tree *leaf) {
227 for (int i = 0; i < leaf->items.size(); i++) {
228 auto &localPart = leaf->items[i];
229 hydro_end_force(&localPart, &COSMO);
230 }
231
232 for (int i = 0; i < leaf->items.size(); i++) {
233 auto &localPart = leaf->items[i];
234
235 localPart.xpart.v_full[0] += localPart.a_hydro[0] * DT_KICK_HYDRO;
236 localPart.xpart.v_full[1] += localPart.a_hydro[1] * DT_KICK_HYDRO;
237
238 hydro_kick_extra(&localPart, &localPart.xpart, dt_kick_therm, 0, 0, 0, 0, &COSMO, &HYDRO_PROPS, nullptr);
239
240 hydro_reset_predicted_values(&localPart, &localPart.xpart, &COSMO, nullptr);
241 }
242};
243
244auto DriftKernel = [](Tree *leaf) {
245 for (int i = 0; i < leaf->items.size(); i++) {
246 auto &localPart = leaf->items[i];
247
248 localPart._x[0] += localPart.xpart.v_full[0] * DT_DRIFT;
249 localPart._x[1] += localPart.xpart.v_full[1] * DT_DRIFT;
250
251 /* Predict velocities (for hydro terms) */
252 localPart.v[0] += localPart.a_hydro[0] * DT_KICK_HYDRO;
253 localPart.v[1] += localPart.a_hydro[1] * DT_KICK_HYDRO;
254
255 hydro_predict_extra(&localPart, &localPart.xpart, DT_DRIFT, dt_therm, 0, &COSMO, &HYDRO_PROPS, nullptr, nullptr);
256
257 /* Compute offsets since last cell construction */
258 for (int k = 0; k < 2; k++) {
259 const float dx = localPart.xpart.v_full[k] * DT_DRIFT;
260 localPart.xpart.x_diff[k] -= dx;
261 localPart.xpart.x_diff_sort[k] -= dx;
262 }
263 }
264};
265
266void dumpData(Tree &q0, int i) {
267 char filePathBuf[128];
268 sprintf(filePathBuf, "./noh2d_%04d.vtk", i);
269 auto fd = open(filePathBuf, DMA_FLAGS, 0644);
270 assert(fd > 0)
271
272 auto fileWriter = DefaultDmaWriter(fd);
273 auto particlesWriter = VtkPointBinaryWriter<decltype(fileWriter), double>(&fileWriter);
274
275 q0.onNode([&](auto *node) {
276 for (auto &p: node->items) particlesWriter.add(p._x[0], p._x[1], 0);
277 }, (u32) 1);
278
279 particlesWriter.addPointDataHeader();
280
281 particlesWriter.addScalarPointDataHeader("id");
282
283 q0.onNode([&](auto *node) {
284 for (auto &p: node->items) particlesWriter.addScalar(p.id);
285 }, (u32) 1);
286
287 particlesWriter.addScalarPointDataHeader("density");
288
289 q0.onNode([&](auto *node) {
290 for (auto &p: node->items) particlesWriter.addScalar(p.rho);
291 }, (u32) 1);
292
293 particlesWriter.addScalarPointDataHeader("pressure");
294
295 q0.onNode([&](auto *node) {
296 for (auto &p: node->items) particlesWriter.addScalar(p.force.pressure);
297 }, (u32) 1);
298
299 particlesWriter.addScalarPointDataHeader("mass");
300
301 q0.onNode([&](auto *node) {
302 for (auto &p: node->items) particlesWriter.addScalar(p.mass);
303 }, (u32) 1);
304
305 particlesWriter.addScalarPointDataHeader("smoothing_length");
306
307 q0.onNode([&](auto *node) {
308 for (auto &p: node->items) particlesWriter.addScalar(p.h);
309 }, (u32) 1);
310
311 particlesWriter.addScalarPointDataHeader("internal_energy");
312
313 q0.onNode([&](auto *node) {
314 for (auto &p: node->items) particlesWriter.addScalar(p.u);
315 }, (u32) 1);
316
317 particlesWriter.finish();
318 close(fd);
319}
320
321int main() {
322 TRACE("Initial conditions", auto tree = createTree();)
323 TRACE("Initial Density", tree.onNode([](auto *node) { node->onNode(ReferenceDensityKernel, 1u); },
325 TRACE("Initial Hydro", tree.onNode([](auto *node) { node->onNode([supernode = node](auto *node) {
326 ReferenceHydroKernel(supernode, node);
327 }, 1u); }, Tree::LeafOrLevelIterator<LEVEL>(&tree));)
328
329 for (int i = 0; i < STEPS; i++) {
330 TRACE("Timestep", {
331 TRACE("Data dump",
332 if (i % 10 == 0 && STEPS_BETWEEN_DUMPS != 0) dumpData(tree, i / STEPS_BETWEEN_DUMPS);
333 )
334
335 TRACE("Kick2 + Drift + Sort", {
336 auto strayParticles = (SpacetreeSorter<Tree, [](auto *node){ Kick2Kernel(node); DriftKernel(node); }>(&tree).run());
337 assert(strayParticles.empty())
338 })
339
340 TRACE("Kick1 + Density", tree.onNode([](auto *node) { node->onNode([](auto *node) {
341 Kick1Kernel(node);
342 ReferenceDensityKernel(node);
343 }, 1u); }, Tree::LeafOrLevelIterator<LEVEL>(&tree));)
344
345 TRACE("Hydro", tree.onNode([](auto *node) { node->onNode([supernode = node](auto *node) {
346 ReferenceHydroKernel(supernode, node);
347 }, 1u); }, Tree::LeafOrLevelIterator<LEVEL>(&tree));)
348 })
349 }
350
351 return 0;
352}
#define assert(...)
Definition LinuxAMD.h:28
auto createTree()
Definition main.cpp:68
#define dt_kick_therm
Definition main.cpp:20
auto Kick2Kernel
Definition main.cpp:226
auto ReferenceHydroKernel
Definition main.cpp:88
auto Kick1Kernel
Definition main.cpp:215
Spacetree< Geo, part, 64, 2, ThreadPoolTaskBackend > Tree
Definition main.cpp:22
#define dt_therm
Definition main.cpp:19
auto GetDefaultParticle
Definition main.cpp:51
auto ReferenceDensityKernel
Definition main.cpp:77
void dumpData(Tree &q0, int i)
Definition main.cpp:266
struct HYDRO_STATS_T HYDRO_STATS
const auto COSMO
Definition main.cpp:42
const auto HYDRO_PROPS
Definition main.cpp:29
int main()
Definition main.cpp:321
auto DriftKernel
Definition main.cpp:244
Geo::Region range
Definition linked_grid.h:51
DmaWriter< 1024 *1024 *1, 4, 512 > DefaultDmaWriter
Definition dma.h:271
constexpr i32 DMA_FLAGS
Definition dma.h:10
const float dx[3]
Definition hydro_iact.h:54
const float const float hi
Definition hydro_iact.h:54
const float const float const float struct part *restrict struct part *restrict pj
Definition hydro_iact.h:55
const float const float const float struct part *restrict struct part *restrict const float const float H
Definition hydro_iact.h:56
const float const float const float struct part *restrict struct part *restrict const float a
Definition hydro_iact.h:55
#define hydro_props_default_viscosity_alpha
void do_density(Tree *leaf)
void runner_do_ghost(Tree *leaf, const hydro_props *hydro_properties, HYDRO_STATS_T *hydro_stats)
#define STEPS
Definition setup_big.h:6
#define N_PART
Definition setup_big.h:5
#define STEPS_BETWEEN_DUMPS
Definition setup_big.h:8
#define DT_KICK_HYDRO
Definition setup_big.h:4
#define DT_DRIFT
Definition setup_big.h:3
Point< dim, fp > Point
Definition geometry.h:242
float h_max_active
Definition main.cpp:26
float h_max
Definition main.cpp:25
LeafOrLevelIterator< Spacetree, LEVEL, IterOrder > LeafOrLevelIterator
Definition spacetree.h:69
void onNode(Callback F, u32 numThreads=0)
Definition spacetree.h:196
InsideRegionLeafsIterator< Spacetree, IterOrder > InsideRegionLeafsIterator
Definition spacetree.h:65
Cosmological parameters.
Definition cosmology.h:34
double a2_inv
Definition cosmology.h:43
Contains all the constants and parameters of the hydro scheme.
Particle fields for the SPH particles.
Definition hydro_part.h:99
float h
Definition hydro_part.h:120
Geo::Point _x
Definition hydro_part.h:108
#define TRACE(name, function_call)
Definition trace.h:15
std::uint32_t u32
Definition type.h:11
std::uint64_t u64
Definition type.h:13