ProblemBase Class Reference

#include <ERF_ProbCommon.H>

Collaboration diagram for ProblemBase:

Public Member Functions
virtual	~ProblemBase ()=default
virtual void	erf_init_dens_hse_dry (amrex::MultiFab &, std::unique_ptr< amrex::MultiFab > &, std::unique_ptr< amrex::MultiFab > &, amrex::Geometry const &, const amrex::Vector< amrex::Real > &, bool, bool)
virtual void	erf_init_const_dens_hse (amrex::MultiFab &)
virtual void	erf_init_const_dens_and_th_hse (amrex::MultiFab &, amrex::MultiFab &, amrex::MultiFab &, amrex::MultiFab &, amrex::MultiFab &, amrex::Real)
virtual void	erf_init_const_dens_and_linear_th_hse (amrex::MultiFab &, amrex::MultiFab &, amrex::MultiFab &, amrex::MultiFab &, amrex::MultiFab &, amrex::Real, std::unique_ptr< amrex::MultiFab > &)
virtual void	erf_init_dens_hse_moist (amrex::MultiFab &, std::unique_ptr< amrex::MultiFab > &, amrex::Geometry const &, const amrex::Vector< amrex::Real > &)
virtual void	init_custom_pert (const amrex::Box &, amrex::Array4< amrex::Real const > const &, amrex::Array4< amrex::Real > const &, amrex::Array4< amrex::Real > const &, amrex::Array4< amrex::Real > const &, amrex::Array4< amrex::Real const > const &, amrex::Array4< amrex::Real const > const &, amrex::Array4< amrex::Real const > const &, amrex::GeometryData const &, amrex::Array4< amrex::Real const > const &, const SolverChoice &, const int)
virtual void	init_custom_pert_vels (const amrex::Box &, const amrex::Box &, const amrex::Box &, amrex::Array4< amrex::Real > const &, amrex::Array4< amrex::Real > const &, amrex::Array4< amrex::Real > const &, amrex::Array4< amrex::Real const > const &, amrex::Array4< amrex::Real const > const &, amrex::GeometryData const &, amrex::Array4< amrex::Real const > const &, amrex::Array4< amrex::Real const > const &, const SolverChoice &, const int)
virtual void	update_rhotheta_sources (const double &, amrex::MultiFab *src, const amrex::Geometry &, std::unique_ptr< amrex::MultiFab > &, const amrex::Vector< amrex::Real > &)
virtual void	update_rhoqt_sources (const double &, amrex::MultiFab *qsrc, const amrex::Geometry &, std::unique_ptr< amrex::MultiFab > &, const amrex::Vector< amrex::Real > &)
virtual void	update_w_subsidence (const double &, amrex::Vector< amrex::Real > &wbar, amrex::Gpu::DeviceVector< amrex::Real > &d_wbar, const amrex::MultiFab &, const amrex::Geometry &geom, std::unique_ptr< amrex::MultiFab > &, const amrex::Vector< amrex::Real > &)
virtual void	update_geostrophic_profile (const double &, amrex::Vector< amrex::Real > &u_geos, amrex::Gpu::DeviceVector< amrex::Real > &d_u_geos, amrex::Vector< amrex::Real > &v_geos, amrex::Gpu::DeviceVector< amrex::Real > &d_v_geos, const amrex::Geometry &geom, std::unique_ptr< amrex::MultiFab > &, const amrex::Vector< amrex::Real > &)
void	init_buildings_surface (const amrex::Geometry &geom, amrex::FArrayBox &buildings_fab, const double &time)
void	init_terrain_surface (const amrex::Geometry &geom, amrex::FArrayBox &terrain_fab, const double &time)
virtual bool	terrain_is_specified ()
void	read_custom_terrain (const std::string &fname, const bool is_usgs, const amrex::Geometry &geom, amrex::FArrayBox &terrain_fab, const double &, bool require_domain_coverage=true, bool zero_outside_table=false)
virtual void	init_custom_terrain (const amrex::Geometry &geom, amrex::FArrayBox &terrain_fab, const double &time)
virtual void	erf_init_rayleigh (amrex::Vector< amrex::Vector< amrex::Real > > &, amrex::Geometry const &, std::unique_ptr< amrex::MultiFab > &, amrex::Real)
void	set_zlevels_stag (const amrex::Vector< amrex::Vector< amrex::Real >> &zlevels_stag)

Static Public Member Functions
static void	make_terrain_surface_periodic (const amrex::Geometry &geom, amrex::FArrayBox &terrain_fab)

Protected Member Functions
void	init_base_parms (amrex::Real rho_0, amrex::Real T_0)
std::string	name ()
amrex::Real	phys_zlo (const amrex::GeometryData &geomdata, int lev) const
	Physical bottom of the domain at AMR level lev. More...
amrex::Real	phys_zhi (const amrex::GeometryData &geomdata, int lev) const
	Physical top of the domain at AMR level lev. More...

Static Protected Member Functions
static bool	use_physical_z ()
	True when heights must come from the vertical levels or z_phys, not from the geometry. More...

