Skip to content

Commit 8165354

Browse files
author
Ben Prather
committed
Annoyances before release
Fix Bondi problem w/Dirichlet bounds in a general way Update bondi.par to demo some more options Be more explicit about B in 1D: only flux_ct works Replace some flag checks with friendlier version
1 parent ca7a009 commit 8165354

16 files changed

Lines changed: 105 additions & 84 deletions

File tree

‎.gitignore‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -25,6 +25,7 @@ log_*.txt
2525
bondi_analytic_*.txt
2626
atmosphere_soln_*.txt
2727
shock_soln_*.txt
28+
parthinput.*
2829

2930
# Editor documents
3031
.project

‎kharma/b_cleanup/b_cleanup.cpp‎

Lines changed: 6 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -390,7 +390,7 @@ TaskStatus B_Cleanup::CornerLaplacian(MeshData<Real>* md, const std::string& p_v
390390
auto rc = md->GetBlockData(i);
391391
auto pmb = rc->GetBlockPointer();
392392
auto dB_block = rc->PackVariables(std::vector<std::string>{"dB"});
393-
if (pmb->boundary_flag[BoundaryFace::inner_x2] == BoundaryFlag::user) {
393+
if (KBoundaries::IsPhysicalBoundary(pmb, BoundaryFace::inner_x2)) {
394394
pmb->par_for("dB_boundary", kb_l.s, kb_l.e, ib_l.s, ib_l.e,
395395
KOKKOS_LAMBDA (const int &k, const int &i) {
396396
dB_block(V1, k, jb.s-1, i) = dB_block(V1, k, jb.s, i);
@@ -399,7 +399,7 @@ TaskStatus B_Cleanup::CornerLaplacian(MeshData<Real>* md, const std::string& p_v
399399
}
400400
);
401401
}
402-
if (pmb->boundary_flag[BoundaryFace::outer_x2] == BoundaryFlag::user) {
402+
if (KBoundaries::IsPhysicalBoundary(pmb, BoundaryFace::outer_x2)) {
403403
pmb->par_for("dB_boundary", kb_l.s, kb_l.e, ib_l.s, ib_l.e,
404404
KOKKOS_LAMBDA (const int &k, const int &i) {
405405
dB_block(V1, k, jb.e+1, i) = dB_block(V1, k, jb.e, i);
@@ -474,14 +474,14 @@ TaskStatus B_Cleanup::CenterLaplacian(MeshData<Real>* md, const std::string& p_v
474474
auto pmb = rc->GetBlockPointer();
475475
auto dB_block = rc->PackVariables(std::vector<std::string>{"dB"});
476476
const IndexRange3 bi2 = KDomain::GetRange(md, IndexDomain::interior, F2);
477-
if (pmb->boundary_flag[BoundaryFace::inner_x2] == BoundaryFlag::user) {
477+
if (KBoundaries::IsPhysicalBoundary(pmb, BoundaryFace::inner_x2)) {
478478
pmb->par_for("dB_boundary", b2.ks, b2.ke, bi2.js, bi2.js, b2.is, b2.ie,
479479
KOKKOS_LAMBDA (const int &k, const int &j, const int &i) {
480480
dB_block(F2, 0, k, j, i) = 0.;
481481
}
482482
);
483483
}
484-
if (pmb->boundary_flag[BoundaryFace::outer_x2] == BoundaryFlag::user) {
484+
if (KBoundaries::IsPhysicalBoundary(pmb, BoundaryFace::outer_x2)) {
485485
pmb->par_for("dB_boundary", b2.ks, b2.ke, bi2.je, bi2.je, b2.is, b2.ie,
486486
KOKKOS_LAMBDA (const int &k, const int &j, const int &i) {
487487
dB_block(F2, 0, k, j, i) = 0.;
@@ -512,7 +512,7 @@ TaskStatus B_Cleanup::CenterLaplacian(MeshData<Real>* md, const std::string& p_v
512512
auto pmb = rc->GetBlockPointer();
513513
auto lap_block = rc->PackVariables(std::vector<std::string>{lap_var});
514514
const IndexRange3 bic = KDomain::GetRange(md, IndexDomain::interior);
515-
if (pmb->boundary_flag[BoundaryFace::inner_x1] == BoundaryFlag::user) {
515+
if (KBoundaries::IsPhysicalBoundary(pmb, BoundaryFace::inner_x2)) {
516516
pmb->par_for("lap_boundary", bc.ks, bc.ke, bc.js, bc.je, bc.is, bic.is,
517517
KOKKOS_LAMBDA (const int &k, const int &j, const int &i) {
518518
lap_block(0, k, j, i) = 0.;
@@ -527,7 +527,7 @@ TaskStatus B_Cleanup::CenterLaplacian(MeshData<Real>* md, const std::string& p_v
527527
auto pmb = rc->GetBlockPointer();
528528
auto lap_block = rc->PackVariables(std::vector<std::string>{lap_var});
529529
const IndexRange3 bic = KDomain::GetRange(md, IndexDomain::interior);
530-
if (pmb->boundary_flag[BoundaryFace::outer_x1] == BoundaryFlag::user) {
530+
if (KBoundaries::IsPhysicalBoundary(pmb, BoundaryFace::outer_x2)) {
531531
pmb->par_for("lap_boundary", bc.ks, bc.ke, bc.js, bc.je, bic.ie, bc.ie,
532532
KOKKOS_LAMBDA (const int &k, const int &j, const int &i) {
533533
lap_block(0, k, j, i) = 0.;

‎kharma/b_ct/b_ct.cpp‎

Lines changed: 6 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -41,6 +41,7 @@
4141

4242
#include <parthenon/parthenon.hpp>
4343
#include <prolong_restrict/pr_ops.hpp>
44+
#include <stdexcept>
4445

4546
using namespace parthenon;
4647
using parthenon::refinement_ops::ProlongateSharedMinMod;
@@ -253,14 +254,14 @@ TaskStatus B_CT::DangerousPtoU(MeshData<Real> *md, IndexDomain domain, bool coar
253254
auto rc = md->GetBlockData(i);
254255
auto pmb = rc->GetBlockPointer();
255256
auto dB_block = rc->PackVariables(std::vector<std::string>{"dB"});
256-
if (pmb->boundary_flag[BoundaryFace::inner_x2] == BoundaryFlag::user) {
257+
if (KBoundaries::IsPhysicalBoundary(pmb, BoundaryFace::inner_x2)) {
257258
pmb->par_for("dB_boundary", bf2.ks, bf2.ke, bf2.is, bf2.ie,
258259
KOKKOS_LAMBDA (const int &k, const int &i) {
259260
dB_block(F2, 0, k, bf2.js, i) = 0.;
260261
}
261262
);
262263
}
263-
if (pmb->boundary_flag[BoundaryFace::outer_x2] == BoundaryFlag::user) {
264+
if (KBoundaries::IsPhysicalBoundary(pmb, BoundaryFace::outer_x2)) {
264265
pmb->par_for("dB_boundary", bf2.ks, bf2.ke, bf2.is, bf2.ie,
265266
KOKKOS_LAMBDA (const int &k, const int &i) {
266267
dB_block(F2, 0, k, bf2.je, i) = 0.;
@@ -286,6 +287,7 @@ TaskStatus B_CT::CalculateEMF(MeshData<Real> *md)
286287
{
287288
auto pmesh = md->GetMeshPointer();
288289
const int ndim = pmesh->ndim;
290+
if (ndim < 2) throw std::runtime_error("Face-centered constrained transport does not support 1D! Use `flux_ct` instead.")
289291

290292
// EMF temporary
291293
auto& emf_pack = md->PackVariables(std::vector<std::string>{"B_CT.emf"});
@@ -510,7 +512,7 @@ TaskStatus B_CT::DerefinePoles(MeshData<Real> *md)
510512
auto bdir = KBoundaries::BoundaryDirection(bface);
511513
auto domain = KBoundaries::BoundaryDomain(bface);
512514
auto binner = KBoundaries::BoundaryIsInner(bface);
513-
if (bdir == X2DIR && pmb->boundary_flag[bface] == BoundaryFlag::user) {
515+
if (bdir == X2DIR && KBoundaries::IsPhysicalBoundary(pmb, bface)) {
514516
// indices
515517
// TODO also get ranges in cells from the beginning rather than using j_p & calculating j_c
516518
IndexRange3 bCC = KDomain::GetRange(rc, IndexDomain::interior, CC);
@@ -667,6 +669,7 @@ double B_CT::MaxDivB(MeshData<Real> *md)
667669
{
668670
auto pmesh = md->GetMeshPointer();
669671
const int ndim = pmesh->ndim;
672+
if (ndim < 2) return 0.;
670673

671674
auto B_U = md->PackVariables(std::vector<std::string>{"cons.fB"});
672675

‎kharma/b_flux_ct/b_flux_ct.cpp‎

Lines changed: 10 additions & 10 deletions
Original file line numberDiff line numberDiff line change
@@ -420,7 +420,7 @@ void ZeroBoundaryFlux(MeshData<Real> *md, IndexDomain domain, bool coarse)
420420
auto& B_F = rc->PackVariablesAndFluxes(std::vector<std::string>{"cons.B"});
421421

422422
if (domain == IndexDomain::inner_x2 &&
423-
pmb->boundary_flag[BoundaryFace::inner_x2] == BoundaryFlag::user) {
423+
KBoundaries::IsPhysicalBoundary(pmb, BoundaryFace::inner_x2)) {
424424
pmb->par_for("fix_flux_b_l", kbs.s, kbs.e, jbf.s, jbf.s, ibs.s, ibs.e,
425425
KOKKOS_LAMBDA (const int &k, const int &j, const int &i) {
426426
B_F.flux(X2DIR, V1, k, j, i) = 0.;
@@ -432,7 +432,7 @@ void ZeroBoundaryFlux(MeshData<Real> *md, IndexDomain domain, bool coarse)
432432
}
433433

434434
if (domain == IndexDomain::outer_x2 &&
435-
pmb->boundary_flag[BoundaryFace::outer_x2] == BoundaryFlag::user) {
435+
KBoundaries::IsPhysicalBoundary(pmb, BoundaryFace::outer_x2)) {
436436
pmb->par_for("fix_flux_b_r", kbs.s, kbs.e, jbf.e, jbf.e, ibs.s, ibs.e,
437437
KOKKOS_LAMBDA (const int &k, const int &j, const int &i) {
438438
B_F.flux(X2DIR, V1, k, j, i) = 0.;
@@ -450,7 +450,7 @@ void ZeroBoundaryFlux(MeshData<Real> *md, IndexDomain domain, bool coarse)
450450
// 2. However, B2 and B3 are normal outflow conditions -- despite the fluxes here, the outflow
451451
// conditions will set them equal to the last zone.
452452
if (domain == IndexDomain::inner_x1 &&
453-
pmb->boundary_flag[BoundaryFace::inner_x1] == BoundaryFlag::user) {
453+
KBoundaries::IsPhysicalBoundary(pmb, BoundaryFace::inner_x1)) {
454454
pmb->par_for("fix_flux_b_in_old", kbs.s, kbs.e, jbs.s, jbs.e, ibf.s, ibf.s,
455455
KOKKOS_LAMBDA (const int &k, const int &j, const int &i) {
456456
B_F.flux(X1DIR, V2, k, j, i) = 0.;
@@ -462,7 +462,7 @@ void ZeroBoundaryFlux(MeshData<Real> *md, IndexDomain domain, bool coarse)
462462
}
463463

464464
if (domain == IndexDomain::outer_x1 &&
465-
pmb->boundary_flag[BoundaryFace::outer_x1] == BoundaryFlag::user) {
465+
KBoundaries::IsPhysicalBoundary(pmb, BoundaryFace::outer_x1)) {
466466
pmb->par_for("fix_flux_b_out_old", kbs.s, kbs.e, jbs.s, jbs.e, ibf.e, ibf.e,
467467
KOKKOS_LAMBDA (const int &k, const int &j, const int &i) {
468468
B_F.flux(X1DIR, V2, k, j, i) = 0.;
@@ -509,7 +509,7 @@ void Bflux0(MeshData<Real> *md, IndexDomain domain, bool coarse)
509509
// Allows nonzero flux across X1 boundary but still keeps divB=0 (turns out effectively to have 0 flux)
510510
// Usable only for Dirichlet conditions
511511
if (domain == IndexDomain::inner_x1 &&
512-
pmb->boundary_flag[BoundaryFace::inner_x1] == BoundaryFlag::user)
512+
KBoundaries::IsPhysicalBoundary(pmb, BoundaryFace::inner_x1))
513513
{
514514
pmb->par_for("fix_flux_b_in", kbs.s, kbs.e, jbs.s, jbs.e, ibf.s, ibf.s, // Hyerin (12/28/22) for 1st & 2nd prescription
515515
KOKKOS_LAMBDA (const int &k, const int &j, const int &i) {
@@ -521,7 +521,7 @@ void Bflux0(MeshData<Real> *md, IndexDomain domain, bool coarse)
521521

522522
}
523523
if (domain == IndexDomain::inner_x2 &&
524-
pmb->boundary_flag[BoundaryFace::inner_x2] == BoundaryFlag::user)
524+
KBoundaries::IsPhysicalBoundary(pmb, BoundaryFace::inner_x2))
525525
{
526526
pmb->par_for("fix_flux_b_in", kbs.s, kbs.e, jbf.s, jbf.s, ibs.s, ibs.e,
527527
KOKKOS_LAMBDA (const int &k, const int &j, const int &i) {
@@ -534,7 +534,7 @@ void Bflux0(MeshData<Real> *md, IndexDomain domain, bool coarse)
534534

535535
// OUTER
536536
if (domain == IndexDomain::outer_x1 &&
537-
pmb->boundary_flag[BoundaryFace::outer_x1] == BoundaryFlag::user)
537+
KBoundaries::IsPhysicalBoundary(pmb, BoundaryFace::outer_x1))
538538
{
539539
pmb->par_for("fix_flux_b_out", kbs.s, kbs.e, jbs.s, jbs.e, ibf.e, ibf.e, // Hyerin (12/28/22) for 1st & 2nd prescription
540540
KOKKOS_LAMBDA (const int &k, const int &j, const int &i) {
@@ -545,7 +545,7 @@ void Bflux0(MeshData<Real> *md, IndexDomain domain, bool coarse)
545545
);
546546
}
547547
if (domain == IndexDomain::outer_x2 &&
548-
pmb->boundary_flag[BoundaryFace::outer_x2] == BoundaryFlag::user)
548+
KBoundaries::IsPhysicalBoundary(pmb, BoundaryFace::outer_x2))
549549
{
550550
pmb->par_for("fix_flux_b_out", kbs.s, kbs.e, jbf.e, jbf.e, ibs.s, ibs.e,
551551
KOKKOS_LAMBDA (const int &k, const int &j, const int &i) {
@@ -563,9 +563,9 @@ IndexRange ValidDivBX1(MeshBlock *pmb)
563563
// intersecting the interior & exterior faces. Don't report these zones, as we expect it.
564564
const IndexRange ibl = pmb->meshblock_data.Get("base")->GetBoundsI(IndexDomain::interior);
565565
bool avoid_inner = (!pmb->packages.Get("B_FluxCT")->Param<bool>("fix_flux_inner_x1") &&
566-
pmb->boundary_flag[BoundaryFace::inner_x1] == BoundaryFlag::user);
566+
KBoundaries::IsPhysicalBoundary(pmb, BoundaryFace::inner_x1));
567567
bool avoid_outer = (!pmb->packages.Get("B_FluxCT")->Param<bool>("fix_flux_outer_x1") &&
568-
pmb->boundary_flag[BoundaryFace::outer_x1] == BoundaryFlag::user);
568+
KBoundaries::IsPhysicalBoundary(pmb, BoundaryFace::outer_x1));
569569
return IndexRange{ibl.s + (avoid_inner), ibl.e + (!avoid_outer)};
570570
}
571571

‎kharma/boundaries/boundaries.cpp‎

Lines changed: 8 additions & 7 deletions
Original file line numberDiff line numberDiff line change
@@ -34,6 +34,7 @@
3434
#include "boundaries.hpp"
3535

3636
#include "bondi.hpp"
37+
#include "boundary_types.hpp"
3738
#include "decs.hpp"
3839
#include "domain.hpp"
3940
#include "kharma.hpp"
@@ -108,7 +109,7 @@ std::shared_ptr<KHARMAPackage> KBoundaries::Initialize(ParameterInput *pin, std:
108109
// Face-centered fields: some duplicate stuff, leaving it separate for now
109110
FC ghost_vars_f = FC({Metadata::FillGhost, Metadata::Face})
110111
- FC({Metadata::GetUserFlag("StartupOnly")});
111-
int nvar_f = 3 * m::max(KHARMA::PackDimension(packages.get(), ghost_vars_f), 1);
112+
int nvar_f = 3 * KHARMA::PackDimension(packages.get(), ghost_vars_f);
112113

113114
// TODO encapsulate this
114115
Metadata m_x1, m_x2, m_x3, m_x1_f, m_x2_f, m_x3_f;
@@ -554,14 +555,14 @@ void KBoundaries::ApplyBoundary(std::shared_ptr<MeshBlockData<Real>> &rc, IndexD
554555
// exactly the same function that Parthenon does. This ensures we're applying
555556
// the same thing, just emulating calling it after X2.
556557
if (params.Get<bool>("fix_corner_inner")) {
557-
if (pmb->boundary_flag[BoundaryFace::inner_x1] == BoundaryFlag::user) {
558+
if (KBoundaries::IsPhysicalBoundary(pmb, BoundaryFace::inner_x1)) {
558559
Flag("FixCorner");
559560
ApplyBoundary(rc, IndexDomain::inner_x1, coarse);
560561
EndFlag();
561562
}
562563
}
563564
if (params.Get<bool>("fix_corner_outer")) {
564-
if (pmb->boundary_flag[BoundaryFace::outer_x1] == BoundaryFlag::user) {
565+
if (KBoundaries::IsPhysicalBoundary(pmb, BoundaryFace::outer_x1)) {
565566
Flag("FixCorner");
566567
ApplyBoundary(rc, IndexDomain::outer_x1, coarse);
567568
EndFlag();
@@ -667,7 +668,7 @@ TaskStatus KBoundaries::FixFlux(MeshData<Real> *md)
667668
if (params.Get<bool>("check_inflow_" + bname)) {
668669
const int m_rho = cons_map["cons.rho"].first;
669670
// ...and if this face of the block corresponds to a global boundary...
670-
if (pmb->boundary_flag[bface] == BoundaryFlag::user) {
671+
if (KBoundaries::IsPhysicalBoundary(pmb, bface)) {
671672
if (binner) {
672673
pmb->par_for(
673674
"zero_inflow_flux_" + bname, b.ks, b.ke, b.js, b.je, b.is, b.ie,
@@ -689,7 +690,7 @@ TaskStatus KBoundaries::FixFlux(MeshData<Real> *md)
689690
// If we should zero flux through this face...
690691
if (params.Get<bool>("zero_flux_" + bname)) {
691692
// ...and if this face of the block corresponds to a global boundary...
692-
if (pmb->boundary_flag[bface] == BoundaryFlag::user) {
693+
if (KBoundaries::IsPhysicalBoundary(pmb, bface)) {
693694
pmb->par_for(
694695
"zero_flux_" + bname, 0, F.GetDim(4) - 1, b.ks, b.ke, b.js, b.je, b.is, b.ie,
695696
KOKKOS_LAMBDA(const int &p, const int &k, const int &j, const int &i) {
@@ -702,7 +703,7 @@ TaskStatus KBoundaries::FixFlux(MeshData<Real> *md)
702703
// If we should replace fluxes with excised versions...
703704
if (params.Get<bool>("excise_flux_" + bname)) {
704705
// ...and if this face of the block corresponds to a global boundary...
705-
if (pmb->boundary_flag[bface] == BoundaryFlag::user) {
706+
if (IsPhysicalBoundary(pmb, bface)) {
706707
if (bdir != 2) throw std::runtime_error("Excised polar fluxes only fully implemented in X2!");
707708

708709
// Pack w/B to match indices with the `Flux.X` below
@@ -891,7 +892,7 @@ void KBoundaries::AddSource(MeshData<Real> *md, MeshData<Real> *mdudt, IndexDoma
891892
// If we should replace fluxes with excised versions...
892893
if (params.Get<bool>("excise_flux_" + bname)) {
893894
// ...and if this face of the block corresponds to a global boundary...
894-
if (pmb->boundary_flag[bface] == BoundaryFlag::user) {
895+
if (IsPhysicalBoundary(pmb, bface)) {
895896
if (bdir != 2) throw std::runtime_error("Excised polar fluxes only fully implemented in X2!");
896897

897898
const IndexRange3 bi = KDomain::GetRange(rc, IndexDomain::interior);

‎kharma/boundaries/dirichlet.cpp‎

Lines changed: 10 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -34,6 +34,7 @@
3434

3535
#include "dirichlet.hpp"
3636

37+
#include "boundary_types.hpp"
3738
#include "domain.hpp"
3839
#include "types.hpp"
3940

@@ -117,9 +118,12 @@ void KBoundaries::FreezeDirichlet(std::shared_ptr<MeshData<Real>> &md)
117118
for (int i=0; i < md->NumBlocks(); i++) {
118119
auto rc = md->GetBlockData(i).get();
119120
auto pmb = rc->GetBlockPointer();
120-
auto domain = BoundaryDomain(bface);
121-
// Set whatever is in that domain as the Dirichlet bound
122-
SetDomainDirichlet(rc, domain, false);
121+
// ...if the face is not internal
122+
if (KBoundaries::IsPhysicalBoundary(pmb, bface)) {
123+
auto domain = BoundaryDomain(bface);
124+
// Set whatever is in that domain as the Dirichlet bound
125+
SetDomainDirichlet(rc, domain, false);
126+
}
123127
}
124128
}
125129
}
@@ -131,8 +135,9 @@ void KBoundaries::FreezeDirichletBlock(MeshBlockData<Real> *rc)
131135
BoundaryFace bface = (BoundaryFace) i;
132136
auto bname = BoundaryName(bface);
133137
auto pmb = rc->GetBlockPointer();
134-
// ...if this boundary is dirichlet...
135-
if (pmb->packages.Get("Boundaries")->Param<std::string>(bname) == "dirichlet") {
138+
// ...if this boundary is dirichlet and physical...
139+
if (pmb->packages.Get("Boundaries")->Param<std::string>(bname) == "dirichlet" &&
140+
KBoundaries::IsPhysicalBoundary(pmb, bface)) {
136141
//std::cout << "Freezing dirichlet " << bname << " on block." << std::endl;
137142
auto domain = BoundaryDomain(bface);
138143
// Set whatever is in that domain as the Dirichlet bound

‎kharma/driver/kharma_driver.cpp‎

Lines changed: 2 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -275,9 +275,7 @@ TaskID KHARMADriver::AddFluxCalculations(TaskID& t_start, TaskList& tl, MeshData
275275
t_calculate_flux3 = tl.AddTask(t_start_fluxes, Flux::GetFlux<RType::mp5, X3DIR>, md);
276276
break;
277277
default:
278-
std::cerr << "Reconstruction type not supported! Main supported reconstructions:" << std::endl
279-
<< "donor_cell, linear_mc, weno5" << std::endl;
280-
throw std::invalid_argument("Unsupported reconstruction algorithm!");
278+
throw std::invalid_argument("Unsupported reconstruction algorithm! Main supported algorithms: linear_mc, weno5, weno5_linear");
281279
}
282280
auto t_calc_fluxes = t_calculate_flux1 | t_calculate_flux2 | t_calculate_flux3;
283281

@@ -379,8 +377,7 @@ TaskID KHARMADriver::AddStateUpdate(TaskID& t_start, TaskList& tl, MeshData<Real
379377
auto& pkgs = pmb0->packages.AllPackages();
380378

381379
// If we're explicitly evolving, UtoP needs a guess (except Kastaun inverter)
382-
if (!pkgs.at("GRMHD")->Param<bool>("implicit") &&
383-
pkgs.at("Inverter")->Param<Inverter::Type>("inverter_type") != Inverter::Type::kastaun) {
380+
if (!pkgs.at("GRMHD")->Param<bool>("implicit")) {
384381
t_copy_prims = tl.AddTask(t_start, Copy<MeshData<Real>>,
385382
std::vector<MetadataFlag>({Metadata::GetUserFlag("HD"), Metadata::GetUserFlag("Primitive")}),
386383
md_sub_step_init, md_update);

‎kharma/flux/fofc.cpp‎

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -84,8 +84,8 @@ TaskStatus Flux::MarkFOFC(MeshData<Real> *guess)
8484
for (int i_block = 0; i_block < guess->NumBlocks(); i_block++) {
8585
auto &rc = guess->GetBlockData(i_block);
8686
auto pmb = rc->GetBlockPointer();
87-
const bool is_inner_x2 = pmb->boundary_flag[BoundaryFace::inner_x2] == BoundaryFlag::user;
88-
const bool is_outer_x2 = pmb->boundary_flag[BoundaryFace::outer_x2] == BoundaryFlag::user;
87+
const bool is_inner_x2 = KBoundaries::IsPhysicalBoundary(pmb, BoundaryFace::inner_x2);
88+
const bool is_outer_x2 = KBoundaries::IsPhysicalBoundary(pmb, BoundaryFace::outer_x2);
8989
if (is_inner_x2 || is_outer_x2) {
9090
auto lfofcflag = rc->PackVariables(std::vector<std::string>{"fofcflag"});
9191
if (is_inner_x2) {

0 commit comments

Comments
 (0)