Skip to content

Commit

Permalink
moving towards enabling EB (#1559)
Browse files Browse the repository at this point in the history
  • Loading branch information
@asalmgren
asalmgren authored Apr 5, 2024
1 parent 8a7fead commit 4a45a7d
Show file tree
Hide file tree
Showing 29 changed files with 771 additions and 686 deletions.
10 changes: 10 additions & 0 deletions Exec/DevTests/EB_Test/amrvis.defaults
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,10 @@
palette /home/almgren/bin/Palette
initialderived volfrac
initialscale 8
expansion 1
numberformat %8.6e
maxpixmapsize 1000000
reservesystemcolors 35
showboxes TRUE
windowheight 650
windowwidth 1100
44 changes: 28 additions & 16 deletions Exec/DevTests/EB_Test/inputs
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2,23 +2,36 @@
max_step = 100
max_step = 0

amr.max_grid_size = 256 256 256

eb2.geometry = terrain

amrex.fpe_trap_invalid = 1

fabarray.mfiter_tile_size = 1024 1024 1024

# PROBLEM SIZE & GEOMETRY
geometry.prob_extent = 16. 16. 16.
amr.n_cell = 16 16 16
geometry.prob_lo = 0. 0. 0.
geometry.prob_hi = 10. 1. 4.

amr.n_cell = 256 8 64 # dx=dy=dz=100 m, Straka et al 1993 / Xue et al 2000

geometry.is_periodic = 0 1 0

geometry.is_periodic = 1 1 0
xlo.type = "Inflow"
xhi.type = "Outflow"

xlo.velocity = 1. 0. 0.
xlo.density = 1.16
xlo.theta = 300.
xlo.scalar = 0.

zlo.type = "SlipWall"
zhi.type = "SlipWall"

# TIME STEP CONTROL
erf.cfl = 0.9 # cfl number for hyperbolic system
erf.no_substepping = 1
erf.fixed_dt = 1.e-5

# DIAGNOSTICS & VERBOSITY
erf.sum_interval = 1 # timesteps between computing mass
Expand All @@ -38,23 +51,22 @@ erf.plot_int_1 = 20 # number of timesteps between plotfiles
erf.plot_vars_1 = density rhoadv_0 x_velocity y_velocity z_velocity pressure scalar volfrac

# SOLVER CHOICE
erf.alpha_T = 0.0
erf.alpha_C = 0.0
erf.use_gravity = false
erf.use_gravity = true
erf.use_coriolis = false
erf.les_type = "None"

erf.les_type = "None"
erf.molec_diff_type = "None"
erf.dynamicViscosity = 0.0

erf.alpha_T = 0.0
erf.alpha_C = 0.0

erf.init_type = "uniform"

# PROBLEM PARAMETERS
prob.rho_0 = 1.0
prob.T_0 = 1.0
prob.A_0 = 1.0
prob.u_0 = 10.0
prob.v_0 = 5.0
prob.rad_0 = 0.125
prob.uRef = 0.0

prob.prob_type = 11
prob.u_0 = 1.0

erf.use_rayleigh_damping = true
prob.dampcoef = 2000. # ==> time scale = 0.0005 s (for testing)
prob.zdamp = 1.0
166 changes: 72 additions & 94 deletions Exec/DevTests/EB_Test/prob.cpp
View file
Original file line number Diff line number Diff line change
@@ -1,4 +1,5 @@
#include "prob.H"
#include "TerrainMetrics.H"

using namespace amrex;

Expand Down Expand Up @@ -66,98 +67,60 @@ Problem::init_custom_pert(
Array4<Real > const& z_vel_pert,
Array4<Real > const& /*r_hse*/,
Array4<Real > const& /*p_hse*/,
Array4<Real const> const& /*z_nd*/,
Array4<Real const> const& /*z_cc*/,
Array4<Real const> const& z_nd,
Array4<Real const> const& z_cc,
GeometryData const& geomdata,
Array4<Real const> const& /*mf_m*/,
Array4<Real const> const& /*mf_u*/,
Array4<Real const> const& /*mf_v*/,
const SolverChoice& sc)
{
const int khi = geomdata.Domain().bigEnd()[2];

const bool use_moisture = (sc.moisture_type != MoistureType::None);

amrex::ParallelFor(bx, [=, parms=parms] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
{
// Geometry
const Real* prob_lo = geomdata.ProbLo();
const Real* prob_hi = geomdata.ProbHi();
const Real* dx = geomdata.CellSize();
const Real x = prob_lo[0] + (i + 0.5) * dx[0];
const Real y = prob_lo[1] + (j + 0.5) * dx[1];
const Real z = prob_lo[2] + (k + 0.5) * dx[2];

// Define a point (xc,yc,zc) at the center of the domain
const Real xc = parms.xc_frac * (prob_lo[0] + prob_hi[0]);
const Real yc = parms.yc_frac * (prob_lo[1] + prob_hi[1]);
const Real zc = parms.zc_frac * (prob_lo[2] + prob_hi[2]);

// Define ellipse parameters
const Real r0 = parms.rad_0 * (prob_hi[0] - prob_lo[0]);

const Real r3d = std::sqrt((x-xc)*(x-xc) + (y-yc)*(y-yc) + (z-zc)*(z-zc));
const Real r2d_xy = std::sqrt((x-xc)*(x-xc) + (y-yc)*(y-yc));
const Real r2d_xz = std::sqrt((x-xc)*(x-xc) + (z-zc)*(z-zc));

if (parms.prob_type == 10)
{
// Set scalar = A_0*exp(-10r^2), where r is distance from center of domain,
// + B_0*sin(x)
state_pert(i, j, k, RhoScalar_comp) = parms.A_0 * exp(-10.*r3d*r3d) + parms.B_0*sin(x);

} else if (parms.prob_type == 11) {
state_pert(i, j, k, RhoScalar_comp) = parms.A_0 * 0.25 * (1.0 + std::cos(PI * std::min(r2d_xy, r0) / r0));
} else if (parms.prob_type == 12) {
state_pert(i, j, k, RhoScalar_comp) = parms.A_0 * 0.25 * (1.0 + std::cos(PI * std::min(r2d_xz, r0) / r0));
} else if (parms.prob_type == 13) {
const Real r0_z = parms.rad_0 * (prob_hi[2] - prob_lo[2]);
//ellipse for mapfac shear validation
const Real r2d_xz_ell = std::sqrt((x-xc)*(x-xc)/(r0*r0) + (z-zc)*(z-zc)/(r0_z*r0_z));
state_pert(i, j, k, RhoScalar_comp) = parms.A_0 * 0.25 * (1.0 + std::cos(PI * std::min(r2d_xz_ell, r0_z) / r0_z));
} else if (parms.prob_type == 14) {
state_pert(i, j, k, RhoScalar_comp) = std::cos(PI*x);
} else {
// Set scalar = A_0 in a ball of radius r0 and 0 elsewhere
if (r3d < r0) {
state_pert(i, j, k, RhoScalar_comp) = parms.A_0;
} else {
state_pert(i, j, k, RhoScalar_comp) = 0.0;
}
}
AMREX_ALWAYS_ASSERT(bx.length()[2] == khi+1);

state_pert(i, j, k, RhoScalar_comp) *= parms.rho_0;
ParallelFor(bx, [=, parms=parms] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
{
// Set scalar = 0 everywhere
state_pert(i, j, k, RhoScalar_comp) = 0.0;

if (use_moisture) {
state_pert(i, j, k, RhoQ1_comp) = 0.0;
state_pert(i, j, k, RhoQ2_comp) = 0.0;
}
});
if (use_moisture) {
state_pert(i, j, k, RhoQ1_comp) = 0.0;
state_pert(i, j, k, RhoQ2_comp) = 0.0;
}
});

// Set the x-velocity
amrex::ParallelFor(xbx, [=, parms=parms] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
{
x_vel_pert(i, j, k) = parms.u_0;
// Set the x-velocity
ParallelFor(xbx, [=, parms=parms] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
{
x_vel_pert(i, j, k) = parms.u_0;
});

const Real* prob_lo = geomdata.ProbLo();
const Real* dx = geomdata.CellSize();
const Real z = prob_lo[2] + (k + 0.5) * dx[2];
// Set the y-velocity
ParallelFor(ybx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
{
y_vel_pert(i, j, k) = 0.0;
});

// Set the x-velocity
x_vel_pert(i, j, k) = parms.u_0 + parms.uRef *
std::log((z + parms.z0)/parms.z0)/
std::log((parms.zRef +parms.z0)/parms.z0);
});
const auto dx = geomdata.CellSize();
GpuArray<Real, AMREX_SPACEDIM> dxInv;
dxInv[0] = 1. / dx[0];
dxInv[1] = 1. / dx[1];
dxInv[2] = 1. / dx[2];

// Set the y-velocity
amrex::ParallelFor(ybx, [=, parms=parms] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
{
y_vel_pert(i, j, k) = parms.v_0;
});
// Set the z-velocity from impenetrable condition
ParallelFor(zbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
{
#ifdef ERF_USE_EB
z_vel_pert(i, j, k) = 0.0;
#else
z_vel_pert(i, j, k) = WFromOmega(i, j, k, 0.0, x_vel_pert, y_vel_pert, z_nd, dxInv);
#endif
});

// Set the z-velocity
amrex::ParallelFor(zbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
{
z_vel_pert(i, j, k) = 0.0;
});
amrex::Gpu::streamSynchronize();
}

/**
Expand All @@ -167,32 +130,47 @@ Problem::init_custom_pert(
*
*/
void
Problem::init_custom_terrain (const amrex::Geometry& /*geom*/,
Problem::init_custom_terrain (const amrex::Geometry& geom,
amrex::FArrayBox& z_phys_nd,
const amrex::Real& /*time*/)
{
// Note that this only sets the terrain value at the ground IF k=0 is in the box
amrex::Print() << "Initializing flat terrain at z=0" << std::endl;

// Bottom of domain
int k0 = 0;

// Domain valid box (z_nd is nodal)
const amrex::Box& domain = geom.Domain();
int domlo_x = domain.smallEnd(0); int domhi_x = domain.bigEnd(0) + 1;

const Real* prob_lo = geom.ProbLo();
const Real* prob_hi = geom.ProbHi();

const Real* dx = geom.CellSize();

// User function parameters
Real a = 0.5;
Real num = 8 * a * a * a;
Real xcen = 0.5 * (prob_lo[0] + prob_hi[0]);

// Grown box with no z range
amrex::Box bx = z_phys_nd.box();
bx.setRange(2,0);

amrex::Array4<amrex::Real> const& z_arr = z_phys_nd.array();

ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int) {
z_arr(i,j,k0) = 2.5;
ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int)
{
// Clip indices for ghost-cells
int ii = amrex::min(amrex::max(i,domlo_x),domhi_x);

// Location of nodes
Real x = (ii * dx[0] - xcen);

// WoA Hill in x-direction
Real height = num / (x*x + 4 * a * a);

// Populate terrain height
z_arr(i,j,k0) = height;

// z_arr(i,j,k0) = 2.5;
});
}

