diff --git a/Docs/sphinx_doc/MOST.rst b/Docs/sphinx_doc/MOST.rst index bc1ed0be4..4498fa57c 100644 --- a/Docs/sphinx_doc/MOST.rst +++ b/Docs/sphinx_doc/MOST.rst @@ -157,6 +157,14 @@ To evaluate the fluxes with MOST, the surface rougness parameter :math:`z_{0}` m If the ``charnock`` method is employed, the :math:`a` constant may be specified with ``erf.most.charnock_constant`` (defaults to 0.0185). If the ``modified_charnock`` method is employed, the depth :math:`d` may be specified with ``erf.most.modified_charnock_depth`` (defaults to 30 m). If the ``wave_coupled`` method is employed, the user must provide wave height and mean wavelength data. +While the MOST methods relevant to air-sea interfaces (``charnock``, ``modified_charnock``, and ``wave_coupled``) dynamically compute :math:`z_{0}`, the ``constant`` case can also allow for spatially varying :math:`z_{0}` if an inputs file is employed. More specifically, one may specify + +:: + + erf.most.roughness_file_name = STRING #Name of file that contains (x,y,z_0) + +in the inputs file and ERF will populate the 2D :math:`z_{0}` array with values contained in the text file. + When computing an average :math:`\overline{\phi}` for the MOST boundary, where :math:`\phi` denotes a generic variable, ERF supports a variety of approaches. Specifically, ``planar averages`` and ``local region averages`` may be computed with or without ``time averaging``. With each averaging methodology, the query point :math:`z` may be determined from the following procedures: specified vertical distance :math:`z_{ref}` from the bottom surface, specified :math:`k_{index}`, or (when employing terrain-fit coordinates) specified normal vector length :math:`z_{ref}`. The available inputs to the MOST boundary and their associated data types are :: diff --git a/Source/BoundaryConditions/ABLMost.H b/Source/BoundaryConditions/ABLMost.H index c486fc6d6..108c22311 100644 --- a/Source/BoundaryConditions/ABLMost.H +++ b/Source/BoundaryConditions/ABLMost.H @@ -6,6 +6,7 @@ #include #include #include +#include #include #include @@ -153,10 +154,20 @@ public: pp.query("most.modified_charnock_depth",depth); } else if (rough_sea_string == "wave_coupled") { rough_type_sea = RoughCalcType::WAVE_COUPLED; + } else if (rough_sea_string == "constant") { + rough_type_sea = RoughCalcType::CONSTANT; } else { amrex::Abort("Undefined MOST roughness type for sea!"); } + // Check if there is a user-specified roughness file to be read + std::string fname; + bool read_z0 = false; + if ( (flux_type == FluxCalcType::MOENG) || + (flux_type == FluxCalcType::ROTATE) ) { + read_z0 = pp.query("most.roughness_file_name", fname); + } + // Size the MOST params for all levels int nlevs = m_geom.size(); z_0.resize(nlevs); @@ -225,11 +236,16 @@ public: // Z0 heights FAB //-------------------------------------------------------- + amrex::Arena* Arena_Used = amrex::The_Arena(); +#ifdef AMREX_USE_GPU + Arena_Used = amrex::The_Pinned_Arena(); +#endif amrex::Box bx = amrex::grow(m_geom[lev].Domain(),ng); bx.setSmall(2,0); bx.setBig(2,0); - z_0[lev].resize(bx,1); + z_0[lev].resize(bx,1,Arena_Used); z_0[lev].setVal(z0_const); + if (read_z0) read_custom_roughness(lev,fname); // 2D MFs for U*, T*, T_surf //-------------------------------------------------------- @@ -340,6 +356,10 @@ public: amrex::MultiFab* z_phys_cc, const PBLHeightEstimator& est); + void + read_custom_roughness (const int& lev, + const std::string& fname); + void update_surf_temp (const amrex::Real& time) { diff --git a/Source/BoundaryConditions/ABLMost.cpp b/Source/BoundaryConditions/ABLMost.cpp index 72c88bef0..95bf5fa1e 100644 --- a/Source/BoundaryConditions/ABLMost.cpp +++ b/Source/BoundaryConditions/ABLMost.cpp @@ -613,3 +613,102 @@ void ABLMost::calc_wstar(const int lev, }); } } + +ABLMost::read_custom_roughness (const int& lev, + const std::string& fname) +{ + // Read the file if we are on the coarsest level + if (lev==0) { + // Only the ioproc reads the file and broadcasts + if (ParallelDescriptor::IOProcessor()) { + Print()<<"Reading MOST roughness file: "<< fname << std::endl; + std::ifstream file(fname); + Gpu::HostVector m_x,m_y,m_z0; + Real value1,value2,value3; + while(file>>value1>>value2>>value3){ + m_x.push_back(value1); + m_y.push_back(value2); + m_z0.push_back(value3); + } + file.close(); + + // Copy data to the GPU + int nnode = m_x.size(); + Gpu::DeviceVector d_x(nnode),d_y(nnode),d_z0(nnode); + Gpu::copy(Gpu::hostToDevice, m_x.begin(), m_x.end(), d_x.begin()); + Gpu::copy(Gpu::hostToDevice, m_y.begin(), m_y.end(), d_y.begin()); + Gpu::copy(Gpu::hostToDevice, m_z0.begin(), m_z0.end(), d_z0.begin()); + Real* xp = d_x.data(); + Real* yp = d_y.data(); + Real* z0p = d_z0.data(); + + // Populate z_phys data + Real tol = 1.0e-4; + auto dx = m_geom[lev].CellSizeArray(); + auto ProbLoArr = m_geom[lev].ProbLoArray(); + int ilo = m_geom[lev].Domain().smallEnd(0); + int jlo = m_geom[lev].Domain().smallEnd(1); + int klo = 0; + int ihi = m_geom[lev].Domain().bigEnd(0); + int jhi = m_geom[lev].Domain().bigEnd(1); + + // Grown box with no z range + Box xybx = z_0[lev].box(); + xybx.setRange(2,0); + + Array4 const& z0_arr = z_0[lev].array(); + 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 + Real x = ProbLoArr[0] + ii * dx[0]; + Real y = ProbLoArr[1] + jj * dx[1]; + int inode = ii + jj * (ihi-ilo+2); // stride is Nx+1 + if (std::sqrt(std::pow(x-xp[inode],2)+std::pow(y-yp[inode],2)) < tol) { + z0_arr(i,j,klo) = z0p[inode]; + } else { + // Unexpected list order, do brute force search + Real z0loc = 0.0; + bool found = false; + for (int n=0; n bc_lo(3,ERFBCType::foextrap); + Vector bc_hi(3,ERFBCType::foextrap); + Vector bcr; bcr.push_back(BCRec(bc_lo.data(),bc_hi.data())); + + // Create ref ratio + IntVect ratio; + for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) { + ratio[idim] = m_geom[lev].Domain().length(idim) / m_geom[0].Domain().length(idim); + } + + // Create interp object and interpolate from the coarsest grid + Interpolater* interp = &cell_cons_interp; + interp->interp(z_0[0] , 0, + z_0[lev], 0, + 1, z_0[lev].box(), + ratio, m_geom[0], m_geom[lev], + bcr, 0, 0, RunOn::Gpu); + } +}