From eb939fc0767dd5def934eaea401b865bdb64f74c Mon Sep 17 00:00:00 2001 From: Ann Almgren Date: Wed, 7 Oct 2026 17:37:27 -0700 Subject: [PATCH] Fix issues #221, #223, #225, #228, #231, #232, #234, #235, #237, #238, #242, #247, #248 #248 ViscBench: fill plotfile components by name rather than by position, so a plotfile written with a different amr.plot_vars ordering still gives the right error norms. #247 ConvertCheckpointGrids: extrapolate into the ghost cells that the periodic FillBoundary does not reach (coarse-fine and non-periodic physical boundaries), which otherwise still held the setVal(10.) sentinel when the interpolation slopes were computed. Also build the target MultiFabs on the source DistributionMapping, since the loops index source and target through a single MFIter. #234 prob_initData: read prob.velocity_factor only for probtype 11, and reset IC.v_x first so prob.velocity_ic can neither be overwritten by it nor stand in for it. #235 Projection::set_outflow_bcs: abort unconditionally on a level that only partially covers an outflow face, instead of only asserting in a debug build. #232 Initialize_bcs: range-check the integer ns.lo_bc/ns.hi_bc path; the values index the six-entry tables in NS_BC.H. #231 Drop the stale derive name diveru from the tutorial inputs files. #228 NavierStokes::Initialize: read ns.set_plot_coveredCell_val, which was documented but never queried. #223 post_timestep_particle: validate particles.timestamp_indices against NUM_STATE at first use. #221 Remove the unused umacFill functor and its stale comment. #237 mlmg_mac_solve: register the coarse/fine ratio on every level above 0, not only when a coarse phi is supplied; the sync solve passes null and MLMG was falling back to a ratio of 2. #238 Remove MacProj::set_outflow_bcs, which only rewrote zeros that mac_project had just set but could still abort with "Cannot yet handle partially refined outflow", together with grids_on_side_of_domain, the then-unused have_divu argument and the mac_proj.do_outflow_bcs key. #242 calc_mut_LES: zero the velocity-gradient FluxBoxes before compVelGrad, which skips covered FABs. #225 set_state_in_checkpoint: correct the abort text -- the flags must mark exactly the missing state types. Co-Authored-By: Claude Opus 5 (1M context) --- Exec/eb_run2d/regtest.2d.bubble | 2 +- Source/MacProj.H | 10 -- Source/MacProj.cpp | 144 +----------------- Source/NS_LES.cpp | 10 ++ Source/NavierStokes.cpp | 20 +++ Source/NavierStokesBase.cpp | 50 +++--- Source/Projection.cpp | 17 ++- Source/prob/prob_init.cpp | 14 +- Tutorials/Bubble/inputs.2d.bubble | 2 +- Tutorials/HotSpot/inputs.2d.average_hotspot | 2 +- Tutorials/HotSpot/inputs.3d.LES_hotspot | 2 +- .../RayleighTaylor/inputs.2d.rayleightaylor | 2 +- .../TaylorGreen/benchmarks/ViscBench.cpp | 11 +- .../ConvertCheckpointGrids.cpp | 35 ++++- 14 files changed, 124 insertions(+), 197 deletions(-) diff --git a/Exec/eb_run2d/regtest.2d.bubble b/Exec/eb_run2d/regtest.2d.bubble index 81045bc65..b6a49593e 100644 --- a/Exec/eb_run2d/regtest.2d.bubble +++ b/Exec/eb_run2d/regtest.2d.bubble @@ -133,4 +133,4 @@ prob.blob_radius = 0.2 prob.density_ic = 2.0 #mag_vort is stepping out of bounds -#amr.derive_plot_vars = mag_vort diveru avg_pressure +#amr.derive_plot_vars = mag_vort avg_pressure diff --git a/Source/MacProj.H b/Source/MacProj.H index 74b89df98..d6dfb15d9 100644 --- a/Source/MacProj.H +++ b/Source/MacProj.H @@ -39,7 +39,6 @@ public: amrex::Real dt, amrex::Real prev_time, const amrex::MultiFab& divu, - int have_divu, const amrex::BCRec& density_math_bc, bool increment_vel_register = true ); @@ -104,14 +103,6 @@ public: void check_div_cond (int level, amrex::MultiFab U_edge[]) const; // - // Boundary conditions. - // - void set_outflow_bcs (int level, - amrex::MultiFab* mac_phi, - const amrex::MultiFab* u_mac, - const amrex::MultiFab& S, - const amrex::MultiFab& divu); - // // Pointers to amr,amrlevel. // amrex::Amr* parent; @@ -132,7 +123,6 @@ public: int finest_level_allocated; static int verbose; - static int do_outflow_bcs; static amrex::Real mac_tol; static amrex::Real mac_abs_tol; static amrex::Real mac_sync_tol; diff --git a/Source/MacProj.cpp b/Source/MacProj.cpp index c69bb01a1..e7134108e 100644 --- a/Source/MacProj.cpp +++ b/Source/MacProj.cpp @@ -3,7 +3,6 @@ #include #include #include -#include #include #ifdef AMREX_USE_EB @@ -25,7 +24,6 @@ int MacProj::verbose; Real MacProj::mac_tol; Real MacProj::mac_abs_tol; Real MacProj::mac_sync_tol; -int MacProj::do_outflow_bcs; int MacProj::check_umac_periodicity; int MacProj::max_order = 4; int MacProj::agglomeration = 1; @@ -49,7 +47,6 @@ MacProj::Initialize () MacProj::mac_tol = 1.0e-12; MacProj::mac_abs_tol = 1.0e-16; MacProj::mac_sync_tol = 1.0e-10; - MacProj::do_outflow_bcs = 1; // // Only check umac periodicity when debugging. Can be overridden on input. // @@ -65,7 +62,6 @@ MacProj::Initialize () pp.query("mac_tol", mac_tol); pp.query("mac_abs_tol", mac_abs_tol); pp.query("mac_sync_tol", mac_sync_tol); - pp.query("do_outflow_bcs", do_outflow_bcs); pp.query("check_umac_periodicity", check_umac_periodicity); pp.query("umac_periodic_test_Tol", umac_periodic_test_Tol); @@ -184,38 +180,6 @@ MacProj::cleanup (int level) // // Projection functions follow ... // -static -bool -grids_on_side_of_domain (const BoxArray& grids, - const Box& domain, - const Orientation& outFace) -{ - const int idir = outFace.coordDir(); - - if (outFace.isLow()) - { - for (int igrid = 0; igrid < grids.size(); igrid++) - { - if (grids[igrid].smallEnd(idir) == domain.smallEnd(idir)) - { - return true; - } - } - } - - if (outFace.isHigh()) - { - for (int igrid = 0; igrid < grids.size(); igrid++) - { - if (grids[igrid].bigEnd(idir) == domain.bigEnd(idir)) - { - return true; - } - } - } - - return false; -} // // Compute the level advance mac projection. @@ -228,7 +192,6 @@ MacProj::mac_project (int level, Real dt, Real time, const MultiFab& divu, - int have_divu, const BCRec& density_math_bc, bool increment_vel_register ) { @@ -262,10 +225,6 @@ MacProj::mac_project (int level, const MultiFab& rhotime = ns.get_rho(time); MultiFab::Copy(S, rhotime, 0, Density, 1, 1); - if (OutFlowBC::HasOutFlowBC(phys_bc) && have_divu && do_outflow_bcs) { - set_outflow_bcs(level, mac_phi, u_mac, S, divu); - } - // // Set up the mac projection // @@ -861,100 +820,6 @@ MacProj::check_div_cond (int level, } } -void -MacProj::set_outflow_bcs (int level, - MultiFab* mac_phi, - const MultiFab* /*u_mac*/, - const MultiFab& /*S*/, - const MultiFab& /*divu*/) -{ - // - // This code is very similar to the outflow BC stuff in the Projection - // class except that here the the phi to be solved for lives on the - // out-directed faces. The projection equation to satisfy is - // - // (1/r)(d/dr)[r/rho dphi/dr] = dv/dr - S - // - bool hasOutFlow; - Orientation outFaces[2*AMREX_SPACEDIM]; - int numOutFlowFaces; - - OutFlowBC::GetOutFlowFaces(hasOutFlow,outFaces,phys_bc,numOutFlowFaces); - - const BoxArray& grids = LevelData[level]->boxArray(); - const Geometry& geom = parent->Geom(level); - const Box& domain = parent->Geom(level).Domain(); - // - // Create 1-wide cc box just outside boundary to hold phi. - // - BoxList ccBoxList, phiBoxList; - // numOutFlowFaces gives the number of outflow faces on the entire - // problem domain - // nOutFlowTouched gives the number of outflow faces a level touches, so - // nOutFlowTouched = numOutFlowFaces for level 0, but - // nOutFlowTouched <= numOutFlowFaces for levels > 0, since - // finer levels may not span the entire problem domain - int nOutFlowTouched = 0; - for (int iface = 0; iface < numOutFlowFaces; iface++) - { - if (grids_on_side_of_domain(grids,geom.Domain(),outFaces[iface])) - { - nOutFlowTouched++; - const int outDir = outFaces[iface].coordDir(); - - Box ccBndBox; - if (outFaces[iface].faceDir() == Orientation::high) - { - ccBndBox = amrex::adjCellHi(domain,outDir,2); - ccBndBox.shift(outDir,-2); - } - else - { - ccBndBox = amrex::adjCellLo(domain,outDir,2); - ccBndBox.shift(outDir,2); - } - ccBoxList.push_back(ccBndBox); - - Box phiBox = amrex::adjCell(domain,outFaces[iface],1); - phiBoxList.push_back(phiBox); - - const Box& valid_ccBndBox = ccBndBox & domain; - const BoxArray uncovered_outflow_ba = amrex::complementIn(valid_ccBndBox,grids); - - if ((!uncovered_outflow_ba.empty()) && - grids.intersects(valid_ccBndBox)) - { - amrex::Error("MacProj: Cannot yet handle partially refined outflow"); - } - } - } - - if ( !ccBoxList.isEmpty() ) - { - BoxArray phiBoxArray(phiBoxList); - phiBoxList.clear(); - - // - // Must do this kind of copy instead of mac_phi->copy(phidat); - // because we're copying onto the ghost cells of the FABs, - // not the valid regions. - // -#ifdef _OPENMP -#pragma omp parallel if (Gpu::notInLaunchRegion()) -#endif - for ( int iface = 0; iface < nOutFlowTouched; ++iface ) - { - for (MFIter mfi(*mac_phi); mfi.isValid(); ++mfi) - { - Box ovlp = (*mac_phi)[mfi].box() & phiBoxArray[iface]; - if (ovlp.ok()) { - (*mac_phi)[mfi].setVal(0,ovlp,0,1); - } - } - } - } -} - // // Structure used by test_umac_periodic(). // @@ -1188,7 +1053,14 @@ MacProj::mlmg_mac_solve (Amr* a_parent, const MultiFab* cphi, const BCRec& a_phy set_mac_solve_bc(mlmg_lobc, mlmg_hibc, a_phys_bc, geom); macproj.setDomainBC(mlmg_lobc, mlmg_hibc); - if (level > 0 && cphi) + // + // cphi is null for the sync solve, which wants a homogeneous Dirichlet + // coarse/fine BC. MLMG still needs to be told the true refinement + // ratio, however, or it assumes 2 and mis-positions the coarse ghost + // value. So register the ratio on every level above 0, with or without + // coarse data. + // + if (level > 0) { macproj.setCoarseFineBC(cphi, a_parent->refRatio(level-1)[0]); } diff --git a/Source/NS_LES.cpp b/Source/NS_LES.cpp index 3f647492c..671ad5ab4 100644 --- a/Source/NS_LES.cpp +++ b/Source/NS_LES.cpp @@ -94,6 +94,16 @@ NavierStokesBase::calc_mut_LES(MultiFab* mu_LES[AMREX_SPACEDIM], const Real time MultiFab** tensorflux = fb.get(); std::array grad_Uvel{AMREX_D_DECL(tensorflux[0], tensorflux[1], tensorflux[2])}; + // + // FluxBoxes does not initialize the face MultiFabs, and compVelGrad skips + // covered FABs (covered tiles when tiling). The loop below evaluates the + // LES formula everywhere, so zero the gradients first to get mu_t = 0 on + // covered faces rather than whatever the arena held. + // + for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) { + grad_Uvel[idim]->setVal(0.0); + } + tensorop.compVelGrad(0,{grad_Uvel},{Uvel},MLLinOp::Location::FaceCenter); // diff --git a/Source/NavierStokes.cpp b/Source/NavierStokes.cpp index 643074fef..218028d61 100644 --- a/Source/NavierStokes.cpp +++ b/Source/NavierStokes.cpp @@ -38,6 +38,14 @@ NavierStokes::Initialize () NavierStokesBase::Initialize(); + { + // + // Documented in RunningProblems.rst; only used with EB. + // + ParmParse pp("ns"); + pp.query("set_plot_coveredCell_val", set_plot_coveredCell_val); + } + // // Set number of state variables. // @@ -95,6 +103,18 @@ NavierStokes::Initialize_bcs () pp.getarr("hi_bc",hi_bc,0,AMREX_SPACEDIM); for (int i = 0; i < AMREX_SPACEDIM; i++) { + // + // phys_bc indexes the six-entry tables in NS_BC.H, so anything + // outside [interior,noslipwall] would read past them. In particular + // PhysBCType also defines inflowoutflow=6, which IAMR does not + // support, and the integer input path is otherwise unchecked. + // + if ( lo_bc[i] < PhysBCType::interior || lo_bc[i] > PhysBCType::noslipwall || + hi_bc[i] < PhysBCType::interior || hi_bc[i] > PhysBCType::noslipwall ) + { + amrex::Abort("NavierStokes::Initialize_bcs: ns.lo_bc/ns.hi_bc must be in [0,5]; see the list of physical BC types in the documentation."); + } + phys_bc.setLo(i,lo_bc[i]); phys_bc.setHi(i,hi_bc[i]); } diff --git a/Source/NavierStokesBase.cpp b/Source/NavierStokesBase.cpp index 6498f5425..3107f2c5a 100644 --- a/Source/NavierStokesBase.cpp +++ b/Source/NavierStokesBase.cpp @@ -55,37 +55,6 @@ struct HomExtDirFill } }; -// -// A dummy function because FillPatch requires something to exist for filling dirichlet boundary conditions, -// even if we know we cannot have an ext_dir BC. -// u_mac BCs are only either periodic (BCType::int_dir) or first order extrapolation (FOEXTRAP). -// -struct umacFill -{ - AMREX_GPU_DEVICE - void operator()( - const amrex::IntVect& /*iv*/, - amrex::Array4 const& /*dummy*/, - const int /*dcomp*/, - const int numcomp, - amrex::GeometryData const& /*geom*/, - const amrex::Real /*time*/, - const amrex::BCRec* bcr, - const int bcomp, - const int /*orig_comp*/) const - { - // Abort if this function is expected to fill an ext_dir BC. - for (int n = bcomp; n < bcomp+numcomp; ++n) { - const amrex::BCRec& bc = bcr[n]; - if ( AMREX_D_TERM( bc.lo(0) == amrex::BCType::ext_dir || bc.hi(0) == amrex::BCType::ext_dir, - || bc.lo(1) == amrex::BCType::ext_dir || bc.hi(1) == amrex::BCType::ext_dir, - || bc.lo(2) == amrex::BCType::ext_dir || bc.hi(2) == amrex::BCType::ext_dir ) ) { - amrex::Abort("NavierStokesBase::umacFill: umac should not have BCType::ext_dir"); - } - } - } -}; - BCRec NavierStokesBase::phys_bc; Projection* NavierStokesBase::projector = nullptr; @@ -2146,7 +2115,7 @@ NavierStokesBase::mac_project (Real time, Vector density_math_bc = fetchBCArray(State_Type,Density,1); - mac_projector->mac_project(level,u_mac,S_old,dt,time,*divu,have_divu, + mac_projector->mac_project(level,u_mac,S_old,dt,time,*divu, density_math_bc[0], increment_vel_register); create_umac_grown(ngrow, divu); @@ -2758,7 +2727,7 @@ NavierStokesBase::set_state_in_checkpoint (Vector& state_in_checkpoint) // Abort if any of the NSB::*_in_checkpoint variables haven't been set by user. // if ( gradp_in_checkpoint<0 || average_in_checkpoint<0 ) - Abort("\n\n Checkpoint file is missing one or more state types. Set both\n ns.gradp_in_checkpoint and ns.avg_in_checkpoint to identify missing\n data. Set to 1 if present in checkpoint, 0 if not present. If unsure,\n try setting both to 0.\n\n If you just activated Time Averaging, you should add \n ns.avg_in_checkpoint=0 ns.gradp_in_checkpoint=1 \n\n"); + Abort("\n\n Checkpoint file is missing one or more state types. Set both\n ns.gradp_in_checkpoint and ns.avg_in_checkpoint to identify missing\n data. Set to 1 if present in checkpoint, 0 if not present. These must\n flag exactly the missing types: the checkpoint is read as a sequential\n stream, so marking a type that is actually present shifts every later\n read and silently restarts from mis-assigned data.\n\n If you just activated Time Averaging, you should add \n ns.avg_in_checkpoint=0 ns.gradp_in_checkpoint=1 \n\n"); // // Tell AmrLevel which types are in the checkpoint, so it knows what to copy. @@ -3987,6 +3956,21 @@ NavierStokesBase::post_timestep_particle (int crse_iteration) n = timestamp_indices.size(); nextras = timestamp_num_extras(); + // + // These index the FillPatched state below. NUM_STATE is not + // yet known when read_particle_params() reads the list, so + // validate here: BaseFab::copy only range-checks the + // component with an AMREX_ASSERT, so a bad entry would read + // outside the FAB in a release build. + // + for (int i = 0; i < n; ++i) + { + if (timestamp_indices[i] < 0 || timestamp_indices[i] >= NUM_STATE) + { + amrex::Abort("NavierStokesBase::post_timestep_particle: particles.timestamp_indices entries must be in [0,NUM_STATE)"); + } + } + int sz = n + nextras; tindices.reserve(sz); diff --git a/Source/Projection.cpp b/Source/Projection.cpp index bb5e43906..88e97e7ca 100644 --- a/Source/Projection.cpp +++ b/Source/Projection.cpp @@ -11,6 +11,8 @@ #include +#include + using namespace amrex; @@ -1820,8 +1822,19 @@ Projection::set_outflow_bcs (int which_call, const Box& valid_state_strip = temp_state_strip & domain; const BoxArray uncovered_outflow_ba = amrex::complementIn(valid_state_strip,Lgrids); - AMREX_ASSERT(uncovered_outflow_ba.empty() || - ! Lgrids.intersects(valid_state_strip)); + // + // A level that touches an outflow face without covering it is not + // supported: it would get no hydrostatic strip of its own, and in + // initialPressureProject its outflow nodes would keep phi = 0 while + // the coarse level gets rho*g*(H-z). Fail loudly rather than only + // in a debug build. + // + if ( !uncovered_outflow_ba.empty() && Lgrids.intersects(valid_state_strip) ) + { + amrex::Abort("Projection::set_outflow_bcs: level " + + std::to_string(lev) + + " only partially covers an outflow face"); + } if ( uncovered_outflow_ba.empty() && fine_level[iface] == -1) { int ii = icount[lev]; diff --git a/Source/prob/prob_init.cpp b/Source/prob/prob_init.cpp index fe3ffb211..58bb6ff72 100644 --- a/Source/prob/prob_init.cpp +++ b/Source/prob/prob_init.cpp @@ -55,7 +55,19 @@ void NavierStokes::prob_initData () pp.query("perturbation_amplitude",IC.pertamp); // for Taylor-Green - pp.query("velocity_factor",IC.v_x); + if (probtype == 11) + { + // + // The Taylor-Green amplitude shares IC.v_x with velocity_ic, so it + // is only read for this probtype -- otherwise a stray + // prob.velocity_factor would silently replace the x-component of + // prob.velocity_ic. Reset v_x first so that, conversely, a stray + // prob.velocity_ic cannot stand in for the mandatory + // prob.velocity_factor (checked in init_TaylorGreen). + // + IC.v_x = 0.0; + pp.query("velocity_factor",IC.v_x); + } pp.query("a", IC.a); pp.query("b", IC.b); pp.query("c", IC.c); diff --git a/Tutorials/Bubble/inputs.2d.bubble b/Tutorials/Bubble/inputs.2d.bubble index c1808d526..39d6e6306 100644 --- a/Tutorials/Bubble/inputs.2d.bubble +++ b/Tutorials/Bubble/inputs.2d.bubble @@ -134,6 +134,6 @@ prob.density_ic = 2.0 #******************************************************************************* # Add vorticity to the variables in the plot files. -amr.derive_plot_vars = mag_vort diveru avg_pressure +amr.derive_plot_vars = mag_vort avg_pressure #******************************************************************************* diff --git a/Tutorials/HotSpot/inputs.2d.average_hotspot b/Tutorials/HotSpot/inputs.2d.average_hotspot index fde3c00dc..5f7551fbb 100644 --- a/Tutorials/HotSpot/inputs.2d.average_hotspot +++ b/Tutorials/HotSpot/inputs.2d.average_hotspot @@ -142,7 +142,7 @@ prob.density_ic = 2.0 #******************************************************************************* # Add vorticity to the variables in the plot files. -amr.derive_plot_vars = mag_vort diveru avg_pressure velocity_average +amr.derive_plot_vars = mag_vort avg_pressure velocity_average amr.plot_vars = x_velocity y_velocity density tracer temp #******************************************************************************* diff --git a/Tutorials/HotSpot/inputs.3d.LES_hotspot b/Tutorials/HotSpot/inputs.3d.LES_hotspot index 7553ce4ed..f2be4abb2 100644 --- a/Tutorials/HotSpot/inputs.3d.LES_hotspot +++ b/Tutorials/HotSpot/inputs.3d.LES_hotspot @@ -138,7 +138,7 @@ amr.blocking_factor = 4 #******************************************************************************* # Add vorticity to the variables in the plot files. -amr.derive_plot_vars = mag_vort diveru avg_pressure +amr.derive_plot_vars = mag_vort avg_pressure #******************************************************************************* diff --git a/Tutorials/RayleighTaylor/inputs.2d.rayleightaylor b/Tutorials/RayleighTaylor/inputs.2d.rayleightaylor index 0936c721b..ad518edd2 100644 --- a/Tutorials/RayleighTaylor/inputs.2d.rayleightaylor +++ b/Tutorials/RayleighTaylor/inputs.2d.rayleightaylor @@ -131,7 +131,7 @@ amr.blocking_factor = 8 #******************************************************************************* # Add vorticity to the variables in the plot files. -amr.derive_plot_vars = mag_vort diveru avg_pressure +amr.derive_plot_vars = mag_vort avg_pressure #******************************************************************************* diff --git a/Tutorials/TaylorGreen/benchmarks/ViscBench.cpp b/Tutorials/TaylorGreen/benchmarks/ViscBench.cpp index 8fa089440..c25e78aa2 100644 --- a/Tutorials/TaylorGreen/benchmarks/ViscBench.cpp +++ b/Tutorials/TaylorGreen/benchmarks/ViscBench.cpp @@ -136,8 +136,14 @@ main (int argc, const int nComp = AMREX_SPACEDIM+1; //amrDataI.NComp(); const Real time = amrDataI.Time(); const int finestLevel = amrDataI.FinestLevel(); - const Vector& derives = amrDataI.PlotVarNames(); + // + // Components the exact solution provides, in the order FORT_VISCBENCH + // fills them. Plotfile data are looked up by these names, not by + // position, so a plotfile written with a different amr.plot_vars + // ordering still gives the right error norms. + // + Vector varNames{ AMREX_D_DECL("x_velocity", "y_velocity", "z_velocity"), "density"}; // // Compute the error @@ -169,7 +175,7 @@ main (int argc, // for (int iComp=0; iComp varNames{ AMREX_D_DECL("x_velocity", "y_velocity", "z_velocity"), "density"}; // if (!errFile.empty()) WritePlotFile(error, amrDataI, errFile, verbose, varNames); diff --git a/Util/ConvertCheckpoint/ConvertCheckpointGrids.cpp b/Util/ConvertCheckpoint/ConvertCheckpointGrids.cpp index aa3760e98..3b2a4e9b2 100644 --- a/Util/ConvertCheckpoint/ConvertCheckpointGrids.cpp +++ b/Util/ConvertCheckpoint/ConvertCheckpointGrids.cpp @@ -635,8 +635,6 @@ static void ConvertData() { falRef_trgt.geom.setPeriodicity({{AMREX_D_DECL(is_periodic_array[0],is_periodic_array[1],is_periodic_array[2])}}); - DistributionMapping dm_trgt{new_grids}; - int ngrow_loc; for (int n = 0; n < falRef_src.state.size(); n++){ @@ -693,8 +691,16 @@ static void ConvertData() { OldData_src -> copy(*(falRef_src.state[n].old_data),0,0,ncomps,0,ngrow_loc); } - MultiFab * NewData_trgt = new MultiFab(new_grids_state,dm_trgt,ncomps,ngrow_loc); - MultiFab * OldData_trgt = new MultiFab(new_grids_state,dm_trgt,ncomps,ngrow_loc); + // + // Box i of new_grids_state is the refinement/coarsening of box i of + // save_grids_state, so the target MultiFabs are built on the source + // DistributionMapping: the loops below index source and target + // through a single MFIter, which is only valid if the two share a + // mapping (an independently built one may order the boxes + // differently under MPI). + // + MultiFab * NewData_trgt = new MultiFab(new_grids_state,dm,ncomps,ngrow_loc); + MultiFab * OldData_trgt = new MultiFab(new_grids_state,dm,ncomps,ngrow_loc); NewData_trgt->setVal(10.); OldData_trgt->setVal(10.); @@ -713,13 +719,28 @@ static void ConvertData() { // The copies above are non-periodic, so ghost cells at periodic // domain boundaries still hold the setVal(10.) sentinel. Fill them // from valid data before computing the interpolation slopes. - // NOTE: ghosts at non-periodic physical boundaries still hold the - // sentinel -- the state descriptors (and hence the real BCs) - // are never restored by this tool. // NewData_src->FillBoundary(cgeom.periodicity()); OldData_src->FillBoundary(cgeom.periodicity()); + // + // Ghost cells not covered by this level's grids -- coarse-fine + // boundaries on levels > 0 and non-periodic physical boundaries -- + // still hold the sentinel. The state descriptors (and hence the + // real BCs) are never restored by this tool, so fill them by + // first-order extrapolation from the adjacent valid data; otherwise + // the limited slopes next to those faces are computed from 10. + // Nodal data (pressure) only reads nodes the box owns, so it needs + // no fill. + // + if (NewData_src->is_cell_centered()) + { + Extrapolater::FirstOrderExtrap(*NewData_src, cgeom, 0, ncomps); + if (has_old) { + Extrapolater::FirstOrderExtrap(*OldData_src, cgeom, 0, ncomps); + } + } + for (MFIter mfi(*NewData_trgt); mfi.isValid(); ++mfi) { FArrayBox& ffab = (*NewData_trgt)[mfi];