Skip to content

Commit ce5ea1c

Browse files
BenWibkingBenWibking
committed
update
1 parent 1aa865d commit ce5ea1c

2 files changed

Lines changed: 121 additions & 63 deletions

File tree

‎src/io/projection.hpp‎

Lines changed: 103 additions & 63 deletions
Original file line numberDiff line numberDiff line change
@@ -14,12 +14,15 @@
1414
#include <vector>
1515

1616
// AMReX headers
17+
#include "AMReX_MFInterpolater.H"
1718
#include "AMReX_MultiFab.H"
1819
#include "AMReX_MultiFabUtil.H"
1920
#include "AMReX_Orientation.H"
21+
#include "AMReX_iMultiFab.H"
2022
#include "AMReX_SPACE.H"
2123
#include "AMReX_VisMF.H"
2224
#include <AMReX.H>
25+
#include <AMReX_BC_TYPES.H>
2326

2427
// YAML headers
2528
#include <yaml-cpp/yaml.h>
@@ -70,69 +73,106 @@ auto ComputePlaneProjection(amrex::Vector<amrex::MultiFab> const &state_new, con
7073
// compute plane-parallel projection of user_f(i, j, k, state) along the given axis.
7174
const BL_PROFILE("quokka::DiagProjection::computePlaneProjection()");
7275

73-
amrex::Vector<amrex::MultiFab> projections(finest_level + 1);
74-
75-
for (int lev = 0; lev <= finest_level; ++lev) {
76-
auto const &level_ba = state_new[lev].boxArray();
77-
amrex::BoxList bl(level_ba.ixType());
78-
for (int i = 0; i < level_ba.size(); ++i) {
79-
bl.push_back(detail::transform_box_to_2D(dir, level_ba[i]));
80-
}
81-
bl.simplify();
82-
amrex::BoxArray ba2d(std::move(bl));
83-
ba2d.removeOverlap();
84-
amrex::DistributionMapping dm2d(ba2d);
85-
86-
projections[lev].define(ba2d, dm2d, 1, 0);
87-
projections[lev].setVal(0.0);
88-
89-
auto const &state = state_new[lev].const_arrays();
90-
auto const &dx = geom[lev].CellSizeArray();
91-
92-
auto plane_local = amrex::ReduceToPlane<ReduceOp, amrex::Real>(static_cast<int>(dir), geom[lev].Domain(), state_new[lev],
93-
[=] AMREX_GPU_DEVICE(int box_no, int i, int j, int k) -> amrex::Real {
94-
return dx[static_cast<int>(dir)] * user_f(i, j, k, state[box_no]);
95-
});
96-
amrex::ParallelDescriptor::ReduceRealSum(plane_local.dataPtr(), static_cast<int>(plane_local.size()));
97-
98-
auto const plane_arr = plane_local.const_array();
99-
auto const proj_arr = projections[lev].arrays();
100-
amrex::ParallelFor(projections[lev], [=] AMREX_GPU_DEVICE(int bx, int i, int j, int k) {
101-
if (dir == amrex::Direction::x) {
102-
proj_arr[bx](i, j, k) = plane_arr(AMREX_D_DECL(0, i, j));
103-
#if AMREX_SPACEDIM >= 2
104-
} else if (dir == amrex::Direction::y) {
105-
proj_arr[bx](i, j, k) = plane_arr(AMREX_D_DECL(i, 0, j));
106-
#endif
107-
#if AMREX_SPACEDIM == 3
108-
} else if (dir == amrex::Direction::z) {
109-
proj_arr[bx](i, j, k) = plane_arr(AMREX_D_DECL(i, j, 0));
110-
#endif
111-
} else {
112-
proj_arr[bx](i, j, k) = 0.0;
113-
}
114-
});
115-
amrex::Gpu::streamSynchronize();
116-
}
117-
118-
// average down
119-
for (int lev = finest_level - 1; lev >= 0; --lev) {
120-
amrex::IntVect rr_2d{AMREX_D_DECL(1, 1, 1)};
121-
if (dir == amrex::Direction::x) {
122-
rr_2d = amrex::IntVect{AMREX_D_DECL(ref_ratio[lev][1], ref_ratio[lev][2], 1)};
123-
#if AMREX_SPACEDIM >= 2
124-
} else if (dir == amrex::Direction::y) {
125-
rr_2d = amrex::IntVect{AMREX_D_DECL(ref_ratio[lev][0], ref_ratio[lev][2], 1)};
126-
#endif
127-
#if AMREX_SPACEDIM == 3
128-
} else if (dir == amrex::Direction::z) {
129-
rr_2d = amrex::IntVect{AMREX_D_DECL(ref_ratio[lev][0], ref_ratio[lev][1], 1)};
130-
#endif
131-
}
132-
amrex::average_down(projections[lev + 1], projections[lev], 0, 1, rr_2d);
133-
}
134-
135-
return projections;}
76+
amrex::Vector<amrex::MultiFab> projections(finest_level + 1);
77+
amrex::Vector<amrex::Geometry> geom2d(finest_level + 1);
78+
79+
for (int lev = 0; lev <= finest_level; ++lev) {
80+
const amrex::Box box2d = detail::transform_box_to_2D(dir, geom[lev].Domain());
81+
const amrex::RealBox domain2d = detail::transform_realbox_to_2D(dir, geom[lev].ProbDomain());
82+
geom2d[lev] = amrex::Geometry(box2d, &domain2d);
83+
}
84+
85+
for (int lev = 0; lev <= finest_level; ++lev) {
86+
auto const &level_ba = state_new[lev].boxArray();
87+
amrex::BoxList bl(level_ba.ixType());
88+
for (int i = 0; i < level_ba.size(); ++i) {
89+
bl.push_back(detail::transform_box_to_2D(dir, level_ba[i]));
90+
}
91+
bl.simplify();
92+
amrex::BoxArray ba2d(std::move(bl));
93+
ba2d.removeOverlap();
94+
amrex::DistributionMapping dm2d(ba2d);
95+
96+
projections[lev].define(ba2d, dm2d, 1, 0);
97+
projections[lev].setVal(0.0);
98+
99+
amrex::iMultiFab mask;
100+
if (lev == finest_level) {
101+
mask.define(state_new[lev].boxArray(), state_new[lev].DistributionMap(), 1, amrex::IntVect(0));
102+
mask.setVal(1);
103+
} else {
104+
mask = amrex::makeFineMask(state_new[lev], state_new[lev + 1], amrex::IntVect(0), ref_ratio[lev], geom[lev].periodicity(), 1, 0);
105+
}
106+
107+
auto const &state = state_new[lev].const_arrays();
108+
auto const &mask_arr = mask.const_arrays();
109+
auto const &dx = geom[lev].CellSizeArray();
110+
111+
auto plane_local = amrex::ReduceToPlane<ReduceOp, amrex::Real>(static_cast<int>(dir), geom[lev].Domain(), state_new[lev],
112+
[=] AMREX_GPU_DEVICE(int box_no, int i, int j, int k) -> amrex::Real {
113+
if (mask_arr[box_no](i, j, k) == 0) {
114+
return 0.0;
115+
}
116+
return dx[static_cast<int>(dir)] * user_f(i, j, k, state[box_no]);
117+
});
118+
amrex::ParallelDescriptor::ReduceRealSum(plane_local.dataPtr(), static_cast<int>(plane_local.size()));
119+
120+
auto const plane_arr = plane_local.const_array();
121+
auto const proj_arr = projections[lev].arrays();
122+
amrex::ParallelFor(projections[lev], [=] AMREX_GPU_DEVICE(int bx, int i, int j, int k) {
123+
if (dir == amrex::Direction::x) {
124+
proj_arr[bx](i, j, k) = plane_arr(AMREX_D_DECL(0, i, j));
125+
#if AMREX_SPACEDIM >= 2
126+
} else if (dir == amrex::Direction::y) {
127+
proj_arr[bx](i, j, k) = plane_arr(AMREX_D_DECL(i, 0, j));
128+
#endif
129+
#if AMREX_SPACEDIM == 3
130+
} else if (dir == amrex::Direction::z) {
131+
proj_arr[bx](i, j, k) = plane_arr(AMREX_D_DECL(i, j, 0));
132+
#endif
133+
} else {
134+
proj_arr[bx](i, j, k) = 0.0;
135+
}
136+
});
137+
amrex::Gpu::streamSynchronize();
138+
}
139+
140+
amrex::Vector<amrex::MultiFab> projections_accum(finest_level + 1);
141+
for (int lev = 0; lev <= finest_level; ++lev) {
142+
projections_accum[lev].define(projections[lev].boxArray(), projections[lev].DistributionMap(), 1, 0);
143+
projections_accum[lev].ParallelCopy(projections[lev], 0, 0, 1, 0, 0);
144+
}
145+
146+
for (int lev = 0; lev < finest_level; ++lev) {
147+
amrex::IntVect rr_2d{AMREX_D_DECL(1, 1, 1)};
148+
if (dir == amrex::Direction::x) {
149+
rr_2d = amrex::IntVect{AMREX_D_DECL(ref_ratio[lev][1], ref_ratio[lev][2], 1)};
150+
#if AMREX_SPACEDIM >= 2
151+
} else if (dir == amrex::Direction::y) {
152+
rr_2d = amrex::IntVect{AMREX_D_DECL(ref_ratio[lev][0], ref_ratio[lev][2], 1)};
153+
#endif
154+
#if AMREX_SPACEDIM == 3
155+
} else if (dir == amrex::Direction::z) {
156+
rr_2d = amrex::IntVect{AMREX_D_DECL(ref_ratio[lev][0], ref_ratio[lev][1], 1)};
157+
#endif
158+
}
159+
160+
amrex::MultiFab coarse_refined(projections[lev + 1].boxArray(), projections[lev + 1].DistributionMap(), 1, 0);
161+
coarse_refined.setVal(0.0);
162+
163+
amrex::Vector<amrex::BCRec> bcs(1);
164+
for (int d = 0; d < AMREX_SPACEDIM; ++d) {
165+
bcs[0].setLo(d, amrex::BCType::int_dir);
166+
bcs[0].setHi(d, amrex::BCType::int_dir);
167+
}
168+
amrex::MFPCInterp interp;
169+
interp.interp(projections_accum[lev], 0, coarse_refined, 0, 1, amrex::IntVect(0), geom2d[lev], geom2d[lev + 1],
170+
geom2d[lev + 1].Domain(), rr_2d, bcs, 0);
171+
amrex::MultiFab::Add(projections_accum[lev + 1], coarse_refined, 0, 0, 1, 0);
172+
}
173+
174+
return projections_accum;
175+
}
136176