Protected Attributes
ProbParmDefaults	base_parms
amrex::Vector< amrex::Vector< amrex::Real > >	m_zlevels_stag
	Staggered vertical levels of AMR level lev (set by set_zlevels_stag) More...

Detailed Description

Class to hold problem-specific routines

Constructor & Destructor Documentation

~ProblemBase()

virtual ProblemBase::~ProblemBase ( )

virtualdefault

Virtual destructor to avoid data leakage with derived class

Member Function Documentation

erf_init_const_dens_and_linear_th_hse()

virtual void ProblemBase::erf_init_const_dens_and_linear_th_hse ( amrex::MultiFab &  , amrex::MultiFab &  , amrex::MultiFab &  , amrex::MultiFab &  , amrex::MultiFab &  , amrex::Real  , std::unique_ptr< amrex::MultiFab > &    )

inlinevirtual

    {
        amrex::Error("Should never call this version of erf_init_const_dens_and_linear_th_hse for "+name()+" problem");
    }

Here is the call graph for this function:

erf_init_const_dens_and_th_hse()

virtual void ProblemBase::erf_init_const_dens_and_th_hse ( amrex::MultiFab &  , amrex::MultiFab &  , amrex::MultiFab &  , amrex::MultiFab &  , amrex::MultiFab &  , amrex::Real    )

inlinevirtual

    {
        amrex::Error("Should never call this version of erf_init_const_dens_and_th_hse for "+name()+" problem");
    }

Here is the call graph for this function:

erf_init_const_dens_hse()

virtual void ProblemBase::erf_init_const_dens_hse ( amrex::MultiFab &  )

inlinevirtual

    {
        amrex::Error("Should never call this version of erf_init_const_dens_hse for "+name()+" problem");
    }

Here is the call graph for this function:

erf_init_dens_hse_dry()

virtual void ProblemBase::erf_init_dens_hse_dry ( amrex::MultiFab &  , std::unique_ptr< amrex::MultiFab > &  , std::unique_ptr< amrex::MultiFab > &  , amrex::Geometry const &  , const amrex::Vector< amrex::Real > &  , bool  , bool    )

inlinevirtual

Function to initialize the hydrostatic reference density

Parameters

[out]	rho_hse	hydrostatic reference density
[in]	z_phys_nd	height coordinate at nodes
[in]	z_phys_cc	height coordinate at cell centers
[in]	geom	container for geometric information

    {
        amrex::Print() << "Hydrostatically balanced density was NOT set"
            << " -- an appropriate init_type should probably have been specified"
            << " (e.g., input_sounding, WRFInput, or Metgrid)"
            << std::endl;
        amrex::Error("Should never call this version of erf_init_dens_hse_dry for "+name()+" problem");
    }

Here is the call graph for this function:

erf_init_dens_hse_moist()

virtual void ProblemBase::erf_init_dens_hse_moist ( amrex::MultiFab &  , std::unique_ptr< amrex::MultiFab > &  , amrex::Geometry const &  , const amrex::Vector< amrex::Real > &    )

inlinevirtual

    {
 
    }

erf_init_rayleigh()

virtual void ProblemBase::erf_init_rayleigh ( amrex::Vector< amrex::Vector< amrex::Real > > &  , amrex::Geometry const &  , std::unique_ptr< amrex::MultiFab > &  , amrex::Real    )

inlinevirtual

Function to define the quantities needed to impose Rayleigh damping

Parameters

[out]	rayleigh_ptrs	= {strength of Rayleigh damping, reference values for xvel/yvel/zvel/theta used to define Rayleigh damping}
[in]	geom	container for geometric information
[in]	z_phys_cc	height coordinate at cell centers

    {
        // Default which does no harm
    }

init_base_parms()

void ProblemBase::init_base_parms ( amrex::Real  rho_0, amrex::Real  T_0  )

inlineprotected

Function to update default base parameters, currently only used for init_type == InitType::Uniform

                                                          {
        base_parms.rho_0 = rho_0;
        base_parms.T_0 = T_0;
    }

init_buildings_surface()

void ProblemBase::init_buildings_surface ( const amrex::Geometry &  geom, amrex::FArrayBox &  buildings_fab, const double &  time  )

inline

Function to perform custom initialization of buildings

Parameters

[in]	geom	container for geometric information
[out]	z_phys_nd	height coordinate at nodes
[in]	time	current time

    {
        // Check if a valid text file exists for the buildings
        std::string fname;
        amrex::ParmParse pp("erf");
        auto valid_fname = pp.queryAdd("buildings_file_name",fname);
 
        if (valid_fname) {
            // Buildings may occupy only a local patch of the physical domain.
            read_custom_terrain(fname,false,geom,buildings_fab,time,
                                /*require_domain_coverage=*/false,
                                /*zero_outside_table=*/true);
        } else {
            init_my_custom_terrain(geom, buildings_fab, time);
        }
 
        make_terrain_surface_periodic(geom, buildings_fab);
    }

Here is the call graph for this function:

init_custom_pert()

