Skip to content

Commit

Permalink
lon and lat for real data cases (#1589)
Browse files Browse the repository at this point in the history 
* Reads in latitude and longitude at mass points for real data cases (init_type of "metgrid" or "real"). These longitude and latitude arrays can be saved in the plotfiles by adding lat_m and lon_m to plot_vars_1 in inputs.

* Updated docs regarding plotfiles to include lat_m and lon_m. Cleaned a stray bit of commented code.

* Fixed whoopsie arising when compiling without netcdf. lat_m and lon_m are only filled for init_type of metgrid or real, both of which will necessitate compiling with netcdf.

---------

Co-authored-by: Aaron M. Lattanzi <[email protected]>
  • Loading branch information
@wiersema1 @AMLattanzi
wiersema1 and AMLattanzi authored Apr 24, 2024
1 parent a750eaa commit 3403574
Show file tree
Hide file tree
Showing 8 changed files with 127 additions and 4 deletions.
8 changes: 8 additions & 0 deletions Docs/sphinx_doc/Plotfiles.rst
View file
Original file line number Diff line number Diff line change
Expand Up @@ -251,6 +251,14 @@ Output Options
| | |
| | |
+-----------------------------+------------------+
| **lat_m** | Latitude at mass |
| | points |
| | |
+-----------------------------+------------------+
| **lon_m** | Longitude at |
| | mass points |
| | |
+-----------------------------+------------------+
| **Kmv** | Vertical |
| | Eddy Diffusivity |
| | of Momentum |
Expand Down
2 changes: 1 addition & 1 deletion Exec/DevTests/MetGrid/inputs
View file
Original file line number Diff line number Diff line change
Expand Up @@ -44,7 +44,7 @@ erf.check_int = 100 # number of timesteps between checkpoints
# PLOTFILES
erf.plot_file_1 = plt # prefix of plotfile name
erf.plot_int_1 = 10 # 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_vars_1 = lat_m lon_m density dens_hse rhoadv_0 x_velocity y_velocity z_velocity pressure temp theta z_phys mapfac pres_hse pert_pres

# SOLVER CHOICE
erf.alpha_T = 1.0
Expand Down
2 changes: 1 addition & 1 deletion Exec/RegTests/WPS_Test/inputs_real_ChisholmView
View file
Original file line number Diff line number Diff line change
Expand Up @@ -39,7 +39,7 @@ erf.check_int = -100 # number of timesteps between checkpoints
# PLOTFILES
erf.plot_file_1 = plt # prefix of plotfile name
erf.plot_int_1 = 1 # number of timesteps between plotfiles
erf.plot_vars_1 = density rhoQ1 x_velocity y_velocity z_velocity pressure temp theta pres_hse z_phys qt qp qv qc
erf.plot_vars_1 = lat_m lon_m density rhoQ1 x_velocity y_velocity z_velocity pressure temp theta pres_hse z_phys qt qp qv qc

# SOLVER CHOICE
erf.alpha_T = 0.0
Expand Down
3 changes: 2 additions & 1 deletion Source/ERF.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -537,6 +537,7 @@ private:
amrex::Vector<amrex::Vector<amrex::FArrayBox>> bdy_data_yhi;

amrex::Real bdy_time_interval;
amrex::Vector<std::unique_ptr<amrex::MultiFab>> lat_m, lon_m;
amrex::Real Latitude;
amrex::Real Longitude;
#endif // ERF_USE_NETCDF
Expand Down Expand Up @@ -792,7 +793,7 @@ private:
"magvel", "divU",
"pres_hse", "dens_hse", "pressure", "pert_pres", "pert_dens", "eq_pot_temp", "num_turb",
"dpdx", "dpdy", "pres_hse_x", "pres_hse_y",
"z_phys", "detJ" , "mapfac",
"z_phys", "detJ" , "mapfac", "lat_m", "lon_m",
// Time averaged velocity
"u_t_avg", "v_t_avg", "w_t_avg", "umag_t_avg",
// eddy viscosity
Expand Down
8 changes: 8 additions & 0 deletions Source/ERF.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -269,6 +269,14 @@ ERF::ERF ()
vel_t_avg.resize(nlevs_max);
t_avg_cnt.resize(nlevs_max);

#ifdef ERF_USE_NETCDF
// Longitude and latitude (only filled for use_real_bcs==True)
if (use_real_bcs) {
lat_m.resize(nlevs_max);
lon_m.resize(nlevs_max);
}
#endif

// Initialize tagging criteria for mesh refinement
refinement_criteria_setup();

Expand Down
39 changes: 39 additions & 0 deletions Source/IO/Plotfile.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -87,6 +87,9 @@ ERF::setPlotVariables (const std::string& pp_plot_var_names, Vector<std::string>
#endif

plot_var_names = tmp_plot_names;

for (int i=0; i<plot_var_names.size(); i++) {
}
}

void
Expand Down Expand Up @@ -749,6 +752,42 @@ ERF::WritePlotFile (int which, Vector<std::string> plot_var_names)
mf_comp ++;
}

