Read z0 from inputs file. (erf-model#1757)

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
 
 ::

Source/BoundaryConditions/ABLMost.H

@@ -137,10 +137,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);
@@ -207,11 +217,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<amrex::RunOn::Device>(z0_const);
+            if (read_z0) read_custom_roughness(lev,fname);
 
             // 2D MFs for U*, T*, T_surf
             //--------------------------------------------------------
@@ -297,6 +312,10 @@ public:
     void
     get_lsm_tsurf (const int& lev);
 
+    void
+    read_custom_roughness (const int& lev,
+                           const std::string& fname);
+
     void
     update_surf_temp (const amrex::Real& time)
     {

Source/BoundaryConditions/ABLMost.cpp

@@ -1,5 +1,7 @@
 #include <ABLMost.H>
 #include <MOSTAverage.H>
+#include <IndexDefines.H>
+#include <AMReX_Interpolater.H>
 
 using namespace amrex;
 
@@ -541,3 +543,103 @@ ABLMost::get_lsm_tsurf (const int& lev)
         });
     }
 }
+
+void
+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<Real> 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<Real> 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<Real> 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<nnode; ++n) {
+                        Real delta=std::sqrt(std::pow(x-xp[n],2)+std::pow(y-yp[n],2));
+                        if (delta < tol) {
+                            found = true;
+                            z0loc = z0p[n];
+                            break;
+                        }
+                    }
+                    AMREX_ASSERT_WITH_MESSAGE(found, "Location read from terrain file does not match the grid!");
+                    amrex::ignore_unused(found);
+                    z0_arr(i,j,klo) = z0loc;
+                }
+            });
+        } // Is ioproc
+
+        int ioproc = ParallelDescriptor::IOProcessorNumber();
+        ParallelDescriptor::Barrier();
+        ParallelDescriptor::Bcast(z_0[lev].dataPtr(),z_0[lev].box().numPts(),ioproc);
+    } else {
+        // Create a BC mapper that uses FOEXTRAP at domain bndry
+        Vector<int> bc_lo(3,ERFBCType::foextrap);
+        Vector<int> bc_hi(3,ERFBCType::foextrap);
+        Vector<BCRec> 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);
+    }
+}