virtual void ProblemBase::init_custom_pert ( const amrex::Box &  , amrex::Array4< amrex::Real const > const &  , amrex::Array4< amrex::Real > const &  , amrex::Array4< amrex::Real > const &  , amrex::Array4< amrex::Real > const &  , amrex::Array4< amrex::Real const > const &  , amrex::Array4< amrex::Real const > const &  , amrex::Array4< amrex::Real const > const &  , amrex::GeometryData const &  , amrex::Array4< amrex::Real const > const &  , const SolverChoice &  , const int    )

inlinevirtual

Function to perform custom initialization of a test problem

Parameters

[in]	bx	cell-centered box on which to initialize scalars
[out]	state	cell-centered variables to be filled in this routine
[out]	r_hse	hydrostatic reference density
[out]	p_hse	hydrostatic reference pressure
[in]	z_nd	height coordinate at nodes
[in]	z_cc	height coordinate at cell centers
[in]	z_sfc_nd	height of the ground at the nodes of the bottom face of each column, at k = 0, over the planar footprint of bx
[in]	mf_m	map factor on cell centers
[in]	mf_u	map factor on x-faces
[in]	mf_v	map factor on y-faces
[in]	sc	SolverChoice structure that carries parameters
[in]	lev	AMR level

    {
        amrex::Print() << "No perturbation to background fields supplied for "
            << name() << " problem" << std::endl;
    }

Here is the call graph for this function:

init_custom_pert_vels()

virtual void ProblemBase::init_custom_pert_vels ( const amrex::Box &  , const amrex::Box &  , const amrex::Box &  , amrex::Array4< amrex::Real > const &  , amrex::Array4< amrex::Real > const &  , amrex::Array4< amrex::Real > const &  , amrex::Array4< amrex::Real const > const &  , amrex::Array4< amrex::Real const > const &  , amrex::GeometryData const &  , amrex::Array4< amrex::Real const > const &  , amrex::Array4< amrex::Real const > const &  , const SolverChoice &  , const int    )

inlinevirtual

Function to perform custom initialization of the velocities

Parameters

[in]	xbx	box on which to initialize x_vel_pert
[in]	ybx	box on which to initialize y_vel_pert
[in]	zbx	box on which to initialize z_vel_pert
[out]	x_vel_pert	x-component of velocity perturbation to be filled in this routine
[out]	y_vel_pert	y-component of velocity perturbation to be filled in this routine
[out]	z_vel_pert	z-component of velocity perturbation to be filled in this routine
[in]	z_nd	height coordinate at nodes
[in]	z_sfc_nd	height of the ground at the nodes of the bottom face of each column, at k = 0 (see init_custom_pert)
[in]	mf_m	map factor on cell centers
[in]	mf_u	map factor on x-faces
[in]	mf_v	map factor on y-faces
[in]	sc	SolverChoice structure that carries parameters
[in]	lev	AMR level

    {
        amrex::Print() << "No perturbation velocities supplied for " << name() << " problem" << std::endl;
    }

Here is the call graph for this function:

init_custom_terrain()

virtual void ProblemBase::init_custom_terrain ( const amrex::Geometry &  geom, amrex::FArrayBox &  terrain_fab, const double &  time  )

inlinevirtual

Function to perform custom initialization of terrain

Note: Terrain functionality can also be used to provide grid stretching.

Parameters

[in]	geom	container for geometric information
[out]	terrain_mf	height coordinate at nodes
[in]	time	current time

    {
        std::string custom_terrain_type = "None";
        amrex::ParmParse pp_prob("prob"); pp_prob.queryAdd("custom_terrain_type",custom_terrain_type);
 
        if (custom_terrain_type != "None")
        {
            amrex::Print() << "Calling custom terrain initialization" << std::endl;
            init_my_custom_terrain(geom,terrain_fab,time);
 
        } else {
            amrex::Print() << "Initializing flat terrain" << std::endl;
            // A flat surface at the bottom of the domain
            terrain_fab.template setVal<amrex::RunOn::Device>(geom.ProbLo(2));
 
            if (SolverChoice::mesh_type == MeshType::VariableDz) {
                SolverChoice::set_mesh_type(MeshType::StretchedDz);
                amrex::Print() << "Resetting mesh type to StretchedDz" << std::endl;
            }
        }
    }

Here is the call graph for this function:

init_terrain_surface()

void ProblemBase::init_terrain_surface ( const amrex::Geometry &  geom, amrex::FArrayBox &  terrain_fab, const double &  time  )

inline

Function to perform custom initialization of terrain

Parameters

[in]	geom	container for geometric information
[out]	z_phys_nd	height coordinate at nodes
[in]	time	current time

    {
        amrex::ParmParse pp("erf");
 
        // NetCDF terrain file takes priority
        std::string fname_nc;
        auto valid_fname_nc = pp.query("terrain_file_name_nc", fname_nc);
 
        std::string fname, fname_usgs;
        auto valid_fname      = pp.query("terrain_file_name",      fname);
        auto valid_fname_USGS = pp.query("terrain_file_name_USGS", fname_usgs);
 
        if (valid_fname_nc) {
#ifdef ERF_USE_NETCDF
            read_terrain_netcdf(fname_nc, geom, terrain_fab, time);
#else
            amrex::Abort("ERF must be compiled with NetCDF support to use terrain_file_name_nc");
#endif
        } else if (valid_fname) {
            read_custom_terrain(fname, false, geom, terrain_fab, time);
 
        } else if (valid_fname_USGS) {
            read_custom_terrain(fname_usgs, true, geom, terrain_fab, time);
 
        } else {
            init_my_custom_terrain (geom, terrain_fab, time);
        }
 
        make_terrain_surface_periodic(geom, terrain_fab);
    }

