Skip to content

Commit

Permalink
Fix flat top with wrfinput and allow for top extrapolation with metgr…
Browse files Browse the repository at this point in the history 
…id. (#1423)
  • Loading branch information
@AMLattanzi
AMLattanzi authored Feb 2, 2024
1 parent 8559185 commit c1ed778
Show file tree
Hide file tree
Showing 4 changed files with 152 additions and 38 deletions.
71 changes: 71 additions & 0 deletions Exec/DevTests/MetGrid/inputs_flowmas_metgrid
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,71 @@
# ------------------ INPUTS TO MAIN PROGRAM -------------------
max_step = 1000

amrex.fpe_trap_invalid = 1
amrex.fpe_trap_zero = 1
amrex.fpe_trap_overflow = 1

# PROBLEM SIZE & GEOMETRY
geometry.prob_extent = 447000 387000 20500
amr.n_cell = 149 129 40

geometry.is_periodic = 0 0 0

xlo.type = "Outflow"
xhi.type = "Outflow"
ylo.type = "Outflow"
yhi.type = "Outflow"
zlo.type = "NoSlipWall"
zhi.type = "SlipWall"

# TIME STEP CONTROL
erf.fixed_dt = 1.0 # fixed time step depending on grid resolution

erf.use_terrain = true
erf.terrain_smoothing = 2

# DIAGNOSTICS & VERBOSITY
erf.sum_interval = -1 # timesteps between computing mass
erf.v = 1 # verbosity in ERF.cpp
amr.v = 1 # verbosity in Amr.cpp

# REFINEMENT / REGRIDDING
amr.max_level = 0 # maximum level number allowed

# CHECKPOINT FILES
erf.check_file = chk # root name of checkpoint file
erf.check_int = 25 # number of timesteps between checkpoints

# PLOTFILES
erf.plot_file_1 = plt # prefix of plotfile name
erf.plot_int_1 = 25 # number of timesteps between plotfiles
erf.plot_vars_1 = density dens_hse rhoadv_0 x_velocity y_velocity z_velocity pressure temp theta z_phys mapfac pres_hse dens_hse pert_pres pert_dens rhoQ1 rhoQ2 rhoQ3 qv qc

# SOLVER CHOICE
erf.alpha_T = 1.0
erf.alpha_C = 1.0
erf.use_gravity = true

erf.molec_diff_type = "None"
erf.les_type = "Smagorinsky"
erf.Cs = 0.1

erf.moisture_model = "Kessler"

erf.use_terrain = true
erf.terrain_smoothing = 0
erf.terrain_z_levels = 0 53.25373067 116.3299934 190.8583428 278.6665766 381.7196783 502.1312033 642.3454199 \
805.2193582 992.889535 1207.388303 1450.761988 1724.488474 2029.43012 2365.936136 \
2733.669456 3131.532666 3557.60783 4009.253851 4483.532727 4978.420557 5494.911024 \
6034.232547 6597.033978 7183.139415 7792.547066 8425.074263 9080.62066 9761.028885 \
10467.85208 11200.42746 11958.47509 12753.28089 13605.36124 14516.05106 15467.59286 \
16453.44037 17477.73597 18503.25481 19533.31103 20500.0

# INITIALIZATION WITH ATM DATA
erf.metgrid_bdy_width = 5
erf.metgrid_bdy_set_width = 1
erf.init_type = "metgrid"
erf.nc_init_file_0 = "met_em.d01.2015-04-01_00_00_00.nc" "met_em.d01.2015-04-01_08_00_00.nc"

#There will be no OpenMP tiling
fabarray.mfiter_tile_size = 1024 1024 1024
17 changes: 9 additions & 8 deletions Exec/DevTests/MetGrid/inputs_fitch_flowmas → .../DevTests/MetGrid/inputs_flowmas_wrfinput
View file
Original file line number Diff line number Diff line change
Expand Up @@ -6,7 +6,7 @@ amrex.fpe_trap_zero = 1
amrex.fpe_trap_overflow = 1

# PROBLEM SIZE & GEOMETRY
geometry.prob_extent = 447000 387000 15000
geometry.prob_extent = 447000 387000 20500
amr.n_cell = 149 129 40

geometry.is_periodic = 0 0 0
Expand Down Expand Up @@ -38,8 +38,8 @@ erf.check_int = 25 # number of timesteps between checkpoints

# PLOTFILES
erf.plot_file_1 = plt # prefix of plotfile name
erf.plot_int_1 = 25 # number of timesteps between plotfiles
erf.plot_vars_1 = density dens_hse rhoadv_0 x_velocity y_velocity z_velocity pressure temp theta z_phys mapfac pres_hse pert_pres
erf.plot_int_1 = 1 # number of timesteps between plotfiles
erf.plot_vars_1 = density dens_hse rhoadv_0 x_velocity y_velocity z_velocity pressure temp theta z_phys mapfac pres_hse dens_hse pert_pres pert_dens rhoQ1 rhoQ2 rhoQ3 qv qc

# SOLVER CHOICE
erf.alpha_T = 1.0
Expand All @@ -50,13 +50,14 @@ erf.molec_diff_type = "None"
erf.les_type = "Smagorinsky"
erf.Cs = 0.1

erf.moisture_model = "NullMoist"
erf.moisture_model = "Kessler"

# INITIALIZATION WITH ATM DATA
erf.metgrid_bdy_width = 7
erf.metgrid_bdy_set_width = 1
erf.init_type = "metgrid"
erf.nc_init_file_0 = "met_em.d01.2015-04-01_00_00_00.nc" "met_em.d01.2015-04-01_08_00_00.nc"
erf.wrfbdy_width = 5
erf.wrfbdy_set_width = 1
erf.init_type = "real"
erf.nc_init_file_0 = "wrfinput_d01"
erf.nc_bdy_file = "wrfbdy_d01"

#There will be no OpenMP tiling
fabarray.mfiter_tile_size = 1024 1024 1024
60 changes: 34 additions & 26 deletions Source/Initialization/ERF_init_from_wrfinput.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -64,7 +64,8 @@ init_state_from_wrfinput (int lev, FArrayBox& state_fab,
const Vector<FArrayBox>& NC_zvel_fab,
const Vector<FArrayBox>& NC_rho_fab,
const Vector<FArrayBox>& NC_rhotheta_fab,
MoistureType moisture_type);
MoistureType moisture_type,
const int& n_qstate);

