Skip to content

Commit

Permalink
WIP more debugging for WPS met_em*.nc initialization and forcing.
Browse files Browse the repository at this point in the history
  • Loading branch information
@wiersema1
wiersema1 committed Oct 30, 2023
1 parent bd57c15 commit df3f0db
Show file tree
Hide file tree
Showing 11 changed files with 64 additions and 82 deletions.
2 changes: 2 additions & 0 deletions Exec/DevTests/MetGrid/inputs
View file
Original file line number Diff line number Diff line change
Expand Up @@ -55,6 +55,8 @@ erf.Cs = 0.1
#erf.terrain_z_levels = 0 130 354 583 816 1054 1549 2068 2615 3193 3803 4450 5142 5892 6709 7603 8591 9702 10967 12442 14230 16610 18711 20752 22133 23960 26579 28493 31236 33699 36068 40000

# 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.2022-06-18_00:00:00.nc" "met_em.d01.2022-06-18_06:00:00.nc" "met_em.d01.2022-06-18_12:00:00.nc" "met_em.d01.2022-06-18_18:00:00.nc"

Expand Down
2 changes: 1 addition & 1 deletion Exec/DevTests/MetGrid/prob.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -12,7 +12,7 @@ struct ProbParm : ProbParmDefaults {
class Problem : public ProblemBase
{
public:
Problem();
Problem(const amrex::Real* problo, const amrex::Real* probhi);

protected:
std::string name() override { return "Metgrid"; }
Expand Down
2 changes: 1 addition & 1 deletion Exec/DevTests/MetGrid/prob.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -8,5 +8,5 @@ amrex_probinit(const amrex_real* problo, const amrex_real* probhi)
return std::make_unique<Problem>(problo, probhi);
}

Problem::Problem()
Problem::Problem(const amrex::Real* problo, const amrex::Real* probhi)
{}
46 changes: 27 additions & 19 deletions Source/BoundaryConditions/BoundaryConditions_metgrid.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -12,7 +12,11 @@ using namespace amrex;
*/

void
ERF::fill_from_metgrid (const Vector<MultiFab*>& mfs, const Real time)
ERF::fill_from_metgrid (const Vector<MultiFab*>& mfs,
const Real time,
bool cons_only,
int icomp_cons,
int ncomp_cons)
{
int lev = 0;

Expand Down Expand Up @@ -54,10 +58,13 @@ ERF::fill_from_metgrid (const Vector<MultiFab*>& mfs, const Real time)
ind_map.push_back( {0} ); // zvel

// Nvars to loop over
Vector<int> comp_var = {NVAR, 1, 1, 1};
Vector<int> comp_var = {ncomp_cons, 1, 1, 1};

// End of vars loop
int var_idx_end = (cons_only) ? Vars::cons + 1 : Vars::NumTypes;

// Loop over all variable types
for (int var_idx = Vars::cons; var_idx < Vars::NumTypes; ++var_idx)
for (int var_idx = Vars::cons; var_idx < var_idx_end; ++var_idx)
{
MultiFab& mf = *mfs[var_idx];

Expand All @@ -69,32 +76,34 @@ ERF::fill_from_metgrid (const Vector<MultiFab*>& mfs, const Real time)
const auto& dom_lo = amrex::lbound(domain);
const auto& dom_hi = amrex::ubound(domain);

// Offset only applys to cons (we may fill a subset of these vars)
int offset = (var_idx == Vars::cons) ? icomp_cons : 0;

// Loop over each component
for (int comp_idx(0); comp_idx < comp_var[var_idx]; ++comp_idx)
{
int width;
int width = metgrid_bdy_set_width;

// Variable can be read from met_em files
//------------------------------------
if (is_read[var_idx][comp_idx])
{
width = metgrid_bdy_width;
int ivar = ind_map[var_idx][comp_idx];
IntVect ng_vect = mf.nGrowVect(); ng_vect[2] = 0;

if (ivar == MetGridBdyVars::U) {
amrex::Print() << "fill_from_metgrid U var_idx=" << var_idx << " comp_idx=" << comp_idx << " ivar=" << ivar << std::endl;
} else if (ivar == MetGridBdyVars::V) {
amrex::Print() << "fill_from_metgrid V var_idx=" << var_idx << " comp_idx=" << comp_idx << " ivar=" << ivar << std::endl;
} else if (ivar == MetGridBdyVars::R) {
amrex::Print() << "fill_from_metgrid R var_idx=" << var_idx << " comp_idx=" << comp_idx << " ivar=" << ivar << std::endl;
} else if (ivar == MetGridBdyVars::T) {
amrex::Print() << "fill_from_metgrid T var_idx=" << var_idx << " comp_idx=" << comp_idx << " ivar=" << ivar << std::endl;
} else if (ivar == MetGridBdyVars::QV) {
amrex::Print() << "fill_from_metgrid QV var_idx=" << var_idx << " comp_idx=" << comp_idx << " ivar=" << ivar << std::endl;
} else {
amrex::Print() << "fill_from_metgrid UNKNOWN" << std::endl;
}
// if (ivar == MetGridBdyVars::U) {
// amrex::Print() << "fill_from_metgrid U var_idx=" << var_idx << " comp_idx=" << comp_idx << " ivar=" << ivar << std::endl;
// } else if (ivar == MetGridBdyVars::V) {
// amrex::Print() << "fill_from_metgrid V var_idx=" << var_idx << " comp_idx=" << comp_idx << " ivar=" << ivar << std::endl;
// } else if (ivar == MetGridBdyVars::R) {
// amrex::Print() << "fill_from_metgrid R var_idx=" << var_idx << " comp_idx=" << comp_idx << " ivar=" << ivar << std::endl;
// } else if (ivar == MetGridBdyVars::T) {
// amrex::Print() << "fill_from_metgrid T var_idx=" << var_idx << " comp_idx=" << comp_idx << " ivar=" << ivar << std::endl;
// } else if (ivar == MetGridBdyVars::QV) {
// amrex::Print() << "fill_from_metgrid QV var_idx=" << var_idx << " comp_idx=" << comp_idx << " ivar=" << ivar << std::endl;
// } else {
// amrex::Print() << "fill_from_metgrid UNKNOWN" << std::endl;
// }

// We have data at fixed time intervals we will call dT
// Then to interpolate, given time, we can define n = (time/dT)
Expand Down Expand Up @@ -173,7 +182,6 @@ ERF::fill_from_metgrid (const Vector<MultiFab*>& mfs, const Real time)
// Grown tilebox so we fill exterior ghost cells as well
Box gbx = mfi.growntilebox(ng_vect);
const Array4<Real>& dest_arr = mf.array(mfi);

Box bx_xlo, bx_xhi, bx_ylo, bx_yhi;
compute_interior_ghost_bxs_xy(gbx, domain, width, 0,
bx_xlo, bx_xhi,
Expand Down
2 changes: 1 addition & 1 deletion Source/BoundaryConditions/ERF_FillPatch.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -251,7 +251,7 @@ ERF::FillIntermediatePatch (int lev, Real time,

// We call this here because it is an ERF routine
if (init_type == "real" && lev==0) fill_from_wrfbdy(mfs,time, cons_only, icomp_cons, ncomp_cons);
if (init_type == "metgrid" && lev==0) fill_from_metgrid(mfs,time);
if (init_type == "metgrid" && lev==0) fill_from_metgrid(mfs,time, cons_only, icomp_cons, ncomp_cons);
#endif

if (m_r2d) fill_from_bndryregs(mfs,time);
Expand Down
9 changes: 6 additions & 3 deletions Source/ERF.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -302,7 +302,10 @@ public:
int icomp_cons = 0,
int ncomp_cons = NVAR);
void fill_from_metgrid (const amrex::Vector<amrex::MultiFab*>& mfs,
amrex::Real time);
amrex::Real time,
bool cons_only = false,
int icomp_cons = 0,
int ncomp_cons = NVAR);
#endif

#ifdef ERF_USE_PARTICLES
Expand Down Expand Up @@ -697,8 +700,8 @@ private:
int wrfbdy_set_width{0};

// NetCDF initialization (met_em) file
int metgrid_bdy_width = 5;
int metgrid_bdy_set_width{1};
int metgrid_bdy_width{0};
int metgrid_bdy_set_width{0};

// Text input_sounding file
static std::string input_sounding_file;
Expand Down
7 changes: 7 additions & 0 deletions Source/ERF.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1015,6 +1015,13 @@ ERF::ReadParameters ()
AMREX_ALWAYS_ASSERT(wrfbdy_set_width >= 0);
AMREX_ALWAYS_ASSERT(wrfbdy_width >= wrfbdy_set_width);

// Query the set and total widths for metgrid_bdy interior ghost cells
pp.query("metgrid_bdy_width", metgrid_bdy_width);
pp.query("metgrid_bdy_set_width", metgrid_bdy_set_width);
AMREX_ALWAYS_ASSERT(metgrid_bdy_width >= 0);
AMREX_ALWAYS_ASSERT(metgrid_bdy_set_width >= 0);
AMREX_ALWAYS_ASSERT(metgrid_bdy_width >= metgrid_bdy_set_width);

// Query the set and total widths for crse-fine interior ghost cells
pp.query("cf_width", cf_width);
pp.query("cf_set_width", cf_set_width);
Expand Down
19 changes: 12 additions & 7 deletions Source/Initialization/ERF_init_from_metgrid.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -435,6 +435,8 @@ ERF::init_from_metgrid(int lev)
// the lateral boundary arrays using the info set aside earlier.
for (int it(0); it < ntimes; it++) {

const Array4<Real const>& R_bcs_arr = fabs_for_bcs[it][MetGridBdyVars::R].const_array();

for (int ivar(MetGridBdyVars::U); ivar < MetGridBdyEnd; ivar++) {

auto xlo_arr = bdy_data_xlo[it][ivar].array();
Expand Down Expand Up @@ -523,18 +525,22 @@ ERF::init_from_metgrid(int lev)
// mixing ratio at west boundary
amrex::ParallelFor(xlo_plane_no_stag, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
xlo_arr(i,j,k,0) = fabs_for_bcs_arr(i,j,k);
xlo_arr(i,j,k,0) *= R_bcs_arr(i,j,k);
});
// mixing ratio at east boundary
amrex::ParallelFor(xhi_plane_no_stag, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
xhi_arr(i,j,k,0) = fabs_for_bcs_arr(i,j,k);
xhi_arr(i,j,k,0) *= R_bcs_arr(i,j,k);
});
// mixing ratio at south boundary
amrex::ParallelFor(ylo_plane_no_stag, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
ylo_arr(i,j,k,0) = fabs_for_bcs_arr(i,j,k);
ylo_arr(i,j,k,0) *= R_bcs_arr(i,j,k);
});
// mixing ratio at north boundary
amrex::ParallelFor(yhi_plane_no_stag, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
yhi_arr(i,j,k,0) = fabs_for_bcs_arr(i,j,k);
yhi_arr(i,j,k,0) *= R_bcs_arr(i,j,k);
});
} // MetGridBdyVars::QV

Expand Down Expand Up @@ -784,7 +790,6 @@ init_state_from_metgrid(const Real l_rdOcp,

ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
rh_to_mxrat(i,j,k,rhum,temp,pres,mxrat);
// mxrat(i,j,k) = 0.0; // TODO: Remove when not needed. Force a dry atmosphere for debugging.
});
}