Referenced by fill_terrain_surface_slab().

Here is the call graph for this function:

Here is the caller graph for this function:

make_terrain_surface_periodic()

static void ProblemBase::make_terrain_surface_periodic ( const amrex::Geometry &  geom, amrex::FArrayBox &  terrain_fab  )

inlinestatic

Give every node in a periodic direction the height of its periodic image inside the domain. The nodes on the high side of a periodic domain are the same nodes as those on the low side, but none of the sources above knows that: the file readers clip a node outside the table to its edge (and look every node up 1e-3 past where it is), and an analytic surface is periodic only to round-off. Left as they are the two copies of a node differ, and so do the metrics of the faces on either side of it – which means the two grids that share a face update it differently, the mass flux through it is no longer the same on both sides, and which copy a fill keeps depends on how the grids are distributed. The low side is the one we keep.

terrain_fab must hold every node of the domain in a periodic direction, since that is where the images are read from. A surface that is far from periodic is still made periodic here, but we say so, since the cells next to the high side then see a step.

    {
        const bool per_x = geom.isPeriodic(0);
        const bool per_y = geom.isPeriodic(1);
        if (!per_x && !per_y) { return; }
 
        const amrex::Box& domain = geom.Domain();
        const int ilo = domain.smallEnd(0); const int nx = domain.length(0);
        const int jlo = domain.smallEnd(1); const int ny = domain.length(1);
 
        const amrex::Box bx = terrain_fab.box();
        const auto blo = amrex::lbound(bx);
        const auto bhi = amrex::ubound(bx);
        amrex::Array4<amrex::Real> const& z_arr = terrain_fab.array();
 
        AMREX_ALWAYS_ASSERT_WITH_MESSAGE(
            (!per_x || (blo.x <= ilo && bhi.x >= ilo+nx-1)) &&
            (!per_y || (blo.y <= jlo && bhi.y >= jlo+ny-1)),
            "make_terrain_surface_periodic: the terrain FArrayBox must hold every node of the domain in a periodic direction");
 
        // How far the surface is from periodic: the largest difference between a node on the
        // high side of the domain and its image on the low side (ghost nodes are left out:
        // they hold whatever the source clipped them to)
        amrex::Box hi_bx = bx & amrex::surroundingNodes(domain);
        hi_bx.setSmall(2, bx.smallEnd(2)); hi_bx.setBig(2, bx.smallEnd(2));
        if (!hi_bx.isEmpty()) {
            amrex::ReduceOps<amrex::ReduceOpMax> reduce_op;
            amrex::ReduceData<amrex::Real> reduce_data(reduce_op);
            reduce_op.eval(hi_bx, reduce_data,
                [=] AMREX_GPU_DEVICE (int i, int j, int k) -> amrex::GpuTuple<amrex::Real>
                {
                    int ii = (per_x) ? ilo + ((i-ilo) % nx + nx) % nx : i;
                    int jj = (per_y) ? jlo + ((j-jlo) % ny + ny) % ny : j;
                    return { (ii != i || jj != j) ? amrex::Math::abs(z_arr(i,j,k) - z_arr(ii,jj,k))
                                                  : amrex::Real(0.0) };
                });
            // Every rank fills its own copy of the whole surface, so there is nothing to reduce
            // across ranks.  Warn once: with a moving mesh we are called several times a step.
            const amrex::Real max_diff = amrex::get<0>(reduce_data.value());
            static bool warned = false;
            if (!warned && max_diff > amrex::Real(0.1) * geom.CellSize(2)) {
                warned = true;
                amrex::Print() << "WARNING: the terrain surface is not periodic in a periodic direction: "
                               << "a node on the high side of the domain differs from its image on the low "
                               << "side by up to " << max_diff << " m.  The low side is used for both, so the "
                               << "cells next to the high side see a step in the surface.\n";
            }
        }
 
        // We only write nodes outside [lo, lo+n) in a periodic direction and only read
        // nodes inside it, so no node is both read and written here
        amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k)
        {
            int ii = (per_x) ? ilo + ((i-ilo) % nx + nx) % nx : i;
            int jj = (per_y) ? jlo + ((j-jlo) % ny + ny) % ny : j;
            if ((ii != i || jj != j) &&
                ii >= blo.x && ii <= bhi.x && jj >= blo.y && jj <= bhi.y) {
                z_arr(i,j,k) = z_arr(ii,jj,k);
            }
        });
    }

Referenced by init_buildings_surface(), and init_terrain_surface().

Here is the call graph for this function:

Here is the caller graph for this function:

name()

std::string ProblemBase::name ( )

