Skip to content

Commit

Permalink
Merge pull request #3 from AMLattanzi/refactor_prob_init
Browse files Browse the repository at this point in the history 
Suppress warnings
  • Loading branch information
@ewquon
ewquon authored Sep 8, 2023
2 parents 0c07671 + 44ce31f commit d3c86a1
Show file tree
Hide file tree
Showing 2 changed files with 49 additions and 148 deletions.
114 changes: 5 additions & 109 deletions Exec/RegTests/DensityCurrent/prob.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -40,7 +40,7 @@ Problem::init_custom_pert(
Array4<Real > const& z_vel,
Array4<Real > const& r_hse,
Array4<Real > const& p_hse,
Array4<Real const> const& z_nd,
Array4<Real const> const& /*z_nd*/,
Array4<Real const> const& z_cc,
#if defined(ERF_USE_MOISTURE)
Array4<Real > const&,
Expand All @@ -60,85 +60,24 @@ Problem::init_custom_pert(

AMREX_ALWAYS_ASSERT(bx.length()[2] == khi+1);

const Real rho_sfc = p_0 / (R_d*parms.T_0);
const Real thetabar = parms.T_0;
const Real dz = geomdata.CellSize()[2];
const Real prob_lo_z = geomdata.ProbLo()[2];

const Real l_x_r = parms.x_r;
//const Real l_x_r = parms.x_r * mf_u(0,0,0); //used to validate constant msf
const Real l_z_r = parms.z_r;
const Real l_x_c = parms.x_c;
const Real l_z_c = parms.z_c;
const Real l_Tpt = parms.T_pert;

#if 0
// These are at cell centers (unstaggered)
Vector<Real> h_r(khi+2);
Vector<Real> h_p(khi+2);

amrex::Gpu::DeviceVector<Real> d_r(khi+2);
amrex::Gpu::DeviceVector<Real> d_p(khi+2);
#endif

const Real rdOcp = sc.rdOcp;

if (z_cc) {

#if 0
// Create a flat box with same horizontal extent but only one cell in vertical
Box b2d = surroundingNodes(bx); // Copy constructor
b2d.setRange(2,0);

ParallelFor(b2d, [=] AMREX_GPU_DEVICE (int i, int j, int)
{
Array1D<Real,0,255> r;;
Array1D<Real,0,255> p;;

init_isentropic_hse_terrain(i,j,rho_sfc,thetabar,&(r(0)),&(p(0)),z_cc,khi);

for (int k = 0; k <= khi; k++) {
r_hse(i,j,k) = r(k);
p_hse(i,j,k) = p(k);
}
r_hse(i,j, -1) = r_hse(i,j,0);
r_hse(i,j,khi+1) = r_hse(i,j,khi);
});
#endif

amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
{
// Geometry (note we must include these here to get the data on device)
const auto prob_lo = geomdata.ProbLo();
const auto dx = geomdata.CellSize();
const auto prob_lo = geomdata.ProbLo();
const auto dx = geomdata.CellSize();

const Real x = prob_lo[0] + (i + 0.5) * dx[0];
const Real z = z_cc(i,j,k);

#if 0
// Temperature that satisfies the EOS given the hydrostatically balanced (r,p)
const Real Tbar_hse = p_hse(i,j,k) / (R_d * r_hse(i,j,k));

Real L = std::sqrt(
std::pow((x - l_x_c)/l_x_r, 2) +
std::pow((z - l_z_c)/l_z_r, 2)
);
Real dT;
if (L > 1.0) {
dT = 0.0;
}
else {
dT = l_Tpt * (std::cos(PI*L) + 1.0)/2.0;
}

// Note: dT is a perturbation in temperature, theta_perturbed is theta PLUS perturbation in theta
Real theta_perturbed = (Tbar_hse+dT)*std::pow(p_0/p_hse(i,j,k), rdOcp);

// This version perturbs rho but not p
state(i, j, k, RhoTheta_comp) = getRhoThetagivenP(p_hse(i,j,k));
state(i, j, k, Rho_comp) = state(i, j, k, RhoTheta_comp) / theta_perturbed;
#endif

Real L = std::sqrt(
std::pow((x - l_x_c)/l_x_r, 2) +
std::pow((z - l_z_c)/l_z_r, 2));
Expand All @@ -165,58 +104,15 @@ Problem::init_custom_pert(
#endif
});
} else {

#if 0
init_isentropic_hse(rho_sfc,thetabar,h_r.data(),h_p.data(),dz,prob_lo_z,khi);

amrex::Gpu::copyAsync(amrex::Gpu::hostToDevice, h_r.begin(), h_r.end(), d_r.begin());
amrex::Gpu::copyAsync(amrex::Gpu::hostToDevice, h_p.begin(), h_p.end(), d_p.begin());

Real* r = d_r.data();
Real* p = d_p.data();
#endif

amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
{
// Geometry (note we must include these here to get the data on device)
const auto prob_lo = geomdata.ProbLo();
const auto dx = geomdata.CellSize();
const auto prob_lo = geomdata.ProbLo();
const auto dx = geomdata.CellSize();

const Real x = prob_lo[0] + (i + 0.5) * dx[0];
const Real z = prob_lo[2] + (k + 0.5) * dx[2];

#if 0
// Temperature that satisfies the EOS given the hydrostatically balanced (r,p)
const Real Tbar_hse = p[k] / (R_d * r[k]);

Real L = std::sqrt(
std::pow((x - l_x_c)/l_x_r, 2) +
std::pow((z - l_z_c)/l_z_r, 2)
);
Real dT;
if (L > 1.0) {
dT = 0.0;
}
else {
dT = l_Tpt * (std::cos(PI*L) + 1.0)/2.0;
}

// Note: dT is a perturbation in temperature, theta_perturbed is theta PLUS perturbation in theta
Real theta_perturbed = (Tbar_hse+dT)*std::pow(p_0/p[k], rdOcp);

// This version perturbs rho but not p
state(i, j, k, RhoTheta_comp) = getRhoThetagivenP(p[k]);
state(i, j, k, Rho_comp) = state(i, j, k, RhoTheta_comp) / theta_perturbed;

if ((i==0) && (j==0))
{
amrex::Print() << "init_custom_pert("<<i<<","<<j<<","<<k<<")"
<< " r["<<k<<"]=" << r[k] << " r_hse=" << r_hse(i,j,k) << " diff: " << r[k] - r_hse(i,j,k)
<< " p["<<k<<"]=" << p[k] << " p_hse=" << p_hse(i,j,k) << " diff: " << p[k] - p_hse(i,j,k)
<< std::endl;
}
#endif

Real L = std::sqrt(
std::pow((x - l_x_c)/l_x_r, 2) +
std::pow((z - l_z_c)/l_z_r, 2));
Expand Down
83 changes: 44 additions & 39 deletions Source/prob_common.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -18,6 +18,12 @@ struct ProbParmDefaults {
class ProblemBase
{
public:

/**
* Virtual destructor to avoid data leakage with derived class
*/
virtual ~ProblemBase() = default;

/**
* Function to initialize the hydrostatic reference density
*
Expand All @@ -27,10 +33,10 @@ public:
* @param[in] geom container for geometric information
*/
virtual void
erf_init_dens_hse (amrex::MultiFab& rho_hse,
std::unique_ptr<amrex::MultiFab>& z_phys_nd,
std::unique_ptr<amrex::MultiFab>& z_phys_cc,
amrex::Geometry const& geom)
erf_init_dens_hse (amrex::MultiFab& /*rho_hse*/,
std::unique_ptr<amrex::MultiFab>& /*z_phys_nd*/,
std::unique_ptr<amrex::MultiFab>& /*z_phys_cc*/,
amrex::Geometry const& /*geom*/)
{
amrex::Print() << "Hydrostatically balanced density was NOT set"
<< " -- an appropriate init_type should probably have been specified"
Expand Down Expand Up @@ -62,34 +68,32 @@ public:
* @param[in] mf_v map factor on y-faces
* @param[in] sc SolverChoice structure that carries parameters
*/

virtual void
init_custom_pert (
const amrex::Box& bx,
const amrex::Box& xbx,
const amrex::Box& ybx,
const amrex::Box& zbx,
amrex::Array4<amrex::Real > const& state,
amrex::Array4<amrex::Real > const& x_vel,
amrex::Array4<amrex::Real > const& y_vel,
amrex::Array4<amrex::Real > const& z_vel,
amrex::Array4<amrex::Real > const& r_hse,
amrex::Array4<amrex::Real > const& p_hse,
amrex::Array4<amrex::Real const> const& z_nd,
amrex::Array4<amrex::Real const> const& z_cc,
init_custom_pert (const amrex::Box& /*bx*/,
const amrex::Box& /*xbx*/,
const amrex::Box& /*ybx*/,
const amrex::Box& /*zbx*/,
amrex::Array4<amrex::Real > const& /*state*/,
amrex::Array4<amrex::Real > const& /*x_vel*/,
amrex::Array4<amrex::Real > const& /*y_vel*/,
amrex::Array4<amrex::Real > const& /*z_vel*/,
amrex::Array4<amrex::Real > const& /*r_hse*/,
amrex::Array4<amrex::Real > const& /*p_hse*/,
amrex::Array4<amrex::Real const> const& /*z_nd*/,
amrex::Array4<amrex::Real const> const& /*z_cc*/,
#if defined(ERF_USE_MOISTURE)
amrex::Array4<amrex::Real > const& qv,
amrex::Array4<amrex::Real > const& qc,
amrex::Array4<amrex::Real > const& qi,
amrex::Array4<amrex::Real > const& /*qv*/,
amrex::Array4<amrex::Real > const& /*qc*/,
amrex::Array4<amrex::Real > const& /*qi*/,
#elif defined(ERF_USE_WARM_NO_PRECIP)
amrex::Array4<amrex::Real > const& qv,
amrex::Array4<amrex::Real > const& qc,
amrex::Array4<amrex::Real > const& /*qv*/,
amrex::Array4<amrex::Real > const& /*qc*/,
#endif
amrex::GeometryData const& geomdata,
amrex::Array4<amrex::Real const> const& mf_m,
amrex::Array4<amrex::Real const> const& mf_u,
amrex::Array4<amrex::Real const> const& mf_v,
const SolverChoice& sc
amrex::GeometryData const& /*geomdata*/,
amrex::Array4<amrex::Real const> const& /*mf_m*/,
amrex::Array4<amrex::Real const> const& /*mf_u*/,
amrex::Array4<amrex::Real const> const& /*mf_v*/,
const SolverChoice& /*sc*/
)
{
amrex::Print() << "No perturbation to background fields supplied for "
Expand All @@ -104,15 +108,15 @@ public:
* @param[in] time current time
*/
virtual void
init_custom_terrain (const amrex::Geometry& geom,
amrex::MultiFab& z_phys_nd,
const amrex::Real& time)
init_custom_terrain (const amrex::Geometry& /*geom*/,
amrex::MultiFab& /*z_phys_nd*/,
const amrex::Real& /*time*/)
{
amrex::Error("Should never call init_custom_terrain for "+name()+" problem");
}

#ifdef ERF_USE_TERRAIN_VELOCITY
virtual amrex::Real compute_terrain_velocity(const amrex::Real time)
virtual amrex::Real compute_terrain_velocity(const amrex::Real /*time*/)
{
amrex::Error("Should never call compute_terrain_velocity for "+name()+" problem");
}
Expand All @@ -129,12 +133,12 @@ public:
* @param[in] geom container for geometric information
*/
virtual void
erf_init_rayleigh (amrex::Vector<amrex::Real>& tau,
amrex::Vector<amrex::Real>& ubar,
amrex::Vector<amrex::Real>& vbar,
amrex::Vector<amrex::Real>& wbar,
amrex::Vector<amrex::Real>& thetabar,
amrex::Geometry const& geom)
erf_init_rayleigh (amrex::Vector<amrex::Real>& /*tau*/,
amrex::Vector<amrex::Real>& /*ubar*/,
amrex::Vector<amrex::Real>& /*vbar*/,
amrex::Vector<amrex::Real>& /*wbar*/,
amrex::Vector<amrex::Real>& /*thetabar*/,
amrex::Geometry const& /*geom*/)
{
amrex::Error("Should never call erf_init_rayleigh for "+name()+" problem");
}
Expand All @@ -143,7 +147,8 @@ public:
* Function to set uniform background density and pressure fields
*/
void
init_uniform (const amrex::Box& bx, amrex::Array4<amrex::Real> const& state);
init_uniform (const amrex::Box& bx,
amrex::Array4<amrex::Real> const& state);

protected:
// Struct to store problem parameters
Expand Down

0 comments on commit d3c86a1

Please sign in to comment.