Skip to content

Commit

Permalink
Compiles but needs verification.
Browse files Browse the repository at this point in the history
  • Loading branch information
@AMLattanzi
AMLattanzi committed Oct 2, 2024
1 parent 53a6ac8 commit 5cf3db2
Show file tree
Hide file tree
Showing 11 changed files with 369 additions and 367 deletions.
2 changes: 1 addition & 1 deletion Exec/DevTests/MetGrid/inputs
View file
Original file line number Diff line number Diff line change
Expand Up @@ -7,7 +7,7 @@ amrex.fpe_trap_overflow = 1

# PROBLEM SIZE & GEOMETRY
geometry.prob_extent = 28000 16000 8000
amr.n_cell = 140 80 100
amr.n_cell = 140 80 100

geometry.is_periodic = 0 0 0

Expand Down
14 changes: 12 additions & 2 deletions Source/BoundaryConditions/ERF_BoundaryConditions_realbdy.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -29,8 +29,14 @@ ERF::fill_from_realbdy (const Vector<MultiFab*>& mfs,
Real oma = 1.0 - alpha;

// Flags for read vars and index mapping
Vector<int> cons_read = {1, 1, 0, 0, 0, 1, 0, 0, 0};
Vector<int> cons_map = {RealBdyVars::R, RealBdyVars::T, 0, 0, 0, RealBdyVars::QV, 0, 0, 0};
Vector<int> cons_read = {0, 1, 0,
0, 0, 1,
0, 0, 0,
0, 0};
Vector<int> cons_map = {Rho_comp, RealBdyVars::T, RhoKE_comp,
RhoQKE_comp, RhoScalar_comp, RealBdyVars::QV,
RhoQ2_comp, RhoQ3_comp, RhoQ4_comp,
RhoQ5_comp, RhoQ6_comp};

Vector<Vector<int>> is_read;
is_read.push_back( cons_read );
Expand Down Expand Up @@ -113,6 +119,7 @@ ERF::fill_from_realbdy (const Vector<MultiFab*>& mfs,
jj = std::min(jj, dom_hi.y);
dest_arr(i,j,k,comp_idx) = oma * bdatxlo_n (ii,jj,k,0)
+ alpha * bdatxlo_np1(ii,jj,k,0);
if (var_idx == Vars::cons) dest_arr(i,j,k,comp_idx) *= dest_arr(i,j,k,Rho_comp);
},
[=] AMREX_GPU_DEVICE (int i, int j, int k)
{
Expand All @@ -121,6 +128,7 @@ ERF::fill_from_realbdy (const Vector<MultiFab*>& mfs,
jj = std::min(jj, dom_hi.y);
dest_arr(i,j,k,comp_idx) = oma * bdatxhi_n (ii,jj,k,0)
+ alpha * bdatxhi_np1(ii,jj,k,0);
if (var_idx == Vars::cons) dest_arr(i,j,k,comp_idx) *= dest_arr(i,j,k,Rho_comp);
});

// y-faces (do not include exterior x ghost cells)
Expand All @@ -130,12 +138,14 @@ ERF::fill_from_realbdy (const Vector<MultiFab*>& mfs,
int jj = std::max(j , dom_lo.y);
dest_arr(i,j,k,comp_idx) = oma * bdatylo_n (i,jj,k,0)
+ alpha * bdatylo_np1(i,jj,k,0);
if (var_idx == Vars::cons) dest_arr(i,j,k,comp_idx) *= dest_arr(i,j,k,Rho_comp);
},
[=] AMREX_GPU_DEVICE (int i, int j, int k)
{
int jj = std::min(j , dom_hi.y);
dest_arr(i,j,k,comp_idx) = oma * bdatyhi_n (i,jj,k,0)
+ alpha * bdatyhi_np1(i,jj,k,0);
if (var_idx == Vars::cons) dest_arr(i,j,k,comp_idx) *= dest_arr(i,j,k,Rho_comp);
});
} // mfi