#if 0
AMREX_GPU_DEVICE
Real
dhdt(int /*i*/, int /*j*/,
const GpuArray<Real,AMREX_SPACEDIM> /*dx*/,
const Real /*time_mt*/, const Real /*delta_t*/)
{
return 0.;
}
#endif
31 changes: 15 additions & 16 deletions Source/Advection/Advection.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -19,18 +19,14 @@ void AdvectionSrcForRho (const amrex::Box& bx,
const amrex::Array4< amrex::Real>& avg_xmom, // These are being defined
const amrex::Array4< amrex::Real>& avg_ymom, // from the rho fluxes
const amrex::Array4< amrex::Real>& avg_zmom,
const amrex::Array4<const amrex::Real>& z_nd,
const amrex::Array4<const amrex::Real>& ax_arr,
const amrex::Array4<const amrex::Real>& ay_arr,
const amrex::Array4<const amrex::Real>& az_arr,
const amrex::Array4<const amrex::Real>& detJ,
const amrex::GpuArray<amrex::Real, AMREX_SPACEDIM>& cellSizeInv,
const amrex::Array4<const amrex::Real>& mf_m,
const amrex::Array4<const amrex::Real>& mf_u,
const amrex::Array4<const amrex::Real>& mf_v,
#ifdef ERF_USE_EB
const amrex::Array4<const amrex::Real>& ax_arr,
const amrex::Array4<const amrex::Real>& ay_arr,
const amrex::Array4<const amrex::Real>& az_arr,
#endif
const bool use_terrain,
const amrex::GpuArray<const amrex::Array4<amrex::Real>, AMREX_SPACEDIM>& flx_arr);