Expand Down Expand Up @@ -946,7 +951,7 @@ calc_rho_p(const int& kmax,
// calculate virtual potential temperature at the surface.
{
#if defined(ERF_USE_MOISTURE) || defined(ERF_USE_WARM_NO_PRECIP)
Real qvf = 1.0+(R_v/R_d-1.0)*Q_vec[0];
Real qvf = 1.0+(R_v/R_d+1.0)*Q_vec[0];
#else
Real qvf = 1.0;
#endif
Expand All @@ -960,7 +965,7 @@ calc_rho_p(const int& kmax,
for (int k=1; k < kmax; k++) {
Real dz = z_vec[k]-z_vec[k-1];
#if defined(ERF_USE_MOISTURE) || defined(ERF_USE_WARM_NO_PRECIP)
Real qvf = 1.0+(R_v/R_d-1.0)*Q_vec[k];
Real qvf = 1.0+(R_v/R_d+1.0)*Q_vec[k];
#else
Real qvf = 1.0;
#endif
Expand All @@ -975,14 +980,14 @@ calc_rho_p(const int& kmax,

// integrate from the top back down to get dry pressure and density.
Pd_vec[kmax-1] = Pm_vec[kmax-1];
Rhod_vec[kmax-1] = 1.0/(R_d/p_0*Theta_vec[kmax-1]*std::pow(Pd_vec[kmax-1]/p_0, -iGamma));
Rhod_vec[kmax-1] = 1.0/(R_d/p_0*theta_v[kmax-1]*std::pow(Pd_vec[kmax-1]/p_0, -iGamma));
for (int k=kmax-2; k >= 0; k--) {
Real dz = z_vec[k+1]-z_vec[k];
Rhod_vec[k] = Rhod_vec[k+1]; // an initial guess.
for (int it=0; it < maxiter; it++) {
Pd_vec[k] = Pd_vec[k+1]+0.5*dz*(Rhod_vec[k]+Rhod_vec[k+1])*CONST_GRAV;
if (Pd_vec[k] < 0.0) Pd_vec[k] = 0.0;
Rhod_vec[k] = 1.0/(R_d/p_0*Theta_vec[k]*std::pow(Pd_vec[k]/p_0, -iGamma));
Rhod_vec[k] = 1.0/(R_d/p_0*theta_v[k]*std::pow(Pd_vec[k]/p_0, -iGamma));
} // it
} // k
}
Expand Down Expand Up @@ -1137,8 +1142,8 @@ init_base_state_from_metgrid(const Real l_rdOcp,
Real RhoQ = r_hse_arr(i,j,k)*new_data(i,j,k,RhoQ_comp);
new_data(i,j,k,RhoQ_comp) = RhoQ;
#endif
// pi_hse_arr(i,j,k) = getExnergivenP(p_hse_arr(i,j,k), l_rdOcp);
pi_hse_arr(i,j,k) = getExnergivenRTh(RhoTheta, l_rdOcp);
pi_hse_arr(i,j,k) = getExnergivenP(p_hse_arr(i,j,k), l_rdOcp);
// pi_hse_arr(i,j,k) = getExnergivenRTh(RhoTheta, l_rdOcp);
});
}

Expand Down
2 changes: 1 addition & 1 deletion Source/TimeIntegration/TI_slow_rhs_fun.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -272,6 +272,7 @@
// Flag for moisture relaxation zone
bool moist_zero = false;
if (init_type=="real" && level==0 && wrfbdy_set_width>0) moist_zero = true;
if (init_type=="metgrid" && level==0 && metgrid_bdy_set_width > 0) moist_zero = true;
#endif

// Moving terrain
Expand Down Expand Up @@ -358,7 +359,6 @@
width = metgrid_bdy_width;
set_width = metgrid_bdy_set_width;
}
amrex::Print() << " DJW[TI_slow_rhs_fun.H]: start_bdy_time \t" << start_bdy_time << std::endl;
wrfbdy_compute_interior_ghost_rhs(init_type, bdy_time_interval, start_bdy_time, new_stage_time, slow_dt,
width-1, set_width, fine_geom,
S_rhs, S_data, bdy_data_xlo, bdy_data_xhi,
Expand Down
53 changes: 5 additions & 48 deletions Source/Utils/InteriorGhostCells.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -128,21 +128,19 @@ wrfbdy_compute_interior_ghost_rhs (const std::string& init_type,

// Time interpolation
Real dT = bdy_time_interval;
Real time_since_start = time;
if (init_type == "real") time_since_start = (time - start_bdy_time) / 1.e10;
Real time_since_start = (time - start_bdy_time) / 1.e10;
if (init_type == "metgrid") time_since_start = time;
int n_time = static_cast<int>( time_since_start / dT);
amrex::Real alpha = (time_since_start - n_time * dT) / dT;
AMREX_ALWAYS_ASSERT( alpha >= 0. && alpha <= 1.0);
amrex::Real oma = 1.0 - alpha;
// amrex::Print() << "time: " << time << " start_bdy_time: " << start_bdy_time << " dT: " << dT << " n_time: " << n_time << " alpha: " << alpha << std::endl;

// Temporary FABs for storage (owned/filled on all ranks)
FArrayBox U_xlo, U_xhi, U_ylo, U_yhi;
FArrayBox V_xlo, V_xhi, V_ylo, V_yhi;
FArrayBox R_xlo, R_xhi, R_ylo, R_yhi;
FArrayBox T_xlo, T_xhi, T_ylo, T_yhi;
#if defined(ERF_USE_MOISTURE) || defined(ERF_USE_WARM_NO_PRECIP)
FArrayBox Q_xlo, Q_xhi, Q_ylo, Q_yhi;
#endif

// Variable index map (WRFBdyVars -> Vars)
Vector<int> var_map = {Vars::xvel, Vars::yvel, Vars::cons, Vars::cons};
Expand All @@ -151,40 +149,19 @@ wrfbdy_compute_interior_ghost_rhs (const std::string& init_type,
// Variable icomp map
Vector<int> comp_map = {0, 0, Rho_comp, RhoTheta_comp};

//#if defined(ERF_USE_MOISTURE)
// var_map.push_back( Vars::cons );
// ivar_map.push_back( IntVar::cons );
// comp_map.push_back( RhoQt_comp );
//#elif defined(ERF_USE_WARM_NO_PRECIP)
// var_map.push_back( Vars::cons );
// ivar_map.push_back( IntVar::cons );
// comp_map.push_back( RhoQv_comp );
//#endif

int BdyEnd, ivarU, ivarV, ivarR, ivarT;
#if defined(ERF_USE_MOISTURE) || defined(ERF_USE_WARM_NO_PRECIP)
int ivarQV;
#endif
if (init_type == "real") {
BdyEnd = WRFBdyVars::NumTypes-3;
ivarU = WRFBdyVars::U;
ivarV = WRFBdyVars::V;
ivarR = WRFBdyVars::R;
ivarT = WRFBdyVars::T;
//#if defined(ERF_USE_MOISTURE) || defined(ERF_USE_WARM_NO_PRECIP)
// BdyEnd += 1;
// ivarQV = WRFBdyVars::QV;
//#endif
} else if (init_type == "metgrid") {
BdyEnd = MetGridBdyVars::NumTypes-1;
ivarU = MetGridBdyVars::U;
ivarV = MetGridBdyVars::V;
ivarR = MetGridBdyVars::R;
ivarT = MetGridBdyVars::T;
//#if defined(ERF_USE_MOISTURE) || defined(ERF_USE_WARM_NO_PRECIP)
// BdyEnd += 1;
// ivarQV = MetGridBdyVars::QV;
//#endif
}

// Size the FABs
Expand Down Expand Up @@ -218,11 +195,6 @@ wrfbdy_compute_interior_ghost_rhs (const std::string& init_type,
} else if (ivar == ivarT){
T_xlo.resize(bx_xlo,1); T_xhi.resize(bx_xhi,1);
T_ylo.resize(bx_ylo,1); T_yhi.resize(bx_yhi,1);
#if defined(ERF_USE_MOISTURE) || defined(ERF_USE_WARM_NO_PRECIP)
} else if (ivar == ivarQV ) {
Q_xlo.resize(bx_xlo,1); Q_xhi.resize(bx_xhi,1);
Q_ylo.resize(bx_ylo,1); Q_yhi.resize(bx_yhi,1);
#endif
} else {
continue;
}
Expand All @@ -241,11 +213,6 @@ wrfbdy_compute_interior_ghost_rhs (const std::string& init_type,
Elixir T_xlo_eli = T_xlo.elixir(); Elixir T_xhi_eli = T_xhi.elixir();
Elixir T_ylo_eli = T_ylo.elixir(); Elixir T_yhi_eli = T_yhi.elixir();

#if defined(ERF_USE_MOISTURE) || defined(ERF_USE_WARM_NO_PRECIP)
Elixir Q_xlo_eli = Q_xlo.elixir(); Elixir Q_xhi_eli = Q_xhi.elixir();
Elixir Q_ylo_eli = Q_ylo.elixir(); Elixir Q_yhi_eli = Q_yhi.elixir();
#endif

// Populate FABs from boundary interpolation
//==========================================================
for (int ivar(ivarU); ivar < BdyEnd; ivar++)
Expand Down Expand Up @@ -280,11 +247,6 @@ wrfbdy_compute_interior_ghost_rhs (const std::string& init_type,
} else if (ivar == ivarT){
arr_xlo = T_xlo.array(); arr_xhi = T_xhi.array();
arr_ylo = T_ylo.array(); arr_yhi = T_yhi.array();
#if defined(ERF_USE_MOISTURE) || defined(ERF_USE_WARM_NO_PRECIP)
} else if (ivar == ivarQV ) {
arr_xlo = Q_xlo.array(); arr_xhi = Q_xhi.array();
arr_ylo = Q_ylo.array(); arr_yhi = Q_yhi.array();
#endif
} else {
continue;
}
Expand Down Expand Up @@ -477,6 +439,8 @@ wrfbdy_compute_interior_ghost_rhs (const std::string& init_type,
}
} // mfi

// return; // DJW debugging return statement to shut off relaxation zone.

// Compute RHS in relaxation region
//==========================================================
for (int ivar(ivarU); ivar < BdyEnd; ivar++)
Expand Down Expand Up @@ -526,13 +490,6 @@ wrfbdy_compute_interior_ghost_rhs (const std::string& init_type,
arr_ylo = T_ylo.array(); arr_yhi = T_yhi.array();
rhs_arr = S_rhs[IntVar::cons].array(mfi);
data_arr = S_data[IntVar::cons].array(mfi);
#if defined(ERF_USE_MOISTURE) || defined(ERF_USE_WARM_NO_PRECIP)
} else if (ivar == ivarQV ) {
arr_xlo = Q_xlo.array(); arr_xhi = Q_xhi.array();
arr_ylo = Q_ylo.array(); arr_yhi = Q_yhi.array();
rhs_arr = S_rhs[IntVar::cons].array(mfi);
data_arr = S_data[IntVar::cons].array(mfi);
#endif
} else {
continue;
}
Expand Down
2 changes: 1 addition & 1 deletion Submodules/AMReX
Submodule AMReX updated 601 files

0 comments on commit df3f0db

Please sign in to comment.