Skip to content

Commit

Permalink
SAM & BndryPlane Clean up (erf-model#1461)
Browse files Browse the repository at this point in the history
* SAM mods from Wei Zhang convo. Also, first pass at fixing the bndry planes issue.

* Fixed Bndry read case.
  • Loading branch information
AMLattanzi authored Feb 28, 2024
1 parent c2b16bb commit 84b658a
Show file tree
Hide file tree
Showing 14 changed files with 150 additions and 150 deletions.
14 changes: 7 additions & 7 deletions Exec/ABL/inputs.read
Original file line number Diff line number Diff line change
Expand Up @@ -9,12 +9,12 @@ fabarray.mfiter_tile_size = 1024 1024 1024
# ORIGINAL PROBLEM SIZE & GEOMETRY
geometry.prob_lo = 0. 0. 0.
geometry.prob_hi = 1024. 1024. 1024.
amr.n_cell = 32 32 32
amr.n_cell = 32 32 32

# THIS PROBLEM SIZE & GEOMETRY
geometry.prob_lo = 256. 256. 0.
geometry.prob_hi = 768. 768. 1024.
amr.n_cell = 16 16 32
geometry.prob_lo = 256. 256. 0.
geometry.prob_hi = 768. 768. 1024.
amr.n_cell = 16 16 32

geometry.is_periodic = 0 1 0

Expand All @@ -30,7 +30,7 @@ erf.fixed_dt = 2.0e-2 # fixed time step depending on grid resolution
# DIAGNOSTICS & VERBOSITY
erf.sum_interval = 1 # timesteps between computing mass
erf.v = 1 # verbosity in ERF.cpp
amr.v = 1 # verbosity in Amr.cpp
amr.v = 1 # verbosity in Amr.cpp

# REFINEMENT / REGRIDDING
amr.max_level = 0 # maximum level number allowed
Expand Down Expand Up @@ -58,7 +58,7 @@ erf.init_type = "uniform"
# PROBLEM PARAMETERS
prob.rho_0 = 1.0
prob.A_0 = 1.0

prob.T_0 = 300.0
prob.U_0 = 10.0
prob.V_0 = 0.0
prob.W_0 = 0.0
Expand All @@ -71,4 +71,4 @@ prob.W_0_Pert_Mag = 0.0

erf.input_bndry_planes = 1
erf.bndry_file = "BndryFiles"
erf.bndry_input_var_names = density temperature velocity
erf.bndry_input_var_names = temperature density velocity
10 changes: 5 additions & 5 deletions Exec/ABL/inputs.write
Original file line number Diff line number Diff line change
Expand Up @@ -7,9 +7,9 @@ amrex.fpe_trap_invalid = 1
fabarray.mfiter_tile_size = 1024 1024 1024

# PROBLEM SIZE & GEOMETRY
geometry.prob_extent = 1024 1024 1024
amr.n_cell = 64 64 64
amr.n_cell = 32 32 32
geometry.prob_lo = 0. 0. 0.
geometry.prob_hi = 1024. 1024. 1024.
amr.n_cell = 32 32 32

geometry.is_periodic = 1 1 0

Expand All @@ -23,7 +23,7 @@ erf.fixed_dt = 2.0e-2 # fixed time step depending on grid resolution
# DIAGNOSTICS & VERBOSITY
erf.sum_interval = 1 # timesteps between computing mass
erf.v = 1 # verbosity in ERF.cpp
amr.v = 1 # verbosity in Amr.cpp
amr.v = 1 # verbosity in Amr.cpp

# REFINEMENT / REGRIDDING
amr.max_level = 0 # maximum level number allowed
Expand Down Expand Up @@ -51,7 +51,7 @@ erf.init_type = "uniform"
# PROBLEM PARAMETERS
prob.rho_0 = 1.0
prob.A_0 = 1.0

prob.T_0 = 300.0
prob.U_0 = 10.0
prob.V_0 = 0.0
prob.W_0 = 0.0
Expand Down
2 changes: 1 addition & 1 deletion Exec/ABL_input_sounding/GNUmakefile
Original file line number Diff line number Diff line change
Expand Up @@ -24,7 +24,7 @@ DEBUG = FALSE
TEST = TRUE
USE_ASSERTION = TRUE

USE_POISSON_SOLVE = TRUE
USE_POISSON_SOLVE = FALSE

# GNU Make
Bpack := ./Make.package
Expand Down
39 changes: 11 additions & 28 deletions Source/BoundaryConditions/BoundaryConditions_bndryreg.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -60,24 +60,15 @@ ERF::fill_from_bndryregs (const Vector<MultiFab*>& mfs, const Real time)
for (MFIter mfi(mf); mfi.isValid(); ++mfi)
{
const Array4<Real>& dest_arr = mf.array(mfi);
Box bx = mfi.validbox();
Box bx = mfi.growntilebox();

// x-faces
{
Box bx_xlo(bx);
bx_xlo.setSmall(0,dom_lo.x-1); bx_xlo.setBig(0,dom_lo.x-1);
bx_xlo.setSmall(1,dom_lo.y); bx_xlo.setBig(1,dom_hi.y);
bx_xlo.setSmall(2,dom_lo.z); bx_xlo.setBig(2,dom_hi.z);
if (var_idx == Vars::xvel) {
bx_xlo.setSmall(0,dom_lo.x); bx_xlo.setBig(0,dom_lo.x);
} else {
bx_xlo.setSmall(0,dom_lo.x-1); bx_xlo.setBig(0,dom_lo.x-1);
}

Box bx_xhi(bx);
bx_xhi.setSmall(1,dom_lo.y ); bx_xhi.setBig(1,dom_hi.y );
bx_xhi.setSmall(2,dom_lo.z ); bx_xhi.setBig(2,dom_hi.z );
bx_xhi.setSmall(0,dom_hi.x+1); bx_xhi.setBig(0,dom_hi.x+1);
Box bx_xlo(bx); bx_xlo.setBig(0,dom_lo.x-1);
if (var_idx == Vars::xvel) bx_xlo.setBig(0,dom_lo.x);

Box bx_xhi(bx); bx_xhi.setSmall(0,dom_hi.x+1);
if (var_idx == Vars::xvel) bx_xhi.setSmall(0,dom_hi.x);

ParallelFor(
bx_xlo, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) {
Expand All @@ -95,19 +86,11 @@ ERF::fill_from_bndryregs (const Vector<MultiFab*>& mfs, const Real time)

// y-faces
{
Box bx_ylo(bx);
bx_ylo.setSmall(0,dom_lo.x); bx_ylo.setBig(0,dom_hi.x);
bx_ylo.setSmall(1,dom_lo.y); bx_ylo.setBig(1,dom_lo.y);
if (var_idx == Vars::yvel) {
bx_ylo.setSmall(1,dom_lo.y); bx_ylo.setBig(1,dom_lo.y);
} else {
bx_ylo.setSmall(1,dom_lo.y-1); bx_ylo.setBig(1,dom_lo.y-1);
}

Box bx_yhi(bx);
bx_yhi.setSmall(0,dom_lo.x ); bx_yhi.setBig(0,dom_hi.x);
bx_yhi.setSmall(2,dom_lo.z ); bx_yhi.setBig(2,dom_hi.z);
bx_yhi.setSmall(1,dom_hi.y+1); bx_yhi.setBig(1,dom_hi.y+1);
Box bx_ylo(bx); bx_ylo.setBig (1,dom_lo.y-1);
if (var_idx == Vars::yvel) bx_ylo.setBig(1,dom_lo.y);

Box bx_yhi(bx); bx_yhi.setSmall(1,dom_hi.y+1);
if (var_idx == Vars::yvel) bx_yhi.setSmall(1,dom_hi.y);

ParallelFor(
bx_ylo, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) {
Expand Down
14 changes: 9 additions & 5 deletions Source/ERF.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -473,6 +473,13 @@ ERF::InitData ()
// Initialize the start time for our CPU-time tracker
startCPUTime = amrex::ParallelDescriptor::second();

// Create the ReadBndryPlanes object so we can read boundary plane data
// m_r2d is used by init_bcs so we must instantiate this class before
if (input_bndry_planes) {
amrex::Print() << "Defining r2d for the first time " << std::endl;
m_r2d = std::make_unique<ReadBndryPlanes>(geom[0], solverChoice.rdOcp);
}

// Map the words in the inputs file to BC types, then translate
// those types into what they mean for each variable
init_bcs();
Expand Down Expand Up @@ -605,14 +612,11 @@ ERF::InitData ()
}

if (input_bndry_planes) {
// Create the ReadBndryPlanes object so we can handle reading of boundary plane data
amrex::Print() << "Defining r2d for the first time " << std::endl;
m_r2d = std::make_unique< ReadBndryPlanes>(geom[0], solverChoice.rdOcp);

// Read the "time.dat" file to know what data is available
m_r2d->read_time_file();

amrex::Real dt_dummy = 1.e200;
// We haven't populated dt yet, set to 0 to ensure assert doesn't crash
amrex::Real dt_dummy = 0.0;
m_r2d->read_input_files(t_new[0],dt_dummy,m_bc_extdir_vals);
}

Expand Down
20 changes: 11 additions & 9 deletions Source/IO/ERF_ReadBndryPlanes.H
Original file line number Diff line number Diff line change
Expand Up @@ -18,21 +18,23 @@ class ReadBndryPlanes
{

public:
explicit ReadBndryPlanes(const amrex::Geometry& geom,
const amrex::Real& rdOcp_in);
explicit ReadBndryPlanes (const amrex::Geometry& geom,
const amrex::Real& rdOcp_in);

void define_level_data(int lev);
void define_level_data (int lev);

void read_time_file();
void read_time_file ();

void read_input_files(amrex::Real time, amrex::Real dt,
amrex::Array<amrex::Array<amrex::Real, AMREX_SPACEDIM*2>, AMREX_SPACEDIM+NVAR_max> m_bc_extdir_vals);
void read_input_files (amrex::Real time,
amrex::Real dt,
amrex::Array<amrex::Array<amrex::Real, AMREX_SPACEDIM*2>, AMREX_SPACEDIM+NVAR_max> m_bc_extdir_vals);

void read_file(int idx, amrex::Vector<std::unique_ptr<PlaneVector>>& data_to_fill,
amrex::Array<amrex::Array<amrex::Real, AMREX_SPACEDIM*2>, AMREX_SPACEDIM+NVAR_max> m_bc_extdir_vals);
void read_file (int idx,
amrex::Vector<std::unique_ptr<PlaneVector>>& data_to_fill,
amrex::Array<amrex::Array<amrex::Real, AMREX_SPACEDIM*2>, AMREX_SPACEDIM+NVAR_max> m_bc_extdir_vals);

// Return the pointer to PlaneVectors at time "time"
amrex::Vector<std::unique_ptr<PlaneVector>>& interp_in_time(const amrex::Real& time);
amrex::Vector<std::unique_ptr<PlaneVector>>& interp_in_time (const amrex::Real& time);

[[nodiscard]] amrex::Real tinterp() const { return m_tinterp; }

Expand Down
52 changes: 34 additions & 18 deletions Source/IO/ERF_ReadBndryPlanes.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -11,8 +11,8 @@ using namespace amrex;
/**
* Return closest index (from lower) of value in vector
*/
AMREX_FORCE_INLINE int
closest_index(const Vector<Real>& vec, const Real value)
AMREX_FORCE_INLINE
int closest_index (const Vector<Real>& vec, const Real value)
{
auto const it = std::upper_bound(vec.begin(), vec.end(), value);
AMREX_ALWAYS_ASSERT(it != vec.end());
Expand All @@ -24,7 +24,8 @@ closest_index(const Vector<Real>& vec, const Real value)
/**
* Return offset vector
*/
AMREX_FORCE_INLINE IntVect offset(const int face_dir, const int normal)
AMREX_FORCE_INLINE
IntVect offset (const int face_dir, const int normal)
{
IntVect offset(IntVect::TheDimensionVector(normal));
if (face_dir == 1) {
Expand All @@ -39,7 +40,7 @@ AMREX_FORCE_INLINE IntVect offset(const int face_dir, const int normal)
* Function in ReadBndryPlanes class for allocating space
* for the boundary plane data ERF will need.
*/
void ReadBndryPlanes::define_level_data(int /*lev*/)
void ReadBndryPlanes::define_level_data (int /*lev*/)
{
amrex::Print() << "ReadBndryPlanes::define_level_data" << std::endl;
// *********************************************************
Expand Down Expand Up @@ -80,7 +81,7 @@ void ReadBndryPlanes::define_level_data(int /*lev*/)
* @param time Constant specifying the time for interpolation
*/
Vector<std::unique_ptr<PlaneVector>>&
ReadBndryPlanes::interp_in_time(const Real& time)
ReadBndryPlanes::interp_in_time (const Real& time)
{
AMREX_ALWAYS_ASSERT(m_tn <= time && time <= m_tnp2);

Expand Down Expand Up @@ -136,7 +137,7 @@ ReadBndryPlanes::interp_in_time(const Real& time)
* @param geom Geometry for the domain
* @param rdOcp_in Real constant for the Rhydberg constant ($R_d$) divided by the specific heat at constant pressure ($c_p$)
*/
ReadBndryPlanes::ReadBndryPlanes(const Geometry& geom, const Real& rdOcp_in)
ReadBndryPlanes::ReadBndryPlanes (const Geometry& geom, const Real& rdOcp_in)
:
m_geom(geom),
m_rdOcp(rdOcp_in)
Expand All @@ -158,7 +159,7 @@ ReadBndryPlanes::ReadBndryPlanes(const Geometry& geom, const Real& rdOcp_in)
is_q1_read = 0;
is_q2_read = 0;
is_KE_read = 0;
is_QKE_read = 0;
is_QKE_read = 0;

if (pp.contains("bndry_input_var_names"))
{
Expand Down Expand Up @@ -194,7 +195,7 @@ ReadBndryPlanes::ReadBndryPlanes(const Geometry& geom, const Real& rdOcp_in)
* Function in ReadBndryPlanes class for reading the external file
* specifying time data and broadcasting this data across MPI ranks.
*/
void ReadBndryPlanes::read_time_file()
void ReadBndryPlanes::read_time_file ()
{
BL_PROFILE("ERF::ReadBndryPlanes::read_time_file");

Expand Down Expand Up @@ -264,8 +265,9 @@ void ReadBndryPlanes::read_time_file()
* @param dt Current timestep
* @param m_bc_extdir_vals Container storing the external dirichlet boundary conditions we are reading from the input files
*/
void ReadBndryPlanes::read_input_files(Real time, Real dt,
Array<Array<Real, AMREX_SPACEDIM*2>,AMREX_SPACEDIM+NVAR_max> m_bc_extdir_vals)
void ReadBndryPlanes::read_input_files (Real time,
Real dt,
Array<Array<Real, AMREX_SPACEDIM*2>,AMREX_SPACEDIM+NVAR_max> m_bc_extdir_vals)
{
BL_PROFILE("ERF::ReadBndryPlanes::read_input_files");

Expand Down Expand Up @@ -337,8 +339,9 @@ void ReadBndryPlanes::read_input_files(Real time, Real dt,
* @param data_to_fill Container for face data on boundaries
* @param m_bc_extdir_vals Container storing the external dirichlet boundary conditions we are reading from the input files
*/
void ReadBndryPlanes::read_file(const int idx, Vector<std::unique_ptr<PlaneVector>>& data_to_fill,
Array<Array<Real, AMREX_SPACEDIM*2>,AMREX_SPACEDIM+NVAR_max> m_bc_extdir_vals)
void ReadBndryPlanes::read_file (const int idx,
Vector<std::unique_ptr<PlaneVector>>& data_to_fill,
Array<Array<Real, AMREX_SPACEDIM*2>,AMREX_SPACEDIM+NVAR_max> m_bc_extdir_vals)
{
const int t_step = m_in_timesteps[idx];
const std::string chkname1 = m_filename + Concatenate("/bndry_output", t_step);
Expand Down Expand Up @@ -374,14 +377,26 @@ void ReadBndryPlanes::read_file(const int idx, Vector<std::unique_ptr<PlaneVecto
if (ori.coordDir() < 2) {
FArrayBox& d = (*data_to_fill[ori])[lev];
const auto& bx = d.box();
Array4<Real> d_arr = d.array();
Array4<Real> d_arr = d.array();
ParallelFor(
bx, ncomp_for_bc, [=] AMREX_GPU_DEVICE(int i, int j, int k, int n) noexcept {
d_arr(i,j,k,n) = 0.;
});
}
}

// Read density for primitive to conserved conversions
std::string filenamer = MultiFabFileFullPrefix(lev, chkname1, level_prefix, "density");
BndryRegister bndry_r(ba, dm, m_in_rad, m_out_rad, m_extent_rad, 1);
bndry_r.setVal(1.0e13);
for (OrientationIter oit; oit != nullptr; ++oit) {
auto ori = oit();
if (ori.coordDir() < 2) {
std::string facenamer = Concatenate(filenamer + '_', ori, 1);
bndry_r[ori].read(facenamer);
}
}

for (int ivar = 0; ivar < m_var_names.size(); ivar++)
{
std::string var_name = m_var_names[ivar];
Expand Down Expand Up @@ -437,8 +452,9 @@ void ReadBndryPlanes::read_file(const int idx, Vector<std::unique_ptr<PlaneVecto
for (MFIter mfi(bndryMF); mfi.isValid(); ++mfi) {

const auto& vbx = mfi.validbox();
const auto& bndry_read_arr = bndry[ori].array(mfi);
const auto& bndry_mf_arr = bndryMF.array(mfi);
const auto& bndry_read_arr = bndry[ori].array(mfi);
const auto& bndry_read_r_arr = bndry_r[ori].array(mfi);
const auto& bndry_mf_arr = bndryMF.array(mfi);

const auto& bx = bbx & vbx;
if (bx.isEmpty()) {
Expand All @@ -455,8 +471,8 @@ void ReadBndryPlanes::read_file(const int idx, Vector<std::unique_ptr<PlaneVecto
if (var_name == "temperature") {
ParallelFor(
bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
Real R1 = bndry_read_arr(i, j, k, n_for_density);
Real R2 = bndry_read_arr(i+v_offset[0],j+v_offset[1],k+v_offset[2],n_for_density);
Real R1 = bndry_read_r_arr(i, j, k, 0);
Real R2 = bndry_read_r_arr(i+v_offset[0],j+v_offset[1],k+v_offset[2],0);
Real T1 = bndry_read_arr(i, j, k, 0);
Real T2 = bndry_read_arr(i+v_offset[0],j+v_offset[1],k+v_offset[2],0);
Real Th1 = getThgivenRandT(R1,T1,rdOcp);
Expand All @@ -467,7 +483,7 @@ void ReadBndryPlanes::read_file(const int idx, Vector<std::unique_ptr<PlaneVecto
var_name == "KE" || var_name == "QKE") {
ParallelFor(
bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
Real R1 = bndry_read_arr(i, j, k, n_for_density);
Real R1 = bndry_read_r_arr(i, j, k, 0);
Real R2 = bndry_read_arr(i+v_offset[0],j+v_offset[1],k+v_offset[2],n_for_density);
bndry_mf_arr(i, j, k, 0) = 0.5 *
( R1 * bndry_read_arr(i, j, k, 0) +
Expand Down
8 changes: 4 additions & 4 deletions Source/IO/ERF_WriteBndryPlanes.H
Original file line number Diff line number Diff line change
Expand Up @@ -14,11 +14,11 @@
class WriteBndryPlanes
{
public:
explicit WriteBndryPlanes(amrex::Vector<amrex::BoxArray>& grids,
amrex::Vector<amrex::Geometry>& geom);
explicit WriteBndryPlanes (amrex::Vector<amrex::BoxArray>& grids,
amrex::Vector<amrex::Geometry>& geom);

void write_planes(int t_step, amrex::Real time,
amrex::Vector<amrex::Vector<amrex::MultiFab>>& vars_new);
void write_planes (int t_step, amrex::Real time,
amrex::Vector<amrex::Vector<amrex::MultiFab>>& vars_new);

private:

Expand Down
11 changes: 6 additions & 5 deletions Source/IO/ERF_WriteBndryPlanes.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -15,7 +15,7 @@ using namespace amrex;
* @param b1 Boundary register containing the source data
* @param b2 Boundary register containing the destination data
*/
void br_shift(OrientationIter oit, const BndryRegister& b1, BndryRegister& b2)
void br_shift (OrientationIter oit, const BndryRegister& b1, BndryRegister& b2)
{
auto ori = oit();
int ncomp = b1[ori].nComp();
Expand All @@ -33,8 +33,9 @@ void br_shift(OrientationIter oit, const BndryRegister& b1, BndryRegister& b2)
Array4<Real const> src_arr = b1[ori][idx].const_array();
int ioff = bx1.smallEnd(0) - bx2.smallEnd(0);
int joff = bx1.smallEnd(1) - bx2.smallEnd(1);
ParallelFor(
bx2, ncomp, [=] AMREX_GPU_DEVICE(int i, int j, int k, int n) noexcept {
ParallelFor(bx2, ncomp,
[=] AMREX_GPU_DEVICE(int i, int j, int k, int n) noexcept
{
dest_arr(i,j,k,n) = src_arr(i+ioff,j+joff,k,n);
});
}
Expand Down Expand Up @@ -114,8 +115,8 @@ WriteBndryPlanes::WriteBndryPlanes(Vector<BoxArray>& grids,
* @param time Current time
* @param vars_new Grid data for all variables across the AMR hierarchy
*/
void WriteBndryPlanes::write_planes(const int t_step, const Real time,
Vector<Vector<MultiFab>>& vars_new)
void WriteBndryPlanes::write_planes (const int t_step, const Real time,
Vector<Vector<MultiFab>>& vars_new)
{
BL_PROFILE("ERF::WriteBndryPlanes::write_planes");

Expand Down
Loading

0 comments on commit 84b658a

Please sign in to comment.