Skip to content

Commit 6365eb1

Browse files
author
fournier2
committed
Adding cubic AMR
1 parent dc37e7e commit 6365eb1

4 files changed

Lines changed: 170 additions & 4 deletions

File tree

‎src/hydro/hydro.cpp‎

Lines changed: 17 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -816,6 +816,23 @@ std::shared_ptr<StateDescriptor> Initialize(ParameterInput *pin) {
816816
"refinement/maxdensity_refine_above");
817817
pkg->AddParam<Real>("refinement/maxdensity_deref_below", deref_below);
818818
pkg->AddParam<Real>("refinement/maxdensity_refine_above", refine_above);
819+
} else if (refine_str == "cubic") {
820+
pkg->CheckRefinementBlock = refinement::other::Cubic;
821+
const auto active = pin->GetOrAddBoolean("refinement", "active", false);
822+
const auto refinement_width =
823+
pin->GetOrAddReal("refinement", "refinement_width", 0.0);
824+
const auto refinement_center_x =
825+
pin->GetOrAddReal("refinement", "refinement_center_x", 0.0);
826+
const auto refinement_center_y =
827+
pin->GetOrAddReal("refinement", "refinement_center_y", 0.0);
828+
const auto refinement_center_z =
829+
pin->GetOrAddReal("refinement", "refinement_center_z", 0.0);
830+
// Adding parameters
831+
pkg->AddParam<>("refinement/active", active);
832+
pkg->AddParam<>("refinement/refinement_width", refinement_width);
833+
pkg->AddParam<>("refinement/refinement_center_x", refinement_center_x);
834+
pkg->AddParam<>("refinement/refinement_center_y", refinement_center_y);
835+
pkg->AddParam<>("refinement/refinement_center_z", refinement_center_z);
819836
} else if (refine_str == "user") {
820837
pkg->CheckRefinementBlock = Hydro::ProblemCheckRefinementBlock;
821838
}

‎src/refinement/other.cpp‎

Lines changed: 65 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -40,5 +40,70 @@ parthenon::AmrTag MaxDensity(MeshBlockData<Real> *rc) {
4040
return parthenon::AmrTag::same;
4141
}
4242

43+
parthenon::AmrTag Cubic(MeshBlockData<Real> *rc) {
44+
45+
auto pmb = rc->GetBlockPointer();
46+
47+
// Check if refinement is active
48+
const bool active = pmb->packages.Get("Hydro")->Param<bool>("refinement/active");
49+
if (!active) {
50+
return parthenon::AmrTag::same; // Skip refinement if inactive
51+
}
52+
53+
const Real refinement_width =
54+
pmb->packages.Get("Hydro")->Param<Real>("refinement/refinement_width");
55+
const Real half_width = refinement_width / 2.0;
56+
57+
// Retrieve center coordinates
58+
const Real center_x =
59+
pmb->packages.Get("Hydro")->Param<Real>("refinement/refinement_center_x");
60+
const Real center_y =
61+
pmb->packages.Get("Hydro")->Param<Real>("refinement/refinement_center_y");
62+
const Real center_z =
63+
pmb->packages.Get("Hydro")->Param<Real>("refinement/refinement_center_z");
64+
65+
// Retrieve bounds
66+
IndexRange ib = pmb->cellbounds.GetBoundsI(IndexDomain::interior);
67+
IndexRange jb = pmb->cellbounds.GetBoundsJ(IndexDomain::interior);
68+
IndexRange kb = pmb->cellbounds.GetBoundsK(IndexDomain::interior);
69+
70+
// Fast bounding box check (block-level)
71+
auto &coords = pmb->coords;
72+
const Real x_min = coords.Xc<1>(ib.s);
73+
const Real x_max = coords.Xc<1>(ib.e);
74+
const Real y_min = coords.Xc<2>(jb.s);
75+
const Real y_max = coords.Xc<2>(jb.e);
76+
const Real z_min = coords.Xc<3>(kb.s);
77+
const Real z_max = coords.Xc<3>(kb.e);
78+
79+
// Check if block is completely outside the cubic region
80+
if (x_min > center_x + half_width || x_max < center_x - half_width ||
81+
y_min > center_y + half_width || y_max < center_y - half_width ||
82+
z_min > center_z + half_width || z_max < center_z - half_width) {
83+
return parthenon::AmrTag::same; // Fully outside, no refinement needed
84+
}
85+
86+
// If the block intersects the cubic region, perform detailed check
87+
bool inside_cubic_region = false;
88+
89+
pmb->par_reduce(
90+
"cubic check refinement", kb.s, kb.e, jb.s, jb.e, ib.s, ib.e,
91+
KOKKOS_LAMBDA(const int k, const int j, const int i, bool &inside) {
92+
const Real x = coords.Xc<1>(i);
93+
const Real y = coords.Xc<2>(j);
94+
const Real z = coords.Xc<3>(k);
95+
96+
if (std::abs(x - center_x) <= half_width &&
97+
std::abs(y - center_y) <= half_width &&
98+
std::abs(z - center_z) <= half_width) {
99+
inside = true;
100+
}
101+
},
102+
103+
Kokkos::LOr<bool>(inside_cubic_region));
104+
105+
return inside_cubic_region ? parthenon::AmrTag::refine : parthenon::AmrTag::same;
106+
}
107+
43108
} // namespace other
44109
} // namespace refinement

