53 auto v_norm = std::sqrt(v[0] * v[0] + v[1] * v[1]);
58 .v = {(float) -v[0], (
float) -v[1], 0},
61 .u = 1.0e-6 / (5.0 / 3 - 1),
63 hydro_first_init_part(&p, &p.xpart);
89 constexpr float a = 1.0f;
90 constexpr float H = 1.0f;
92 if (leaf->items.empty()) [[unlikely]]
return;
94 auto activeRegion = Tree::Geo::Region(leaf->range.low - 0.001, leaf->range.high + 0.001);
95 auto *activeRegionParent = leaf->getRegionParent(activeRegion,
true);
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];
102 auto hig =
hi * kernel_gamma;
103 auto hig2 = hig * hig;
105 if (!activeCell->range.overlaps({pi->_x - hig, pi->_x + hig}))
continue;
107 for (
auto activeIdx = 0; activeIdx < activeCell->items.size(); activeIdx++) {
108 auto *
pj = &activeCell->items[activeIdx];
112 for (
int k = 0; k < Geo::D::N; k++) {
117 auto hjg =
pj->
h * kernel_gamma;
118 auto hjg2 = hjg * hjg;
121 if (r2 < hig2 | r2 < hjg2) {
124 runner_iact_nonsym_force(r2,
dx,
hi,
pj->
h, pi,
pj,
a,
H);
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];
134 auto hig =
hi * kernel_gamma;
135 auto hig2 = hig * hig;
137 if (!activeCell->range.overlaps({pi->_x - hig, pi->_x + hig}))
continue;
139 for (
auto activeIdx = 0; activeIdx < activeCell->items.size(); activeIdx++) {
140 auto *
pj = &activeCell->items[activeIdx];
144 for (
int k = 0; k < Geo::D::N; k++) {
149 auto hjg =
pj->
h * kernel_gamma;
150 auto hjg2 = hjg * hjg;
155 if (r2 < hig2 | r2 < hjg2) {
158 runner_iact_force(r2,
dx,
hi,
pj->
h, pi,
pj,
a,
H);
164 auto selfIact = [](
Tree *leaf) {
165 for (
auto localIdx = 0; localIdx < leaf->items.size(); localIdx++) {
166 auto *pi = &leaf->items[localIdx];
168 auto hig =
hi * kernel_gamma;
169 auto hig2 = hig * hig;
171 for (
auto activeIdx = localIdx + 1; activeIdx < leaf->items.size(); activeIdx++) {
172 auto *
pj = &leaf->items[activeIdx];
176 for (
int k = 0; k < Geo::D::N; k++) {
183 auto hjg =
pj->
h * kernel_gamma;
184 auto hjg2 = hjg * hjg;
187 if (r2 < hig2 | r2 < hjg2) {
190 runner_iact_force(r2,
dx,
hi,
pj->
h, pi,
pj,
a,
H);
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);
204 auto isActiveCellLocalCell = leaf == activeCell;
205 if (isActiveCellLocalCell) [[unlikely]] {
210 auto isLocalLowerAddressThanActive = (
u64) leaf < (
u64) activeCell;
211 if (isLocalLowerAddressThanActive) symIact(leaf, activeCell);
227 for (
int i = 0; i < leaf->items.size(); i++) {
228 auto &localPart = leaf->items[i];
229 hydro_end_force(&localPart, &
COSMO);
232 for (
int i = 0; i < leaf->items.size(); i++) {
233 auto &localPart = leaf->items[i];
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;
240 hydro_reset_predicted_values(&localPart, &localPart.xpart, &
COSMO,
nullptr);
245 for (
int i = 0; i < leaf->items.size(); i++) {
246 auto &localPart = leaf->items[i];
248 localPart._x[0] += localPart.xpart.v_full[0] *
DT_DRIFT;
249 localPart._x[1] += localPart.xpart.v_full[1] *
DT_DRIFT;
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;
267 char filePathBuf[128];
268 sprintf(filePathBuf,
"./noh2d_%04d.vtk", i);
269 auto fd = open(filePathBuf,
DMA_FLAGS, 0644);
275 q0.
onNode([&](
auto *node) {
276 for (
auto &p: node->items) particlesWriter.add(p._x[0], p._x[1], 0);
279 particlesWriter.addPointDataHeader();
281 particlesWriter.addScalarPointDataHeader(
"id");
283 q0.
onNode([&](
auto *node) {
284 for (
auto &p: node->items) particlesWriter.addScalar(p.id);
287 particlesWriter.addScalarPointDataHeader(
"density");
289 q0.
onNode([&](
auto *node) {
290 for (
auto &p: node->items) particlesWriter.addScalar(p.rho);
293 particlesWriter.addScalarPointDataHeader(
"pressure");
295 q0.
onNode([&](
auto *node) {
296 for (
auto &p: node->items) particlesWriter.addScalar(p.force.pressure);
299 particlesWriter.addScalarPointDataHeader(
"mass");
301 q0.
onNode([&](
auto *node) {
302 for (
auto &p: node->items) particlesWriter.addScalar(p.mass);
305 particlesWriter.addScalarPointDataHeader(
"smoothing_length");
307 q0.
onNode([&](
auto *node) {
308 for (
auto &p: node->items) particlesWriter.addScalar(p.h);
311 particlesWriter.addScalarPointDataHeader(
"internal_energy");
313 q0.
onNode([&](
auto *node) {
314 for (
auto &p: node->items) particlesWriter.addScalar(p.u);
317 particlesWriter.finish();
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);
329 for (
int i = 0;
i <
STEPS;
i++) {
335 TRACE(
"Kick2 + Drift + Sort", {
337 assert(strayParticles.empty())
340 TRACE(
"Kick1 + Density",
tree.onNode([](
auto *node) { node->onNode([](auto *node) {
342 ReferenceDensityKernel(node);
345 TRACE(
"Hydro",
tree.onNode([](
auto *node) { node->onNode([supernode = node](
auto *node) {
auto ReferenceHydroKernel
Spacetree< Geo, part, 64, 2, ThreadPoolTaskBackend > Tree
auto ReferenceDensityKernel
void dumpData(Tree &q0, int i)
struct HYDRO_STATS_T HYDRO_STATS