#ifdef ERF_USE_NETCDF
if (use_real_bcs) {
if (containerHasElement(plot_var_names, "lat_m")) {
#ifdef _OPENMP
#pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
#endif
for ( MFIter mfi(mf[lev],TilingIfNotGPU()); mfi.isValid(); ++mfi)
{
const Box& bx = mfi.tilebox();
const Array4<Real>& derdat = mf[lev].array(mfi);
const Array4<Real>& data = lat_m[lev]->array(mfi);
ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
derdat(i, j, k, mf_comp) = data(i,j,0);
});
}
mf_comp ++;
} // lat_m
if (containerHasElement(plot_var_names, "lon_m")) {
#ifdef _OPENMP
#pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
#endif
for ( MFIter mfi(mf[lev],TilingIfNotGPU()); mfi.isValid(); ++mfi)
{
const Box& bx = mfi.tilebox();
const Array4<Real>& derdat = mf[lev].array(mfi);
const Array4<Real>& data = lon_m[lev]->array(mfi);
ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
derdat(i, j, k, mf_comp) = data(i,j,0);
});
}
mf_comp ++;
} // lon_m
} // use_real_bcs
#endif


if (solverChoice.time_avg_vel) {
if (containerHasElement(plot_var_names, "u_t_avg")) {
#ifdef _OPENMP
Expand Down
30 changes: 29 additions & 1 deletion Source/Initialization/ERF_init_from_metgrid.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -140,7 +140,7 @@ ERF::init_from_metgrid (int lev)
// This defines all the z(i,j,k) values given z(i,j,0) from above.
init_terrain_grid(lev, geom[lev], *z_phys, zlevels_stag);

// Copy SST and LANDMASK data into MF and iMF data structures
// Copy LATITUDE, LONGITUDE, SST and LANDMASK data into MF and iMF data structures
auto& ba = lev_new[Vars::cons].boxArray();
auto& dm = lev_new[Vars::cons].DistributionMap();
auto ngv = lev_new[Vars::cons].nGrowVect(); ngv[2] = 0;
Expand Down Expand Up @@ -193,6 +193,34 @@ ERF::init_from_metgrid (int lev)
} else {
for (int it = 0; it < ntimes; ++it) lmask_lev[lev][it] = nullptr;
}
lat_m[lev] = std::make_unique<MultiFab>(ba2d,dm,1,ngv);
for ( MFIter mfi(*(lat_m[lev]), TilingIfNotGPU()); mfi.isValid(); ++mfi ) {
Box gtbx = mfi.growntilebox();
FArrayBox& dst = (*(lat_m[lev]))[mfi];
FArrayBox& src = NC_LAT_fab[0];
const Array4< Real>& dst_arr = dst.array();
const Array4<const Real>& src_arr = src.const_array();
ParallelFor(gtbx, [=] AMREX_GPU_DEVICE (int i, int j, int) noexcept
{
int li = amrex::min(amrex::max(i, i_lo), i_hi);
int lj = amrex::min(amrex::max(j, j_lo), j_hi);
dst_arr(i,j,0) = src_arr(li,lj,0);
});
}
lon_m[lev] = std::make_unique<MultiFab>(ba2d,dm,1,ngv);
for ( MFIter mfi(*(lon_m[lev]), TilingIfNotGPU()); mfi.isValid(); ++mfi ) {
Box gtbx = mfi.growntilebox();
FArrayBox& dst = (*(lon_m[lev]))[mfi];
FArrayBox& src = NC_LON_fab[0];
const Array4< Real>& dst_arr = dst.array();
const Array4<const Real>& src_arr = src.const_array();
ParallelFor(gtbx, [=] AMREX_GPU_DEVICE (int i, int j, int) noexcept
{
int li = amrex::min(amrex::max(i, i_lo), i_hi);
int lj = amrex::min(amrex::max(j, j_lo), j_hi);
dst_arr(i,j,0) = src_arr(li,lj,0);
});
}

for (int it = 0; it < ntimes; it++) {
// Verify that the grid size and resolution from met_em file matches that in geom (from ERF inputs file).
Expand Down
39 changes: 39 additions & 0 deletions Source/Initialization/ERF_init_from_wrfinput.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -185,6 +185,45 @@ ERF::init_from_wrfinput (int lev)
solverChoice.moisture_type, n_qstate);
} // mf