‎src/refinement/refinement.hpp‎

Lines changed: 2 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -20,7 +20,8 @@ AmrTag VelocityGradient(MeshBlockData<Real> *rc);
2020
} // namespace gradient
2121
namespace other {
2222
parthenon::AmrTag MaxDensity(MeshBlockData<Real> *rc);
23-
}
23+
parthenon::AmrTag Cubic(MeshBlockData<Real> *rc);
24+
} // namespace other
2425

2526
} // namespace refinement
2627

‎src/tracers/tracers.cpp‎

Lines changed: 86 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -336,9 +336,19 @@ std::shared_ptr<StateDescriptor> Initialize(ParameterInput *pin) {
336336
// packages.
337337
const bool mhd = pin->GetString("hydro", "fluid") == "glmmhd";
338338

339-
PARTHENON_REQUIRE_THROWS(
340-
pin->GetString("parthenon/mesh", "refinement") != "adaptive",
341-
"Tracers/swarms currently only supported on non-adaptive meshes.");
339+
// Check if tracers/swarms are compatible with the mesh refinement type
340+
const std::string mesh_refinement = pin->GetString("parthenon/mesh", "refinement");
341+
const bool is_adaptive = (mesh_refinement == "adaptive");
342+
343+
bool is_cubic_refinement = false;
344+
if (is_adaptive) {
345+
const std::string refinement_type = pin->GetString("refinement", "type");
346+
is_cubic_refinement = (refinement_type == "cubic");
347+
}
348+
349+
PARTHENON_REQUIRE_THROWS(!is_adaptive || is_cubic_refinement,
350+
"Tracers/swarms currently only supported on non-adaptive "
351+
"meshes or with cubic adaptive refinement.");
342352

343353
if (mhd) {
344354
tracers_pkg->AddSwarmValue("B_x", swarm_name, real_swarmvalue_metadata);
@@ -347,6 +357,12 @@ std::shared_ptr<StateDescriptor> Initialize(ParameterInput *pin) {
347357
tracers_pkg->AddSwarmValue("rot_B_x", swarm_name, real_swarmvalue_metadata);
348358
tracers_pkg->AddSwarmValue("rot_B_y", swarm_name, real_swarmvalue_metadata);
349359
tracers_pkg->AddSwarmValue("rot_B_z", swarm_name, real_swarmvalue_metadata);
360+
tracers_pkg->AddSwarmValue("tens_B_x", swarm_name, real_swarmvalue_metadata);
361+
tracers_pkg->AddSwarmValue("tens_B_y", swarm_name, real_swarmvalue_metadata);
362+
tracers_pkg->AddSwarmValue("tens_B_z", swarm_name, real_swarmvalue_metadata);
363+
tracers_pkg->AddSwarmValue("grad_B2_x", swarm_name, real_swarmvalue_metadata);
364+
tracers_pkg->AddSwarmValue("grad_B2_y", swarm_name, real_swarmvalue_metadata);
365+
tracers_pkg->AddSwarmValue("grad_B2_z", swarm_name, real_swarmvalue_metadata);
350366
}
351367

352368
auto nscalars = pin->GetOrAddInteger("hydro", "nscalars", 0);
@@ -1010,13 +1026,25 @@ TaskStatus FillTracers(MeshData<Real> *md, parthenon::SimTime &tm) {
10101026
auto rot_B_x = vel_x.Get();
10111027
auto rot_B_y = vel_x.Get();
10121028
auto rot_B_z = vel_x.Get();
1029+
auto tens_B_x = vel_x.Get();
1030+
auto tens_B_y = vel_x.Get();
1031+
auto tens_B_z = vel_x.Get();
1032+
auto grad_B2_x = vel_x.Get();
1033+
auto grad_B2_y = vel_x.Get();
1034+
auto grad_B2_z = vel_x.Get();
10131035
if (mhd) {
10141036
B_x = swarm->Get<Real>("B_x").Get();
10151037
B_y = swarm->Get<Real>("B_y").Get();
10161038
B_z = swarm->Get<Real>("B_z").Get();
10171039
rot_B_x = swarm->Get<Real>("rot_B_x").Get();
10181040
rot_B_y = swarm->Get<Real>("rot_B_y").Get();
10191041
rot_B_z = swarm->Get<Real>("rot_B_z").Get();
1042+
tens_B_x = swarm->Get<Real>("tens_B_x").Get();
1043+
tens_B_y = swarm->Get<Real>("tens_B_y").Get();
1044+
tens_B_z = swarm->Get<Real>("tens_B_z").Get();
1045+
grad_B2_x = swarm->Get<Real>("grad_B2_x").Get();
1046+
grad_B2_y = swarm->Get<Real>("grad_B2_y").Get();
1047+
grad_B2_z = swarm->Get<Real>("grad_B2_z").Get();
10201048
}
10211049

10221050
auto &density = swarm->Get<Real>("density").Get();
@@ -1136,6 +1164,61 @@ TaskStatus FillTracers(MeshData<Real> *md, parthenon::SimTime &tm) {
11361164
rot_B_x(n) = rotBx;
11371165
rot_B_y(n) = rotBy;
11381166
rot_B_z(n) = rotBz;
1167+
1168+
// --- Magnetic pressure gradient: -grad(B^2 / 2) ---
1169+
const Real dB2_dx =
1170+
((prim_pack(b, IB1, k, j, i + 1) * prim_pack(b, IB1, k, j, i + 1) +
1171+
prim_pack(b, IB2, k, j, i + 1) * prim_pack(b, IB2, k, j, i + 1) +
1172+
((ndim == 3) ? prim_pack(b, IB3, k, j, i + 1) *
1173+
prim_pack(b, IB3, k, j, i + 1)
1174+
: 0.0)) -
1175+
(prim_pack(b, IB1, k, j, i - 1) * prim_pack(b, IB1, k, j, i - 1) +
1176+
prim_pack(b, IB2, k, j, i - 1) * prim_pack(b, IB2, k, j, i - 1) +
1177+
((ndim == 3) ? prim_pack(b, IB3, k, j, i - 1) *
1178+
prim_pack(b, IB3, k, j, i - 1)
1179+
: 0.0))) /
1180+
(2.0 * dx);
1181+
1182+
const Real dB2_dy =
1183+
((prim_pack(b, IB1, k, j + 1, i) * prim_pack(b, IB1, k, j + 1, i) +
1184+
prim_pack(b, IB2, k, j + 1, i) * prim_pack(b, IB2, k, j + 1, i) +
1185+
((ndim == 3) ? prim_pack(b, IB3, k, j + 1, i) *
1186+
prim_pack(b, IB3, k, j + 1, i)
1187+
: 0.0)) -
1188+
(prim_pack(b, IB1, k, j - 1, i) * prim_pack(b, IB1, k, j - 1, i) +
1189+
prim_pack(b, IB2, k, j - 1, i) * prim_pack(b, IB2, k, j - 1, i) +
1190+
((ndim == 3) ? prim_pack(b, IB3, k, j - 1, i) *
1191+
prim_pack(b, IB3, k, j - 1, i)
1192+
: 0.0))) /
1193+
(2.0 * dy);
1194+
1195+
const Real dB2_dz = (ndim == 3) ? ((prim_pack(b, IB1, k + 1, j, i) *
1196+
prim_pack(b, IB1, k + 1, j, i) +
1197+
prim_pack(b, IB2, k + 1, j, i) *
1198+
prim_pack(b, IB2, k + 1, j, i) +
1199+
prim_pack(b, IB3, k + 1, j, i) *
1200+
prim_pack(b, IB3, k + 1, j, i)) -
1201+
(prim_pack(b, IB1, k - 1, j, i) *
1202+
prim_pack(b, IB1, k - 1, j, i) +
1203+
prim_pack(b, IB2, k - 1, j, i) *
1204+
prim_pack(b, IB2, k - 1, j, i) +
1205+
prim_pack(b, IB3, k - 1, j, i) *
1206+
prim_pack(b, IB3, k - 1, j, i))) /
1207+
(2.0 * dz)
1208+
: 0.0;
1209+
1210+
grad_B2_x(n) = dB2_dx;
1211+
grad_B2_y(n) = dB2_dy;
1212+
grad_B2_z(n) = dB2_dz;
1213+
1214+
// --- Magnetic tension: (B · ∇)B ---
1215+
const Real tens_x = Bx * dBx_dx + By * dBx_dy + Bz * dBx_dz;
1216+
const Real tens_y = Bx * dBy_dx + By * dBy_dy + Bz * dBy_dz;
1217+
const Real tens_z = Bx * dBz_dx + By * dBz_dy + Bz * dBz_dz;
1218+
1219+
tens_B_x(n) = tens_x;
1220+
tens_B_y(n) = tens_y;
1221+
tens_B_z(n) = tens_z;
11391222
}
11401223

11411224
// Central differences using neighbors

0 commit comments

Comments
 (0)