Skip to content

Commit

Permalink
Remove mixed-type coupling (#1319)
Browse files Browse the repository at this point in the history
  • Loading branch information
@asalmgren
asalmgren authored Nov 24, 2023
1 parent 162677b commit 40e64ed
Show file tree
Hide file tree
Showing 7 changed files with 83 additions and 140 deletions.
12 changes: 4 additions & 8 deletions Docs/sphinx_doc/MeshRefinement.rst
View file
Original file line number Diff line number Diff line change
Expand Up @@ -7,8 +7,7 @@
Mesh Refinement
===============

ERF allows both static and dynamic mesh refinement, as well as three different options for
how the coarse and fine data.
ERF allows both static and dynamic mesh refinement, as well as the choice of one-way or two-way coupling.

Note that any tagged region will be covered by one or more boxes. The user may
specify the refinement criteria and/or region to be covered, but not the decomposition of the region into
Expand Down Expand Up @@ -124,11 +123,11 @@ computed by dividing the variable named rhotheta by the variable named density.
Coupling Types
--------------

ERF supports one-way, two-way, and "mixed" coupling between levels; this is a run-time input
ERF supports one-way and two-way coupling between levels; this is a run-time input

::

erf.coupling_type = "OneWay" or "TwoWay" or "Mixed"
erf.coupling_type = "OneWay" or "TwoWay"

By one-way coupling, we mean that between each pair of refinement levels,
the coarse level communicates data to the fine level to serve as boundary conditions
Expand Down Expand Up @@ -157,10 +156,7 @@ the fine level also communicates data back to the coarse level in two ways:
- A "reflux" operation is performed for all cell-centered data; this updates values on the coarser
level outside of regions covered by the finer level.

We define "mixed" coupling as using the two-way coupling algorithm for all cell-centered quantities except for
:math:`\rho` and :math:`\rho \theta.`

We note that all three coupling schemes are conservative for mass because the fluxes for the continuity
We note that both coupling schemes conserve mass because the fluxes for the continuity
equation are the momenta themselves, which are interpolated on faces at the coarse-fine interface.
Other advected quantities which are advanced in conservation form will lose conservation with one-way coupling.
Two-way coupling ensures conservation of the advective contribution to all scalar updates but
Expand Down
6 changes: 1 addition & 5 deletions Source/DataStructs/DataStruct.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -21,7 +21,7 @@ enum struct ABLDriverType {
};

enum struct CouplingType {
OneWay, TwoWay, Mixed
OneWay, TwoWay
};

enum struct TerrainType {
Expand Down Expand Up @@ -198,8 +198,6 @@ struct SolverChoice {
coupling_type = CouplingType::TwoWay;
} else if (coupling_type_string == "OneWay") {
coupling_type = CouplingType::OneWay;
} else if (coupling_type_string == "Mixed") {
coupling_type = CouplingType::Mixed;
} else {
amrex::Abort("Dont know this coupling_type");
}
Expand All @@ -219,8 +217,6 @@ struct SolverChoice {
amrex::Print() << "Using two-way coupling " << std::endl;
} else if (coupling_type == CouplingType::OneWay) {
amrex::Print() << "Using one-way coupling " << std::endl;
} else if (coupling_type == CouplingType::Mixed) {
amrex::Print() << "Using mixed coupling " << std::endl;
}

if (terrain_type == TerrainType::Static) {
Expand Down
2 changes: 2 additions & 0 deletions Source/ERF.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -354,6 +354,8 @@ private:

void init_uniform (int lev);

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

void initialize_integrator (int lev, amrex::MultiFab& cons_mf, amrex::MultiFab& vel_mf);

// Compute a vector of new MultiFabs by copying from valid region and filling ghost cells -
Expand Down
43 changes: 16 additions & 27 deletions Source/ERF.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -318,12 +318,10 @@ ERF::post_timestep (int nstep, Real time, Real dt_lev0)
particleData.Redistribute();
#endif

if (solverChoice.coupling_type == CouplingType::TwoWay ||
solverChoice.coupling_type == CouplingType::Mixed)
if (solverChoice.coupling_type == CouplingType::TwoWay)
{
// If we doing mixed rather than two-way coupling, we only reflux the "slow" variables (not rho and rhotheta)
int src_comp_reflux = (solverChoice.coupling_type == CouplingType::TwoWay) ? 0 : 2;
int num_comp_reflux = NVAR - src_comp_reflux;
int src_comp_reflux = 0;
int num_comp_reflux = NVAR;
bool use_terrain = solverChoice.use_terrain;
for (int lev = finest_level-1; lev >= 0; lev--)
{
Expand Down Expand Up @@ -534,10 +532,9 @@ ERF::InitData ()
}


if (solverChoice.coupling_type == CouplingType::TwoWay ||
solverChoice.coupling_type == CouplingType::Mixed) {
int src_comp_reflux = (solverChoice.coupling_type == CouplingType::TwoWay) ? 0 : 2;
int num_comp_reflux = NVAR - src_comp_reflux;
if (solverChoice.coupling_type == CouplingType::TwoWay) {
int src_comp_reflux = 0;
int num_comp_reflux = NVAR;
AverageDown(src_comp_reflux, num_comp_reflux);
}

Expand Down Expand Up @@ -570,8 +567,7 @@ ERF::InitData ()
}

// Initialize flux registers (whether we start from scratch or restart)
if (solverChoice.coupling_type == CouplingType::TwoWay ||
solverChoice.coupling_type == CouplingType::Mixed) {
if (solverChoice.coupling_type == CouplingType::TwoWay) {
advflux_reg[0] = nullptr;
for (int lev = 1; lev <= finest_level; lev++)
{
Expand Down Expand Up @@ -1138,10 +1134,9 @@ ERF::MakeHorizontalAverages ()
int lev = 0;

// First, average down all levels (if doing two-way coupling)
if (solverChoice.coupling_type == CouplingType::TwoWay ||
solverChoice.coupling_type == CouplingType::Mixed) {
int src_comp_reflux = (solverChoice.coupling_type == CouplingType::TwoWay) ? 0 : 2;
int num_comp_reflux = NVAR - src_comp_reflux;
if (solverChoice.coupling_type == CouplingType::TwoWay) {
int src_comp_reflux = 0;
int num_comp_reflux = NVAR;
AverageDown(src_comp_reflux, num_comp_reflux);
}

Expand Down Expand Up @@ -1273,8 +1268,7 @@ ERF::MakeDiagnosticAverage (Vector<Real>& h_havg, MultiFab& S, int n)
void
ERF::AverageDown (int scomp, int ncomp)
{
AMREX_ALWAYS_ASSERT(solverChoice.coupling_type == CouplingType::TwoWay ||
solverChoice.coupling_type == CouplingType::Mixed);
AMREX_ALWAYS_ASSERT(solverChoice.coupling_type == CouplingType::TwoWay);
for (int lev = finest_level-1; lev >= 0; --lev)
{
AverageDownTo(lev, scomp, ncomp);
Expand All @@ -1285,8 +1279,7 @@ ERF::AverageDown (int scomp, int ncomp)
void
ERF::AverageDownTo (int crse_lev, int scomp, int ncomp) // NOLINT
{
AMREX_ALWAYS_ASSERT(solverChoice.coupling_type == CouplingType::TwoWay ||
solverChoice.coupling_type == CouplingType::Mixed);
AMREX_ALWAYS_ASSERT(solverChoice.coupling_type == CouplingType::TwoWay);

for (int var_idx = 0; var_idx < Vars::NumTypes; ++var_idx) {
const BoxArray& ba(vars_new[crse_lev][var_idx].boxArray());
Expand Down Expand Up @@ -1353,8 +1346,7 @@ ERF::AverageDownTo (int crse_lev, int scomp, int ncomp) // NOLINT
void
ERF::Construct_ERFFillPatchers (int lev)
{
AMREX_ALWAYS_ASSERT(solverChoice.coupling_type == CouplingType::OneWay ||
solverChoice.coupling_type == CouplingType::Mixed);
AMREX_ALWAYS_ASSERT(solverChoice.coupling_type == CouplingType::OneWay);

auto& fine_new = vars_new[lev];
auto& crse_new = vars_new[lev-1];
Expand All @@ -1363,8 +1355,7 @@ ERF::Construct_ERFFillPatchers (int lev)
auto& dm_fine = fine_new[Vars::cons].DistributionMap();
auto& dm_crse = crse_new[Vars::cons].DistributionMap();

// NOTE: if "Mixed", then crse-fine set/relaxation only done on Rho/RhoTheta
int ncomp = (solverChoice.coupling_type == CouplingType::OneWay) ? NVAR : 2;
int ncomp = NVAR;

FPr_c.emplace_back(ba_fine, dm_fine, geom[lev] ,
ba_crse, dm_crse, geom[lev-1],
Expand All @@ -1383,8 +1374,7 @@ ERF::Construct_ERFFillPatchers (int lev)
void
ERF::Define_ERFFillPatchers (int lev)
{
AMREX_ALWAYS_ASSERT(solverChoice.coupling_type == CouplingType::OneWay ||
solverChoice.coupling_type == CouplingType::Mixed);
AMREX_ALWAYS_ASSERT(solverChoice.coupling_type == CouplingType::OneWay);

auto& fine_new = vars_new[lev];
auto& crse_new = vars_new[lev-1];
Expand All @@ -1393,8 +1383,7 @@ ERF::Define_ERFFillPatchers (int lev)
auto& dm_fine = fine_new[Vars::cons].DistributionMap();
auto& dm_crse = crse_new[Vars::cons].DistributionMap();

// NOTE: if "Mixed", then crse-fine set/relaxation only done on Rho/RhoTheta
int ncomp = (solverChoice.coupling_type == CouplingType::OneWay) ? NVAR : 2;
int ncomp = NVAR;

FPr_c[lev-1].Define(ba_fine, dm_fine, geom[lev] ,
ba_crse, dm_crse, geom[lev-1],
Expand Down
145 changes: 56 additions & 89 deletions Source/ERF_make_new_level.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -42,6 +42,8 @@ void ERF::MakeNewLevelFromScratch (int lev, Real time, const BoxArray& ba,
auto& lev_new = vars_new[lev];
auto& lev_old = vars_old[lev];

init_stuff(lev, ba, dm);

lev_new[Vars::cons].define(ba, dm, Cons::NumVars, ngrow_state);
lev_old[Vars::cons].define(ba, dm, Cons::NumVars, ngrow_state);

Expand All @@ -54,18 +56,6 @@ void ERF::MakeNewLevelFromScratch (int lev, Real time, const BoxArray& ba,
lev_new[Vars::zvel].define(convert(ba, IntVect(0,0,1)), dm, 1, IntVect(ngrow_vels,ngrow_vels,0));
lev_old[Vars::zvel].define(convert(ba, IntVect(0,0,1)), dm, 1, IntVect(ngrow_vels,ngrow_vels,0));

// ********************************************************************************************
// These are just used for scratch in the time integrator but we might as well define them here
// ********************************************************************************************
rU_old[lev].define(convert(ba, IntVect(1,0,0)), dm, 1, ngrow_vels);
rU_new[lev].define(convert(ba, IntVect(1,0,0)), dm, 1, ngrow_vels);

rV_old[lev].define(convert(ba, IntVect(0,1,0)), dm, 1, ngrow_vels);
rV_new[lev].define(convert(ba, IntVect(0,1,0)), dm, 1, ngrow_vels);

rW_old[lev].define(convert(ba, IntVect(0,0,1)), dm, 1, ngrow_vels);
rW_new[lev].define(convert(ba, IntVect(0,0,1)), dm, 1, ngrow_vels);

//********************************************************************************************
// Microphysics
// *******************************************************************************************
Expand All @@ -76,29 +66,6 @@ void ERF::MakeNewLevelFromScratch (int lev, Real time, const BoxArray& ba,
// Build the data structures for calculating diffusive/turbulent terms
update_arrays(lev, ba, dm);

// ********************************************************************************************
// Map factors
// ********************************************************************************************
BoxList bl2d = ba.boxList();
for (auto& b : bl2d) {
b.setRange(2,0);
}
BoxArray ba2d(std::move(bl2d));

mapfac_m[lev] = std::make_unique<MultiFab>(ba2d,dm,1,3);
mapfac_u[lev] = std::make_unique<MultiFab>(convert(ba2d,IntVect(1,0,0)),dm,1,3);
mapfac_v[lev] = std::make_unique<MultiFab>(convert(ba2d,IntVect(0,1,0)),dm,1,3);
if (solverChoice.test_mapfactor) {
mapfac_m[lev]->setVal(0.5);
mapfac_u[lev]->setVal(0.5);
mapfac_v[lev]->setVal(0.5);
}
else {
mapfac_m[lev]->setVal(1.);
mapfac_u[lev]->setVal(1.);
mapfac_v[lev]->setVal(1.);
}

// ********************************************************************************************
// Base state holds r_0, pres_0, pi_0 (in that order)
// ********************************************************************************************
Expand Down Expand Up @@ -148,15 +115,6 @@ 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;

// ********************************************************************************************
// Define Theta_prim storage if using MOST BC
// ********************************************************************************************
if (phys_bc_type[Orientation(Direction::z,Orientation::low)] == ERF_BC::MOST) {
Theta_prim[lev] = std::make_unique<MultiFab>(ba,dm,1,IntVect(ngrow_state,ngrow_state,0));
} else {
Theta_prim[lev] = nullptr;
}

// ********************************************************************************************
// Initialize the integrator class
// ********************************************************************************************
Expand Down Expand Up @@ -217,17 +175,7 @@ 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());

// ********************************************************************************************
// These are just used for scratch in the time integrator but we might as well define them here
// ********************************************************************************************
rU_old[lev].define(convert(ba, IntVect(1,0,0)), dm, 1, crse_new[Vars::xvel].nGrowVect());
rU_new[lev].define(convert(ba, IntVect(1,0,0)), dm, 1, crse_new[Vars::xvel].nGrowVect());

rV_old[lev].define(convert(ba, IntVect(0,1,0)), dm, 1, crse_new[Vars::yvel].nGrowVect());
rV_new[lev].define(convert(ba, IntVect(0,1,0)), dm, 1, crse_new[Vars::yvel].nGrowVect());

rW_old[lev].define(convert(ba, IntVect(0,0,1)), dm, 1, crse_new[Vars::zvel].nGrowVect());
rW_new[lev].define(convert(ba, IntVect(0,0,1)), dm, 1, crse_new[Vars::zvel].nGrowVect());
init_stuff(lev, ba, dm);

t_new[lev] = time;
t_old[lev] = time - 1.e200;
Expand Down Expand Up @@ -269,17 +217,7 @@ 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));

// ********************************************************************************************
// These are just used for scratch in the time integrator but we might as well define them here
// ********************************************************************************************
rU_old[lev].define(convert(ba, IntVect(1,0,0)), dm, 1, ngrow_vels);
rU_new[lev].define(convert(ba, IntVect(1,0,0)), dm, 1, ngrow_vels);

rV_old[lev].define(convert(ba, IntVect(0,1,0)), dm, 1, ngrow_vels);
rV_new[lev].define(convert(ba, IntVect(0,1,0)), dm, 1, ngrow_vels);

rW_old[lev].define(convert(ba, IntVect(0,0,1)), dm, 1, ngrow_vels);
rW_new[lev].define(convert(ba, IntVect(0,0,1)), dm, 1, ngrow_vels);
init_stuff(lev,ba,dm);

// ********************************************************************************************
// This will fill the temporary MultiFabs with data from vars_new
Expand Down Expand Up @@ -319,29 +257,6 @@ ERF::RemakeLevel (int lev, Real time, const BoxArray& ba, const DistributionMapp
base_state[lev].setVal(0.);
initHSE(lev);

// ********************************************************************************************
// Map factors
// ********************************************************************************************
BoxList bl2d = ba.boxList();
for (auto& b : bl2d) {
b.setRange(2,0);
}
BoxArray ba2d(std::move(bl2d));

mapfac_m[lev] = std::make_unique<MultiFab>(ba2d,dm,1,3);
mapfac_u[lev] = std::make_unique<MultiFab>(convert(ba2d,IntVect(1,0,0)),dm,1,3);
mapfac_v[lev] = std::make_unique<MultiFab>(convert(ba2d,IntVect(0,1,0)),dm,1,3);
if (solverChoice.test_mapfactor) {
mapfac_m[lev]->setVal(0.5);
mapfac_u[lev]->setVal(0.5);
mapfac_v[lev]->setVal(0.5);
}
else {
mapfac_m[lev]->setVal(1.);
mapfac_u[lev]->setVal(1.);
mapfac_v[lev]->setVal(1.);
}

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

Expand Down Expand Up @@ -465,6 +380,58 @@ ERF::initialize_integrator (int lev, MultiFab& cons_mf, MultiFab& vel_mf)
z_phys_nd[lev], detJ_cc[lev]);
}

void ERF::init_stuff(int lev, const BoxArray& ba, const DistributionMapping& dm)
{
// The number of ghost cells for density must be 1 greater than that for velocity
// so that we can go back in forth betwen velocity and momentum on all faces
int ngrow_state = ComputeGhostCells(solverChoice.advChoice, solverChoice.use_NumDiff) + 1;
int ngrow_vels = ComputeGhostCells(solverChoice.advChoice, solverChoice.use_NumDiff);

// ********************************************************************************************
// These are just used for scratch in the time integrator but we might as well define them here
// ********************************************************************************************
rU_old[lev].define(convert(ba, IntVect(1,0,0)), dm, 1, ngrow_vels);
rU_new[lev].define(convert(ba, IntVect(1,0,0)), dm, 1, ngrow_vels);

rV_old[lev].define(convert(ba, IntVect(0,1,0)), dm, 1, ngrow_vels);
rV_new[lev].define(convert(ba, IntVect(0,1,0)), dm, 1, ngrow_vels);

rW_old[lev].define(convert(ba, IntVect(0,0,1)), dm, 1, ngrow_vels);
rW_new[lev].define(convert(ba, IntVect(0,0,1)), dm, 1, ngrow_vels);

// ********************************************************************************************
// Define Theta_prim storage if using MOST BC
// ********************************************************************************************
if (phys_bc_type[Orientation(Direction::z,Orientation::low)] == ERF_BC::MOST) {
Theta_prim[lev] = std::make_unique<MultiFab>(ba,dm,1,IntVect(ngrow_state,ngrow_state,0));
} else {
Theta_prim[lev] = nullptr;
}

// ********************************************************************************************
// Map factors
// ********************************************************************************************
BoxList bl2d = ba.boxList();
for (auto& b : bl2d) {
b.setRange(2,0);
}
BoxArray ba2d(std::move(bl2d));

mapfac_m[lev] = std::make_unique<MultiFab>(ba2d,dm,1,3);
mapfac_u[lev] = std::make_unique<MultiFab>(convert(ba2d,IntVect(1,0,0)),dm,1,3);
mapfac_v[lev] = std::make_unique<MultiFab>(convert(ba2d,IntVect(0,1,0)),dm,1,3);
if (solverChoice.test_mapfactor) {
mapfac_m[lev]->setVal(0.5);
mapfac_u[lev]->setVal(0.5);
mapfac_v[lev]->setVal(0.5);
}
else {
mapfac_m[lev]->setVal(1.);
mapfac_u[lev]->setVal(1.);
mapfac_v[lev]->setVal(1.);
}
}

//
// Delete level data (overrides the pure virtual function in AmrCore)
//
Expand Down
Loading

0 comments on commit 40e64ed

Please sign in to comment.