inlineprotected

                      {
        std::string prob_name;
        amrex::ParmParse pp("erf");
        pp.get("prob_name",prob_name); return prob_name;
    }

Referenced by erf_init_const_dens_and_linear_th_hse(), erf_init_const_dens_and_th_hse(), erf_init_const_dens_hse(), erf_init_dens_hse_dry(), init_custom_pert(), init_custom_pert_vels(), update_geostrophic_profile(), update_rhoqt_sources(), update_rhotheta_sources(), and update_w_subsidence().

Here is the call graph for this function:

Here is the caller graph for this function:

phys_zhi()

amrex::Real ProblemBase::phys_zhi ( const amrex::GeometryData &  geomdata, int  lev  ) const

inlineprotected

Physical top of the domain at AMR level lev.

                                                                                        {
        return use_physical_z() ? m_zlevels_stag[lev].back() : geomdata.ProbHi(2);
    }

Here is the call graph for this function:

phys_zlo()

amrex::Real ProblemBase::phys_zlo ( const amrex::GeometryData &  geomdata, int  lev  ) const

inlineprotected

Physical bottom of the domain at AMR level lev.

                                                                                        {
        return use_physical_z() ? m_zlevels_stag[lev].front() : geomdata.ProbLo(2);
    }

Here is the call graph for this function:

read_custom_terrain()

void ProblemBase::read_custom_terrain ( const std::string &  fname, const bool  is_usgs, const amrex::Geometry &  geom, amrex::FArrayBox &  terrain_fab, const double &  , bool  require_domain_coverage = true, bool  zero_outside_table = false  )

