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];