void
init_msfs_from_wrfinput (int lev, FArrayBox& msfu_fab,
Expand All @@ -73,7 +74,8 @@ init_msfs_from_wrfinput (int lev, FArrayBox& msfu_fab,
const Vector<FArrayBox>& NC_MSFV_fab,
const Vector<FArrayBox>& NC_MSFM_fab);
void
init_terrain_from_wrfinput (int lev, const Box& domain, FArrayBox& z_phys,
init_terrain_from_wrfinput (int lev, const Real& z_top,
const Box& domain, FArrayBox& z_phys,
const Vector<FArrayBox>& NC_PH_fab,
const Vector<FArrayBox>& NC_PHB_fab);

Expand Down Expand Up @@ -132,7 +134,7 @@ ERF::init_from_wrfinput (int lev)
}

auto& lev_new = vars_new[lev];

int n_qstate = micro.Get_Qstate_Size();
#ifdef _OPENMP
#pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
#endif
Expand All @@ -149,7 +151,8 @@ ERF::init_from_wrfinput (int lev)
init_state_from_wrfinput(lev, cons_fab, xvel_fab, yvel_fab, zvel_fab,
NC_QVAPOR_fab, NC_QCLOUD_fab, NC_QRAIN_fab,
NC_xvel_fab, NC_yvel_fab, NC_zvel_fab,
NC_rho_fab, NC_rhoth_fab, solverChoice.moisture_type);
NC_rho_fab, NC_rhoth_fab,
solverChoice.moisture_type, n_qstate);
} // mf

#ifdef _OPENMP
Expand All @@ -168,14 +171,14 @@ ERF::init_from_wrfinput (int lev)
} // mf

const Box& domain = geom[lev].Domain();

