Skip to content

Commit

Permalink
fix for different regridding options
Browse files Browse the repository at this point in the history
  • Loading branch information
@asalmgren
asalmgren committed Nov 28, 2023
1 parent 40e64ed commit a6029a9
Show file tree
Hide file tree
Showing 8 changed files with 161 additions and 61 deletions.
3 changes: 2 additions & 1 deletion Exec/inputs.mesh_refinement
View file
Original file line number Diff line number Diff line change
Expand Up @@ -5,9 +5,10 @@
# REFINEMENT / REGRIDDING
amr.max_level = 1 # maximum level number allowed
amr.ref_ratio_vect = 2 2 1 3 3 1 # refinement ratio written as lev0_x lev0_y lev0_z lev1_x lev1_y lev1_z ...
amr.regrid_int = 2 2 2 2 # how often to regrid
amr.n_error_buf = 2 2 2 2 # number of buffer cells in error est

erf.regrid_int = 2 2 2 2 # how often to regrid

# How to specify a region for static refinement
erf.refinement_indicators = box1
erf.box1.in_box_lo = .15 .25
Expand Down
2 changes: 1 addition & 1 deletion Source/BoundaryConditions/ERF_FillPatch.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -73,7 +73,7 @@ ERF::FillPatch (int lev, Real time, const Vector<MultiFab*>& mfs, bool fillset)
} // var_idx