/** Compute advection tendency for all scalars other than density and potential temperature */
Expand All @@ -46,7 +42,6 @@ void AdvectionSrcForScalars (const amrex::Box& bx,
const amrex::Array4<const amrex::Real>& mf_m,
const AdvType horiz_adv_type, const AdvType vert_adv_type,
const amrex::Real horiz_upw_frac, const amrex::Real vert_upw_frac,
const bool use_terrain,
const amrex::GpuArray<const amrex::Array4<amrex::Real>, AMREX_SPACEDIM>& flx_arr,
const amrex::Box& domain,
const amrex::BCRec* bc_ptr_h);
Expand All @@ -60,7 +55,10 @@ void AdvectionSrcForMom (const amrex::Box& bxx, const amrex::Box& bxy, const amr
const amrex::Array4<const amrex::Real>& w ,
const amrex::Array4<const amrex::Real>& rho_u , const amrex::Array4<const amrex::Real>& rho_v,
const amrex::Array4<const amrex::Real>& Omega ,
const amrex::Array4<const amrex::Real>& z_nd , const amrex::Array4<const amrex::Real>& detJ,
const amrex::Array4<const amrex::Real>& ax,
const amrex::Array4<const amrex::Real>& ay,
const amrex::Array4<const amrex::Real>& az,
const amrex::Array4<const amrex::Real>& detJ,
const amrex::GpuArray<amrex::Real, AMREX_SPACEDIM>& cellSizeInv,
const amrex::Array4<const amrex::Real>& mf_m,
const amrex::Array4<const amrex::Real>& mf_u,
Expand Down Expand Up @@ -90,10 +88,10 @@ AdvectionSrcForOpenBC_Tangent_Xmom (const amrex::Box& bxx,
const amrex::Array4<const amrex::Real>& rho_u,
const amrex::Array4<const amrex::Real>& rho_v,
const amrex::Array4<const amrex::Real>& Omega,
const amrex::Array4<const amrex::Real>& z_nd,
const amrex::Array4<const amrex::Real>& ax,
const amrex::Array4<const amrex::Real>& az,
const amrex::Array4<const amrex::Real>& detJ,
const amrex::GpuArray<amrex::Real, AMREX_SPACEDIM>& cellSizeInv,
const bool use_terrain,
const bool do_lo=false);

/** Compute advection tendencies for v momentum tangent to open bc*/
Expand All @@ -105,10 +103,10 @@ AdvectionSrcForOpenBC_Tangent_Ymom (const amrex::Box& bxy,
const amrex::Array4<const amrex::Real>& rho_u,
const amrex::Array4<const amrex::Real>& rho_v,
const amrex::Array4<const amrex::Real>& Omega,
const amrex::Array4<const amrex::Real>& z_nd,
const amrex::Array4<const amrex::Real>& ay,
const amrex::Array4<const amrex::Real>& az,
const amrex::Array4<const amrex::Real>& detJ,
const amrex::GpuArray<amrex::Real, AMREX_SPACEDIM>& cellSizeInv,
const bool use_terrain,
const bool do_lo=false);

/** Compute advection tendencies for w momentum tangent to open bc*/
Expand All @@ -120,10 +118,12 @@ AdvectionSrcForOpenBC_Tangent_Zmom (const amrex::Box& bxz,
const amrex::Array4<const amrex::Real>& rho_u,
const amrex::Array4<const amrex::Real>& rho_v,
const amrex::Array4<const amrex::Real>& Omega,
const amrex::Array4<const amrex::Real>& z_nd,
const amrex::Array4<const amrex::Real>& ax,
const amrex::Array4<const amrex::Real>& ay,
const amrex::Array4<const amrex::Real>& az,
const amrex::Array4<const amrex::Real>& detJ,
const amrex::GpuArray<amrex::Real, AMREX_SPACEDIM>& cellSizeInv,
const bool use_terrain, const int domhi_z,
const int domhi_z,
const bool do_lo=false);

/** Compute advection tendencies for cons tangential to BC (2nd order)*/
Expand All @@ -139,7 +139,6 @@ AdvectionSrcForOpenBC_Tangent_Cons (const amrex::Box& bx,
const amrex::Array4<const amrex::Real>& avg_zmom,
const amrex::Array4<const amrex::Real>& detJ,
const amrex::GpuArray<amrex::Real, AMREX_SPACEDIM>& cellSizeInv,
const bool use_terrain,
const bool do_lo=false);

/** Compute advection tendencies in normal dir for vars tangent to open bc */
Expand Down
Loading

0 comments on commit 4a45a7d

Please sign in to comment.