inline

    {
        // Read, checked and broadcast once per file; other levels, regrids, the EB build and
        // the station slabs reuse the table
        const auto& table = erf_terrain_file::ReadHeightTable(fname, is_usgs,
                                                              geom.ProbLoArray(), geom.ProbHiArray(),
                                                              require_domain_coverage);
        const int nx = table.nx;
        const int ny = table.ny;
        const amrex::Vector<amrex::Real>& m_xterrain = table.x;
        const amrex::Vector<amrex::Real>& m_yterrain = table.y;
        const amrex::Vector<amrex::Real>& m_zterrain = table.z;
 
        const int nz = nx * ny;
        const int nx_vals = is_usgs ? (nx * ny) : nx;
 
        // Copy data to the GPU
        amrex::Gpu::DeviceVector<amrex::Real> d_xterrain(nx_vals),d_yterrain(ny),d_zterrain(nz);
        amrex::Gpu::copy(amrex::Gpu::hostToDevice, m_xterrain.begin(), m_xterrain.end(), d_xterrain.begin());
        amrex::Gpu::copy(amrex::Gpu::hostToDevice, m_yterrain.begin(), m_yterrain.end(), d_yterrain.begin());
        amrex::Gpu::copy(amrex::Gpu::hostToDevice, m_zterrain.begin(), m_zterrain.end(), d_zterrain.begin());
 
        amrex::Real* d_xt = d_xterrain.data();
        amrex::Real* d_yt = d_yterrain.data();
        amrex::Real* d_zt = d_zterrain.data();
 
        auto dx = geom.CellSizeArray();
        auto ProbLoArr = geom.ProbLoArray();
 
        int ilo = geom.Domain().smallEnd(0);
        int jlo = geom.Domain().smallEnd(1);
        int klo = geom.Domain().smallEnd(2);
        int ihi = geom.Domain().bigEnd(0) + 1;
        int jhi = geom.Domain().bigEnd(1) + 1;
 
        amrex::Box zbx = terrain_fab.box();
        amrex::Array4<amrex::Real> const& z_arr = terrain_fab.array();
 
        amrex::ParallelFor(zbx, [=] AMREX_GPU_DEVICE (int i, int j, int /*k*/)
        {
            // Clip indices for ghost cells
            int ii = amrex::min(amrex::max(i,ilo),ihi);
            int jj = amrex::min(amrex::max(j,jlo),jhi);
 
            // Location of nodes
            amrex::Real x_node = ProbLoArr[0] + ii * dx[0];
            amrex::Real y_node = ProbLoArr[1] + jj * dx[1];
            amrex::Real x = x_node + amrex::Real(1e-3);
            amrex::Real y = y_node + amrex::Real(1e-3);
 
            // A local building table describes no buildings beyond its footprint.
            // Use the unshifted node position so a node on the last table point
            // remains inside, despite the interpolation offset above.
            if (zero_outside_table) {
                const bool out_x = nx > 1 && (x_node < d_xt[0] || x_node > d_xt[nx-1]);
                const bool out_y = ny > 1 && (y_node < d_yt[0] || y_node > d_yt[ny-1]);
                if (out_x || out_y) {
                    z_arr(i,j,klo) = amrex::Real(0.0);
                    return;
                }
            }
 
            int ind11, ind12, ind21, ind22;
            amrex::Real x1, x2, y1, y2;
 
            int iindex_terrain=-1;
            int jindex_terrain=-1;
            // *******************************************************************
            // NOTE: usgs-format is contiguous in x
            // *******************************************************************
            if (is_usgs) {
 
                // NOTE: the searches below must distinguish "the node is below the first
                //       table value" (bracket index -1) from "no table value exceeds the
                //       node" (search never triggered), so we can not use the initial
                //       value of -1 as the not-found flag
                bool jfound = false;
                for (int it = 0; it < ny && !jfound; it++) {
                    if (d_yt[it] > y) {
                        jindex_terrain = it-1; jfound = true;
                    }
                }
                if (!jfound) {
                    jindex_terrain = ny-1;
                }
                // Keep the 2x2 stencil inside the table
                jindex_terrain = amrex::min(amrex::max(jindex_terrain,0), ny-2);
 
                int gstart = (jindex_terrain  )*nx;
                int gend   = (jindex_terrain+1)*nx-1;
                bool ifound = false;
                for (int it = gstart; it <= gend && !ifound; it++) {
                    if (d_xt[it] > x) {
                        iindex_terrain = it-gstart-1; ifound = true;
                    }
                }
                if (!ifound) {
                    iindex_terrain = nx-1;
                }
                // Keep the 2x2 stencil inside the table
                iindex_terrain = amrex::min(amrex::max(iindex_terrain,0), nx-2);
 
                // Define the four values to interpolate between
                ind11 = jindex_terrain*nx + iindex_terrain; // (x1,y1)
                ind12 = ind11+nx;                           // (x1,y2)
                ind21 = ind11+1;                            // (x2,y1)
                ind22 = ind12+1;                            // (x2,y2)
 
                x1 = d_xt[ind11];
                x2 = d_xt[ind21];
                y1 = d_yt[jindex_terrain];
                y2 = d_yt[jindex_terrain+1];
 
                amrex::Real denom = (x2-x1)*(y2-y1);
 
                if (denom == amrex::Real(0.0)) {
                    // Degenerate cell (repeated x or y in the table) -- use the corner value
                    z_arr(i,j,klo) = d_zt[ind11];
                } else {
                    amrex::Real  w_11 = (x2-x)*(y2-y)/denom; // (x1,y1)
                    amrex::Real  w_12 = (x2-x)*(y-y1)/denom; // (x1,y2)
                    amrex::Real  w_21 = (x-x1)*(y2-y)/denom; // (x2,y1)
                    amrex::Real  w_22 = (x-x1)*(y-y1)/denom; // (x2,y2)
 
                    // Do bilinear interpolation
                    z_arr(i,j,klo) = w_11*d_zt[ind11] + w_12*d_zt[ind12] + w_21*d_zt[ind21] + w_22*d_zt[ind22];
                }
 
            } else {
 
                bool jfound = false;
                for (int it = 0; it < ny && !jfound; it++) {
                    if (d_yt[it] > y) {
                        jindex_terrain = it-1; jfound = true;
                    }
                }
                if (!jfound) {
                    // No table value exceeds y -- fall back to the last row
                    jindex_terrain = ny-1;
                }
                jindex_terrain = amrex::max(jindex_terrain,0);
 
                bool ifound = false;
                for (int it = 0; it < nx && !ifound; it++) {
                    if (d_xt[it] > x) {
                        iindex_terrain = it-1; ifound = true;
                    }
                }
                if (!ifound) {
                    // No table value exceeds x -- fall back to the last column
                    iindex_terrain = nx-1;
                }
                iindex_terrain = amrex::max(iindex_terrain,0);
 
                // NOTE: these must be clamped since the fallbacks above can leave
                //       iindex_terrain = nx-1 (jindex_terrain = ny-1)
                int ip1 = amrex::min(iindex_terrain+1,nx-1);
                int jp1 = amrex::min(jindex_terrain+1,ny-1);
 
                // Define the four values to interpolate between
                x1 = d_xt[iindex_terrain];
                x2 = d_xt[ip1];
                y1 = d_yt[jindex_terrain];
                y2 = d_yt[jp1];
 
#if 1
                // *******************************************************************
                // NOTE: this format is contiguous in y to match the AMR-Wind read
                // *******************************************************************
                ind11 = iindex_terrain * ny + jindex_terrain;  // (x1,y1)
                ind21 = ip1            * ny + jindex_terrain;  // (x2,y1)
 
                ind12 = iindex_terrain * ny + jp1;             // (x1,y2)
                ind22 = ip1            * ny + jp1;             // (x2,y2)
#else
                // *******************************************************************
                // NOTE: this format is contiguous in x as an alternative
                // *******************************************************************
 
                ind11 = jindex_terrain * nx + iindex_terrain;  // (x1,y1)
                ind12 = jp1            * nx + iindex_terrain;  // (x1,y2)
 
                ind21 = jindex_terrain * nx + ip1;             // (x2,y1)
                ind22 = jp1            * nx + ip1;             // (x2,y2)
#endif
 
                // Only interpolate in a direction if we have a non-degenerate interval there
                bool interp_x = (ip1 != iindex_terrain) && (x2 != x1);
                bool interp_y = (jp1 != jindex_terrain) && (y2 != y1);
 
                if (!interp_x && !interp_y)
                {
                    z_arr(i,j,klo) = d_zt[ind11];
                }
                else if (interp_x && !interp_y)
                {
                    amrex::Real  w_11 = (x2-x); // (x1,y1)
                    amrex::Real  w_21 = (x-x1); // (x2,y1)
                    amrex::Real denom = (x2-x1);
                    z_arr(i,j,klo) = (w_11*d_zt[ind11] + w_21*d_zt[ind21])/denom;
                }
                else if (!interp_x && interp_y)
                {
                    amrex::Real  w_11 = (y2-y); // (x1,y1)
                    amrex::Real  w_12 = (y-y1); // (x1,y2)
                    amrex::Real denom = (y2-y1);
                    z_arr(i,j,klo) = (w_11*d_zt[ind11] + w_12*d_zt[ind12])/denom;
                }
                else
                {
                    amrex::Real  w_11 = (x2-x)*(y2-y); // (x1,y1)
                    amrex::Real  w_21 = (x-x1)*(y2-y); // (x2,y1)
                    amrex::Real  w_12 = (x2-x)*(y-y1); // (x1,y2)
                    amrex::Real  w_22 = (x-x1)*(y-y1); // (x2,y2)
                    amrex::Real denom = (x2-x1)*(y2-y1);
                    z_arr(i,j,klo) = (w_11*d_zt[ind11] + w_12*d_zt[ind12] + w_21*d_zt[ind21] + w_22*d_zt[ind22]) / denom;
                }
            } // usgs?
        });
 
        // The kernel reads the device vectors above: finish it before they are freed
        amrex::Gpu::streamSynchronize();
    }