Expand Down
11 changes: 4 additions & 7 deletions Source/ERF_IndexDefines.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -85,9 +85,8 @@ namespace RealBdyVars {
enum {
U = 0,
V = 1,
R = 2,
T = 3,
QV,
T = 2,
QV = 3,
NumTypes
};
}
Expand All @@ -96,8 +95,7 @@ namespace WRFBdyVars {
enum {
U = 0,
V = 1,
R = 2,
T = 3,
T = 2,
QV , // water vapor
MU , // bdy perturbation dry air mass in column (we will get mub from the initial data)
PC , // p_s - p_top = dry hydrostatic pressure difference between the surface and the model top
Expand All @@ -109,8 +107,7 @@ namespace MetGridBdyVars {
enum {
U = 0,
V = 1,
R = 2,
T = 3,
T = 2,
QV,
NumTypes
};
Expand Down
324 changes: 147 additions & 177 deletions Source/IO/ERF_ReadFromWRFBdy.cpp
View file

Large diffs are not rendered by default.

3 changes: 0 additions & 3 deletions Source/IO/ERF_ReadFromWRFInput.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -158,8 +158,5 @@ read_from_wrfinput (int lev,

const Real theta_ref = 300.0;
NC_rhotheta_fab.template plus<RunOn::Device>(theta_ref);

// Now multiply by rho to get (rho theta) instead of theta
NC_rhotheta_fab.template mult<RunOn::Device>(NC_rho_fab,0,0,1);
}
#endif // ERF_USE_NETCDF
5 changes: 2 additions & 3 deletions Source/Initialization/ERF_Metgrid_utils.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -97,8 +97,7 @@ init_base_state_from_metgrid (const bool use_moisture,
amrex::FArrayBox& z_phys_cc_fab,
const amrex::Vector<amrex::FArrayBox>& NC_ght_fab,
const amrex::Vector<amrex::FArrayBox>& NC_psfc_fab,
amrex::Vector<amrex::Vector<amrex::FArrayBox>>& fabs_for_bcs,
const amrex::Array4<const int>& mask_c_arr);
amrex::Vector<amrex::Vector<amrex::FArrayBox>>& fabs_for_bcs);

AMREX_FORCE_INLINE
AMREX_GPU_DEVICE
Expand All @@ -117,7 +116,7 @@ calc_rho_p (const int& kmax,
amrex::Real* Pm_vec)
{
const int maxiter = 10;
const amrex::Real tol = 1.0e-12;
const amrex::Real tol = 1.0e-10;

// Calculate or use moist pressure at the surface.
amrex::Real Psurf;
Expand Down
50 changes: 10 additions & 40 deletions Source/Initialization/ERF_init_from_metgrid.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -340,8 +340,6 @@ ERF::init_from_metgrid (int lev)
FArrayBox& cons_fab = lev_new[Vars::cons][mfi];
FArrayBox& z_phys_nd_fab = (*z_phys)[mfi];

const Array4<const int>& mask_c_arr = mask_c->const_array(mfi);

// Fill base state data using origin data (initialization and BC arrays)
// p_hse calculate dry pressure
// r_hse calculate dry density
Expand All @@ -352,7 +350,7 @@ ERF::init_from_metgrid (int lev)
flag_psfc,
cons_fab, r_hse_fab, p_hse_fab, pi_hse_fab,
z_phys_nd_fab, NC_ght_fab, NC_psfc_fab,
fabs_for_bcs, mask_c_arr);
fabs_for_bcs);
} // mf

// FillBoundary to populate the internal halo cells
Expand All @@ -378,7 +376,7 @@ ERF::init_from_metgrid (int lev)
// Otherwise, we make the total width match the set width.
if (real_width-1 <= real_set_width) real_width = real_set_width;
Print() << "Running with specification width: " << real_set_width
<< " and relaxation width: " << real_width - real_set_width << std::endl;
<< " and relaxation width: " << real_width - real_set_width << std::endl;