137177
void WriteProjection(amrex::Direction dir, std::unordered_map<std::string, amrex::Vector<amrex::MultiFab>> const &proj,
138178
amrex::Vector<amrex::Geometry> const &geom, amrex::Vector<amrex::IntVect> const &ref_ratio, amrex::Real time, int istep,

‎src/problems/ShockCloud/testShockCloud.cpp‎

Lines changed: 18 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -646,6 +646,24 @@ auto QuokkaSimulation<ShockCloud>::ComputeProjections(const amrex::Direction dir
646646
return (H_mass_fraction * rho_wind) / m_H;
647647
});
648648

649+
// diagnostic: difference between total and partials
650+
proj["nH_residual"] = quokka::diagnostics::ComputePlaneProjection<amrex::ReduceOpSum>(
651+
state_new_cc_, finestLevel(), geom, ref_ratio, dir, [=] AMREX_GPU_DEVICE(int i, int j, int k, amrex::Array4<const Real> const &state) noexcept {
652+
Real const rho = state(i, j, k, HydroSystem<ShockCloud>::density_index);
653+
Real const rho_cloud = state(i, j, k, HydroSystem<ShockCloud>::scalar0_index + 1);
654+
Real const rho_wind = state(i, j, k, HydroSystem<ShockCloud>::scalar0_index + 2);
655+
Real const residual = rho - (rho_cloud + rho_wind);
656+
return (H_mass_fraction * residual) / m_H;
657+
});
658+
659+
if (amrex::ParallelDescriptor::IOProcessor()) {
660+
for (int lev = 0; lev <= finestLevel(); ++lev) {
661+
const amrex::Real max_residual = proj["nH_residual"][lev].norminf(0, 0, true);
662+
amrex::Print() << "ShockCloud diag dir=" << quokka::diagnostics::detail::direction_to_string(dir) << " lev=" << lev
663+
<< " max|nH_residual|=" << max_residual << "\n";
664+
}
665+
}
666+
649667
return proj;
650668
}
651669

0 commit comments

Comments
 (0)