const Real& z_top = geom[lev].ProbHi(2);
if (solverChoice.use_terrain)
{
std::unique_ptr<MultiFab>& z_phys = z_phys_nd[lev];
for ( MFIter mfi(lev_new[Vars::cons], TilingIfNotGPU()); mfi.isValid(); ++mfi )
{
FArrayBox& z_phys_nd_fab = (*z_phys)[mfi];
init_terrain_from_wrfinput(lev, domain, z_phys_nd_fab, NC_PH_fab, NC_PHB_fab);
init_terrain_from_wrfinput(lev, z_top, domain, z_phys_nd_fab, NC_PH_fab, NC_PHB_fab);
} // mf

make_J (geom[lev],*z_phys_nd[lev],* detJ_cc[lev]);
Expand Down Expand Up @@ -265,7 +268,8 @@ init_state_from_wrfinput (int lev,
const Vector<FArrayBox>& NC_zvel_fab,
const Vector<FArrayBox>& NC_rho_fab,
const Vector<FArrayBox>& NC_rhotheta_fab,
MoistureType moisture_type)
MoistureType moisture_type,
const int& n_qstate)
{
int nboxes = NC_xvel_fab.size();
for (int idx = 0; idx < nboxes; idx++)
Expand Down Expand Up @@ -299,8 +303,10 @@ init_state_from_wrfinput (int lev,
state_fab.template copy<RunOn::Device>(NC_QCLOUD_fab[idx], 0, RhoQ2_comp, 1);
state_fab.template mult<RunOn::Device>(NC_rho_fab[idx] , 0, RhoQ2_comp, 1);

state_fab.template copy<RunOn::Device>(NC_QRAIN_fab[idx], 0, RhoQ3_comp, 1);
state_fab.template mult<RunOn::Device>(NC_rho_fab[idx] , 0, RhoQ3_comp, 1);
if (n_qstate >= 3) {
state_fab.template copy<RunOn::Device>(NC_QRAIN_fab[idx], 0, RhoQ3_comp, 1);
state_fab.template mult<RunOn::Device>(NC_rho_fab[idx] , 0, RhoQ3_comp, 1);
}
}
} // idx
}
Expand Down Expand Up @@ -389,7 +395,8 @@ init_base_state_from_wrfinput (int lev, const Box& valid_bx, const Real l_rdOcp,
* @param NC_PHB_fab Vector of FArrayBox objects storing WRF terrain coordinate data (PHB)
*/
void
init_terrain_from_wrfinput (int lev, const Box& domain, FArrayBox& z_phys,
init_terrain_from_wrfinput (int lev, const Real& z_top,
const Box& domain, FArrayBox& z_phys,
const Vector<FArrayBox>& NC_PH_fab,
const Vector<FArrayBox>& NC_PHB_fab)
{
Expand Down Expand Up @@ -440,27 +447,28 @@ init_terrain_from_wrfinput (int lev, const Box& domain, FArrayBox& z_phys,
nc_ph_arr (ii,jj-1,khi-1) + nc_ph_arr (ii-1,jj-1,khi-1) +
nc_phb_arr(ii,jj ,khi-1) + nc_phb_arr(ii-1,jj ,khi-1) +
nc_phb_arr(ii,jj-1,khi-1) + nc_phb_arr(ii-1,jj-1,khi-1) ) / CONST_GRAV;
z_arr(i, j, k) = 2.0 * z_khi - z_khim1;
} else {
z_arr(i, j, k) = 2.0 * z_top - z_khim1;
} else if (k == khi) {
z_arr(i, j, k) = z_top;
} else {
z_arr(i, j, k) = 0.25 * ( nc_ph_arr (ii,jj ,k) + nc_ph_arr (ii-1,jj ,k) +
nc_ph_arr (ii,jj-1,k) + nc_ph_arr (ii-1,jj-1,k) +
nc_phb_arr(ii,jj ,k) + nc_phb_arr(ii-1,jj ,k) +
nc_phb_arr(ii,jj-1,k) + nc_phb_arr(ii-1,jj-1,k) ) / CONST_GRAV;

//
// Fill values outside the domain on the fly (we will need these to make detJ in ghost cells)
//
if (i == domlo.x) z_arr(i-1,j,k) = z_arr(i,j,k);
if (i == domhi.x) z_arr(i+1,j,k) = z_arr(i,j,k);
if (j == domlo.y) z_arr(i,j-1,k) = z_arr(i,j,k);
if (j == domhi.y) z_arr(i,j+1,k) = z_arr(i,j,k);

if (i == domlo.x && j == domlo.y) z_arr(i-1,j-1,k) = z_arr(i,j,k);
if (i == domlo.x && j == domhi.y) z_arr(i-1,j+1,k) = z_arr(i,j,k);
if (i == domhi.x && j == domlo.y) z_arr(i+1,j-1,k) = z_arr(i,j,k);
if (i == domhi.x && j == domhi.y) z_arr(i+1,j+1,k) = z_arr(i,j,k);

} // k

//
// Fill values outside the domain on the fly (we will need these to make detJ in ghost cells)
//
if (i == domlo.x) z_arr(i-1,j,k) = z_arr(i,j,k);
if (i == domhi.x) z_arr(i+1,j,k) = z_arr(i,j,k);
if (j == domlo.y) z_arr(i,j-1,k) = z_arr(i,j,k);
if (j == domhi.y) z_arr(i,j+1,k) = z_arr(i,j,k);

if (i == domlo.x && j == domlo.y) z_arr(i-1,j-1,k) = z_arr(i,j,k);
if (i == domlo.x && j == domhi.y) z_arr(i-1,j+1,k) = z_arr(i,j,k);
if (i == domhi.x && j == domlo.y) z_arr(i+1,j-1,k) = z_arr(i,j,k);
if (i == domhi.x && j == domhi.y) z_arr(i+1,j+1,k) = z_arr(i,j,k);
});
} // idx
}
Expand Down
42 changes: 38 additions & 4 deletions Source/Initialization/Metgrid_utils.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -189,8 +189,8 @@ interpolate_column_metgrid (const int& i,
amrex::Real z0, z1;
int klow = -1;
int khi0 = -1;
amrex::Real dzlow = 10000.0;
amrex::Real dzhi0 = -10000.0;
amrex::Real dzlow = 1.0e12;
amrex::Real dzhi0 = -1.0e12;
for (int kk = 0; kk < kmax_orig; kk++) {
amrex::Real orig_z_stag = 0.0;
if (stag == 'M') {
Expand Down Expand Up @@ -218,6 +218,7 @@ interpolate_column_metgrid (const int& i,
orig_z_stag = 0.5*(orig_z(i,j,kk)+orig_z(i,j-1,kk));
}
}

amrex::Real dz = z - orig_z_stag;
if ((dz < 0.0) && (dz > dzhi0)) {
dzhi0 = dz;
Expand All @@ -231,10 +232,10 @@ interpolate_column_metgrid (const int& i,
}
} // kk

// extrapolate below the bottom surface
if (klow == -1) {
// extrapolate
int khi1 = -1;
amrex::Real dzhi1 = -10000.0;
amrex::Real dzhi1 = -1.0e12;
for (int kk = 0; kk < kmax_orig; kk++) {
amrex::Real orig_z_stag = 0.0;
if (stag == 'M') {
Expand Down Expand Up @@ -272,6 +273,39 @@ interpolate_column_metgrid (const int& i,
amrex::Real y0 = orig_data(i,j,khi0,src_comp);
amrex::Real y1 = orig_data(i,j,khi1,src_comp);
return ( y0-(y1-y0)/(z1-z0)*(z0-z) );

// Extrapolate above the top surface
} else if (khi0 == -1) {
khi0 = klow - 1;
int khi1 = klow;
if (stag == 'M') {
z0 = orig_z(i,j,khi0);
}
else if (stag == 'X') {
if (i == 0) {
z0 = orig_z(i,j,khi0);
}
else if (i == imax_orig) {
z0 = orig_z(imax_orig-1,j,khi0);
}
else {
z0 = 0.5*(orig_z(i,j,khi0)+orig_z(i-1,j,khi0));
}
}
else if (stag == 'Y') {
if (j == 0) {
z0 = orig_z(i,j,khi0);
}
else if (j == jmax_orig) {
z0 = orig_z(i,jmax_orig-1,khi0);
}
else {
z0 = 0.5*(orig_z(i,j,khi0)+orig_z(i,j-1,khi0));
}
}
amrex::Real y0 = orig_data(i,j,khi0,src_comp);
amrex::Real y1 = orig_data(i,j,khi1,src_comp);
return ( y0+(y1-y0)/(z1-z0)*(z-z0) );
} else {
// interpolate
amrex::Real y0 = orig_data(i,j,klow,src_comp);
Expand Down

0 comments on commit c1ed778

Please sign in to comment.