Referenced by init_buildings_surface(), and init_terrain_surface().

Here is the call graph for this function:

Here is the caller graph for this function:

set_zlevels_stag()

void ProblemBase::set_zlevels_stag ( const amrex::Vector< amrex::Vector< amrex::Real >> &  zlevels_stag )

inline

Store the staggered vertical levels of every AMR level. On a mesh other than ConstantDz, geometry.prob_lo[2], prob_hi[2] and CellSize(2) only bound and step the computational coordinate; the physical heights are these levels (and z_phys_cc/nd).

                                                                                    {
        m_zlevels_stag = zlevels_stag;
    }

terrain_is_specified()

virtual bool ProblemBase::terrain_is_specified ( )

inlinevirtual

Function to report whether the user has actually specified a terrain surface, i.e. whether init_terrain_surface will do anything other than return flat terrain.

This mirrors the source precedence in init_terrain_surface, and is deliberately kept next to it so that the two cannot drift apart. A problem that overrides the terrain unconditionally, rather than keying off prob.custom_terrain_type, should override this to return true.

    {
        amrex::ParmParse pp("erf");
        std::string fname;
        if (pp.query("terrain_file_name_nc",   fname) ||
            pp.query("terrain_file_name",      fname) ||
            pp.query("terrain_file_name_USGS", fname)) {
            return true;
        }
 
        std::string custom_terrain_type = "None";
        amrex::ParmParse pp_prob("prob");
        pp_prob.query("custom_terrain_type", custom_terrain_type);
        return (custom_terrain_type != "None");
    }

Here is the call graph for this function:

update_geostrophic_profile()

virtual void ProblemBase::update_geostrophic_profile ( const double &  , amrex::Vector< amrex::Real > &  u_geos, amrex::Gpu::DeviceVector< amrex::Real > &  d_u_geos, amrex::Vector< amrex::Real > &  v_geos, amrex::Gpu::DeviceVector< amrex::Real > &  d_v_geos, const amrex::Geometry &  geom, std::unique_ptr< amrex::MultiFab > &  , const amrex::Vector< amrex::Real > &    )

inlinevirtual

Function to update user-specified geostrophic wind profile.

Parameters

[in]	time	current time
[out]	u_geos	geostrophic wind profile
[out]	v_geos	geostrophic wind profile
[in]	geom	container for geometric information
[in]	z_phys_cc	height coordinate at cell centers
[in]	zlevels_stag	nominal heights of the z faces of this level

    {
        if (u_geos.empty()) return;
 
        amrex::Warning("Geostrophic wind profile not defined for "+name()+" problem");
 
        const int khi              = geom.Domain().bigEnd()[2];
        // const amrex::Real* prob_lo = geom.ProbLo();
        // const auto dx              = geom.CellSize();
        for (int k = 0; k <= khi; k++)
        {
            // const amrex::Real z_cc = prob_lo[2] + (k+myhalf)* dx[2];
            // set RHS term of RhoTheta equation based on time, z_cc here
            u_geos[k] = zero;
            v_geos[k] = zero;
        }
 
        // Copy from host version to device version
        amrex::Gpu::copy(amrex::Gpu::hostToDevice, u_geos.begin(), u_geos.end(), d_u_geos.begin());
        amrex::Gpu::copy(amrex::Gpu::hostToDevice, v_geos.begin(), v_geos.end(), d_v_geos.begin());
    }

Here is the call graph for this function:

update_rhoqt_sources()