// Set up boxes for lateral boundary arrays.
bdy_data_xlo.resize(ntimes);
Expand Down Expand Up @@ -431,11 +429,6 @@ ERF::init_from_metgrid (int lev)
bdy_data_xhi[it].push_back(FArrayBox(xhi_plane_y_stag, 1));
bdy_data_ylo[it].push_back(FArrayBox(ylo_plane_y_stag, 1));
bdy_data_yhi[it].push_back(FArrayBox(yhi_plane_y_stag, 1));
} else if (ivar == MetGridBdyVars::R) {
bdy_data_xlo[it].push_back(FArrayBox(xlo_plane_no_stag, 1));
bdy_data_xhi[it].push_back(FArrayBox(xhi_plane_no_stag, 1));
bdy_data_ylo[it].push_back(FArrayBox(ylo_plane_no_stag, 1));
bdy_data_yhi[it].push_back(FArrayBox(yhi_plane_no_stag, 1));
} else if (ivar == MetGridBdyVars::T) {
bdy_data_xlo[it].push_back(FArrayBox(xlo_plane_no_stag, 1));
bdy_data_xhi[it].push_back(FArrayBox(xhi_plane_no_stag, 1));
Expand All @@ -459,12 +452,8 @@ ERF::init_from_metgrid (int lev)
// we only need the whole domain processed at initialization and the lateral boundaries
// at subsequent times. We can optimize this later if needed. For now, we need to fill
// the lateral boundary arrays using the info set aside earlier.
bool multiply_rho = false;
amrex::Box xlo_plane, xhi_plane, ylo_plane, yhi_plane;
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 All @@ -474,50 +463,38 @@ ERF::init_from_metgrid (int lev)
const Array4<Real const>& fabs_for_bcs_arr = fabs_for_bcs[it][ivar].const_array();

if (ivar == MetGridBdyVars::U) {
multiply_rho = false;
xlo_plane = xlo_plane_x_stag; xhi_plane = xhi_plane_x_stag;
ylo_plane = ylo_plane_x_stag; yhi_plane = yhi_plane_x_stag;
} else if (ivar == MetGridBdyVars::V) {
multiply_rho = false;
xlo_plane = xlo_plane_y_stag; xhi_plane = xhi_plane_y_stag;
ylo_plane = ylo_plane_y_stag; yhi_plane = yhi_plane_y_stag;
} else if (ivar == MetGridBdyVars::R) {
multiply_rho = false;
xlo_plane = xlo_plane_no_stag; xhi_plane = xhi_plane_no_stag;
ylo_plane = ylo_plane_no_stag; yhi_plane = yhi_plane_no_stag;
} else if (ivar == MetGridBdyVars::T) {
multiply_rho = false;
xlo_plane = xlo_plane_no_stag; xhi_plane = xhi_plane_no_stag;
ylo_plane = ylo_plane_no_stag; yhi_plane = yhi_plane_no_stag;
} else if (ivar == MetGridBdyVars::QV) {
multiply_rho = false;
xlo_plane = xlo_plane_no_stag; xhi_plane = xhi_plane_no_stag;
ylo_plane = ylo_plane_no_stag; yhi_plane = yhi_plane_no_stag;
} // MetGridBdyVars::QV

// west boundary
ParallelFor(xlo_plane, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
{
Real Factor = (multiply_rho) ? R_bcs_arr(i,j,k) : 1.0;
xlo_arr(i,j,k,0) = fabs_for_bcs_arr(i,j,k)*Factor;
xlo_arr(i,j,k,0) = fabs_for_bcs_arr(i,j,k);
});
// xvel at east boundary
ParallelFor(xhi_plane, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
{
Real Factor = (multiply_rho) ? R_bcs_arr(i,j,k) : 1.0;
xhi_arr(i,j,k,0) = fabs_for_bcs_arr(i,j,k)*Factor;
xhi_arr(i,j,k,0) = fabs_for_bcs_arr(i,j,k);
});
// xvel at south boundary
ParallelFor(ylo_plane, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
{
Real Factor = (multiply_rho) ? R_bcs_arr(i,j,k) : 1.0;
ylo_arr(i,j,k,0) = fabs_for_bcs_arr(i,j,k)*Factor;
ylo_arr(i,j,k,0) = fabs_for_bcs_arr(i,j,k);
});
// xvel at north boundary
ParallelFor(yhi_plane, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
{
Real Factor = (multiply_rho) ? R_bcs_arr(i,j,k) : 1.0;
yhi_arr(i,j,k,0) = fabs_for_bcs_arr(i,j,k)*Factor;
yhi_arr(i,j,k,0) = fabs_for_bcs_arr(i,j,k);
});

} // ivar
Expand Down Expand Up @@ -803,8 +780,7 @@ init_base_state_from_metgrid (const bool use_moisture,
FArrayBox& z_phys_cc_fab,
const Vector<FArrayBox>& /*NC_ght_fab*/,
const Vector<FArrayBox>& NC_psfc_fab,
Vector<Vector<FArrayBox>>& fabs_for_bcs,
const amrex::Array4<const int>& mask_c_arr)
Vector<Vector<FArrayBox>>& fabs_for_bcs)
{
int RhoQ_comp = RhoQ1_comp;
int kmax = amrex::ubound(valid_bx).z;
Expand Down Expand Up @@ -885,11 +861,11 @@ init_base_state_from_metgrid (const bool use_moisture,

ParallelFor(valid_bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
{
// Multiply by Rho to get conserved vars
// Multiply by Rho to get conserved vars
Real Qv = 0.0;
new_data(i,j,k,Rho_comp) = r_hse_arr(i,j,k);
new_data(i,j,k,Rho_comp) = r_hse_arr(i,j,k);
new_data(i,j,k,RhoTheta_comp) *= r_hse_arr(i,j,k);
if (use_moisture){
if (use_moisture) {
Qv = new_data(i,j,k,RhoQ_comp);
new_data(i,j,k,RhoQ_comp) *= r_hse_arr(i,j,k);
}
Expand Down Expand Up @@ -963,7 +939,6 @@ init_base_state_from_metgrid (const bool use_moisture,
valid_bx2d.setRange(2,0);
auto const orig_psfc = NC_psfc_fab[it].const_array();
auto const new_z = z_phys_cc_fab.const_array();
auto r_arr = fabs_for_bcs[it][MetGridBdyVars::R].array();
auto Theta_arr = fabs_for_bcs[it][MetGridBdyVars::T].array();
auto Q_arr = (use_moisture ) ? fabs_for_bcs[it][MetGridBdyVars::QV].array() : Array4<Real>{};
auto p_hse_arr = p_hse_bcs_fab.array();
Expand All @@ -983,11 +958,6 @@ init_base_state_from_metgrid (const bool use_moisture,

for (int k=0; k<=kmax; k++) {
p_hse_arr(i,j,k) = Pm_vec[k];
if (mask_c_arr(i,j,k)) {
r_arr(i,j,k) = Rhom_vec[k];
if (use_moisture) Q_arr(i,j,k) = Rhom_vec[k]*Q_vec[k];
Theta_arr(i,j,k) = Rhom_vec[k]*Thetad_vec[k];
}
} // k
});
} // it
Expand Down
Loading

0 comments on commit 5cf3db2

Please sign in to comment.