From 34035744de457457923333ae017d2cd7775c1398 Mon Sep 17 00:00:00 2001 From: wiersema1 <46633497+wiersema1@users.noreply.github.com> Date: Wed, 24 Apr 2024 15:20:59 -0700 Subject: [PATCH] lon and lat for real data cases (#1589) * 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 <103702284+AMLattanzi@users.noreply.github.com> --- Docs/sphinx_doc/Plotfiles.rst | 8 ++++ Exec/DevTests/MetGrid/inputs | 2 +- .../WPS_Test/inputs_real_ChisholmView | 2 +- Source/ERF.H | 3 +- Source/ERF.cpp | 8 ++++ Source/IO/Plotfile.cpp | 39 +++++++++++++++++++ .../Initialization/ERF_init_from_metgrid.cpp | 30 +++++++++++++- .../Initialization/ERF_init_from_wrfinput.cpp | 39 +++++++++++++++++++ 8 files changed, 127 insertions(+), 4 deletions(-) diff --git a/Docs/sphinx_doc/Plotfiles.rst b/Docs/sphinx_doc/Plotfiles.rst index 146249f4d..a13332a1d 100644 --- a/Docs/sphinx_doc/Plotfiles.rst +++ b/Docs/sphinx_doc/Plotfiles.rst @@ -251,6 +251,14 @@ Output Options | | | | | | +-----------------------------+------------------+ +| **lat_m** | Latitude at mass | +| | points | +| | | ++-----------------------------+------------------+ +| **lon_m** | Longitude at | +| | mass points | +| | | ++-----------------------------+------------------+ | **Kmv** | Vertical | | | Eddy Diffusivity | | | of Momentum | diff --git a/Exec/DevTests/MetGrid/inputs b/Exec/DevTests/MetGrid/inputs index ccf0ee7dc..6219d7d40 100644 --- a/Exec/DevTests/MetGrid/inputs +++ b/Exec/DevTests/MetGrid/inputs @@ -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 diff --git a/Exec/RegTests/WPS_Test/inputs_real_ChisholmView b/Exec/RegTests/WPS_Test/inputs_real_ChisholmView index dfd5efb4f..de034cab9 100644 --- a/Exec/RegTests/WPS_Test/inputs_real_ChisholmView +++ b/Exec/RegTests/WPS_Test/inputs_real_ChisholmView @@ -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 diff --git a/Source/ERF.H b/Source/ERF.H index 9e78dccf6..36052a29b 100644 --- a/Source/ERF.H +++ b/Source/ERF.H @@ -537,6 +537,7 @@ private: amrex::Vector> bdy_data_yhi; amrex::Real bdy_time_interval; + amrex::Vector> lat_m, lon_m; amrex::Real Latitude; amrex::Real Longitude; #endif // ERF_USE_NETCDF @@ -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 diff --git a/Source/ERF.cpp b/Source/ERF.cpp index 1044f3335..aec7f0ac3 100644 --- a/Source/ERF.cpp +++ b/Source/ERF.cpp @@ -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(); diff --git a/Source/IO/Plotfile.cpp b/Source/IO/Plotfile.cpp index c78dab745..e486d5314 100644 --- a/Source/IO/Plotfile.cpp +++ b/Source/IO/Plotfile.cpp @@ -87,6 +87,9 @@ ERF::setPlotVariables (const std::string& pp_plot_var_names, Vector #endif plot_var_names = tmp_plot_names; + + for (int i=0; i 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& derdat = mf[lev].array(mfi); + const Array4& 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& derdat = mf[lev].array(mfi); + const Array4& 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 diff --git a/Source/Initialization/ERF_init_from_metgrid.cpp b/Source/Initialization/ERF_init_from_metgrid.cpp index 2a0a6a505..71dcc1748 100644 --- a/Source/Initialization/ERF_init_from_metgrid.cpp +++ b/Source/Initialization/ERF_init_from_metgrid.cpp @@ -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; @@ -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(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& 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(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& 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). diff --git a/Source/Initialization/ERF_init_from_wrfinput.cpp b/Source/Initialization/ERF_init_from_wrfinput.cpp index 9e59f2adf..fb35b79b1 100644 --- a/Source/Initialization/ERF_init_from_wrfinput.cpp +++ b/Source/Initialization/ERF_init_from_wrfinput.cpp @@ -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(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& 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(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& 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