virtual void ProblemBase::update_rhoqt_sources ( const double &  , amrex::MultiFab *  qsrc, const amrex::Geometry &  , std::unique_ptr< amrex::MultiFab > &  , const amrex::Vector< amrex::Real > &    )

inlinevirtual

Function to update user-specified moisture source terms that can vary with time and height.

Parameters

[in]	time	current time
[out]	rhoqt_source	moisture forcing profile
[in]	geom	container for geometric information
[in]	z_phys_cc	height coordinate at cell centers
[in]	zlevels_stag	nominal heights of the z faces of this level

    {
        if (qsrc->empty()) return;
 
        amrex::Warning("Moisture forcing not defined for "+name()+" problem");
        for ( amrex::MFIter mfi(*qsrc, amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi )
        {
            const auto &box = mfi.tilebox();
            const amrex::Array4<amrex::Real>& qsrc_arr = qsrc->array(mfi);
            // src is a spatial function if erf.spatial_moisture_forcing = true
            // otherwise, qsrc_arr is defined only in over Z (box x and y dimensions are 1)
            ParallelFor(box, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
                // const amrex::Real z_cc = prob_lo[2] + (k+myhalf)* dx[2];
                // set RHS term of RhoQ1 equation based on time, z_cc here
                qsrc_arr(i, j, k) = zero;
            });
        }
    }

Here is the call graph for this function:

update_rhotheta_sources()

virtual void ProblemBase::update_rhotheta_sources ( const double &  , amrex::MultiFab *  src, const amrex::Geometry &  , std::unique_ptr< amrex::MultiFab > &  , const amrex::Vector< amrex::Real > &    )

inlinevirtual

Function to update user-specified temperature source terms that can vary with time and height.

Parameters

[in]	time	current time
[out]	rhotheta_source	forcing profile
[in]	geom	container for geometric information
[in]	z_phys_cc	height coordinate at cell centers
[in]	zlevels_stag	nominal heights of the z faces of this level

    {
        if (src->empty()) return;
 
        amrex::Warning("Temperature forcing not defined for "+name()+" problem");
        for ( amrex::MFIter mfi(*src, amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi )
        {
            const auto &box = mfi.tilebox();
            const amrex::Array4<amrex::Real>& src_arr = src->array(mfi);
            // src is a spatial function if erf.spatial_rhotheta_forcing = true
            // otherwise, qsrc_arr is defined only in over Z (box x and y dimensions are 1)
            ParallelFor(box, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
                // const amrex::Real z_cc = prob_lo[2] + (k+myhalf)* dx[2];
                // set RHS term of RhoTheta equation based on time, z_cc here
                src_arr(i, j, k) = zero;
            });
        }
    }

Here is the call graph for this function:

update_w_subsidence()

virtual void ProblemBase::update_w_subsidence ( const double &  , amrex::Vector< amrex::Real > &  wbar, amrex::Gpu::DeviceVector< amrex::Real > &  d_wbar, const amrex::MultiFab &  , const amrex::Geometry &  geom, std::unique_ptr< amrex::MultiFab > &  , const amrex::Vector< amrex::Real > &    )

inlinevirtual

Function to update the vertical velocity profile, used to add subsidence source terms for x-mom, y-mom, rho*theta, rho*Q1, and rho*Q2.

TODO: Currently, this is only called by InitData, so there is no time dependence.

Parameters

[in]	time	current time
[out]	wbar	w vel forcing profile
[in]	geom	container for geometric information
[in]	z_phys_cc	height coordinate at cell centers
[in]	zlevels_stag	nominal heights of the z faces of this level

    {
        if (wbar.empty()) return;
 
        amrex::Warning("Moisture forcing not defined for "+name()+" problem");
 
        const int khi              = geom.Domain().bigEnd()[2];
        // const amrex::Real* prob_lo = geom.ProbLo();
        // const auto dx              = geom.CellSize();
        for (int k = 0; k <= khi; k++)
        {
            // const amrex::Real z_cc = prob_lo[2] + (k+myhalf)* dx[2];
            // set vertical velocity profile based on time, z_cc here
            wbar[k] = zero;
        }
 
        // Copy from host version to device version
        amrex::Gpu::copy(amrex::Gpu::hostToDevice, wbar.begin(), wbar.end(), d_wbar.begin());
    }

Here is the call graph for this function:

use_physical_z()

static bool ProblemBase::use_physical_z ( )

inlinestaticprotected

True when heights must come from the vertical levels or z_phys, not from the geometry.

{ return SolverChoice::mesh_type != MeshType::ConstantDz; }

Referenced by phys_zhi(), and phys_zlo().

Here is the caller graph for this function:

Member Data Documentation

base_parms

ProbParmDefaults ProblemBase::base_parms

protected

Referenced by init_base_parms().

m_zlevels_stag

amrex::Vector<amrex::Vector<amrex::Real> > ProblemBase::m_zlevels_stag

protected

Staggered vertical levels of AMR level lev (set by set_zlevels_stag)

Referenced by phys_zhi(), phys_zlo(), and set_zlevels_stag().

---

The documentation for this class was generated from the following file:

• Source/ERF_ProbCommon.H