// Coarse-Fine set region
if (lev>0 && solverChoice.coupling_type == CouplingType::OneWay && cf_set_width>0 && fillset) {
if (lev>0 && cf_set_width>0 && fillset) {
FPr_c[lev-1].FillSet(*mfs[Vars::cons], time, null_bc, domain_bcs_type);
FPr_u[lev-1].FillSet(*mfs[Vars::xvel], time, null_bc, domain_bcs_type);
FPr_v[lev-1].FillSet(*mfs[Vars::yvel], time, null_bc, domain_bcs_type);
Expand Down
2 changes: 1 addition & 1 deletion Source/ERF.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -336,7 +336,7 @@ private:
// set covered coarse cells to be the average of overlying fine cells
void AverageDown (int scomp, int ncomp);

void update_arrays (int lev, const amrex::BoxArray& ba, const amrex::DistributionMapping& dm);
void update_diffusive_arrays (int lev, const amrex::BoxArray& ba, const amrex::DistributionMapping& dm);

void update_terrain_arrays (int lev, amrex::Real time);

Expand Down
43 changes: 25 additions & 18 deletions Source/ERF.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -681,21 +681,24 @@ ERF::InitData ()
// Fill ghost cells/faces
for (int lev = 0; lev <= finest_level; ++lev)
{
//
// NOTE: we must set up the FillPatcher object before calling FillPatch at a fine level
//
if (solverChoice.coupling_type != CouplingType::TwoWay && cf_width>0 && lev>0) {
if (lev > 0 && cf_width > 0) {
Construct_ERFFillPatchers(lev);
Register_ERFFillPatchers(lev);
}

auto& lev_new = vars_new[lev];

//
// We don't use the FillPatcher in this call because
// we don't need to fill the interior data at this point.
//
bool fillset = false;
auto& lev_new = vars_new[lev];
FillPatch(lev, t_new[lev],
{&lev_new[Vars::cons],&lev_new[Vars::xvel],&lev_new[Vars::yvel],&lev_new[Vars::zvel]});
{&lev_new[Vars::cons],&lev_new[Vars::xvel],&lev_new[Vars::yvel],&lev_new[Vars::zvel]},
fillset);

// We need to fill the ghost cell values of the base state in case it wasn't
// done in the initialization
//
// We fill the ghost cell values of the base state in case it wasn't done in the initialization
//
base_state[lev].FillBoundary(geom[lev].periodicity());

// For moving terrain only
Expand Down Expand Up @@ -1102,9 +1105,9 @@ ERF::ReadParameters ()
// Query the set and total widths for crse-fine interior ghost cells
pp.query("cf_width", cf_width);
pp.query("cf_set_width", cf_set_width);
AMREX_ALWAYS_ASSERT(cf_width >= 0);
AMREX_ALWAYS_ASSERT(cf_set_width >= 0);
AMREX_ALWAYS_ASSERT(cf_width >= cf_set_width);
if (cf_width < 0 || cf_set_width < 0 || cf_width < cf_set_width) {
amrex::Abort("You must set cf_width >= cf_set_width >= 0");
}

// AmrMesh iterate on grids?
bool iterate(true);
Expand All @@ -1124,6 +1127,10 @@ ERF::ReadParameters ()
if (verbose > 0) {
solverChoice.display();
}

if (solverChoice.coupling_type == CouplingType::TwoWay && cf_width > 0) {
amrex::Abort("For two-way coupling you must set cf_width = 0");
}
}

// Create horizontal average quantities for 5 variables:
Expand Down Expand Up @@ -1402,12 +1409,12 @@ ERF::Define_ERFFillPatchers (int lev)
void
ERF::Register_ERFFillPatchers (int lev)
{
auto& lev_new = vars_new[lev-1];
auto& lev_old = vars_old[lev-1];
FPr_c[lev-1].RegisterCoarseData({&lev_old[Vars::cons], &lev_new[Vars::cons]}, {t_old[lev-1], t_new[lev-1]});
FPr_u[lev-1].RegisterCoarseData({&lev_old[Vars::xvel], &lev_new[Vars::xvel]}, {t_old[lev-1], t_new[lev-1]});
FPr_v[lev-1].RegisterCoarseData({&lev_old[Vars::yvel], &lev_new[Vars::yvel]}, {t_old[lev-1], t_new[lev-1]});
FPr_w[lev-1].RegisterCoarseData({&lev_old[Vars::zvel], &lev_new[Vars::zvel]}, {t_old[lev-1], t_new[lev-1]});
auto& lev_new = vars_new[lev];
auto& lev_old = vars_old[lev];
FPr_c[lev].RegisterCoarseData({&lev_old[Vars::cons], &lev_new[Vars::cons]}, {t_old[lev], t_new[lev]});
FPr_u[lev].RegisterCoarseData({&lev_old[Vars::xvel], &lev_new[Vars::xvel]}, {t_old[lev], t_new[lev]});
FPr_v[lev].RegisterCoarseData({&lev_old[Vars::yvel], &lev_new[Vars::yvel]}, {t_old[lev], t_new[lev]});
FPr_w[lev].RegisterCoarseData({&lev_old[Vars::zvel], &lev_new[Vars::zvel]}, {t_old[lev], t_new[lev]});
}

#ifdef ERF_USE_MULTIBLOCK
Expand Down
122 changes: 114 additions & 8 deletions Source/ERF_make_new_level.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -63,8 +63,10 @@ void ERF::MakeNewLevelFromScratch (int lev, Real time, const BoxArray& ba,
qmoist[lev].define(ba, dm, 6, ngrow_state); // qv, qc, qi, qr, qs, qg
#endif

// ********************************************************************************************
// Build the data structures for calculating diffusive/turbulent terms
update_arrays(lev, ba, dm);
// ********************************************************************************************
update_diffusive_arrays(lev, ba, dm);

// ********************************************************************************************
// Base state holds r_0, pres_0, pi_0 (in that order)
Expand Down Expand Up @@ -115,6 +117,21 @@ void ERF::MakeNewLevelFromScratch (int lev, Real time, const BoxArray& ba,
sst_lev[lev].resize(1); sst_lev[lev][0] = nullptr;
lmask_lev[lev].resize(1); lmask_lev[lev][0] = nullptr;

// ********************************************************************************************
// Create the ERFFillPatcher object
// ********************************************************************************************
if (cf_width > 0 && lev > 0) {
Construct_ERFFillPatchers(lev);
Define_ERFFillPatchers(lev);
}

// ********************************************************************************************
// Initialize the boundary conditions
// ********************************************************************************************
physbcs[lev] = std::make_unique<ERFPhysBCFunct> (lev, geom[lev], domain_bcs_type, domain_bcs_type_d,
solverChoice.terrain_type, m_bc_extdir_vals, m_bc_neumann_vals,
z_phys_nd[lev], detJ_cc[lev]);

// ********************************************************************************************
// Initialize the integrator class
// ********************************************************************************************
Expand Down Expand Up @@ -175,21 +192,73 @@ ERF::MakeNewLevelFromCoarse (int lev, Real time, const BoxArray& ba,
lev_new[Vars::zvel].define(convert(ba, IntVect(0,0,1)), dm, 1, crse_new[Vars::zvel].nGrowVect());
lev_old[Vars::zvel].define(convert(ba, IntVect(0,0,1)), dm, 1, crse_new[Vars::zvel].nGrowVect());

//********************************************************************************************
// Microphysics
// *******************************************************************************************
#if defined(ERF_USE_MOISTURE)
int ngrow_state = ComputeGhostCells(solverChoice.advChoice, solverChoice.use_NumDiff) + 1;
qmoist[lev].define(ba, dm, 6, ngrow_state); // qv, qc, qi, qr, qs, qg
#endif

init_stuff(lev, ba, dm);

t_new[lev] = time;
t_old[lev] = time - 1.e200;

// ********************************************************************************************
// Build the data structures for terrain-related quantities
// ********************************************************************************************
update_terrain_arrays(lev, time);

// ********************************************************************************************
// Update the base state at this level
// ********************************************************************************************
base_state[lev].define(ba,dm,3,1);
base_state[lev].setVal(0.);
initHSE(lev);

// ********************************************************************************************
// Build the data structures for calculating diffusive/turbulent terms
update_arrays(lev, ba, dm);
// ********************************************************************************************
update_diffusive_arrays(lev, ba, dm);

// ********************************************************************************************
// Initialize flux register at new level
// ********************************************************************************************
if (solverChoice.coupling_type == CouplingType::TwoWay) {
advflux_reg.resize(lev+1);
advflux_reg[lev] = new YAFluxRegister(ba, grids[lev-1],
dm, dmap[lev-1],
geom[lev], geom[lev-1],
ref_ratio[lev-1], lev, NVAR);
}

// *****************************************************************************************************
// Initialize the boundary conditions (after initializing the terrain but before calling FillCoarsePatch
// *****************************************************************************************************
physbcs[lev] = std::make_unique<ERFPhysBCFunct> (lev, geom[lev], domain_bcs_type, domain_bcs_type_d,
solverChoice.terrain_type, m_bc_extdir_vals, m_bc_neumann_vals,
z_phys_nd[lev], detJ_cc[lev]);

// ********************************************************************************************
// Fill data at the new level by interpolation from the coarser level
// ********************************************************************************************
FillCoarsePatch(lev, time, {&lev_new[Vars::cons],&lev_new[Vars::xvel],
&lev_new[Vars::yvel],&lev_new[Vars::zvel]});

// ********************************************************************************************
// Initialize the integrator class
// ********************************************************************************************
dt_mri_ratio[lev] = dt_mri_ratio[lev-1];
initialize_integrator(lev, lev_new[Vars::cons], lev_new[Vars::xvel]);

// ********************************************************************************************
// If we are making a new level then the FillPatcher for this level hasn't been allocated yet
// ********************************************************************************************
if (cf_width > 0) {
Construct_ERFFillPatchers(lev);
Define_ERFFillPatchers(lev);
}
}

// Remake an existing level using provided BoxArray and DistributionMapping and
Expand Down Expand Up @@ -217,8 +286,29 @@ ERF::RemakeLevel (int lev, Real time, const BoxArray& ba, const DistributionMapp
temp_lev_new[Vars::zvel].define(convert(ba, IntVect(0,0,1)), dm, 1, IntVect(ngrow_vels,ngrow_vels,0));
temp_lev_old[Vars::zvel].define(convert(ba, IntVect(0,0,1)), dm, 1, IntVect(ngrow_vels,ngrow_vels,0));

// ********************************************************************************************
// Remake the flux registers
// ********************************************************************************************
if (solverChoice.coupling_type == CouplingType::TwoWay) {
delete advflux_reg[lev];
advflux_reg[lev] = new YAFluxRegister(ba, grids[lev-1],
dm, dmap[lev-1],
geom[lev], geom[lev-1],
ref_ratio[lev-1], lev, NVAR);
}

//********************************************************************************************
// Microphysics
// *******************************************************************************************
#if defined(ERF_USE_MOISTURE)
qmoist[lev].define(ba, dm, 6, ngrow_state); // qv, qc, qi, qr, qs, qg
#endif

init_stuff(lev,ba,dm);

BoxArray ba_old(vars_new[lev][Vars::cons].boxArray());
DistributionMapping dm_old(vars_new[lev][Vars::cons].DistributionMap());

// ********************************************************************************************
// This will fill the temporary MultiFabs with data from vars_new
// ********************************************************************************************
Expand All @@ -245,7 +335,7 @@ ERF::RemakeLevel (int lev, Real time, const BoxArray& ba, const DistributionMapp
t_old[lev] = time - 1.e200;

// Build the data structures for calculating diffusive/turbulent terms
update_arrays(lev, ba, dm);
update_diffusive_arrays(lev, ba, dm);

// Build the data structures for terrain-related quantities
update_terrain_arrays(lev, time);
Expand All @@ -257,12 +347,32 @@ ERF::RemakeLevel (int lev, Real time, const BoxArray& ba, const DistributionMapp
base_state[lev].setVal(0.);
initHSE(lev);

// ********************************************************************************************
// Initialize the boundary conditions
// ********************************************************************************************
physbcs[lev] = std::make_unique<ERFPhysBCFunct> (lev, geom[lev], domain_bcs_type, domain_bcs_type_d,
solverChoice.terrain_type, m_bc_extdir_vals, m_bc_neumann_vals,
z_phys_nd[lev], detJ_cc[lev]);

// ********************************************************************************************
// Initialize the integrator class
// ********************************************************************************************

initialize_integrator(lev, vars_new[lev][Vars::cons], vars_new[lev][Vars::xvel]);

// We need to re-define the FillPatcher if the grids have changed
if (lev > 0 && cf_width > 0) {
bool ba_changed = (ba != ba_old);
bool dm_changed = (dm != dm_old);
if (ba_changed || dm_changed) {
Define_ERFFillPatchers(lev);
}
}
}


void
ERF::update_arrays (int lev, const BoxArray& ba, const DistributionMapping& dm)
ERF::update_diffusive_arrays (int lev, const BoxArray& ba, const DistributionMapping& dm)
{
// ********************************************************************************************
// Diffusive terms
Expand Down Expand Up @@ -374,10 +484,6 @@ ERF::initialize_integrator (int lev, MultiFab& cons_mf, MultiFab& vel_mf)
mri_integrator_mem[lev]->setNoSubstepping(solverChoice.no_substepping);
mri_integrator_mem[lev]->setIncompressible(solverChoice.incompressible);
mri_integrator_mem[lev]->setForceFirstStageSingleSubstep(solverChoice.force_stage1_single_substep);

physbcs[lev] = std::make_unique<ERFPhysBCFunct> (lev, geom[lev], domain_bcs_type, domain_bcs_type_d,
solverChoice.terrain_type, m_bc_extdir_vals, m_bc_neumann_vals,
z_phys_nd[lev], detJ_cc[lev]);
}

void ERF::init_stuff(int lev, const BoxArray& ba, const DistributionMapping& dm)
Expand Down
9 changes: 6 additions & 3 deletions Source/IO/Plotfile.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -104,11 +104,14 @@ ERF::WritePlotFile (int which, Vector<std::string> plot_var_names)
if (ncomp_mf == 0)
return;

// We fillpatch here because some of the derived quantities require derivatives
// which require ghost cells to be filled
// We Fillpatch here because some of the derived quantities require derivatives
// which require ghost cells to be filled. We do not need to call FillPatcher
// because we don't need to set interior fine points.
for (int lev = 0; lev <= finest_level; ++lev) {
bool fillset = false;
FillPatch(lev, t_new[lev], {&vars_new[lev][Vars::cons], &vars_new[lev][Vars::xvel],
&vars_new[lev][Vars::yvel], &vars_new[lev][Vars::zvel]});
&vars_new[lev][Vars::yvel], &vars_new[lev][Vars::zvel]},
fillset);
}

Vector<MultiFab> mf(finest_level+1);
Expand Down
35 changes: 9 additions & 26 deletions Source/TimeIntegration/ERF_TimeStep.cpp
View file
Original file line number Diff line number Diff line change
@@ -1,5 +1,4 @@
#include <ERF.H>
#include <TileNoZ.H>
#include <Utils.H>

using namespace amrex;
Expand All @@ -25,26 +24,15 @@ ERF::timeStep (int lev, Real time, int iteration)
// regrid changes level "lev+1" so we don't regrid on max_level
// also make sure we don't regrid fine levels again if
// it was taken care of during a coarser regrid
if (lev < max_level && istep[lev] > last_regrid_step[lev])
if (lev < max_level)
{
if (istep[lev] % regrid_int == 0)
if ( (istep[lev] % regrid_int == 0) && (istep[lev] > last_regrid_step[lev]) )
{
// regrid could add newly refine levels (if finest_level < max_level)
// so we save the previous finest level index
int old_finest = finest_level;
regrid(lev, time);

// NOTE: Define & Register index lev backwards (so we add 1 here)
// Redefine & register the ERFFillpatcher objects
if (solverChoice.coupling_type != CouplingType::TwoWay && cf_width > 0)
{
Define_ERFFillPatchers(lev+1);
Register_ERFFillPatchers(lev+1);
if (lev < max_level-1) {
Define_ERFFillPatchers(lev+2);
Register_ERFFillPatchers(lev+2);
}
}
regrid(lev, time);

// mark that we have regridded this level already
for (int k = lev; k <= finest_level; ++k) {
Expand All @@ -55,8 +43,8 @@ ERF::timeStep (int lev, Real time, int iteration)
for (int k = old_finest+1; k <= finest_level; ++k) {
dt[k] = dt[k-1] / MaxRefRatio(k-1);
}
}
}
} // if
} // lev
}

// Update what we call "old" and "new" time
Expand All @@ -72,9 +60,10 @@ ERF::timeStep (int lev, Real time, int iteration)
// Advance a single level for a single time step
Advance(lev, time, dt[lev], iteration, nsubsteps[lev]);

#ifdef ERF_USE_PARTICLES
particleData.Redistribute();
#endif
// We store the old and new data at level "lev" once we have finished the full advance at this level
if (finest_level > 0 && lev < finest_level && cf_set_width > 0) {
Register_ERFFillPatchers(lev);
}

++istep[lev];

Expand All @@ -92,12 +81,6 @@ ERF::timeStep (int lev, Real time, int iteration)
Real strt_time_for_fine = time + (i-1)*dt[lev+1];
timeStep(lev+1, strt_time_for_fine, i);
}

if (solverChoice.coupling_type == CouplingType::TwoWay) {
int scomp = 0;
int ncomp = NVAR;
AverageDownTo(lev, scomp, ncomp); // average lev+1 down to lev
}
}
}

Expand Down
Loading

0 comments on commit a6029a9

Please sign in to comment.