From 693ffc077095a65e90591f30e9f274a927199cfd Mon Sep 17 00:00:00 2001 From: Ann Almgren Date: Wed, 19 Jun 2024 17:22:38 -0700 Subject: [PATCH 1/9] This commit enables vertical refinement but there may still be some fillpatching issues (#1651) --- Source/BoundaryConditions/ERF_FillPatch.cpp | 70 ++++++++--- Source/BoundaryConditions/ERF_PhysBCFunct.H | 5 +- Source/BoundaryConditions/ERF_PhysBCFunct.cpp | 5 +- Source/ERF.H | 2 +- Source/ERF.cpp | 4 +- Source/ERF_Tagging.cpp | 3 +- Source/IO/Plotfile.cpp | 118 ++++++++++-------- Source/Initialization/ERF_init1d.cpp | 2 + 8 files changed, 133 insertions(+), 76 deletions(-) diff --git a/Source/BoundaryConditions/ERF_FillPatch.cpp b/Source/BoundaryConditions/ERF_FillPatch.cpp index dc630079b..e6c7c7fc9 100644 --- a/Source/BoundaryConditions/ERF_FillPatch.cpp +++ b/Source/BoundaryConditions/ERF_FillPatch.cpp @@ -108,7 +108,7 @@ ERF::FillPatch (int lev, Real time, FillPatchSingleLevel(*mfs_vel[Vars::zvel], time, fmf, ftime, icomp, icomp, ncomp_w, geom[lev], *physbcs_w_no_terrain[lev], BCVars::zvel_bc); (*physbcs_w[lev])(*mfs_vel[Vars::zvel],*mfs_vel[Vars::xvel],*mfs_vel[Vars::yvel], - ngvect_vels,time,BCVars::xvel_bc,BCVars::yvel_bc,BCVars::zvel_bc); + ngvect_vels,time,BCVars::zvel_bc); } // !cons_only } else { @@ -122,8 +122,9 @@ ERF::FillPatch (int lev, Real time, mapper = &cell_cons_interp; FillPatchTwoLevels(mf_c, time, cmf, ctime, fmf, ftime, 0, 0, mf_c.nComp(), geom[lev-1], geom[lev], - null_bc, BCVars::cons_bc, null_bc, BCVars::cons_bc, refRatio(lev-1), - mapper, domain_bcs_type, BCVars::cons_bc); + *physbcs_cons[lev-1], BCVars::cons_bc, + *physbcs_cons[lev ], BCVars::cons_bc, + refRatio(lev-1), mapper, domain_bcs_type, BCVars::cons_bc); if (!cons_only) { mapper = &face_cons_linear_interp; @@ -133,24 +134,34 @@ ERF::FillPatch (int lev, Real time, cmf = {&vars_old[lev-1][Vars::xvel], &vars_new[lev-1][Vars::xvel]}; FillPatchTwoLevels(mf_u, time, cmf, ctime, fmf, ftime, 0, 0, 1, geom[lev-1], geom[lev], - null_bc, BCVars::xvel_bc, null_bc, BCVars::xvel_bc, refRatio(lev-1), - mapper, domain_bcs_type, BCVars::xvel_bc); + *physbcs_u[lev-1], BCVars::xvel_bc, + *physbcs_u[lev ], BCVars::xvel_bc, + refRatio(lev-1), mapper, domain_bcs_type, BCVars::xvel_bc); MultiFab& mf_v = *mfs_vel[Vars::yvel]; fmf = {&vars_old[lev ][Vars::yvel], &vars_new[lev ][Vars::yvel]}; cmf = {&vars_old[lev-1][Vars::yvel], &vars_new[lev-1][Vars::yvel]}; FillPatchTwoLevels(mf_v, time, cmf, ctime, fmf, ftime, 0, 0, 1, geom[lev-1], geom[lev], - null_bc, BCVars::yvel_bc, null_bc, BCVars::yvel_bc, refRatio(lev-1), - mapper, domain_bcs_type, BCVars::yvel_bc); + *physbcs_v[lev-1], BCVars::yvel_bc, + *physbcs_v[lev ], BCVars::yvel_bc, + refRatio(lev-1), mapper, domain_bcs_type, BCVars::yvel_bc); + + // We note there is an issue here -- we use the no-terrain version to fill the physical + // bcs of the coarse data used to interpolate. We later fix the fine data with the + // correct bcs but if there was an error due to the interpolation with the wrong bcs, + // we will not necessarily be able to fix that. MultiFab& mf_w = *mfs_vel[Vars::zvel]; fmf = {&vars_old[lev ][Vars::zvel], &vars_new[lev ][Vars::zvel]}; cmf = {&vars_old[lev-1][Vars::zvel], &vars_new[lev-1][Vars::zvel]}; FillPatchTwoLevels(mf_w, time, cmf, ctime, fmf, ftime, 0, 0, 1, geom[lev-1], geom[lev], - null_bc, BCVars::zvel_bc, null_bc, BCVars::zvel_bc, refRatio(lev-1), - mapper, domain_bcs_type, BCVars::zvel_bc); + *physbcs_w_no_terrain[lev-1], BCVars::zvel_bc, + *physbcs_w_no_terrain[lev ], BCVars::zvel_bc, + refRatio(lev-1), mapper, domain_bcs_type, BCVars::zvel_bc); + (*physbcs_w[lev])(*mfs_vel[Vars::zvel],*mfs_vel[Vars::xvel],*mfs_vel[Vars::yvel], + ngvect_vels,time,BCVars::zvel_bc); } // !cons_only } // lev > 0 @@ -175,7 +186,7 @@ ERF::FillPatch (int lev, Real time, (*physbcs_u[lev])(*mfs_vel[Vars::xvel],0,1,ngvect_vels,time,BCVars::xvel_bc); (*physbcs_v[lev])(*mfs_vel[Vars::yvel],0,1,ngvect_vels,time,BCVars::yvel_bc); (*physbcs_w[lev])(*mfs_vel[Vars::zvel],*mfs_vel[Vars::xvel],*mfs_vel[Vars::yvel], - ngvect_vels,time,BCVars::xvel_bc,BCVars::yvel_bc,BCVars::zvel_bc); + ngvect_vels,time,BCVars::zvel_bc); } } @@ -287,7 +298,7 @@ ERF::FillIntermediatePatch (int lev, Real time, Geom(lev).Domain(), domain_bcs_type); } - // We now start working on VELOCITY + // We now start working on conserved quantities + VELOCITY for (int var_idx = 0; var_idx < Vars::NumTypes; ++var_idx) { if (cons_only && var_idx != Vars::cons) continue; @@ -331,7 +342,7 @@ ERF::FillIntermediatePatch (int lev, Real time, if (lev == 0) { - // This fills fine-fine ghost values of VELOCITY + // This fills fine-fine ghost values of cons and VELOCITY (not momentum) mf.FillBoundary(icomp,ncomp,ngvect,geom[lev].periodicity()); } else @@ -348,10 +359,33 @@ ERF::FillIntermediatePatch (int lev, Real time, Vector cmf = {&vars_old[lev-1][var_idx], &vars_new[lev-1][var_idx]}; Vector ctime = {t_old[lev-1], t_new[lev-1]}; - FillPatchTwoLevels(mf, time, cmf, ctime, fmf, {time}, - icomp, icomp, ncomp, geom[lev-1], geom[lev], - null_bc, 0, null_bc, 0, refRatio(lev-1), - mapper, domain_bcs_type, bccomp); + if (var_idx == Vars::cons) { + FillPatchTwoLevels(mf, time, cmf, ctime, fmf, {time}, + icomp, icomp, ncomp, geom[lev-1], geom[lev], + *physbcs_cons[lev-1], BCVars::cons_bc, + *physbcs_cons[lev ], BCVars::cons_bc, + refRatio(lev-1), mapper, domain_bcs_type, bccomp); + } else if (var_idx == Vars::xvel) { + FillPatchTwoLevels(mf, time, cmf, ctime, fmf, {time}, + icomp, icomp, ncomp, geom[lev-1], geom[lev], + *physbcs_u[lev-1], BCVars::xvel_bc, + *physbcs_u[lev ], BCVars::xvel_bc, + refRatio(lev-1), mapper, domain_bcs_type, bccomp); + } else if (var_idx == Vars::yvel) { + FillPatchTwoLevels(mf, time, cmf, ctime, fmf, {time}, + icomp, icomp, ncomp, geom[lev-1], geom[lev], + *physbcs_v[lev-1], BCVars::yvel_bc, + *physbcs_v[lev ], BCVars::yvel_bc, + refRatio(lev-1), mapper, domain_bcs_type, bccomp); + } else if (var_idx == Vars::zvel) { + FillPatchTwoLevels(mf, time, cmf, ctime, fmf, {time}, + icomp, icomp, ncomp, geom[lev-1], geom[lev], + *physbcs_w_no_terrain[lev-1], BCVars::zvel_bc, + *physbcs_w_no_terrain[lev ], BCVars::zvel_bc, + refRatio(lev-1), mapper, domain_bcs_type, bccomp); + (*physbcs_w[lev])(*mfs_vel[Vars::zvel],*mfs_vel[Vars::xvel],*mfs_vel[Vars::yvel], + ngvect,time,BCVars::zvel_bc); + } } // lev > 0 } // var_idx @@ -376,7 +410,7 @@ ERF::FillIntermediatePatch (int lev, Real time, (*physbcs_u[lev])(*mfs_vel[Vars::xvel],0,1,ngvect_vels,time,BCVars::xvel_bc); (*physbcs_v[lev])(*mfs_vel[Vars::yvel],0,1,ngvect_vels,time,BCVars::yvel_bc); (*physbcs_w[lev])(*mfs_vel[Vars::zvel],*mfs_vel[Vars::xvel],*mfs_vel[Vars::yvel], - ngvect_vels,time,BCVars::xvel_bc,BCVars::yvel_bc,BCVars::zvel_bc); + ngvect_vels,time,BCVars::zvel_bc); } // *************************************************************************** @@ -533,7 +567,7 @@ ERF::FillCoarsePatch (int lev, Real time) ( *physbcs_u[lev])(vars_new[lev][Vars::xvel],0,1 ,ngvect_vels,time,BCVars::xvel_bc); ( *physbcs_v[lev])(vars_new[lev][Vars::yvel],0,1 ,ngvect_vels,time,BCVars::yvel_bc); ( *physbcs_w[lev])(vars_new[lev][Vars::zvel],vars_new[lev][Vars::xvel],vars_new[lev][Vars::yvel], - ngvect_vels,time,BCVars::xvel_bc,BCVars::yvel_bc,BCVars::zvel_bc); + ngvect_vels,time,BCVars::zvel_bc); // *************************************************************************** // Since lev > 0 here we don't worry about m_r2d or wrfbdy data diff --git a/Source/BoundaryConditions/ERF_PhysBCFunct.H b/Source/BoundaryConditions/ERF_PhysBCFunct.H index 82eeedbbb..c0000105b 100644 --- a/Source/BoundaryConditions/ERF_PhysBCFunct.H +++ b/Source/BoundaryConditions/ERF_PhysBCFunct.H @@ -214,8 +214,7 @@ public: * @param[in] use_real_bcs if true then we fill boundary conditions for interior locations */ void operator() (amrex::MultiFab& mf, amrex::MultiFab& xvel, amrex::MultiFab& yvel, - amrex::IntVect const& nghost, const amrex::Real time, - int bccomp_u, int bccomp_v, int bccomp_w); + amrex::IntVect const& nghost, const amrex::Real time, int bccomp); void impose_lateral_zvel_bcs (const amrex::Array4& dest_arr, const amrex::Array4& xvel_arr, @@ -223,7 +222,7 @@ public: const amrex::Box& bx, const amrex::Box& domain, const amrex::Array4& z_nd, const amrex::GpuArray dxInv, - int bccomp_w); + int bccomp); void impose_vertical_zvel_bcs (const amrex::Array4& dest_arr, const amrex::Array4& xvel_arr, const amrex::Array4& yvel_arr, diff --git a/Source/BoundaryConditions/ERF_PhysBCFunct.cpp b/Source/BoundaryConditions/ERF_PhysBCFunct.cpp index 1e5e12d49..88b623da9 100644 --- a/Source/BoundaryConditions/ERF_PhysBCFunct.cpp +++ b/Source/BoundaryConditions/ERF_PhysBCFunct.cpp @@ -193,10 +193,13 @@ void ERFPhysBCFunct_v::operator() (MultiFab& mf, int /*icomp*/, int /*ncomp*/, void ERFPhysBCFunct_w::operator() (MultiFab& mf, MultiFab& xvel, MultiFab& yvel, IntVect const& nghost, const Real /*time*/, - const int bccomp_u, const int bccomp_v, const int bccomp_w) + const int bccomp_w) { BL_PROFILE("ERFPhysBCFunct_w::()"); + int bccomp_u = BCVars::xvel_bc; + int bccomp_v = BCVars::yvel_bc; + if (m_geom.isAllPeriodic()) return; const auto& domain = m_geom.Domain(); diff --git a/Source/ERF.H b/Source/ERF.H index e7122ed5a..a1ad20c05 100644 --- a/Source/ERF.H +++ b/Source/ERF.H @@ -318,7 +318,7 @@ public: // calls AmrCore constructor -> AmrMesh constructor ERF (const amrex::RealBox& rb, int max_level_in, const amrex::Vector& n_cell_in, int coord, - const amrex::Vector& ref_ratios, + const amrex::Vector& ref_ratio, const amrex::Array& is_per, std::string prefix); diff --git a/Source/ERF.cpp b/Source/ERF.cpp index 7b69f8428..a884317f8 100644 --- a/Source/ERF.cpp +++ b/Source/ERF.cpp @@ -1773,10 +1773,10 @@ ERF::Define_ERFFillPatchers (int lev) // constructor used when ERF is created by a multiblock driver ERF::ERF (const RealBox& rb, int max_level_in, const Vector& n_cell_in, int coord, - const Vector& ref_ratios, + const Vector& ref_ratio, const Array& is_per, std::string prefix) - : AmrCore(rb, max_level_in, n_cell_in, coord, ref_ratios, is_per) + : AmrCore(rb, max_level_in, n_cell_in, coord, ref_ratio, is_per) { SetParmParsePrefix(prefix); diff --git a/Source/ERF_Tagging.cpp b/Source/ERF_Tagging.cpp index e8af57f23..3dc2cff81 100644 --- a/Source/ERF_Tagging.cpp +++ b/Source/ERF_Tagging.cpp @@ -156,7 +156,8 @@ ERF::refinement_criteria_setup () int khi = static_cast((box_hi[2] - plo[2])/dx[2]-1); Box bx(IntVect(ilo,jlo,klo),IntVect(ihi,jhi,khi)); if ( (ilo%ref_ratio[lev_for_box-1][0] != 0) || ((ihi+1)%ref_ratio[lev_for_box-1][0] != 0) || - (jlo%ref_ratio[lev_for_box-1][1] != 0) || ((jhi+1)%ref_ratio[lev_for_box-1][1] != 0) ) + (jlo%ref_ratio[lev_for_box-1][1] != 0) || ((jhi+1)%ref_ratio[lev_for_box-1][1] != 0) || + (klo%ref_ratio[lev_for_box-1][2] != 0) || ((khi+1)%ref_ratio[lev_for_box-1][2] != 0) ) amrex::Error("Fine box is not legit with this ref_ratio"); boxes_at_level[lev_for_box].push_back(bx); Print() << "Saving in 'boxes at level' as " << bx << std::endl; diff --git a/Source/IO/Plotfile.cpp b/Source/IO/Plotfile.cpp index c9e0ff65e..bd547e6f3 100644 --- a/Source/IO/Plotfile.cpp +++ b/Source/IO/Plotfile.cpp @@ -1318,72 +1318,90 @@ ERF::WritePlotFile (int which, Vector plot_var_names) } else { // multilevel - Vector r2(finest_level); - Vector g2(finest_level+1); - Vector mf2(finest_level+1); + if (plotfile_type == "amrex") { - mf2[0].define(grids[0], dmap[0], ncomp_mf, 0); + if (ref_ratio[0][2] == 1) { - // Copy level 0 as is - MultiFab::Copy(mf2[0],mf[0],0,0,mf[0].nComp(),0); + Vector r2(finest_level); + Vector g2(finest_level+1); + Vector mf2(finest_level+1); - // Define a new multi-level array of Geometry's so that we pass the new "domain" at lev > 0 - Array periodicity = - {Geom()[0].isPeriodic(0),Geom()[0].isPeriodic(1),Geom()[0].isPeriodic(2)}; - g2[0].define(Geom()[0].Domain(),&(Geom()[0].ProbDomain()),0,periodicity.data()); + mf2[0].define(grids[0], dmap[0], ncomp_mf, 0); - if (plotfile_type == "amrex") { - r2[0] = IntVect(1,1,ref_ratio[0][0]); - for (int lev = 1; lev <= finest_level; ++lev) { - if (lev > 1) { - r2[lev-1][0] = 1; - r2[lev-1][1] = 1; - r2[lev-1][2] = r2[lev-2][2] * ref_ratio[lev-1][0]; - } + // Copy level 0 as is + MultiFab::Copy(mf2[0],mf[0],0,0,mf[0].nComp(),0); - mf2[lev].define(refine(grids[lev],r2[lev-1]), dmap[lev], ncomp_mf, 0); + // Define a new multi-level array of Geometry's so that we pass the new "domain" at lev > 0 + Array periodicity = + {Geom()[0].isPeriodic(0),Geom()[0].isPeriodic(1),Geom()[0].isPeriodic(2)}; + g2[0].define(Geom()[0].Domain(),&(Geom()[0].ProbDomain()),0,periodicity.data()); - // Set the new problem domain - Box d2(Geom()[lev].Domain()); - d2.refine(r2[lev-1]); + r2[0] = IntVect(1,1,ref_ratio[0][0]); + for (int lev = 1; lev <= finest_level; ++lev) { + if (lev > 1) { + r2[lev-1][0] = 1; + r2[lev-1][1] = 1; + r2[lev-1][2] = r2[lev-2][2] * ref_ratio[lev-1][0]; + } - g2[lev].define(d2,&(Geom()[lev].ProbDomain()),0,periodicity.data()); - } + mf2[lev].define(refine(grids[lev],r2[lev-1]), dmap[lev], ncomp_mf, 0); - // Do piecewise interpolation of mf into mf2 - for (int lev = 1; lev <= finest_level; ++lev) { - Interpolater* mapper_c = &pc_interp; - InterpFromCoarseLevel(mf2[lev], t_new[lev], mf[lev], - 0, 0, ncomp_mf, - geom[lev], g2[lev], - null_bc_for_fill, 0, null_bc_for_fill, 0, - r2[lev-1], mapper_c, domain_bcs_type, 0); - } + // Set the new problem domain + Box d2(Geom()[lev].Domain()); + d2.refine(r2[lev-1]); - // Define an effective ref_ratio which is isotropic to be passed into WriteMultiLevelPlotfile - Vector rr(finest_level); - for (int lev = 0; lev < finest_level; ++lev) { - rr[lev] = IntVect(ref_ratio[lev][0],ref_ratio[lev][1],ref_ratio[lev][0]); - } + g2[lev].define(d2,&(Geom()[lev].ProbDomain()),0,periodicity.data()); + } - Print() << "Writing plotfile " << plotfilename << "\n"; - if (solverChoice.use_terrain) { - WriteMultiLevelPlotfileWithTerrain(plotfilename, finest_level+1, - GetVecOfConstPtrs(mf), - GetVecOfConstPtrs(mf_nd), - varnames, - t_new[0], istep); - } else { - WriteMultiLevelPlotfile(plotfilename, finest_level+1, - GetVecOfConstPtrs(mf2), varnames, - g2, t_new[0], istep, rr); - } + // Do piecewise interpolation of mf into mf2 + for (int lev = 1; lev <= finest_level; ++lev) { + Interpolater* mapper_c = &pc_interp; + InterpFromCoarseLevel(mf2[lev], t_new[lev], mf[lev], + 0, 0, ncomp_mf, + geom[lev], g2[lev], + null_bc_for_fill, 0, null_bc_for_fill, 0, + r2[lev-1], mapper_c, domain_bcs_type, 0); + } + + // Define an effective ref_ratio which is isotropic to be passed into WriteMultiLevelPlotfile + Vector rr(finest_level); + for (int lev = 0; lev < finest_level; ++lev) { + rr[lev] = IntVect(ref_ratio[lev][0],ref_ratio[lev][1],ref_ratio[lev][0]); + } + + Print() << "Writing plotfile " << plotfilename << "\n"; + if (solverChoice.use_terrain) { + WriteMultiLevelPlotfileWithTerrain(plotfilename, finest_level+1, + GetVecOfConstPtrs(mf), + GetVecOfConstPtrs(mf_nd), + varnames, + t_new[0], istep); + } else { + WriteMultiLevelPlotfile(plotfilename, finest_level+1, + GetVecOfConstPtrs(mf2), varnames, + g2, t_new[0], istep, rr); + } + + } else if (ref_ratio[0][2] != 1) { + if (solverChoice.use_terrain) { + WriteMultiLevelPlotfileWithTerrain(plotfilename, finest_level+1, + GetVecOfConstPtrs(mf), + GetVecOfConstPtrs(mf_nd), + varnames, + t_new[0], istep); + } else { + WriteMultiLevelPlotfile(plotfilename, finest_level+1, + GetVecOfConstPtrs(mf), varnames, + geom, t_new[0], istep, ref_ratio); + } + } // ref_ratio test writeJobInfo(plotfilename); #ifdef ERF_USE_PARTICLES particleData.Checkpoint(plotfilename); #endif + #ifdef ERF_USE_NETCDF } else if (plotfile_type == "netcdf" || plotfile_type == "NetCDF") { for (int lev = 0; lev <= finest_level; ++lev) { diff --git a/Source/Initialization/ERF_init1d.cpp b/Source/Initialization/ERF_init1d.cpp index c41ab360a..f86cfb99d 100644 --- a/Source/Initialization/ERF_init1d.cpp +++ b/Source/Initialization/ERF_init1d.cpp @@ -207,6 +207,7 @@ ERF::initHSE (int lev) MultiFab p_hse (base_state[lev], make_alias, 1, 1); // p_0 is second component MultiFab pi_hse(base_state[lev], make_alias, 2, 1); // pi_0 is third component +#if 0 if (lev > 0) { // Interp all three components: rho, p, pi @@ -222,6 +223,7 @@ ERF::initHSE (int lev) null_bc, 0, null_bc, 0, refRatio(lev-1), mapper, domain_bcs_type, bccomp); } +#endif // Initial r_hse may or may not be in HSE -- defined in prob.cpp if (solverChoice.use_moist_background){ From 584e73f9e63906d1c9d18f8c807f0a15fb99b4a6 Mon Sep 17 00:00:00 2001 From: "Aaron M. Lattanzi" <103702284+AMLattanzi@users.noreply.github.com> Date: Thu, 20 Jun 2024 12:17:39 -0700 Subject: [PATCH 2/9] Add average of Kturb for output. (#1652) --- Source/ERF.H | 1 + Source/IO/ERF_Write1DProfiles.cpp | 185 ++++++++++++++++-------------- 2 files changed, 98 insertions(+), 88 deletions(-) diff --git a/Source/ERF.H b/Source/ERF.H index a1ad20c05..d92977148 100644 --- a/Source/ERF.H +++ b/Source/ERF.H @@ -161,6 +161,7 @@ public: void derive_diag_profiles (amrex::Gpu::HostVector& h_avg_u , amrex::Gpu::HostVector& h_avg_v , amrex::Gpu::HostVector& h_avg_w , amrex::Gpu::HostVector& h_avg_rho, amrex::Gpu::HostVector& h_avg_th , amrex::Gpu::HostVector& h_avg_ksgs, + amrex::Gpu::HostVector& h_avg_kturb, amrex::Gpu::HostVector& h_avg_qv , amrex::Gpu::HostVector& h_avg_qc , amrex::Gpu::HostVector& h_avg_qr , amrex::Gpu::HostVector& h_avg_wqv , amrex::Gpu::HostVector& h_avg_wqc , diff --git a/Source/IO/ERF_Write1DProfiles.cpp b/Source/IO/ERF_Write1DProfiles.cpp index 4b6d9db3e..d8c4ec0de 100644 --- a/Source/IO/ERF_Write1DProfiles.cpp +++ b/Source/IO/ERF_Write1DProfiles.cpp @@ -27,7 +27,7 @@ ERF::write_1D_profiles (Real time) { // Define the 1d arrays we will need Gpu::HostVector h_avg_u, h_avg_v, h_avg_w; - Gpu::HostVector h_avg_rho, h_avg_th, h_avg_ksgs; + Gpu::HostVector h_avg_rho, h_avg_th, h_avg_ksgs, h_avg_kturb; Gpu::HostVector h_avg_qv, h_avg_qc, h_avg_qr, h_avg_wqv, h_avg_wqc, h_avg_wqr, h_avg_qi, h_avg_qs, h_avg_qg; Gpu::HostVector h_avg_uth, h_avg_vth, h_avg_wth, h_avg_thth; Gpu::HostVector h_avg_uu, h_avg_uv, h_avg_uw, h_avg_vv, h_avg_vw, h_avg_ww; @@ -38,7 +38,7 @@ ERF::write_1D_profiles (Real time) if (NumDataLogs() > 1) { derive_diag_profiles(h_avg_u, h_avg_v, h_avg_w, - h_avg_rho, h_avg_th, h_avg_ksgs, + h_avg_rho, h_avg_th, h_avg_ksgs, h_avg_kturb, h_avg_qv, h_avg_qc, h_avg_qr, h_avg_wqv, h_avg_wqc, h_avg_wqr, h_avg_qi, h_avg_qs, h_avg_qg, h_avg_uu, h_avg_uv, h_avg_uw, h_avg_vv, h_avg_vw, h_avg_ww, h_avg_uth, h_avg_vth, h_avg_wth, h_avg_thth, @@ -70,10 +70,11 @@ ERF::write_1D_profiles (Real time) } data_log1 << std::setw(datwidth) << std::setprecision(timeprecision) << time << " " << std::setw(datwidth) << std::setprecision(datprecision) << z << " " - << h_avg_u[k] << " " << h_avg_v[k] << " " << h_avg_w[k] << " " - << h_avg_rho[k] << " " << h_avg_th[k] << " " << h_avg_ksgs[k] << " " - << h_avg_qv[k] << " " << h_avg_qc[k] << " " << h_avg_qr[k] << " " - << h_avg_qi[k] << " " << h_avg_qs[k] << " " << h_avg_qg[k] + << h_avg_u[k] << " " << h_avg_v[k] << " " << h_avg_w[k] << " " + << h_avg_rho[k] << " " << h_avg_th[k] << " " << h_avg_ksgs[k] << " " + << h_avg_kturb[k] << " " << h_avg_qv[k] << " " << h_avg_qc[k] << " " + << h_avg_qr[k] << " " << h_avg_qi[k] << " " << h_avg_qs[k] << " " + << h_avg_qg[k] << std::endl; } // loop over z } // if good @@ -183,7 +184,7 @@ ERF::write_1D_profiles (Real time) void ERF::derive_diag_profiles(Gpu::HostVector& h_avg_u , Gpu::HostVector& h_avg_v , Gpu::HostVector& h_avg_w, Gpu::HostVector& h_avg_rho , Gpu::HostVector& h_avg_th , Gpu::HostVector& h_avg_ksgs, - Gpu::HostVector& h_avg_qv , Gpu::HostVector& h_avg_qc , Gpu::HostVector& h_avg_qr, + Gpu::HostVector& h_avg_kturb, Gpu::HostVector& h_avg_qv , Gpu::HostVector& h_avg_qc , Gpu::HostVector& h_avg_qr, Gpu::HostVector& h_avg_wqv , Gpu::HostVector& h_avg_wqc, Gpu::HostVector& h_avg_wqr, Gpu::HostVector& h_avg_qi , Gpu::HostVector& h_avg_qs , Gpu::HostVector& h_avg_qg, Gpu::HostVector& h_avg_uu , Gpu::HostVector& h_avg_uv , Gpu::HostVector& h_avg_uw, @@ -198,14 +199,15 @@ void ERF::derive_diag_profiles(Gpu::HostVector& h_avg_u , Gpu::HostVecto // We assume that this is always called at level 0 int lev = 0; - bool l_use_KE = (solverChoice.turbChoice[lev].les_type == LESType::Deardorff); - bool l_use_QKE = solverChoice.turbChoice[lev].use_QKE && solverChoice.turbChoice[lev].advect_QKE; + bool l_use_Turb = (solverChoice.turbChoice[lev].les_type != LESType::None); + bool l_use_KE = (solverChoice.turbChoice[lev].les_type == LESType::Deardorff); + bool l_use_QKE = solverChoice.turbChoice[lev].use_QKE && solverChoice.turbChoice[lev].advect_QKE; - // This will hold rho, theta, ksgs, uu, uv, uw, vv, vw, ww, uth, vth, wth, - // 0 1 2 3 4 5 6 7 8 9 10 11 + // This will hold rho, theta, ksgs, kturb, uu, uv, uw, vv, vw, ww, uth, vth, wth, + // 0 1 2 3 4 5 6 7 8 9 10 11 12 // thth, uiuiu, uiuiv, uiuiw, p, pu, pv, pw, qv, qc, qr, wqv, wqc, wqr, qi, qs, qg - // 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 - MultiFab mf_out(grids[lev], dmap[lev], 29, 0); + // 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 + MultiFab mf_out(grids[lev], dmap[lev], 30, 0); MultiFab mf_vels(grids[lev], dmap[lev], AMREX_SPACEDIM, 0); @@ -214,7 +216,7 @@ void ERF::derive_diag_profiles(Gpu::HostVector& h_avg_u , Gpu::HostVecto MultiFab w_cc(mf_vels, make_alias, 2, 1); // w at cell centers average_face_to_cellcenter(mf_vels,0, - Array{&vars_new[lev][Vars::xvel],&vars_new[lev][Vars::yvel],&vars_new[lev][Vars::zvel]}); + Array{&vars_new[lev][Vars::xvel],&vars_new[lev][Vars::yvel],&vars_new[lev][Vars::zvel]}); int zdir = 2; auto domain = geom[0].Domain(); @@ -248,6 +250,8 @@ void ERF::derive_diag_profiles(Gpu::HostVector& h_avg_u , Gpu::HostVecto const Array4& w_cc_arr = w_cc.array(mfi); const Array4& cons_arr = mf_cons.array(mfi); const Array4& p0_arr = p_hse.array(mfi); + const Array4& eta_arr = (l_use_Turb) ? eddyDiffs_lev[lev]->const_array(mfi) : + Array4{}; ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { @@ -255,42 +259,46 @@ void ERF::derive_diag_profiles(Gpu::HostVector& h_avg_u , Gpu::HostVecto fab_arr(i, j, k, 0) = cons_arr(i,j,k,Rho_comp); fab_arr(i, j, k, 1) = theta; Real ksgs = 0.0; - if (l_use_KE) + if (l_use_KE) { ksgs = cons_arr(i,j,k,RhoKE_comp) / cons_arr(i,j,k,Rho_comp); - else if (l_use_QKE) + } else if (l_use_QKE) { ksgs = cons_arr(i,j,k,RhoQKE_comp) / cons_arr(i,j,k,Rho_comp); + } fab_arr(i, j, k, 2) = ksgs; - fab_arr(i, j, k, 3) = u_cc_arr(i,j,k) * u_cc_arr(i,j,k); // u*u - fab_arr(i, j, k, 4) = u_cc_arr(i,j,k) * v_cc_arr(i,j,k); // u*v - fab_arr(i, j, k, 5) = u_cc_arr(i,j,k) * w_cc_arr(i,j,k); // u*w - fab_arr(i, j, k, 6) = v_cc_arr(i,j,k) * v_cc_arr(i,j,k); // v*v - fab_arr(i, j, k, 7) = v_cc_arr(i,j,k) * w_cc_arr(i,j,k); // v*w - fab_arr(i, j, k, 8) = w_cc_arr(i,j,k) * w_cc_arr(i,j,k); // w*w - fab_arr(i, j, k, 9) = u_cc_arr(i,j,k) * theta; // u*th - fab_arr(i, j, k,10) = v_cc_arr(i,j,k) * theta; // v*th - fab_arr(i, j, k,11) = w_cc_arr(i,j,k) * theta; // w*th - fab_arr(i, j, k,12) = theta * theta; // th*th + Real kturb = 0.0; + if (l_use_Turb) kturb = eta_arr(i,j,k,EddyDiff::Mom_h); + fab_arr(i, j, k, 3) = kturb; + fab_arr(i, j, k, 4) = u_cc_arr(i,j,k) * u_cc_arr(i,j,k); // u*u + fab_arr(i, j, k, 5) = u_cc_arr(i,j,k) * v_cc_arr(i,j,k); // u*v + fab_arr(i, j, k, 6) = u_cc_arr(i,j,k) * w_cc_arr(i,j,k); // u*w + fab_arr(i, j, k, 7) = v_cc_arr(i,j,k) * v_cc_arr(i,j,k); // v*v + fab_arr(i, j, k, 8) = v_cc_arr(i,j,k) * w_cc_arr(i,j,k); // v*w + fab_arr(i, j, k, 9) = w_cc_arr(i,j,k) * w_cc_arr(i,j,k); // w*w + fab_arr(i, j, k,10) = u_cc_arr(i,j,k) * theta; // u*th + fab_arr(i, j, k,11) = v_cc_arr(i,j,k) * theta; // v*th + fab_arr(i, j, k,12) = w_cc_arr(i,j,k) * theta; // w*th + fab_arr(i, j, k,13) = theta * theta; // th*th Real uiui = fab_arr(i,j,k,3) + fab_arr(i,j,k,6) + fab_arr(i,j,k,8); - fab_arr(i, j, k,13) = uiui * u_cc_arr(i,j,k); // (ui*ui)*u - fab_arr(i, j, k,14) = uiui * v_cc_arr(i,j,k); // (ui*ui)*v - fab_arr(i, j, k,15) = uiui * w_cc_arr(i,j,k); // (ui*ui)*w + fab_arr(i, j, k,14) = uiui * u_cc_arr(i,j,k); // (ui*ui)*u + fab_arr(i, j, k,15) = uiui * v_cc_arr(i,j,k); // (ui*ui)*v + fab_arr(i, j, k,16) = uiui * w_cc_arr(i,j,k); // (ui*ui)*w if (!use_moisture) { Real p = getPgivenRTh(cons_arr(i, j, k, RhoTheta_comp)); p -= p0_arr(i,j,k); - fab_arr(i, j, k,16) = p; // p' - fab_arr(i, j, k,17) = p * u_cc_arr(i,j,k); // p'*u - fab_arr(i, j, k,18) = p * v_cc_arr(i,j,k); // p'*v - fab_arr(i, j, k,19) = p * w_cc_arr(i,j,k); // p'*w - fab_arr(i, j, k,20) = 0.; // qv - fab_arr(i, j, k,21) = 0.; // qc - fab_arr(i, j, k,22) = 0.; // qr - fab_arr(i, j, k,23) = 0.; // w*qv - fab_arr(i, j, k,24) = 0.; // w*qc - fab_arr(i, j, k,25) = 0.; // w*qr - fab_arr(i, j, k,26) = 0.; // qi - fab_arr(i, j, k,27) = 0.; // qs - fab_arr(i, j, k,28) = 0.; // qg + fab_arr(i, j, k,17) = p; // p' + fab_arr(i, j, k,18) = p * u_cc_arr(i,j,k); // p'*u + fab_arr(i, j, k,19) = p * v_cc_arr(i,j,k); // p'*v + fab_arr(i, j, k,20) = p * w_cc_arr(i,j,k); // p'*w + fab_arr(i, j, k,21) = 0.; // qv + fab_arr(i, j, k,22) = 0.; // qc + fab_arr(i, j, k,23) = 0.; // qr + fab_arr(i, j, k,24) = 0.; // w*qv + fab_arr(i, j, k,25) = 0.; // w*qc + fab_arr(i, j, k,26) = 0.; // w*qr + fab_arr(i, j, k,27) = 0.; // qi + fab_arr(i, j, k,28) = 0.; // qs + fab_arr(i, j, k,29) = 0.; // qg } }); } // mfi @@ -321,63 +329,64 @@ void ERF::derive_diag_profiles(Gpu::HostVector& h_avg_u , Gpu::HostVecto Real p = getPgivenRTh(cons_arr(i, j, k, RhoTheta_comp), qv_arr(i,j,k)); p -= p0_arr(i,j,k); - fab_arr(i, j, k,16) = p; // p' - fab_arr(i, j, k,17) = p * u_cc_arr(i,j,k); // p'*u - fab_arr(i, j, k,18) = p * v_cc_arr(i,j,k); // p'*v - fab_arr(i, j, k,19) = p * w_cc_arr(i,j,k); // p'*w - fab_arr(i, j, k,20) = cons_arr(i,j,k,RhoQ1_comp) / cons_arr(i,j,k,Rho_comp); // qv - fab_arr(i, j, k,21) = cons_arr(i,j,k,RhoQ2_comp) / cons_arr(i,j,k,Rho_comp); // qc - fab_arr(i, j, k,22) = cons_arr(i,j,k,RhoQr_comp) / cons_arr(i,j,k,Rho_comp); // qr - fab_arr(i, j, k,23) = w_cc_arr(i,j,k) * cons_arr(i,j,k,RhoQ1_comp) / cons_arr(i,j,k,Rho_comp); // w*qv - fab_arr(i, j, k,24) = w_cc_arr(i,j,k) * cons_arr(i,j,k,RhoQ2_comp) / cons_arr(i,j,k,Rho_comp); // w*qc - fab_arr(i, j, k,25) = w_cc_arr(i,j,k) * cons_arr(i,j,k,RhoQr_comp) / cons_arr(i,j,k,Rho_comp); // w*qr + fab_arr(i, j, k,17) = p; // p' + fab_arr(i, j, k,18) = p * u_cc_arr(i,j,k); // p'*u + fab_arr(i, j, k,19) = p * v_cc_arr(i,j,k); // p'*v + fab_arr(i, j, k,20) = p * w_cc_arr(i,j,k); // p'*w + fab_arr(i, j, k,21) = cons_arr(i,j,k,RhoQ1_comp) / cons_arr(i,j,k,Rho_comp); // qv + fab_arr(i, j, k,22) = cons_arr(i,j,k,RhoQ2_comp) / cons_arr(i,j,k,Rho_comp); // qc + fab_arr(i, j, k,23) = cons_arr(i,j,k,RhoQr_comp) / cons_arr(i,j,k,Rho_comp); // qr + fab_arr(i, j, k,24) = w_cc_arr(i,j,k) * cons_arr(i,j,k,RhoQ1_comp) / cons_arr(i,j,k,Rho_comp); // w*qv + fab_arr(i, j, k,25) = w_cc_arr(i,j,k) * cons_arr(i,j,k,RhoQ2_comp) / cons_arr(i,j,k,Rho_comp); // w*qc + fab_arr(i, j, k,26) = w_cc_arr(i,j,k) * cons_arr(i,j,k,RhoQr_comp) / cons_arr(i,j,k,Rho_comp); // w*qr if (n_qstate > 3) { - fab_arr(i, j, k,26) = cons_arr(i,j,k,RhoQ3_comp) / cons_arr(i,j,k,Rho_comp); // qi - fab_arr(i, j, k,27) = cons_arr(i,j,k,RhoQ5_comp) / cons_arr(i,j,k,Rho_comp); // qs - fab_arr(i, j, k,28) = cons_arr(i,j,k,RhoQ6_comp) / cons_arr(i,j,k,Rho_comp); // qg + fab_arr(i, j, k,27) = cons_arr(i,j,k,RhoQ3_comp) / cons_arr(i,j,k,Rho_comp); // qi + fab_arr(i, j, k,28) = cons_arr(i,j,k,RhoQ5_comp) / cons_arr(i,j,k,Rho_comp); // qs + fab_arr(i, j, k,29) = cons_arr(i,j,k,RhoQ6_comp) / cons_arr(i,j,k,Rho_comp); // qg } else { - fab_arr(i, j, k,26) = 0.0; // qi - fab_arr(i, j, k,27) = 0.0; // qs - fab_arr(i, j, k,28) = 0.0; // qg + fab_arr(i, j, k,27) = 0.0; // qi + fab_arr(i, j, k,28) = 0.0; // qs + fab_arr(i, j, k,29) = 0.0; // qg } }); } // mfi } // use_moisture - h_avg_rho = sumToLine(mf_out, 0,1,domain,zdir); - h_avg_th = sumToLine(mf_out, 1,1,domain,zdir); - h_avg_ksgs = sumToLine(mf_out, 2,1,domain,zdir); - h_avg_uu = sumToLine(mf_out, 3,1,domain,zdir); - h_avg_uv = sumToLine(mf_out, 4,1,domain,zdir); - h_avg_uw = sumToLine(mf_out, 5,1,domain,zdir); - h_avg_vv = sumToLine(mf_out, 6,1,domain,zdir); - h_avg_vw = sumToLine(mf_out, 7,1,domain,zdir); - h_avg_ww = sumToLine(mf_out, 8,1,domain,zdir); - h_avg_uth = sumToLine(mf_out, 9,1,domain,zdir); - h_avg_vth = sumToLine(mf_out,10,1,domain,zdir); - h_avg_wth = sumToLine(mf_out,11,1,domain,zdir); - h_avg_thth = sumToLine(mf_out,12,1,domain,zdir); - h_avg_uiuiu= sumToLine(mf_out,13,1,domain,zdir); - h_avg_uiuiv= sumToLine(mf_out,14,1,domain,zdir); - h_avg_uiuiw= sumToLine(mf_out,15,1,domain,zdir); - h_avg_p = sumToLine(mf_out,16,1,domain,zdir); - h_avg_pu = sumToLine(mf_out,17,1,domain,zdir); - h_avg_pv = sumToLine(mf_out,18,1,domain,zdir); - h_avg_pw = sumToLine(mf_out,19,1,domain,zdir); - h_avg_qv = sumToLine(mf_out,20,1,domain,zdir); - h_avg_qc = sumToLine(mf_out,21,1,domain,zdir); - h_avg_qr = sumToLine(mf_out,22,1,domain,zdir); - h_avg_wqv = sumToLine(mf_out,23,1,domain,zdir); - h_avg_wqc = sumToLine(mf_out,24,1,domain,zdir); - h_avg_wqr = sumToLine(mf_out,25,1,domain,zdir); - h_avg_qi = sumToLine(mf_out,26,1,domain,zdir); - h_avg_qs = sumToLine(mf_out,27,1,domain,zdir); - h_avg_qg = sumToLine(mf_out,28,1,domain,zdir); + h_avg_rho = sumToLine(mf_out, 0,1,domain,zdir); + h_avg_th = sumToLine(mf_out, 1,1,domain,zdir); + h_avg_ksgs = sumToLine(mf_out, 2,1,domain,zdir); + h_avg_kturb = sumToLine(mf_out, 3,1,domain,zdir); + h_avg_uu = sumToLine(mf_out, 4,1,domain,zdir); + h_avg_uv = sumToLine(mf_out, 5,1,domain,zdir); + h_avg_uw = sumToLine(mf_out, 6,1,domain,zdir); + h_avg_vv = sumToLine(mf_out, 7,1,domain,zdir); + h_avg_vw = sumToLine(mf_out, 8,1,domain,zdir); + h_avg_ww = sumToLine(mf_out, 9,1,domain,zdir); + h_avg_uth = sumToLine(mf_out,10,1,domain,zdir); + h_avg_vth = sumToLine(mf_out,11,1,domain,zdir); + h_avg_wth = sumToLine(mf_out,12,1,domain,zdir); + h_avg_thth = sumToLine(mf_out,13,1,domain,zdir); + h_avg_uiuiu = sumToLine(mf_out,14,1,domain,zdir); + h_avg_uiuiv = sumToLine(mf_out,15,1,domain,zdir); + h_avg_uiuiw = sumToLine(mf_out,16,1,domain,zdir); + h_avg_p = sumToLine(mf_out,17,1,domain,zdir); + h_avg_pu = sumToLine(mf_out,18,1,domain,zdir); + h_avg_pv = sumToLine(mf_out,19,1,domain,zdir); + h_avg_pw = sumToLine(mf_out,20,1,domain,zdir); + h_avg_qv = sumToLine(mf_out,21,1,domain,zdir); + h_avg_qc = sumToLine(mf_out,22,1,domain,zdir); + h_avg_qr = sumToLine(mf_out,23,1,domain,zdir); + h_avg_wqv = sumToLine(mf_out,24,1,domain,zdir); + h_avg_wqc = sumToLine(mf_out,25,1,domain,zdir); + h_avg_wqr = sumToLine(mf_out,26,1,domain,zdir); + h_avg_qi = sumToLine(mf_out,27,1,domain,zdir); + h_avg_qs = sumToLine(mf_out,28,1,domain,zdir); + h_avg_qg = sumToLine(mf_out,29,1,domain,zdir); // Divide by the total number of cells we are averaging over int h_avg_u_size = static_cast(h_avg_u.size()); for (int k = 0; k < h_avg_u_size; ++k) { - h_avg_rho[k] /= area_z; h_avg_ksgs[k] /= area_z; + h_avg_rho[k] /= area_z; h_avg_ksgs[k] /= area_z; h_avg_kturb[k] /= area_z; h_avg_th[k] /= area_z; h_avg_thth[k] /= area_z; h_avg_uu[k] /= area_z; h_avg_uv[k] /= area_z; h_avg_uw[k] /= area_z; h_avg_vv[k] /= area_z; h_avg_vw[k] /= area_z; h_avg_ww[k] /= area_z; From 2851437de54787174878d0923b2c2c9f33d44fcf Mon Sep 17 00:00:00 2001 From: japham0 <156237070+japham0@users.noreply.github.com> Date: Thu, 20 Jun 2024 19:14:24 -0700 Subject: [PATCH 3/9] first pass at coupling with WW3 (#1653) * first pass at coupling with WW3 * Move the ifdef so we don't have to worry about finding mpi.h --------- Co-authored-by: AMLattanzi --- CMake/BuildERFExe.cmake | 1 + Exec/ABL/inputs_mpmd | 68 +++++++++++++++ Exec/Make.ERF | 4 + Source/ERF.H | 7 ++ Source/ERF.cpp | 13 ++- Source/ERF_make_new_arrays.cpp | 26 ++++++ Source/ERF_read_waves.cpp | 107 ++++++++++++++++++++++++ Source/Make.package | 4 + Source/TimeIntegration/ERF_TimeStep.cpp | 4 + Source/main.cpp | 19 ++++- 10 files changed, 251 insertions(+), 2 deletions(-) create mode 100644 Exec/ABL/inputs_mpmd create mode 100644 Source/ERF_read_waves.cpp diff --git a/CMake/BuildERFExe.cmake b/CMake/BuildERFExe.cmake index 412c1d926..5f76b5d72 100644 --- a/CMake/BuildERFExe.cmake +++ b/CMake/BuildERFExe.cmake @@ -107,6 +107,7 @@ function(build_erf_lib erf_lib_name) ${SRC_DIR}/ERF.cpp ${SRC_DIR}/ERF_make_new_arrays.cpp ${SRC_DIR}/ERF_make_new_level.cpp + ${SRC_DIR}/ERF_read_waves.cpp ${SRC_DIR}/ERF_Tagging.cpp ${SRC_DIR}/Advection/AdvectionSrcForMom.cpp ${SRC_DIR}/Advection/AdvectionSrcForState.cpp diff --git a/Exec/ABL/inputs_mpmd b/Exec/ABL/inputs_mpmd new file mode 100644 index 000000000..d03df316a --- /dev/null +++ b/Exec/ABL/inputs_mpmd @@ -0,0 +1,68 @@ +# ------------------ INPUTS TO MAIN PROGRAM ------------------- +max_step = 5 + +amrex.fpe_trap_invalid = 0 + +fabarray.mfiter_tile_size = 1024 1024 1024 + +# PROBLEM SIZE & GEOMETRY +geometry.prob_extent = 191000 91000 1024 +amr.n_cell = 191 91 4 + +geometry.is_periodic = 1 1 0 + +zlo.type = "NoSlipWall" +zhi.type = "SlipWall" + +# TIME STEP CONTROL +erf.fixed_dt = 10 # fixed time step depending on grid resolution +erf.max_step=1 +# 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 = 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 rhoadv_0 x_velocity y_velocity z_velocity pressure temp theta + +# SOLVER CHOICE +erf.alpha_T = 0.0 +erf.alpha_C = 1.0 +erf.use_gravity = false + +erf.molec_diff_type = "None" +erf.les_type = "Smagorinsky" +erf.Cs = 0.1 + +erf.init_type = "uniform" + +# MOST BOUNDARY (DEFAULT IS ADIABATIC FOR THETA) +zlo.type = "Most" +# erf.most.z0 = 0.1 Don't set z0 +erf.most.zref = 128.0 # set 1/2 of the vertical dz: 1024 / 4 = 256 -> 256/2 = 128 +erf.most.roughness_type_sea = wave_coupled + +# PROBLEM PARAMETERS +prob.rho_0 = 1.0 +prob.A_0 = 1.0 + +prob.U_0 = 10.0 +prob.V_0 = 0.0 +prob.W_0 = 0.0 +prob.T_0 = 300.0 + +# Higher values of perturbations lead to instability +# Instability seems to be coming from BC +prob.U_0_Pert_Mag = 0.08 +prob.V_0_Pert_Mag = 0.08 # +prob.W_0_Pert_Mag = 0.0 + diff --git a/Exec/Make.ERF b/Exec/Make.ERF index 7efb8b1fa..8f81e8092 100644 --- a/Exec/Make.ERF +++ b/Exec/Make.ERF @@ -155,6 +155,10 @@ ifeq ($(USE_WINDFARM), TRUE) INCLUDE_LOCATIONS += $(ERF_WINDFARM_EWP_DIR) endif +ifeq ($(USE_WW3_COUPLING), TRUE) + DEFINES += -DERF_USE_WW3_COUPLING +endif + ERF_LSM_DIR = $(ERF_SOURCE_DIR)/LandSurfaceModel include $(ERF_LSM_DIR)/Make.package VPATH_LOCATIONS += $(ERF_LSM_DIR) diff --git a/Source/ERF.H b/Source/ERF.H index d92977148..4cb59e790 100644 --- a/Source/ERF.H +++ b/Source/ERF.H @@ -238,6 +238,10 @@ public: // compute dt from CFL considerations amrex::Real estTimeStep (int lev, long& dt_fast_ratio) const; +#ifdef ERF_USE_WW3_COUPLING + void read_waves (int lev); +#endif + // Interface for advancing the data at one level by one "slow" timestep void advance_dycore (int level, amrex::Vector& state_old, @@ -716,6 +720,9 @@ private: // Wave coupling data amrex::Vector> Hwave; amrex::Vector> Lwave; + amrex::Vector> Hwave_onegrid; + amrex::Vector> Lwave_onegrid; + bool finished_wave = false; // array of flux registers amrex::Vector advflux_reg; diff --git a/Source/ERF.cpp b/Source/ERF.cpp index a884317f8..26faab742 100644 --- a/Source/ERF.cpp +++ b/Source/ERF.cpp @@ -276,7 +276,13 @@ ERF::ERF () Hwave[lev] = nullptr; Lwave[lev] = nullptr; } - + Hwave_onegrid.resize(nlevs_max); + Lwave_onegrid.resize(nlevs_max); + for (int lev = 0; lev < max_level; ++lev) + { + Hwave_onegrid[lev] = nullptr; + Lwave_onegrid[lev] = nullptr; + } // Theta prim for MOST Theta_prim.resize(nlevs_max); @@ -900,6 +906,11 @@ ERF::InitData () } } +#ifdef ERF_USE_WW3_COUPLING + int lev = 0; + read_waves(lev); +#endif + // Configure ABLMost params if used MostWall boundary condition // NOTE: we must set up the MOST routine after calling FillPatch // in order to have lateral ghost cells filled (MOST + terrain interp). diff --git a/Source/ERF_make_new_arrays.cpp b/Source/ERF_make_new_arrays.cpp index 08ccf435a..7e2e5d561 100644 --- a/Source/ERF_make_new_arrays.cpp +++ b/Source/ERF_make_new_arrays.cpp @@ -245,6 +245,32 @@ ERF::init_stuff (int lev, const BoxArray& ba, const DistributionMapping& dm, } #endif + +#ifdef ERF_USE_WW3_COUPLING + // create a new BoxArray and DistributionMapping for a MultiFab with 1 box + BoxArray ba_onegrid(geom[lev].Domain()); + BoxList bl2d_onegrid = ba_onegrid.boxList(); + for (auto& b : bl2d_onegrid) { + b.setRange(2,0); + } + BoxArray ba2d_onegrid(std::move(bl2d_onegrid)); + Vector pmap; + pmap.resize(1); + pmap[0]=0; + DistributionMapping dm_onegrid(ba2d_onegrid); + dm_onegrid.define(pmap); + + Hwave_onegrid[lev] = std::make_unique(ba2d_onegrid,dm_onegrid,1,IntVect(1,1,0)); + Lwave_onegrid[lev] = std::make_unique(ba2d_onegrid,dm_onegrid,1,IntVect(1,1,0)); + Hwave[lev] = std::make_unique(ba2d,dm,1,IntVect(3,3,0)); + Lwave[lev] = std::make_unique(ba2d,dm,1,IntVect(3,3,0)); + std::cout< +#include +#include +#include + +using namespace amrex; + +void +ERF::read_waves (int lev) +{ + for ( MFIter mfi(*Hwave_onegrid[lev],false); mfi.isValid(); ++mfi) + { + //auto my_H_ptr = Hwave[lev]->array(mfi); + //auto my_L_ptr = Lwave[lev]->array(mfi); + const auto & bx = mfi.validbox(); + + amrex::Print() << " HERE " << bx << std::endl; + // How to declare my_H_ptr directly? + amrex::Array4 my_H_arr = Hwave_onegrid[lev]->array(mfi); + amrex::Array4 my_L_arr = Lwave_onegrid[lev]->array(mfi); + + Real* my_H_ptr = my_H_arr.dataPtr(); + Real* my_L_ptr = my_L_arr.dataPtr(); + + int rank_offset = amrex::MPMD::MyProc() - amrex::ParallelDescriptor::MyProc(); + int this_root, other_root; + if (rank_offset == 0) { // First program + this_root = 0; + other_root = amrex::ParallelDescriptor::NProcs(); + } else { + this_root = rank_offset; + other_root = 0; + } + + + int nx=2147483647; // sanity check + int ny=2147483647; // sanity check + + //JUST RECV + if (amrex::MPMD::MyProc() == this_root) { + if (rank_offset == 0) // the first program + { + MPI_Recv(&nx, 1, MPI_INT, other_root, 1, MPI_COMM_WORLD, MPI_STATUS_IGNORE); + MPI_Recv(&ny, 1, MPI_INT, other_root, 7, MPI_COMM_WORLD, MPI_STATUS_IGNORE); + } + else // the second program + { + MPI_Recv(&nx, 1, MPI_INT, other_root, 0, MPI_COMM_WORLD, MPI_STATUS_IGNORE); + MPI_Recv(&ny, 1, MPI_INT, other_root, 6, MPI_COMM_WORLD, MPI_STATUS_IGNORE); + } + //This may not be necessary + ParallelDescriptor::Bcast(&nx, 1); + ParallelDescriptor::Bcast(&ny, 1); + } + if((nx)*(ny) > 0) { + int nsealm = (nx)*ny; + Print()<ParallelCopy(*Hwave_onegrid[lev]); + Hwave[lev]->FillBoundary(geom[lev].periodicity()); + Lwave[lev]->ParallelCopy(*Lwave_onegrid[lev]); + Lwave[lev]->FillBoundary(geom[lev].periodicity()); + amrex::Print() << "HWAVE BOX " << (*Hwave[lev])[0].box() << std::endl; + + print_state(*Hwave[lev],IntVect(103,-3,0),0,IntVect(3,3,0)); + print_state(*Hwave[lev],IntVect(103,88,0),0,IntVect(3,3,0)); +} +#endif + diff --git a/Source/Make.package b/Source/Make.package index 119eb03f0..b62909334 100644 --- a/Source/Make.package +++ b/Source/Make.package @@ -14,3 +14,7 @@ CEXE_sources += ERF_make_new_level.cpp CEXE_sources += ERF_make_new_arrays.cpp CEXE_sources += Derive.cpp CEXE_headers += Derive.H + +ifeq ($(USE_WW3_COUPLING),TRUE) +CEXE_sources += ERF_read_waves.cpp +endif diff --git a/Source/TimeIntegration/ERF_TimeStep.cpp b/Source/TimeIntegration/ERF_TimeStep.cpp index 2212ffcbf..4084498b7 100644 --- a/Source/TimeIntegration/ERF_TimeStep.cpp +++ b/Source/TimeIntegration/ERF_TimeStep.cpp @@ -63,6 +63,10 @@ ERF::timeStep (int lev, Real time, int /*iteration*/) << " with dt = " << dt[lev] << std::endl; } +#ifdef ERF_USE_WW3_COUPLING + read_waves(lev); +#endif + // Advance a single level for a single time step Advance(lev, time, dt[lev], istep[lev], nsubsteps[lev]); diff --git a/Source/main.cpp b/Source/main.cpp index a79f97c75..2f473b85f 100644 --- a/Source/main.cpp +++ b/Source/main.cpp @@ -11,6 +11,11 @@ #include #endif +#ifdef ERF_USE_WW3_COUPLING +#include +#include +#endif + std::string inputs_name; using namespace amrex; @@ -55,6 +60,7 @@ void add_par () { */ int main (int argc, char* argv[]) { + #ifdef AMREX_USE_MPI MPI_Init(&argc, &argv); #endif @@ -96,7 +102,12 @@ int main (int argc, char* argv[]) } } } +#ifdef ERF_USE_WW3_COUPLING + MPI_Comm comm = amrex::MPMD::Initialize(argc, argv); + amrex::Initialize(argc,argv,true,comm,add_par); +#else amrex::Initialize(argc,argv,true,MPI_COMM_WORLD,add_par); +#endif // Save the inputs file name for later. if (!strchr(argv[1], '=')) { @@ -228,9 +239,15 @@ int main (int argc, char* argv[]) // destroy timer for profiling BL_PROFILE_VAR_STOP(pmain); - +#ifdef ERF_USE_WW3_COUPLING + MPI_Barrier(MPI_COMM_WORLD); +#endif amrex::Finalize(); #ifdef AMREX_USE_MPI +#ifdef ERF_USE_WW3_COUPLING + amrex::MPMD::Finalize(); +#else MPI_Finalize(); #endif +#endif } From 966cdb7910261e077f82cd8f511eb13c42f81ca9 Mon Sep 17 00:00:00 2001 From: Weiqun Zhang Date: Thu, 20 Jun 2024 21:56:12 -0700 Subject: [PATCH 4/9] MRISplitIntegrator no longer derives from amrex::IntegratorBase. (#1655) --- Source/TimeIntegration/ERF_MRI.H | 12 +++++------- 1 file changed, 5 insertions(+), 7 deletions(-) diff --git a/Source/TimeIntegration/ERF_MRI.H b/Source/TimeIntegration/ERF_MRI.H index e0e33afb1..48225bf29 100644 --- a/Source/TimeIntegration/ERF_MRI.H +++ b/Source/TimeIntegration/ERF_MRI.H @@ -12,7 +12,7 @@ #include template -class MRISplitIntegrator : public amrex::IntegratorBase +class MRISplitIntegrator { private: /** @@ -90,12 +90,12 @@ public: initialize_data(S_data); } - void initialize (const T& S_data) override + void initialize (const T& S_data) { initialize_data(S_data); } - ~MRISplitIntegrator () override = default; + ~MRISplitIntegrator () = default; // Declare a default move constructor so we ensure the destructor is // not called when we return an object of this class by value @@ -184,7 +184,7 @@ public: return rhs; } - amrex::Real advance (T& S_old, T& S_new, amrex::Real time, const amrex::Real time_step) override + amrex::Real advance (T& S_old, T& S_new, amrex::Real time, const amrex::Real time_step) { BL_PROFILE_REGION("MRI_advance"); using namespace amrex; @@ -348,9 +348,7 @@ public: return timestep; } - void time_interpolate (const T& /* S_new */, const T& /*S_old*/, amrex::Real /*timestep_fraction*/, T& /*data*/) override {} - - void map_data (std::function Map) override + void map_data (std::function Map) { for (auto& F : T_store) { Map(*F); From 5ce79bb568662a1a4b2898b7074869f2fddf98bc Mon Sep 17 00:00:00 2001 From: japham0 <156237070+japham0@users.noreply.github.com> Date: Fri, 21 Jun 2024 16:51:13 -0700 Subject: [PATCH 5/9] Updates to read land/sea mask (#1656) * changes to read land/sea mask * changes to read land/sea mask --- Source/ERF_read_waves.cpp | 51 ++++++++++++++++----------------------- 1 file changed, 21 insertions(+), 30 deletions(-) diff --git a/Source/ERF_read_waves.cpp b/Source/ERF_read_waves.cpp index 9af9410b4..ec9032992 100644 --- a/Source/ERF_read_waves.cpp +++ b/Source/ERF_read_waves.cpp @@ -12,12 +12,10 @@ ERF::read_waves (int lev) { for ( MFIter mfi(*Hwave_onegrid[lev],false); mfi.isValid(); ++mfi) { - //auto my_H_ptr = Hwave[lev]->array(mfi); - //auto my_L_ptr = Lwave[lev]->array(mfi); + const auto & bx = mfi.validbox(); - amrex::Print() << " HERE " << bx << std::endl; - // How to declare my_H_ptr directly? + amrex::Print() << " HERE " << bx << std::endl; amrex::Array4 my_H_arr = Hwave_onegrid[lev]->array(mfi); amrex::Array4 my_L_arr = Lwave_onegrid[lev]->array(mfi); @@ -35,17 +33,17 @@ ERF::read_waves (int lev) } - int nx=2147483647; // sanity check + int nx=2147483647; int ny=2147483647; // sanity check - //JUST RECV + //JUST RECEIVED if (amrex::MPMD::MyProc() == this_root) { - if (rank_offset == 0) // the first program + if (rank_offset == 0) // First program { MPI_Recv(&nx, 1, MPI_INT, other_root, 1, MPI_COMM_WORLD, MPI_STATUS_IGNORE); MPI_Recv(&ny, 1, MPI_INT, other_root, 7, MPI_COMM_WORLD, MPI_STATUS_IGNORE); } - else // the second program + else // Second program { MPI_Recv(&nx, 1, MPI_INT, other_root, 0, MPI_COMM_WORLD, MPI_STATUS_IGNORE); MPI_Recv(&ny, 1, MPI_INT, other_root, 6, MPI_COMM_WORLD, MPI_STATUS_IGNORE); @@ -54,17 +52,9 @@ ERF::read_waves (int lev) ParallelDescriptor::Bcast(&nx, 1); ParallelDescriptor::Bcast(&ny, 1); } + if((nx)*(ny) > 0) { int nsealm = (nx)*ny; - Print()<ParallelCopy(*Hwave_onegrid[lev]); @@ -99,9 +80,19 @@ ERF::read_waves (int lev) Lwave[lev]->ParallelCopy(*Lwave_onegrid[lev]); Lwave[lev]->FillBoundary(geom[lev].periodicity()); amrex::Print() << "HWAVE BOX " << (*Hwave[lev])[0].box() << std::endl; + for (MFIter mfi(*Hwave[lev],TilingIfNotGPU()); mfi.isValid(); ++mfi) { + Box bx = mfi.tilebox(); + const Array4& Hwave_arr = Hwave[lev]->const_array(mfi); + const Array4& Lmask_arr = lmask_lev[lev][0]->array(mfi); + ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k){ + if (Hwave_arr(i,j,k)<0) { + Lmask_arr(i,j,k) = 1; + } else { + Lmask_arr(i,j,k) = 0; + } + }); + } - print_state(*Hwave[lev],IntVect(103,-3,0),0,IntVect(3,3,0)); - print_state(*Hwave[lev],IntVect(103,88,0),0,IntVect(3,3,0)); } #endif From a0ea7699865abb63bbe44bc504ba444e8d131b76 Mon Sep 17 00:00:00 2001 From: "Dustin (Ting-Hsuan) Ma" Date: Mon, 24 Jun 2024 09:25:09 -0700 Subject: [PATCH 6/9] Turbulent inflow generation technique (#1657) * Box perturbation setup start * perturbation box implementation start * added geometry into CalcTurbPert() * Using index of box to created boxes * Index Perturbation box initialization working * added Perturbation box initialization in init_stuff(), next create class for PerturbationRegion * created 3 functions init_PerturbationRegion(), calc_TurbPert_amplitude(), and calc_TurbPert_updateTime. Now PB boxes initialize at t=0s * cleaned up some portions of the code. Now start implemented t>0s PB * start of struct TurbPert implementation * created TurbPert structure and compiles/runs fines. Now t>=0 uses same function to create turbulent perturbations. * name change to TurbulentPerturbation structure * made small modifications to perturbation funcion call. Start playing with momentum perturbation. * Averaging paraview out put through impermeability enforcement. Added in random number generator in Perturbation amplitude. * Cleaned up some code, and added comments for tomorrow start. Tomorrow: Implement update-time and amplitude * distribution mapping debug w/ Aaron * debugging mf global access * mf now is seen outside of struct * moving on PB implementation to using vectors for per box unique value. * Replaced all perturbation boxes from boxArray and vector to MultiFab. Current version works with np=1, but not np>1. * fixed parallelization problem. Succesfully accessed velocity at initialization. Moving on to calculate update time/think about how to calculate perturbation amplitude on the fly. * Created elapsed time array for mf in tpi. Start writing perturbation amplitude calculations. * mirgrated over to ssh github check * PB u_mean calc per bx using box union (volume averaged). PB items uses vectors. PB frequency triggeres amplitude calculation. PB apply_tpi applied random number perturbation. * Added in Ma&senocak2023 PB amp formulation. Made PB velocity average thread safe. Created some initial constant values for PB. * debugging ubx output (start) * Bug fixed: object.makeSlab() will modify the object, while makeSlab(object,...) will make a copy then edit copy. * Generalized slab average formulation. Can begin some kind of small test to perturb flow. General setup complete * (June 18) TurbulentPerturbation struct now own all computation for PB. Slab average implemented. Start runtime tests. * deleted trailing white space to pass github check * Modified some names, made sure CMake compile. * Changes in attemping to pass workout test * TurbPertStruct.H doesn't like ParallelFor within it? * TurbPertStruct.H ParallelFor loops now using pass-by-value so GPU does not access illegal memory * Update TurbPertStruct.H cuda pass-by- workout test * quick fix for CUDA error * attemping to fix oneAPI SYCL workflow error * Inflow perturbation 1st PR. Cleared all github workflow error * flipped order of solverChoice.pert_type to avoid level error. ready for 1st PR --------- Co-authored-by: Ting-Hsuan (Dustin) Ma Co-authored-by: Aaron M. Lattanzi <103702284+AMLattanzi@users.noreply.github.com> --- CMake/BuildERFExe.cmake | 1 + Exec/ABL/GNUmakefile | 2 +- Source/DataStructs/DataStruct.H | 27 +- Source/DataStructs/Make.package | 1 + Source/DataStructs/TurbPertStruct.H | 328 ++++++++++++++++++ Source/ERF.H | 9 + Source/ERF.cpp | 9 + Source/ERF_make_new_arrays.cpp | 13 + Source/Initialization/ERF_init_TurbPert.cpp | 116 +++++++ Source/Initialization/ERF_init_custom.cpp | 1 - Source/Initialization/Make.package | 1 + Source/SourceTerms/ERF_make_sources.cpp | 15 +- Source/SourceTerms/Src_headers.H | 4 +- Source/TimeIntegration/ERF_Advance.cpp | 10 +- Source/TimeIntegration/ERF_advance_dycore.cpp | 3 + Source/TimeIntegration/TI_slow_rhs_fun.H | 6 +- Source/prob_common.H | 10 + 17 files changed, 547 insertions(+), 9 deletions(-) create mode 100644 Source/DataStructs/TurbPertStruct.H create mode 100644 Source/Initialization/ERF_init_TurbPert.cpp diff --git a/CMake/BuildERFExe.cmake b/CMake/BuildERFExe.cmake index 5f76b5d72..34f7640d4 100644 --- a/CMake/BuildERFExe.cmake +++ b/CMake/BuildERFExe.cmake @@ -140,6 +140,7 @@ function(build_erf_lib erf_lib_name) ${SRC_DIR}/Initialization/ERF_init_from_metgrid.cpp ${SRC_DIR}/Initialization/ERF_init_uniform.cpp ${SRC_DIR}/Initialization/ERF_init1d.cpp + ${SRC_DIR}/Initialization/ERF_init_TurbPert.cpp ${SRC_DIR}/Initialization/ERF_input_sponge.cpp ${SRC_DIR}/IO/Checkpoint.cpp ${SRC_DIR}/IO/ERF_ReadBndryPlanes.cpp diff --git a/Exec/ABL/GNUmakefile b/Exec/ABL/GNUmakefile index c85179014..f26adcb3a 100644 --- a/Exec/ABL/GNUmakefile +++ b/Exec/ABL/GNUmakefile @@ -19,7 +19,7 @@ USE_HIP = FALSE USE_SYCL = FALSE # Debugging -DEBUG = FALSE +DEBUG = TRUE TEST = TRUE USE_ASSERTION = TRUE diff --git a/Source/DataStructs/DataStruct.H b/Source/DataStructs/DataStruct.H index f3ebcb0b4..44a027b15 100644 --- a/Source/DataStructs/DataStruct.H +++ b/Source/DataStructs/DataStruct.H @@ -15,6 +15,7 @@ #include #include #include +#include enum struct ABLDriverType { None, PressureGradient, GeostrophicWind @@ -60,6 +61,11 @@ enum Sponge { ubar_sponge, vbar_sponge, nvars_sponge }; +enum struct PerturbationType { + BPM, CPM, None + // Box perturbation method, Cell perturbation method, None +}; + /** * Container holding many of the algorithmic options and parameters */ @@ -224,12 +230,28 @@ struct SolverChoice { abl_driver_type = ABLDriverType::PressureGradient; } else if (!abl_driver_type_string.compare("GeostrophicWind")) { abl_driver_type = ABLDriverType::GeostrophicWind; - } else if (!abl_driver_type_string.compare("None")){ + } else if (!abl_driver_type_string.compare("None")) { abl_driver_type = ABLDriverType::None; // No ABL driver for simulating classical fluid dynamics problems } else { amrex::Error("Don't know this abl_driver_type"); } + // Which type of inflow turbulent generation + static std::string turb_pert_type_string = "None"; + pp.query("inlet_perturbation_type",turb_pert_type_string); + if (!turb_pert_type_string.compare("BoxPerturbation")) { + pert_type = PerturbationType::BPM; + amrex::Print() << "Box perturbation method selected\n"; + } else if (!turb_pert_type_string.compare("CellPerturbation")) { + pert_type = PerturbationType::CPM; + amrex::Print() << "Cell perturbation method selected\n"; + } else if (!turb_pert_type_string.compare("None")) { + pert_type = PerturbationType::None; + amrex::Print() << "No perturbation method selected\n"; + } else { + amrex::Abort("Dont know this inlet_perturbation_type"); + } + amrex::Vector abl_pressure_grad_in = {0.0, 0.0, 0.0}; pp.queryarr("abl_pressure_grad",abl_pressure_grad_in); for(int i = 0; i < AMREX_SPACEDIM; ++i) abl_pressure_grad[i] = abl_pressure_grad_in[i]; @@ -487,6 +509,9 @@ struct SolverChoice { // User wishes to output time averaged velocity fields bool time_avg_vel = false; + // Type of perturbation + PerturbationType pert_type; + // Numerical diffusion bool use_NumDiff{false}; amrex::Real NumDiffCoeff{0.}; diff --git a/Source/DataStructs/Make.package b/Source/DataStructs/Make.package index b82b7b0ad..c1bc1767f 100644 --- a/Source/DataStructs/Make.package +++ b/Source/DataStructs/Make.package @@ -3,3 +3,4 @@ CEXE_headers += DiffStruct.H CEXE_headers += AdvStruct.H CEXE_headers += SpongeStruct.H CEXE_headers += TurbStruct.H +CEXE_headers += TurbPertStruct.H diff --git a/Source/DataStructs/TurbPertStruct.H b/Source/DataStructs/TurbPertStruct.H new file mode 100644 index 000000000..5d5021248 --- /dev/null +++ b/Source/DataStructs/TurbPertStruct.H @@ -0,0 +1,328 @@ +#ifndef _TURB_PERT_STRUCT_H_ +#define _TURB_PERT_STRUCT_H_ + +#include +#include + +/** + * Container holding quantities related to turbulent perturbation parameters + */ + +struct TurbulentPerturbation { + + public: + + ~TurbulentPerturbation () { + amrex::Print() << "~TurbulentPerturbation destructor\n"; + } + +//#define DIAGONAL_REF +#define LENGTH_REF + +/* void init_tpi_const(const amrex::Real T_0) + { + // Bulk Richardson number + tpi_nonDim = 0.042; + + // Ambient Temperature + tpi_Tinfty = T_0; + }*/ + + // Initializing Perturbation Region + // Note: Box will always be initialized at cell centered + void init_tpi (const int lev, + const amrex::IntVect& nx, + amrex::GpuArray const dx) + + { + // TODO: Use user Param to define this (Currently HARDCODED) + // Perturbation box region setup + int pertBox_offset = 0; + amrex::IntVect regionLo(0+pertBox_offset, 0, 0); + amrex::IntVect regionHi(7+pertBox_offset, nx[1], nx[2]); + amrex::IntVect boxSize(8,8,8); + //amrex::IntVect boxSize(4,4,4); // HACK + + // TODO: This will need to be changed later on + tpi_Lpb = 8*dx[0]; + tpi_Wpb = 8*dx[1]; + tpi_Hpb = 8*dx[2]; + + #ifdef DIAGONAL_REF + tpi_lref = sqrt(boxSize(0)*boxSize(0)*dx[0]*dx[0] + boxSize(1)*boxSize(1)*dx[1]*dx[1]); + #elif defined(LENGTH_REF) + tpi_lref = 8*dx[0]; + #endif + + // Creating structure box array for conserved quantity + amrex::Box tmp_bx(regionLo, regionHi, amrex::IntVect(0,0,0)); + amrex::BoxArray tmp_ba(tmp_bx); + tmp_ba.maxSize(boxSize); + pb_ba.push_back(tmp_ba); + + // Initializing mean magnitude vector + pb_mag.resize(pb_ba[lev].size(), 0.); + + // Set size of vector and initialize + pb_freq.resize(pb_ba[lev].size(), -1.0); + pb_local_etime.resize(pb_ba[lev].size(), 0.0); + pb_amp.resize(pb_ba[lev].size(), 0.0); + + // Function check point message + amrex::Print() << "Turbulent perturbation region initialized with ba " << pb_ba[lev].size() << " boxes\n"; + + // Function method check point message + amrex::Print() << "Using "; + #ifdef DIAGONAL_REF + amrex::Print() << "box diagonal "; + #elif defined(LENGTH_REF) + amrex::Print() << "box length "; + #endif + amrex::Print() << "distance as reference for turbulent inlet perturbation: tpi_lref = " << tpi_lref << "\n"; + } + + // Perturbation update frequency check. + // This function will trigger calc_tpi_meanMag() and calc_tpi_amp(), when + // the local elasped time is greater than the perturbation frequency + void calc_tpi_update (const int lev, + const amrex::Real dt, + amrex::MultiFab& mf_xvel, + amrex::MultiFab& mf_yvel, + amrex::MultiFab& mf_cons) + { + for (int boxIdx = 0; boxIdx < pb_ba[lev].size(); boxIdx++) { + if ( pb_freq[boxIdx] <= pb_local_etime[boxIdx] ) { + // Compute mean velocity magnitude within perturbation box + calc_tpi_meanMag(boxIdx, lev, mf_xvel, mf_yvel, mf_cons); + + // Ma and Senocak (2015) Eq. 15 + if (pb_mag[boxIdx] !=0.) { + pb_freq[boxIdx] = tpi_lref / pb_mag[boxIdx]; + } + + // Trigger amplitude calculation per perturbation box + calc_tpi_amp(boxIdx); + + // Reset local elapsed time + pb_local_etime[boxIdx] = 0.; + } else { + // Increase by timestep + pb_local_etime[boxIdx] += dt; + } // if + } // for + } + +// TODO: Test the difference between these two +//#define USE_SLAB_AVERAGE +#define USE_VOLUME_AVERAGE + + // Perturbation box mean velocity magnitude calculation + // This is pulled into the strucutre to also utilize during runtime + void calc_tpi_meanMag (const int boxIdx, + const int lev, + amrex::MultiFab& mf_xvel, + amrex::MultiFab& mf_yvel, + amrex::MultiFab& mf_cons) + + { + // Creating local copy of PB box array and magnitude + const amrex::BoxArray m_pb_ba = pb_ba[lev]; + amrex::Real* m_pb_mag = get_pb_mag(); + + // Storage of averages per PB + // Index: 0=u (vol/slab_lo), 1=v (vol/slab_lo) + // 2=u (slab_hi), 3=v (slab_hi) + amrex::Gpu::DeviceVector tmp_d(4,0.); + amrex::Real* tmp = tmp_d.data(); + + // Averaging u components + for (amrex::MFIter mfi(mf_xvel, TileNoZ()) ; mfi.isValid(); ++mfi) { + const amrex::Box &vbx = mfi.validbox(); + amrex::Box pbx = convert(m_pb_ba[boxIdx], vbx.ixType()); + amrex::Box ubx = pbx & vbx; + + // Operation over box union + if (ubx.ok()) { + const amrex::Array4 &xvel_arry = mf_xvel.array(mfi); + + #ifdef USE_VOLUME_AVERAGE + int npts = ubx.numPts(); + ParallelFor(amrex::Gpu::KernelInfo().setReduction(true), ubx, [=] + AMREX_GPU_DEVICE(int i, int j, int k, amrex::Gpu::Handler const& handler) noexcept { + amrex::Gpu::deviceReduceSum(&tmp[0], xvel_arry(i,j,k), handler); + }); + tmp[0] /= (amrex::Real) npts; + #endif + + #ifdef USE_SLAB_AVERAGE + amrex::Box ubxSlab_lo = makeSlab(ubx,2,ubx.smallEnd(2)); + amrex::Box ubxSlab_hi = makeSlab(ubx,2,ubx.bigEnd(2)); + int npts_lo = ubxSlab_lo.numPts(); + int npts_hi = ubxSlab_hi.numPts(); + + // Average u in the low slab + ParallelFor(amrex::Gpu::KernelInfo().setReduction(true), ubxSlab_lo, [=] + AMREX_GPU_DEVICE(int i, int j, int k, amrex::Gpu::Handler const& handler) noexcept { + amrex::Gpu::deviceReduceSum(&tmp[0], xvel_arry(i,j,k), handler); + }); + tmp[0] /= (amrex::Real) npts_lo; + + // Average u in the high slab + ParallelFor(amrex::Gpu::KernelInfo().setReduction(true), ubxSlab_hi, [=] + AMREX_GPU_DEVICE(int i, int j, int k, amrex::Gpu::Handler const& handler) noexcept { + amrex::Gpu::deviceReduceSum(&tmp[2], xvel_arry(i,j,k), handler); + }); + tmp[2] /= (amrex::Real) npts_hi; + #endif + } // if + } // MFIter + + // Averaging v components + for (amrex::MFIter mfi(mf_yvel, TileNoZ()); mfi.isValid(); ++mfi) { + amrex::Box const& vbx = mfi.validbox(); + amrex::Box pbx = convert(m_pb_ba[boxIdx], vbx.ixType()); + amrex::Box ubx = pbx & vbx; + + // Operation over box union + if (ubx.ok()) { + const amrex::Array4 &yvel_arry = mf_yvel.array(mfi); + + #ifdef USE_VOLUME_AVERAGE + int npts = ubx.numPts(); + ParallelFor(amrex::Gpu::KernelInfo().setReduction(true), ubx, [=] + AMREX_GPU_DEVICE(int i, int j, int k, amrex::Gpu::Handler const& handler) noexcept { + amrex::Gpu::deviceReduceSum(&tmp[1], yvel_arry(i,j,k), handler); + }); + tmp[1] /= (amrex::Real) npts; + #endif + + #ifdef USE_SLAB_AVERAGE + amrex::Box ubxSlab_lo = makeSlab(ubx,2,ubx.smallEnd(2)); + amrex::Box ubxSlab_hi = makeSlab(ubx,2,ubx.bigEnd(2)); + int npts_lo = ubxSlab_lo.numPts(); + int npts_hi = ubxSlab_hi.numPts(); + + // Average v in the low slab + ParallelFor(amrex::Gpu::KernelInfo().setReduction(true), ubxSlab_lo, [=] + AMREX_GPU_DEVICE(int i, int j, int k, amrex::Gpu::Handler const& handler) noexcept { + amrex::Gpu::deviceReduceSum(&tmp[1], yvel_arry(i,j,k), handler); + }); + tmp[1] /= (amrex::Real) npts_lo; + + // Average v in the high slab + ParallelFor(amrex::Gpu::KernelInfo().setReduction(true), ubxSlab_hi, [=] + AMREX_GPU_DEVICE(int i, int j, int k, amrex::Gpu::Handler const& handler) noexcept { + amrex::Gpu::deviceReduceSum(&tmp[3], yvel_arry(i,j,k), handler); + }); + tmp[3] /= (amrex::Real) npts_hi; + #endif + } // if + } // MFIter + + // Computing the average magnitude within PB + for (amrex::MFIter mfi(mf_cons, TileNoZ()); mfi.isValid(); ++mfi) { + amrex::Box const& vbx = mfi.validbox(); + amrex::Box pbx = convert(m_pb_ba[boxIdx], vbx.ixType()); + amrex::Box ubx = pbx & vbx; + if (ubx.ok()) { + #ifdef USE_SLAB_AVERAGE + m_pb_mag[boxIdx] = 0.5*(sqrt(tmp[0]*tmp[0] + tmp[1]*tmp[1]) + sqrt(tmp[2]*tmp[2] + tmp[3]*tmp[3])); + #endif + + #ifdef USE_VOLUME_AVERAGE + m_pb_mag[boxIdx] = sqrt(tmp[0]*tmp[0] + tmp[1]*tmp[1]); + #endif + } // if + } // MFIter + + } + +//#define ECKERT_FORMULATION +#define BULK_RICHARDSON_FORMULATION + + // Perturbation amplitude calcluation + void calc_tpi_amp (const int boxIdx) + { + amrex::Real g = CONST_GRAV; + amrex::Real beta = 1./tpi_Tinfty; + amrex::Real Um = pb_mag[boxIdx]; + + // Reset the perturbation amplitude + pb_amp[boxIdx] = 0.; + + #ifdef BULK_RICHARDSON_FORMULATION + // Ma and Senocak (2023) Eq. 17 + pb_amp[boxIdx] = tpi_nonDim * Um * Um * Um / (g * beta * tpi_lref * tpi_Hpb); + #endif + + } + +//#define INDEX_PERTURB + + // Applying perturbation amplitude onto source term (Umphrey and Senocak 2016) + // Random perturbation is applied as a white noise on the temperature source term + // it becomes colored through the filtering of temperature transport equation from + // the eddy diffusivity. + void apply_tpi (const int lev, + const amrex::Box& vbx, // box union from upper level + const int comp, // Component to modify + const amrex::IndexType m_ixtype, // IntVect type of src_arr + amrex::Array4 const& src_arr)// Source Array to apply perturbation + { + for (int boxIdx = 0; boxIdx < pb_ba[lev].size(); boxIdx++) { + amrex::Box pbx = convert(pb_ba[lev][boxIdx], m_ixtype); + amrex::Box ubx = pbx & vbx; + if (ubx.ok()) { + amrex::Real amp_copy = pb_amp[boxIdx]; // This exposes a copy so the GPU can capture by value + ParallelForRNG(ubx, [=] AMREX_GPU_DEVICE(int i, int j, int k, const amrex::RandomEngine& engine) noexcept { + // Random number generation between 0 and 1 + amrex::Real rand_double = amrex::Random(engine); + + // Rescale between -1 and 1 and multiple with amplitude + src_arr(i,j,k,comp) = (rand_double*2.0 - 1.0) * amp_copy; + + // For box region debug only + #ifdef INDEX_PERTURB + src_arr(i,j,k,comp) = (amrex::Real) (boxIdx + 1.); + #endif + }); + } + } + } + + + // Quality of life function definitions + amrex::Real* get_pb_mag() { return pb_mag.data(); } + + void debug () + { + for (int i = 0; i < pb_mag.size(); i++) { + amrex::Print() << "[" << i<< "] pb_mag=" << pb_mag[i] + << " | pb_freq=" << pb_freq[i] + << " | pb_local_etime=" << pb_local_etime[i] <<"\n"; + } + } + + // Public constant constants + amrex::Real tpi_Tinfty; // [K] + amrex::Real tpi_nonDim; + + // Public data members + amrex::Vector pb_ba; // PB box array + amrex::Vector pb_mag; // BP mean magnitude [m/s] + + private: + + // Private constant constants + amrex::Real tpi_Hpb; // PB height [m] + amrex::Real tpi_Lpb; // PB length [m] + amrex::Real tpi_Wpb; // PB width [m] + amrex::Real tpi_lref; // PB reference length [m] + + // Storage Arrays + amrex::Vector pb_freq; // PB frequency [1/s] + amrex::Vector pb_local_etime; // PB local elapsed time [s] + amrex::Vector pb_amp; // PB perturbation amplitude Ri:[K] Ec:[K/s] + +}; +#endif diff --git a/Source/ERF.H b/Source/ERF.H index 4cb59e790..926d1d183 100644 --- a/Source/ERF.H +++ b/Source/ERF.H @@ -29,6 +29,7 @@ #include #include +#include #include #include #include @@ -423,6 +424,11 @@ private: amrex::Vector& lev_new, amrex::Vector& lev_old); + // Initialize the Turbulent perturbation + void turbPert_constants (const int lev); + void turbPert_update (const int lev, const amrex::Real dt); + void turbPert_amplitude (const int lev); + // Initialize the integrator object void initialize_integrator (int lev, amrex::MultiFab& cons_mf, amrex::MultiFab& vel_mf); @@ -842,6 +848,9 @@ private: // algorithm choices static SolverChoice solverChoice; + // Turbulent perturbation structure + TurbulentPerturbation turbPert; + #ifdef ERF_USE_PARTICLES // Particle container with all particle species ParticleData particleData; diff --git a/Source/ERF.cpp b/Source/ERF.cpp index 26faab742..7fff9ac0e 100644 --- a/Source/ERF.cpp +++ b/Source/ERF.cpp @@ -1216,6 +1216,15 @@ ERF::init_only (int lev, Real time) input_sponge(lev); } + // Initialize turbulent perturbation + if (solverChoice.pert_type == PerturbationType::BPM) { + if (lev == 0) { + turbPert_constants(lev); + turbPert_update(lev, 0.); + turbPert_amplitude(lev); + } + } + } // read in some parameters from inputs file diff --git a/Source/ERF_make_new_arrays.cpp b/Source/ERF_make_new_arrays.cpp index 7e2e5d561..710abb39e 100644 --- a/Source/ERF_make_new_arrays.cpp +++ b/Source/ERF_make_new_arrays.cpp @@ -304,6 +304,16 @@ ERF::init_stuff (int lev, const BoxArray& ba, const DistributionMapping& dm, } #endif + //********************************************************* + // Turbulent perturbation region initialization + //********************************************************* + // TODO: Test perturbation on multiple levels + if (solverChoice.pert_type == PerturbationType::BPM) { + if (lev == 0) { + turbPert.init_tpi(lev, geom[lev].Domain().bigEnd(), geom[lev].CellSizeArray()); + } + } + // // Define the land mask here and set it to all land // NOTE: the logic below will BREAK if we have any grids not touching the bottom boundary @@ -476,3 +486,6 @@ ERF::initialize_bcs (int lev) (lev, geom[lev], domain_bcs_type, domain_bcs_type_d, m_bc_extdir_vals, m_bc_neumann_vals, use_real_bcs); } + + + diff --git a/Source/Initialization/ERF_init_TurbPert.cpp b/Source/Initialization/ERF_init_TurbPert.cpp new file mode 100644 index 000000000..beb2c924d --- /dev/null +++ b/Source/Initialization/ERF_init_TurbPert.cpp @@ -0,0 +1,116 @@ +/** + * \ERF_init_TurbPert.cpp + */ + +//DUSTIN MA: May 28th, 2024 + +#define DEBUG_PERTBOX_MSG + +#include +#include +#include +#include + +using namespace amrex; + +void +ERF::turbPert_constants(const int lev) +{ + prob->init_turbPert_const(turbPert); +} + +void +ERF::turbPert_update (const int lev, const Real local_dt) +{ + // Grabing data from velocity field + auto& lev_new = vars_new[lev]; + + // Accessing local data + MultiFab cons_data(lev_new[Vars::cons].boxArray(), lev_new[Vars::cons].DistributionMap(), + lev_new[Vars::cons].nComp() , lev_new[Vars::cons].nGrow()); + MultiFab xvel_data(lev_new[Vars::xvel].boxArray(), lev_new[Vars::xvel].DistributionMap(), 1, lev_new[Vars::xvel].nGrowVect()); + MultiFab yvel_data(lev_new[Vars::yvel].boxArray(), lev_new[Vars::yvel].DistributionMap(), 1, lev_new[Vars::yvel].nGrowVect()); + MultiFab::Copy (cons_data, lev_new[Vars::cons], 0, 0, 1, lev_new[Vars::xvel].nGrowVect()); + MultiFab::Copy (xvel_data, lev_new[Vars::xvel], 0, 0, 1, lev_new[Vars::xvel].nGrowVect()); + MultiFab::Copy (yvel_data, lev_new[Vars::yvel], 0, 0, 1, lev_new[Vars::yvel].nGrowVect()); + + // Computing perturbation update time + turbPert.calc_tpi_update(lev, local_dt, xvel_data, yvel_data, cons_data); + + #ifdef DEBUG_PERTBOX_MSG + turbPert.debug(); + #endif + + Print() << "Turbulent perturbation update time and amplitude initialized\n"; +} + +// Calculate the perturbation region amplitude. +// This function heavily emmulates the ERF::init_custom () +void +ERF::turbPert_amplitude (int lev) +{ + auto& lev_new = vars_new[lev]; + + MultiFab r_hse(base_state[lev], make_alias, 0, 1); // r_0 is first component + MultiFab p_hse(base_state[lev], make_alias, 1, 1); // p_0 is second component + + MultiFab cons_pert(lev_new[Vars::cons].boxArray(), lev_new[Vars::cons].DistributionMap(), + lev_new[Vars::cons].nComp() , lev_new[Vars::cons].nGrow()); + MultiFab xvel_pert(lev_new[Vars::xvel].boxArray(), lev_new[Vars::xvel].DistributionMap(), 1, lev_new[Vars::xvel].nGrowVect()); + MultiFab yvel_pert(lev_new[Vars::yvel].boxArray(), lev_new[Vars::yvel].DistributionMap(), 1, lev_new[Vars::yvel].nGrowVect()); + MultiFab zvel_pert(lev_new[Vars::zvel].boxArray(), lev_new[Vars::zvel].DistributionMap(), 1, lev_new[Vars::zvel].nGrowVect()); + + // Only storing for conserved quantity for now + auto m_ixtype = cons_pert.boxArray().ixType(); + + // Default perturbations to zero + cons_pert.setVal(0.); + xvel_pert.setVal(0.); + yvel_pert.setVal(0.); + zvel_pert.setVal(0.); + + int fix_random_seed = 0; + ParmParse pp("erf"); pp.query("fix_random_seed", fix_random_seed); + // Note that the value of 1024UL is not significant -- the point here is just to set the + // same seed for all MPI processes for the purpose of regression testing + if (fix_random_seed) { + Print() << "Fixing the random seed" << std::endl; + InitRandom(1024UL); + } + +#ifdef _OPENMP +#pragma omp parallel if (amrex::Gpu::notInLaunchRegion()) +#endif + for (MFIter mfi(lev_new[Vars::cons], TileNoZ()); mfi.isValid(); ++mfi) { + const Box &bx = mfi.validbox(); + const Box &xbx = mfi.tilebox(IntVect(1,0,0)); + const Box &ybx = mfi.tilebox(IntVect(0,1,0)); + const Box &zbx = mfi.tilebox(IntVect(0,0,1)); + + // Perturbation on to different components + const auto &cons_pert_arr = cons_pert.array(mfi); + const auto &xvel_pert_arr = xvel_pert.array(mfi); + const auto &yvel_pert_arr = yvel_pert.array(mfi); + const auto &zvel_pert_arr = zvel_pert.array(mfi); + + Array4 cons_arr = lev_new[Vars::cons].const_array(mfi); + Array4 z_nd_arr = (solverChoice.use_terrain) ? z_phys_nd[lev]->const_array(mfi) : Array4{}; + Array4 z_cc_arr = (solverChoice.use_terrain) ? z_phys_cc[lev]->const_array(mfi) : Array4{}; + + Array4 mf_m = mapfac_m[lev]->array(mfi); + Array4 mf_u = mapfac_m[lev]->array(mfi); + Array4 mf_v = mapfac_m[lev]->array(mfi); + + Array4 r_hse_arr = r_hse.array(mfi); + Array4 p_hse_arr = p_hse.array(mfi); + + turbPert.apply_tpi(lev, bx, RhoTheta_comp, m_ixtype, cons_pert_arr); // Using boxArray + prob->init_custom_pert(bx, xbx, ybx, zbx, cons_arr, cons_pert_arr, + xvel_pert_arr, yvel_pert_arr, zvel_pert_arr, + r_hse_arr, p_hse_arr, z_nd_arr, z_cc_arr, + geom[lev].data(), mf_m, mf_u, mf_v, + solverChoice); + } // mfi + // This initialization adds the initial perturbation direction on the RhoTheta field + MultiFab::Add(lev_new[Vars::cons], cons_pert, RhoTheta_comp, RhoTheta_comp, 1, cons_pert.nGrow()); +} diff --git a/Source/Initialization/ERF_init_custom.cpp b/Source/Initialization/ERF_init_custom.cpp index 98e935d9f..61f3368e8 100644 --- a/Source/Initialization/ERF_init_custom.cpp +++ b/Source/Initialization/ERF_init_custom.cpp @@ -61,7 +61,6 @@ ERF::init_custom (int lev) const Box &ybx = mfi.tilebox(IntVect(0,1,0)); const Box &zbx = mfi.tilebox(IntVect(0,0,1)); - const auto &cons_pert_arr = cons_pert.array(mfi); const auto &xvel_pert_arr = xvel_pert.array(mfi); const auto &yvel_pert_arr = yvel_pert.array(mfi); diff --git a/Source/Initialization/Make.package b/Source/Initialization/Make.package index 685889786..21444ddce 100644 --- a/Source/Initialization/Make.package +++ b/Source/Initialization/Make.package @@ -5,6 +5,7 @@ CEXE_sources += ERF_init_from_input_sounding.cpp CEXE_sources += ERF_input_sponge.cpp CEXE_sources += ERF_init_bcs.cpp CEXE_sources += ERF_init1d.cpp +CEXE_sources += ERF_init_TurbPert.cpp ifeq ($(USE_WINDFARM),TRUE) CEXE_sources += ERF_init_windfarm.cpp diff --git a/Source/SourceTerms/ERF_make_sources.cpp b/Source/SourceTerms/ERF_make_sources.cpp index bb1e38fa6..69fbc5002 100644 --- a/Source/SourceTerms/ERF_make_sources.cpp +++ b/Source/SourceTerms/ERF_make_sources.cpp @@ -44,12 +44,13 @@ void make_sources (int level, const Real* dptr_rhotheta_src, const Real* dptr_rhoqt_src, const Real* dptr_wbar_sub, - const Vector d_rayleigh_ptrs_at_lev) + const Vector d_rayleigh_ptrs_at_lev, + TurbulentPerturbation& turbPert) { BL_PROFILE_REGION("erf_make_sources()"); // ***************************************************************************** - // Initialize source to zero since we re-compute it ever RK stage + // Initialize source to zero since we re-compute it every RK stage // ***************************************************************************** source.setVal(0.0); @@ -303,6 +304,16 @@ void make_sources (int level, if(!(solverChoice.spongeChoice.sponge_type == "input_sponge")){ ApplySpongeZoneBCsForCC(solverChoice.spongeChoice, geom, bx, cell_src, cell_data); } + + // ************************************************************************************* + // Add perturbation + // ************************************************************************************* + if (solverChoice.pert_type == PerturbationType::BPM) { + auto m_ixtype = S_data[IntVars::cons].boxArray().ixType(); + + // Apply stored values onto cell_src + turbPert.apply_tpi(level, bx, RhoTheta_comp, m_ixtype, cell_src); + } } // mfi } // OMP } diff --git a/Source/SourceTerms/Src_headers.H b/Source/SourceTerms/Src_headers.H index da30c721c..5cbb0c7a2 100644 --- a/Source/SourceTerms/Src_headers.H +++ b/Source/SourceTerms/Src_headers.H @@ -3,6 +3,7 @@ #include #include "DataStruct.H" +#include "TurbPertStruct.H" #ifdef ERF_USE_EB #include @@ -36,7 +37,8 @@ void make_sources (int level, int nrk, const amrex::Real* dptr_rhotheta_src, const amrex::Real* dptr_rhoqt_src, const amrex::Real* dptr_wbar_sub, - const amrex::Vector d_rayleigh_ptrs_at_lev); + const amrex::Vector d_rayleigh_ptrs_at_lev, + TurbulentPerturbation& turbPert); void make_mom_sources (int level, int nrk, amrex::Real dt, diff --git a/Source/TimeIntegration/ERF_Advance.cpp b/Source/TimeIntegration/ERF_Advance.cpp index 60fd99139..02448f2bf 100644 --- a/Source/TimeIntegration/ERF_Advance.cpp +++ b/Source/TimeIntegration/ERF_Advance.cpp @@ -35,7 +35,15 @@ ERF::Advance (int lev, Real time, Real dt_lev, int iteration, int /*ncycle*/) MultiFab& V_new = vars_new[lev][Vars::yvel]; MultiFab& W_new = vars_new[lev][Vars::zvel]; - + // TODO: Can test on multiple levels later + // Only apply to level 0 + // DUSTIN MA + if (solverChoice.pert_type == PerturbationType::BPM) { + if (lev == 0) { + turbPert.calc_tpi_update(lev, dt_lev, U_old, V_old, S_old); + turbPert.debug(); + } + } // configure ABLMost params if used MostWall boundary condition if (phys_bc_type[Orientation(Direction::z,Orientation::low)] == ERF_BC::MOST) { if (m_most) { diff --git a/Source/TimeIntegration/ERF_advance_dycore.cpp b/Source/TimeIntegration/ERF_advance_dycore.cpp index d4691bf38..4e59b5c9d 100644 --- a/Source/TimeIntegration/ERF_advance_dycore.cpp +++ b/Source/TimeIntegration/ERF_advance_dycore.cpp @@ -66,6 +66,9 @@ void ERF::advance_dycore(int level, Real* dptr_rhoqt_src = solverChoice.custom_moisture_forcing ? d_rhoqt_src[level].data() : nullptr; Real* dptr_wbar_sub = solverChoice.custom_w_subsidence ? d_w_subsid[level].data() : nullptr; + // Turbulent Perturbation Pointer + //Real* dptr_rhotheta_src = solverChoice.pert_type ? d_rhotheta_src[level].data() : nullptr; + Vector d_rayleigh_ptrs_at_lev; d_rayleigh_ptrs_at_lev.resize(Rayleigh::nvars); bool rayleigh_damp_any = (solverChoice.rayleigh_damp_U ||solverChoice.rayleigh_damp_V || diff --git a/Source/TimeIntegration/TI_slow_rhs_fun.H b/Source/TimeIntegration/TI_slow_rhs_fun.H index 54a709392..0c50a2d22 100644 --- a/Source/TimeIntegration/TI_slow_rhs_fun.H +++ b/Source/TimeIntegration/TI_slow_rhs_fun.H @@ -44,7 +44,8 @@ fine_geom, solverChoice, mapfac_u[level], mapfac_v[level], dptr_rhotheta_src, dptr_rhoqt_src, - dptr_wbar_sub, d_rayleigh_ptrs_at_lev); + dptr_wbar_sub, d_rayleigh_ptrs_at_lev, + turbPert); // Moving terrain if ( solverChoice.use_terrain && (solverChoice.terrain_type == TerrainType::Moving) ) @@ -413,7 +414,8 @@ fine_geom, solverChoice, mapfac_u[level], mapfac_v[level], dptr_rhotheta_src, dptr_rhoqt_src, - dptr_wbar_sub, d_rayleigh_ptrs_at_lev); + dptr_wbar_sub, d_rayleigh_ptrs_at_lev, + turbPert); int n_qstate = micro->Get_Qstate_Size(); make_mom_sources(level, nrk, slow_dt, S_data, S_prim, diff --git a/Source/prob_common.H b/Source/prob_common.H index 0e35a1869..acba961b7 100644 --- a/Source/prob_common.H +++ b/Source/prob_common.H @@ -365,6 +365,16 @@ public: }); } + /** + * Function to set constant values for turbulent pertubation calculation + */ + void + init_turbPert_const(TurbulentPerturbation& turbPert) + { + turbPert.tpi_Tinfty = base_parms.T_0; + turbPert.tpi_nonDim = 0.042; + } + protected: // Struct to store problem parameters ProbParmDefaults base_parms; From e571ffe625c319da2506ba2a7e88cdb5a5ee1a8c Mon Sep 17 00:00:00 2001 From: "Aaron M. Lattanzi" <103702284+AMLattanzi@users.noreply.github.com> Date: Tue, 25 Jun 2024 13:09:50 -0700 Subject: [PATCH 7/9] Terrain From File (#1658) * first pass at terrain initialization from file. * Uncomment time for the new function call. --- Exec/DevTests/MovingTerrain/prob.cpp | 94 +++++----- Exec/RegTests/ParticlesOverWoA/prob.cpp | 106 ++++++----- Exec/RegTests/Terrain2d_Cylinder/prob.cpp | 106 ++++++----- Exec/RegTests/Terrain3d_Hemisphere/prob.cpp | 96 +++++----- Exec/RegTests/WitchOfAgnesi/prob.cpp | 194 ++++++++++---------- Exec/SpongeTest/prob.cpp | 106 ++++++----- Source/prob_common.H | 80 ++++++++ 7 files changed, 457 insertions(+), 325 deletions(-) diff --git a/Exec/DevTests/MovingTerrain/prob.cpp b/Exec/DevTests/MovingTerrain/prob.cpp index c251cda5b..82b7a346b 100644 --- a/Exec/DevTests/MovingTerrain/prob.cpp +++ b/Exec/DevTests/MovingTerrain/prob.cpp @@ -155,47 +155,57 @@ Problem::init_custom_terrain (const Geometry& geom, MultiFab& z_phys_nd, const Real& time) { - // Domain cell size and real bounds - auto dx = geom.CellSizeArray(); - - // Domain valid box (z_nd is nodal) - const amrex::Box& domain = geom.Domain(); - int domlo_x = domain.smallEnd(0); int domhi_x = domain.bigEnd(0) + 1; - int domlo_z = domain.smallEnd(2); - - // Number of ghost cells - int ngrow = z_phys_nd.nGrow(); - - // Populate bottom plane - int k0 = domlo_z; - - Real Ampl = parms.Ampl; - Real wavelength = 100.; - Real kp = 2.0 * PI / wavelength; - Real g = CONST_GRAV; - Real omega = std::sqrt(g * kp); - - for (MFIter mfi(z_phys_nd,amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi) - { - // Grown box with no z range - amrex::Box xybx = mfi.growntilebox(ngrow); - xybx.setRange(2,0); - - amrex::Array4 const& z_arr = z_phys_nd.array(mfi); - - ParallelFor(xybx, [=] AMREX_GPU_DEVICE (int i, int j, int) { - - // Clip indices for ghost-cells - int ii = amrex::min(amrex::max(i,domlo_x),domhi_x); - - // Location of nodes - Real x = ii * dx[0]; - - // Wave height - Real height = Ampl * std::sin(kp * x - omega * time); - - // Populate terrain height - z_arr(i,j,k0) = height; - }); + // Check if a valid csv file exists for the terrain + std::string fname; + amrex::ParmParse pp("erf"); + auto valid_fname = pp.query("terrain_file_name",fname); + if (valid_fname) { + this->read_custom_terrain(fname,geom,z_phys_nd,time); + } else { + + // Domain cell size and real bounds + auto dx = geom.CellSizeArray(); + + // Domain valid box (z_nd is nodal) + const amrex::Box& domain = geom.Domain(); + int domlo_x = domain.smallEnd(0); int domhi_x = domain.bigEnd(0) + 1; + int domlo_z = domain.smallEnd(2); + + // Number of ghost cells + int ngrow = z_phys_nd.nGrow(); + + // Populate bottom plane + int k0 = domlo_z; + + Real Ampl = parms.Ampl; + Real wavelength = 100.; + Real kp = 2.0 * PI / wavelength; + Real g = CONST_GRAV; + Real omega = std::sqrt(g * kp); + + for (MFIter mfi(z_phys_nd,amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi) + { + // Grown box with no z range + amrex::Box xybx = mfi.growntilebox(ngrow); + xybx.setRange(2,0); + + amrex::Array4 const& z_arr = z_phys_nd.array(mfi); + + ParallelFor(xybx, [=] AMREX_GPU_DEVICE (int i, int j, int) + { + + // Clip indices for ghost-cells + int ii = amrex::min(amrex::max(i,domlo_x),domhi_x); + + // Location of nodes + Real x = ii * dx[0]; + + // Wave height + Real height = Ampl * std::sin(kp * x - omega * time); + + // Populate terrain height + z_arr(i,j,k0) = height; + }); + } } } diff --git a/Exec/RegTests/ParticlesOverWoA/prob.cpp b/Exec/RegTests/ParticlesOverWoA/prob.cpp index 422823dde..b65838449 100644 --- a/Exec/RegTests/ParticlesOverWoA/prob.cpp +++ b/Exec/RegTests/ParticlesOverWoA/prob.cpp @@ -126,54 +126,64 @@ void Problem::init_custom_terrain( const Geometry& geom, MultiFab& z_phys_nd, - const Real& /*time*/) + const Real& time) { - // Domain cell size and real bounds - auto dx = geom.CellSizeArray(); - auto ProbLoArr = geom.ProbLoArray(); - auto ProbHiArr = geom.ProbHiArray(); - - // Domain valid box (z_nd is nodal) - const amrex::Box& domain = geom.Domain(); - int domlo_x = domain.smallEnd(0); int domhi_x = domain.bigEnd(0) + 1; - // int domlo_y = domain.smallEnd(1); int domhi_y = domain.bigEnd(1) + 1; - int domlo_z = domain.smallEnd(2); - - // User function parameters - Real a = 0.5; - Real num = 8 * a * a * a; - Real xcen = 0.5 * (ProbLoArr[0] + ProbHiArr[0]); - // Real ycen = 0.5 * (ProbLoArr[1] + ProbHiArr[1]); - - // Number of ghost cells - int ngrow = z_phys_nd.nGrow(); - - // Populate bottom plane - int k0 = domlo_z; - - for ( amrex::MFIter mfi(z_phys_nd,amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi ) - { - // Grown box with no z range - amrex::Box xybx = mfi.growntilebox(ngrow); - xybx.setRange(2,0); - - amrex::Array4 const& z_arr = z_phys_nd.array(mfi); - - ParallelFor(xybx, [=] AMREX_GPU_DEVICE (int i, int j, int) { - - // Clip indices for ghost-cells - int ii = amrex::min(amrex::max(i,domlo_x),domhi_x); - // int jj = amrex::min(amrex::max(j,domlo_y),domhi_y); - - // Location of nodes - Real x = (ii * dx[0] - xcen); - // Real y = (jj * dx[1] - ycen); - - // WoA Hill in x-direction - Real height = num / (x*x + 4 * a * a); - - // Populate terrain height - z_arr(i,j,k0) = height; - }); + // Check if a valid csv file exists for the terrain + std::string fname; + amrex::ParmParse pp("erf"); + auto valid_fname = pp.query("terrain_file_name",fname); + if (valid_fname) { + this->read_custom_terrain(fname,geom,z_phys_nd,time); + } else { + + // Domain cell size and real bounds + auto dx = geom.CellSizeArray(); + auto ProbLoArr = geom.ProbLoArray(); + auto ProbHiArr = geom.ProbHiArray(); + + // Domain valid box (z_nd is nodal) + const amrex::Box& domain = geom.Domain(); + int domlo_x = domain.smallEnd(0); int domhi_x = domain.bigEnd(0) + 1; + // int domlo_y = domain.smallEnd(1); int domhi_y = domain.bigEnd(1) + 1; + int domlo_z = domain.smallEnd(2); + + // User function parameters + Real a = 0.5; + Real num = 8 * a * a * a; + Real xcen = 0.5 * (ProbLoArr[0] + ProbHiArr[0]); + // Real ycen = 0.5 * (ProbLoArr[1] + ProbHiArr[1]); + + // Number of ghost cells + int ngrow = z_phys_nd.nGrow(); + + // Populate bottom plane + int k0 = domlo_z; + + for ( amrex::MFIter mfi(z_phys_nd,amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi ) + { + // Grown box with no z range + amrex::Box xybx = mfi.growntilebox(ngrow); + xybx.setRange(2,0); + + amrex::Array4 const& z_arr = z_phys_nd.array(mfi); + + ParallelFor(xybx, [=] AMREX_GPU_DEVICE (int i, int j, int) + { + + // Clip indices for ghost-cells + int ii = amrex::min(amrex::max(i,domlo_x),domhi_x); + // int jj = amrex::min(amrex::max(j,domlo_y),domhi_y); + + // Location of nodes + Real x = (ii * dx[0] - xcen); + // Real y = (jj * dx[1] - ycen); + + // WoA Hill in x-direction + Real height = num / (x*x + 4 * a * a); + + // Populate terrain height + z_arr(i,j,k0) = height; + }); + } } } diff --git a/Exec/RegTests/Terrain2d_Cylinder/prob.cpp b/Exec/RegTests/Terrain2d_Cylinder/prob.cpp index 948d8c9b2..f99f0e224 100644 --- a/Exec/RegTests/Terrain2d_Cylinder/prob.cpp +++ b/Exec/RegTests/Terrain2d_Cylinder/prob.cpp @@ -87,54 +87,64 @@ void Problem::init_custom_terrain( const Geometry& geom, MultiFab& z_phys_nd, - const Real& /*time*/) + const Real& time) { - // Domain cell size and real bounds - auto dx = geom.CellSizeArray(); - auto ProbLoArr = geom.ProbLoArray(); - auto ProbHiArr = geom.ProbHiArray(); - - // Domain valid box (z_nd is nodal) - const amrex::Box& domain = geom.Domain(); - int domlo_x = domain.smallEnd(0); int domhi_x = domain.bigEnd(0) + 1; - // int domlo_y = domain.smallEnd(1); int domhi_y = domain.bigEnd(1) + 1; - int domlo_z = domain.smallEnd(2); - - // User function parameters - Real a = 0.5; - Real num = 8 * a * a * a; - Real xcen = 0.5 * (ProbLoArr[0] + ProbHiArr[0]); - Real ycen = 0.5 * (ProbLoArr[1] + ProbHiArr[1]); - - // Number of ghost cells - int ngrow = z_phys_nd.nGrow(); - - // Populate bottom plane - int k0 = domlo_z; - - for ( amrex::MFIter mfi(z_phys_nd,amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi ) - { - // Grown box with no z range - amrex::Box xybx = mfi.growntilebox(ngrow); - xybx.setRange(2,0); - - amrex::Array4 const& z_arr = z_phys_nd.array(mfi); - - ParallelFor(xybx, [=] AMREX_GPU_DEVICE (int i, int j, int) { - - // Clip indices for ghost-cells - int ii = amrex::min(amrex::max(i,domlo_x),domhi_x); - - // Location of nodes - Real x = (ii * dx[0] - xcen); - - if(fabs(x)read_custom_terrain(fname,geom,z_phys_nd,time); + } else { + + // Domain cell size and real bounds + auto dx = geom.CellSizeArray(); + auto ProbLoArr = geom.ProbLoArray(); + auto ProbHiArr = geom.ProbHiArray(); + + // Domain valid box (z_nd is nodal) + const amrex::Box& domain = geom.Domain(); + int domlo_x = domain.smallEnd(0); int domhi_x = domain.bigEnd(0) + 1; + // int domlo_y = domain.smallEnd(1); int domhi_y = domain.bigEnd(1) + 1; + int domlo_z = domain.smallEnd(2); + + // User function parameters + Real a = 0.5; + Real num = 8 * a * a * a; + Real xcen = 0.5 * (ProbLoArr[0] + ProbHiArr[0]); + Real ycen = 0.5 * (ProbLoArr[1] + ProbHiArr[1]); + + // Number of ghost cells + int ngrow = z_phys_nd.nGrow(); + + // Populate bottom plane + int k0 = domlo_z; + + for ( amrex::MFIter mfi(z_phys_nd,amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi ) + { + // Grown box with no z range + amrex::Box xybx = mfi.growntilebox(ngrow); + xybx.setRange(2,0); + + amrex::Array4 const& z_arr = z_phys_nd.array(mfi); + + ParallelFor(xybx, [=] AMREX_GPU_DEVICE (int i, int j, int) + { + + // Clip indices for ghost-cells + int ii = amrex::min(amrex::max(i,domlo_x),domhi_x); + + // Location of nodes + Real x = (ii * dx[0] - xcen); + + if(fabs(x)read_custom_terrain(fname,geom,z_phys_nd,time); + } else { + + // Domain cell size and real bounds + auto dx = geom.CellSizeArray(); + auto ProbLoArr = geom.ProbLoArray(); + auto ProbHiArr = geom.ProbHiArray(); + + // Domain valid box (z_nd is nodal) + const amrex::Box& domain = geom.Domain(); + int domlo_x = domain.smallEnd(0); int domhi_x = domain.bigEnd(0) + 1; + int domlo_y = domain.smallEnd(1); int domhi_y = domain.bigEnd(1) + 1; + int domlo_z = domain.smallEnd(2); + + // User function parameters + Real a = 0.5; + Real xcen = 0.5 * (ProbLoArr[0] + ProbHiArr[0]); + Real ycen = 0.5 * (ProbLoArr[1] + ProbHiArr[1]); + + // Number of ghost cells + int ngrow = z_phys_nd.nGrow(); + + // Populate bottom plane + int k0 = domlo_z; + + for ( amrex::MFIter mfi(z_phys_nd,amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi ) + { + // Grown box with no z range + amrex::Box xybx = mfi.growntilebox(ngrow); + xybx.setRange(2,0); - amrex::Array4 const& z_arr = z_phys_nd.array(mfi); + amrex::Array4 const& z_arr = z_phys_nd.array(mfi); - ParallelFor(xybx, [=] AMREX_GPU_DEVICE (int i, int j, int) { + ParallelFor(xybx, [=] AMREX_GPU_DEVICE (int i, int j, int) + { - // Clip indices for ghost-cells - int ii = amrex::min(amrex::max(i,domlo_x),domhi_x); - int jj = amrex::min(amrex::max(j,domlo_y),domhi_y); + // Clip indices for ghost-cells + int ii = amrex::min(amrex::max(i,domlo_x),domhi_x); + int jj = amrex::min(amrex::max(j,domlo_y),domhi_y); - // Location of nodes - Real x = (ii * dx[0] - xcen); - Real y = (jj * dx[1] - ycen); + // Location of nodes + Real x = (ii * dx[0] - xcen); + Real y = (jj * dx[1] - ycen); - if(std::pow(x*x + y*y, 0.5) < a){ - z_arr(i,j,k0) = std::pow(a*a - x*x - y*y , 0.5); - } - else{ - z_arr(i,j,k0) = 0.0; - } + if(std::pow(x*x + y*y, 0.5) < a){ + z_arr(i,j,k0) = std::pow(a*a - x*x - y*y , 0.5); + } + else{ + z_arr(i,j,k0) = 0.0; + } - }); + }); + } } } diff --git a/Exec/RegTests/WitchOfAgnesi/prob.cpp b/Exec/RegTests/WitchOfAgnesi/prob.cpp index 45c64c93d..82e89f3a1 100644 --- a/Exec/RegTests/WitchOfAgnesi/prob.cpp +++ b/Exec/RegTests/WitchOfAgnesi/prob.cpp @@ -4,14 +4,14 @@ using namespace amrex; std::unique_ptr -amrex_probinit( +amrex_probinit ( const amrex_real* /*problo*/, const amrex_real* /*probhi*/) { return std::make_unique(); } -Problem::Problem() +Problem::Problem () { // Parse params amrex::ParmParse pp("prob"); @@ -25,7 +25,7 @@ Problem::Problem() } void -Problem::init_custom_pert( +Problem::init_custom_pert ( const Box& bx, const Box& xbx, const Box& ybx, @@ -35,119 +35,121 @@ Problem::init_custom_pert( Array4 const& x_vel_pert, Array4 const& y_vel_pert, Array4 const& z_vel_pert, - Array4 const& r_hse, - Array4 const& p_hse, + Array4 const& /*r_hse*/, + Array4 const& /*p_hse*/, Array4 const& z_nd, - Array4 const& z_cc, + Array4 const& /*z_cc*/, GeometryData const& geomdata, Array4 const& /*mf_m*/, Array4 const& /*mf_u*/, Array4 const& /*mf_v*/, const SolverChoice& sc) { - const int khi = geomdata.Domain().bigEnd()[2]; + const int khi = geomdata.Domain().bigEnd()[2]; const bool use_moisture = (sc.moisture_type != MoistureType::None); - AMREX_ALWAYS_ASSERT(bx.length()[2] == khi+1); + AMREX_ALWAYS_ASSERT(bx.length()[2] == khi+1); - amrex::ParallelFor(bx, [=, parms=parms] AMREX_GPU_DEVICE(int i, int j, int k) noexcept - { - // Geometry (note we must include these here to get the data on device) - const auto prob_lo = geomdata.ProbLo(); - const auto dx = geomdata.CellSize(); - const Real x = prob_lo[0] + (i + 0.5) * dx[0]; - const Real z = z_cc(i,j,k); + amrex::ParallelFor(bx, [=, parms=parms] AMREX_GPU_DEVICE(int i, int j, int k) noexcept + { + // Set scalar = 0 everywhere + state_pert(i, j, k, RhoScalar_comp) = 0.0; - // Set scalar = 0 everywhere - state_pert(i, j, k, RhoScalar_comp) = 0.0; + if (use_moisture) { + state_pert(i, j, k, RhoQ1_comp) = 0.0; + state_pert(i, j, k, RhoQ2_comp) = 0.0; + } + }); - if (use_moisture) { - state_pert(i, j, k, RhoQ1_comp) = 0.0; - state_pert(i, j, k, RhoQ2_comp) = 0.0; - } - }); - - // Set the x-velocity - amrex::ParallelFor(xbx, [=, parms=parms] AMREX_GPU_DEVICE(int i, int j, int k) noexcept - { - x_vel_pert(i, j, k) = parms.U_0; - }); - - // Set the y-velocity - amrex::ParallelFor(ybx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept - { - y_vel_pert(i, j, k) = 0.0; - }); - - const auto dx = geomdata.CellSize(); - amrex::GpuArray dxInv; - dxInv[0] = 1. / dx[0]; - dxInv[1] = 1. / dx[1]; - dxInv[2] = 1. / dx[2]; - - // Set the z-velocity from impenetrable condition - amrex::ParallelFor(zbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept - { - z_vel_pert(i, j, k) = WFromOmega(i, j, k, 0.0, x_vel_pert, y_vel_pert, z_nd, dxInv); - }); - - amrex::Gpu::streamSynchronize(); + // Set the x-velocity + amrex::ParallelFor(xbx, [=, parms=parms] AMREX_GPU_DEVICE(int i, int j, int k) noexcept + { + x_vel_pert(i, j, k) = parms.U_0; + }); + + // Set the y-velocity + amrex::ParallelFor(ybx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept + { + y_vel_pert(i, j, k) = 0.0; + }); + + const auto dx = geomdata.CellSize(); + amrex::GpuArray dxInv; + dxInv[0] = 1. / dx[0]; + dxInv[1] = 1. / dx[1]; + dxInv[2] = 1. / dx[2]; + + // Set the z-velocity from impenetrable condition + amrex::ParallelFor(zbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept + { + z_vel_pert(i, j, k) = WFromOmega(i, j, k, 0.0, x_vel_pert, y_vel_pert, z_nd, dxInv); + }); + amrex::Gpu::streamSynchronize(); } void -Problem::init_custom_terrain( +Problem::init_custom_terrain ( const Geometry& geom, MultiFab& z_phys_nd, - const Real& /*time*/) + const Real& time) { - // Domain cell size and real bounds - auto dx = geom.CellSizeArray(); - auto ProbLoArr = geom.ProbLoArray(); - auto ProbHiArr = geom.ProbHiArray(); - - // Domain valid box (z_nd is nodal) - const amrex::Box& domain = geom.Domain(); - int domlo_x = domain.smallEnd(0); int domhi_x = domain.bigEnd(0) + 1; - // int domlo_y = domain.smallEnd(1); int domhi_y = domain.bigEnd(1) + 1; - int domlo_z = domain.smallEnd(2); - - // User function parameters - Real a = 0.5; - Real num = 8 * a * a * a; - Real xcen = 0.5 * (ProbLoArr[0] + ProbHiArr[0]); - // Real ycen = 0.5 * (ProbLoArr[1] + ProbHiArr[1]); - - // Number of ghost cells - int ngrow = z_phys_nd.nGrow(); - - // Populate bottom plane - int k0 = domlo_z; - - for ( amrex::MFIter mfi(z_phys_nd,amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi ) - { - // Grown box with no z range - amrex::Box xybx = mfi.growntilebox(ngrow); - xybx.setRange(2,0); - - amrex::Array4 const& z_arr = z_phys_nd.array(mfi); - - ParallelFor(xybx, [=] AMREX_GPU_DEVICE (int i, int j, int) { - - // Clip indices for ghost-cells - int ii = amrex::min(amrex::max(i,domlo_x),domhi_x); - // int jj = amrex::min(amrex::max(j,domlo_y),domhi_y); - - // Location of nodes - Real x = (ii * dx[0] - xcen); - // Real y = (jj * dx[1] - ycen); - - // WoA Hill in x-direction - Real height = num / (x*x + 4 * a * a); - // Populate terrain height - z_arr(i,j,k0) = height; - }); + // Check if a valid csv file exists for the terrain + std::string fname; + amrex::ParmParse pp("erf"); + auto valid_fname = pp.query("terrain_file_name",fname); + if (valid_fname) { + this->read_custom_terrain(fname,geom,z_phys_nd,time); + } else { + // Domain cell size and real bounds + auto dx = geom.CellSizeArray(); + auto ProbLoArr = geom.ProbLoArray(); + auto ProbHiArr = geom.ProbHiArray(); + + // Domain valid box (z_nd is nodal) + const amrex::Box& domain = geom.Domain(); + int domlo_x = domain.smallEnd(0); int domhi_x = domain.bigEnd(0) + 1; + // int domlo_y = domain.smallEnd(1); int domhi_y = domain.bigEnd(1) + 1; + int domlo_z = domain.smallEnd(2); + + // User function parameters + Real a = 0.5; + Real num = 8 * a * a * a; + Real xcen = 0.5 * (ProbLoArr[0] + ProbHiArr[0]); + // Real ycen = 0.5 * (ProbLoArr[1] + ProbHiArr[1]); + + // Number of ghost cells + int ngrow = z_phys_nd.nGrow(); + + // Populate bottom plane + int k0 = domlo_z; + + for ( amrex::MFIter mfi(z_phys_nd,amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi ) + { + // Grown box with no z range + amrex::Box xybx = mfi.growntilebox(ngrow); + xybx.setRange(2,0); + + amrex::Array4 const& z_arr = z_phys_nd.array(mfi); + + ParallelFor(xybx, [=] AMREX_GPU_DEVICE (int i, int j, int) + { + // Clip indices for ghost-cells + int ii = amrex::min(amrex::max(i,domlo_x),domhi_x); + // int jj = amrex::min(amrex::max(j,domlo_y),domhi_y); + + // Location of nodes + Real x = (ii * dx[0] - xcen); + // Real y = (jj * dx[1] - ycen); + + // WoA Hill in x-direction + Real height = num / (x*x + 4 * a * a); + + // Populate terrain height + z_arr(i,j,k0) = height; + }); + } } } diff --git a/Exec/SpongeTest/prob.cpp b/Exec/SpongeTest/prob.cpp index 469ba545a..2cb3b36e2 100644 --- a/Exec/SpongeTest/prob.cpp +++ b/Exec/SpongeTest/prob.cpp @@ -87,54 +87,64 @@ void Problem::init_custom_terrain( const Geometry& geom, MultiFab& z_phys_nd, - const Real& /*time*/) + const Real& time) { - // Domain cell size and real bounds - auto dx = geom.CellSizeArray(); - auto ProbLoArr = geom.ProbLoArray(); - auto ProbHiArr = geom.ProbHiArray(); - - // Domain valid box (z_nd is nodal) - const amrex::Box& domain = geom.Domain(); - int domlo_x = domain.smallEnd(0); int domhi_x = domain.bigEnd(0) + 1; - // int domlo_y = domain.smallEnd(1); int domhi_y = domain.bigEnd(1) + 1; - int domlo_z = domain.smallEnd(2); - - // User function parameters - Real a = 0.5; - Real num = 8 * a * a * a; - Real xcen = 0.5 * (ProbLoArr[0] + ProbHiArr[0]); - Real ycen = 0.5 * (ProbLoArr[1] + ProbHiArr[1]); - - // Number of ghost cells - int ngrow = z_phys_nd.nGrow(); - - // Populate bottom plane - int k0 = domlo_z; - - for ( amrex::MFIter mfi(z_phys_nd,amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi ) - { - // Grown box with no z range - amrex::Box xybx = mfi.growntilebox(ngrow); - xybx.setRange(2,0); - - amrex::Array4 const& z_arr = z_phys_nd.array(mfi); - - ParallelFor(xybx, [=] AMREX_GPU_DEVICE (int i, int j, int) { - - // Clip indices for ghost-cells - int ii = amrex::min(amrex::max(i,domlo_x),domhi_x); - - // Location of nodes - Real x = (ii * dx[0] - xcen); - - if(fabs(x)read_custom_terrain(fname,geom,z_phys_nd,time); + } else { + + // Domain cell size and real bounds + auto dx = geom.CellSizeArray(); + auto ProbLoArr = geom.ProbLoArray(); + auto ProbHiArr = geom.ProbHiArray(); + + // Domain valid box (z_nd is nodal) + const amrex::Box& domain = geom.Domain(); + int domlo_x = domain.smallEnd(0); int domhi_x = domain.bigEnd(0) + 1; + // int domlo_y = domain.smallEnd(1); int domhi_y = domain.bigEnd(1) + 1; + int domlo_z = domain.smallEnd(2); + + // User function parameters + Real a = 0.5; + Real num = 8 * a * a * a; + Real xcen = 0.5 * (ProbLoArr[0] + ProbHiArr[0]); + Real ycen = 0.5 * (ProbLoArr[1] + ProbHiArr[1]); + + // Number of ghost cells + int ngrow = z_phys_nd.nGrow(); + + // Populate bottom plane + int k0 = domlo_z; + + for ( amrex::MFIter mfi(z_phys_nd,amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi ) + { + // Grown box with no z range + amrex::Box xybx = mfi.growntilebox(ngrow); + xybx.setRange(2,0); + + amrex::Array4 const& z_arr = z_phys_nd.array(mfi); + + ParallelFor(xybx, [=] AMREX_GPU_DEVICE (int i, int j, int) + { + + // Clip indices for ghost-cells + int ii = amrex::min(amrex::max(i,domlo_x),domhi_x); + + // Location of nodes + Real x = (ii * dx[0] - xcen); + + if(fabs(x) m_xterrain,m_yterrain,m_zterrain; + amrex::Real value1,value2,value3; + while(file>>value1>>value2>>value3){ + m_xterrain.push_back(value1); + m_yterrain.push_back(value2); + m_zterrain.push_back(value3); + } + file.close(); + + // Copy data to the GPU + int ncell = m_xterrain.size(); + amrex::Gpu::DeviceVector d_xterrain(ncell),d_yterrain(ncell),d_zterrain(ncell); + amrex::Gpu::copy(amrex::Gpu::hostToDevice, m_xterrain.begin(), m_xterrain.end(), d_xterrain.begin()); + amrex::Gpu::copy(amrex::Gpu::hostToDevice, m_yterrain.begin(), m_yterrain.end(), d_yterrain.begin()); + amrex::Gpu::copy(amrex::Gpu::hostToDevice, m_zterrain.begin(), m_zterrain.end(), d_zterrain.begin()); + amrex::Real* d_xt = d_xterrain.data(); + amrex::Real* d_yt = d_yterrain.data(); + amrex::Real* d_zt = d_zterrain.data(); + + // Populate z_phys data + amrex::Real tol = 1.0e-4; + int ngrow = z_phys_nd.nGrow(); + auto dx = geom.CellSizeArray(); + auto ProbLoArr = geom.ProbLoArray(); + int ilo = geom.Domain().smallEnd(0); + int jlo = geom.Domain().smallEnd(1); + int klo = geom.Domain().smallEnd(2); + int ihi = geom.Domain().bigEnd(0); + int jhi = geom.Domain().bigEnd(1); + for (amrex::MFIter mfi(z_phys_nd,amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi ) + { + // Grown box with no z range + amrex::Box xybx = mfi.growntilebox(ngrow); + xybx.setRange(2,0); + + amrex::Array4 const& z_arr = z_phys_nd.array(mfi); + amrex::ParallelFor(xybx, [=] AMREX_GPU_DEVICE (int i, int j, int /*k*/) + { + // Clip indices for ghost-cells + int ii = amrex::min(amrex::max(i,ilo),ihi); + int jj = amrex::min(amrex::max(j,jlo),jhi); + + // Location of nodes + bool found = false; + amrex::Real x = ProbLoArr[0] + ii * dx[0]; + amrex::Real y = ProbLoArr[1] + jj * dx[1]; + amrex::Real zloc = 0.0; + for(int n=0; n Date: Tue, 25 Jun 2024 14:30:18 -0700 Subject: [PATCH 8/9] Limit surface roughness to on land value. (#1659) --- Source/BoundaryConditions/MOSTStress.H | 15 +++++++++------ 1 file changed, 9 insertions(+), 6 deletions(-) diff --git a/Source/BoundaryConditions/MOSTStress.H b/Source/BoundaryConditions/MOSTStress.H index 272213257..b15bd2410 100644 --- a/Source/BoundaryConditions/MOSTStress.H +++ b/Source/BoundaryConditions/MOSTStress.H @@ -256,8 +256,8 @@ struct adiabatic_wave_coupled u_star_arr(i,j,k) = mdata.kappa * umm_arr(i,j,k) / std::log(mdata.zref / z0_arr(i,j,k)); do { ustar = u_star_arr(i,j,k); - z0 = std::max(1200.0 * Hwave_arr(i,j,k) * std::pow( Hwave_arr(i,j,k)/(Lwave_arr(i,j,k)+eps), 4.5 ) - + 0.11 * eta_arr(ie,je,k,EddyDiff::Mom_v) / ustar, z0_eps); + z0 = std::min( std::max(1200.0 * Hwave_arr(i,j,k) * std::pow( Hwave_arr(i,j,k)/(Lwave_arr(i,j,k)+eps), 4.5 ) + + 0.11 * eta_arr(ie,je,k,EddyDiff::Mom_v) / ustar, z0_eps), z0_max ); u_star_arr(i,j,k) = mdata.kappa * umm_arr(i,j,k) / std::log(mdata.zref / z0); ++iter; } while ((std::abs(u_star_arr(i,j,k) - ustar) > tol) && iter <= max_iters); @@ -273,6 +273,7 @@ private: const amrex::Real tol = 1.0e-5; const amrex::Real eps = 1e-15; const amrex::Real z0_eps = 1.0e-6; + const amrex::Real z0_max = 0.1; }; @@ -515,8 +516,8 @@ struct surface_flux_wave_coupled u_star_arr(i,j,k) = mdata.kappa * umm_arr(i,j,k) / std::log(mdata.zref / z0_arr(i,j,k)); do { ustar = u_star_arr(i,j,k); - z0 = std::max(1200.0 * Hwave_arr(i,j,k) * std::pow( Hwave_arr(i,j,k)/(Lwave_arr(i,j,k)+eps), 4.5 ) - + 0.11 * eta_arr(ie,je,k,EddyDiff::Mom_v) / ustar, z0_eps); + z0 = std::min( std::max(1200.0 * Hwave_arr(i,j,k) * std::pow( Hwave_arr(i,j,k)/(Lwave_arr(i,j,k)+eps), 4.5 ) + + 0.11 * eta_arr(ie,je,k,EddyDiff::Mom_v) / ustar, z0_eps), z0_max ); Olen = -ustar * ustar * ustar * tm_arr(i,j,k) / (mdata.kappa * mdata.gravity * mdata.surf_temp_flux); zeta = mdata.zref / Olen; @@ -539,6 +540,7 @@ private: const amrex::Real tol = 1.0e-5; const amrex::Real eps = 1e-15; const amrex::Real z0_eps = 1.0e-6; + const amrex::Real z0_max = 0.1; }; @@ -788,8 +790,8 @@ struct surface_temp_wave_coupled u_star_arr(i,j,k) = mdata.kappa * umm_arr(i,j,k) / std::log(mdata.zref / z0_arr(i,j,k)); do { ustar = u_star_arr(i,j,k); - z0 = std::max(1200.0 * Hwave_arr(i,j,k) * std::pow( Hwave_arr(i,j,k)/(Lwave_arr(i,j,k)+eps), 4.5 ) - + 0.11 * eta_arr(ie,je,k,EddyDiff::Mom_v) / ustar, z0_eps); + z0 = std::min( std::max(1200.0 * Hwave_arr(i,j,k) * std::pow( Hwave_arr(i,j,k)/(Lwave_arr(i,j,k)+eps), 4.5 ) + + 0.11 * eta_arr(ie,je,k,EddyDiff::Mom_v) / ustar, z0_eps), z0_max ); tflux = -(tm_arr(i,j,k) - t_surf_arr(i,j,k)) * ustar * mdata.kappa / (std::log(mdata.zref / z0) - psi_h); Olen = -ustar * ustar * ustar * tm_arr(i,j,k) / @@ -813,6 +815,7 @@ private: const amrex::Real tol = 1.0e-5; const amrex::Real eps = 1e-15; const amrex::Real z0_eps = 1.0e-6; + const amrex::Real z0_max = 0.1; }; From aa98170cc61e8e53eda49bfe1d1c1e339c2222ea Mon Sep 17 00:00:00 2001 From: "Jean M. Sexton" Date: Wed, 26 Jun 2024 19:21:53 -0400 Subject: [PATCH 9/9] Add WW3 submodule to point to mpmd branch (#1660) --- .gitmodules | 3 +++ Submodules/WW3 | 1 + 2 files changed, 4 insertions(+) create mode 160000 Submodules/WW3 diff --git a/.gitmodules b/.gitmodules index 7a95f042f..92e06008c 100644 --- a/.gitmodules +++ b/.gitmodules @@ -7,3 +7,6 @@ [submodule "Submodules/RRTMGP"] path = Submodules/RRTMGP url = https://github.com/E3SM-Project/rte-rrtmgp +[submodule "Submodules/WW3"] + path = Submodules/WW3 + url = https://github.com/erf-model/WW3 diff --git a/Submodules/WW3 b/Submodules/WW3 new file mode 160000 index 000000000..f5fa58228 --- /dev/null +++ b/Submodules/WW3 @@ -0,0 +1 @@ +Subproject commit f5fa5822849e33976cdc935a91b304c60bc2a976