auto& ba = lev_new[Vars::cons].boxArray();
auto& dm = lev_new[Vars::cons].DistributionMap();
auto ngv = lev_new[Vars::cons].nGrowVect(); ngv[2] = 0;
BoxList bl2d = ba.boxList();
for (auto& b : bl2d) {
b.setRange(2,0);
}
BoxArray ba2d(std::move(bl2d));
int i_lo = geom[lev].Domain().smallEnd(0); int i_hi = geom[lev].Domain().bigEnd(0);
int j_lo = geom[lev].Domain().smallEnd(1); int j_hi = geom[lev].Domain().bigEnd(1);
lat_m[lev] = std::make_unique<MultiFab>(ba2d,dm,1,ngv);
for ( MFIter mfi(*(lat_m[lev]), TilingIfNotGPU()); mfi.isValid(); ++mfi ) {
Box gtbx = mfi.growntilebox();
FArrayBox& dst = (*(lat_m[lev]))[mfi];
FArrayBox& src = NC_LAT_fab[0];
const Array4< Real>& dst_arr = dst.array();
const Array4<const Real>& src_arr = src.const_array();
ParallelFor(gtbx, [=] AMREX_GPU_DEVICE (int i, int j, int) noexcept
{
int li = amrex::min(amrex::max(i, i_lo), i_hi);
int lj = amrex::min(amrex::max(j, j_lo), j_hi);
dst_arr(i,j,0) = src_arr(li,lj,0);
});
}
lon_m[lev] = std::make_unique<MultiFab>(ba2d,dm,1,ngv);
for ( MFIter mfi(*(lon_m[lev]), TilingIfNotGPU()); mfi.isValid(); ++mfi ) {
Box gtbx = mfi.growntilebox();
FArrayBox& dst = (*(lon_m[lev]))[mfi];
FArrayBox& src = NC_LON_fab[0];
const Array4< Real>& dst_arr = dst.array();
const Array4<const Real>& src_arr = src.const_array();
ParallelFor(gtbx, [=] AMREX_GPU_DEVICE (int i, int j, int) noexcept
{
int li = amrex::min(amrex::max(i, i_lo), i_hi);
int lj = amrex::min(amrex::max(j, j_lo), j_hi);
dst_arr(i,j,0) = src_arr(li,lj,0);
});
}

#ifdef _OPENMP
#pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
#endif
Expand Down

0 comments on commit 3403574

Please sign in to comment.