Replace manual interp with utils function

Exawind · Nov 13, 2024 · 3d33c9f · 3d33c9f
1 parent e843c6c
commit 3d33c9f
Show file tree

Hide file tree

Showing 4 changed files with 16 additions and 52 deletions.
diff --git a/amr-wind/equation_systems/icns/source_terms/BodyForce.cpp b/amr-wind/equation_systems/icns/source_terms/BodyForce.cpp
@@ -121,10 +121,9 @@ void BodyForce::operator()(
     const auto& nph_time = 0.5 * (m_time.current_time() + m_time.new_time());
     const auto& problo = m_mesh.Geom(lev).ProbLoArray();
     const auto& dx = m_mesh.Geom(lev).CellSizeArray();
-    const int lp1 = lev + 1;
-    const int nh_max = (int)m_prof_x.size() - 2;
 
     const amrex::Real* force_ht = m_ht.data();
+    const amrex::Real* force_ht_end = m_ht.data() + m_ht.size();
     const amrex::Real* force_x = m_prof_x.data();
     const amrex::Real* force_y = m_prof_y.data();
 
@@ -134,19 +133,10 @@ void BodyForce::operator()(
             bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
                 amrex::IntVect iv(i, j, k);
                 const amrex::Real ht = problo[2] + (iv[2] + 0.5) * dx[2];
-                const int il = amrex::min(k / lp1, nh_max);
-                const int ir = il + 1;
-                amrex::Real fx;
-                amrex::Real fy;
-
-                fx = force_x[il] + ((force_x[ir] - force_x[il]) /
-                                    (force_ht[ir] - force_ht[il])) *
-                                       (ht - force_ht[il]);
-
-                fy = force_y[il] + ((force_y[ir] - force_y[il]) /
-                                    (force_ht[ir] - force_ht[il])) *
-                                       (ht - force_ht[il]);
-
+                const amrex::Real fx = amr_wind::interp::linear(
+                    force_ht, force_ht_end, force_x, ht);
+                const amrex::Real fy = amr_wind::interp::linear(
+                    force_ht, force_ht_end, force_y, ht);
                 src_term(i, j, k, 0) += fx;
                 src_term(i, j, k, 1) += fy;
             });

diff --git a/amr-wind/equation_systems/icns/source_terms/HurricaneForcing.cpp b/amr-wind/equation_systems/icns/source_terms/HurricaneForcing.cpp
@@ -4,6 +4,7 @@
 #include "amr-wind/core/vs/vstraits.H"
 #include "amr-wind/wind_energy/ABL.H"
 #include "amr-wind/core/FieldUtils.H"
+#include "amr-wind/utilities/linear_interpolation.H"
 
 #include "AMReX_ParmParse.H"
 #include "AMReX_Gpu.H"
@@ -61,32 +62,20 @@ void HurricaneForcing::operator()(
     const amrex::Real f = m_coriolis_factor;
 
     // Mean velocity profile used to compute background hurricane forcing term
-    //
-    // Assumes that the velocity profile is at the cell-centers of the finest
-    // level grid. For finer meshes, it will extrapolate beyond the Level0
-    // cell-centers for the lo/hi cells.
 
     const int idir = m_axis;
-    const int nh_max = static_cast<int>(m_vel_ht.size()) - 2;
-    const int lp1 = lev + 1;
     const amrex::Real* heights = m_vel_ht.data();
+    const amrex::Real* heights_end = m_vel_ht.data() + m_vel_ht.size();
     const amrex::Real* vals = m_vel_vals.data();
 
     amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
         amrex::IntVect iv(i, j, k);
         const amrex::Real ht = problo[idir] + (iv[idir] + 0.5) * dx[idir];
 
-        const int il = amrex::min(k / lp1, nh_max);
-        const int ir = il + 1;
-
         const amrex::Real umean =
-            vals[3 * il + 0] + ((vals[3 * ir + 0] - vals[3 * il + 0]) /
-                                (heights[ir] - heights[il])) *
-                                   (ht - heights[il]);
+            amr_wind::interp::linear(heights, heights_end, vals, ht, 3, 0);
         const amrex::Real vmean =
-            vals[3 * il + 1] + ((vals[3 * ir + 1] - vals[3 * il + 1]) /
-                                (heights[ir] - heights[il])) *
-                                   (ht - heights[il]);
+            amr_wind::interp::linear(heights, heights_end, vals, ht, 3, 1);
 
         // Gradient velocities varying with height
         const amrex::Real V_z = V * (Vzh - ht) / Vzh;

diff --git a/amr-wind/equation_systems/temperature/source_terms/BodyForce.cpp b/amr-wind/equation_systems/temperature/source_terms/BodyForce.cpp
@@ -1,6 +1,7 @@
 #include "amr-wind/equation_systems/temperature/source_terms/BodyForce.H"
 #include "amr-wind/CFDSim.H"
 #include "amr-wind/utilities/trig_ops.H"
+#include "amr-wind/utilities/linear_interpolation.H"
 
 #include "AMReX_ParmParse.H"
 #include "AMReX_Gpu.H"
@@ -72,10 +73,9 @@ void BodyForce::operator()(
 
     const auto& problo = m_mesh.Geom(lev).ProbLoArray();
     const auto& dx = m_mesh.Geom(lev).CellSizeArray();
-    const int lp1 = lev + 1;
-    const int nh_max = (int)m_prof_theta.size() - 2;
 
     const amrex::Real* force_ht = m_ht.data();
+    const amrex::Real* force_ht_end = m_ht.data() + m_ht.size();
     const amrex::Real* force_theta = m_prof_theta.data();
 
     if (m_type == "height_varying" || m_type == "height-varying") {
@@ -84,14 +84,8 @@ void BodyForce::operator()(
             bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
                 amrex::IntVect iv(i, j, k);
                 const amrex::Real ht = problo[2] + (iv[2] + 0.5) * dx[2];
-                const int il = amrex::min(k / lp1, nh_max);
-                const int ir = il + 1;
-                amrex::Real ftheta;
-
-                ftheta =
-                    force_theta[il] + ((force_theta[ir] - force_theta[il]) /
-                                       (force_ht[ir] - force_ht[il])) *
-                                          (ht - force_ht[il]);
+                const amrex::Real ftheta = amr_wind::interp::linear(
+                    force_ht, force_ht_end, force_theta, ht);
 
                 src_term(i, j, k, 0) += ftheta;
             });

diff --git a/amr-wind/equation_systems/temperature/source_terms/HurricaneTempForcing.cpp b/amr-wind/equation_systems/temperature/source_terms/HurricaneTempForcing.cpp
@@ -4,6 +4,7 @@
 #include "amr-wind/core/vs/vstraits.H"
 #include "amr-wind/wind_energy/ABL.H"
 #include "amr-wind/core/FieldUtils.H"
+#include "amr-wind/utilities/linear_interpolation.H"
 
 #include "AMReX_ParmParse.H"
 #include "AMReX_Gpu.H"
@@ -44,34 +45,24 @@ void HurricaneTempForcing::operator()(
     const amrex::Real dTzh = m_dTzh;
 
     // Mean velocity profile used to compute background hurricane forcing term
-    //
-    // Assumes that the velocity profile is at the cell-centers of the finest
-    // level grid. For finer meshes, it will extrapolate beyond the Level0
-    // cell-centers for the lo/hi cells.
 
     const int idir = m_axis;
-    const int nh_max = static_cast<int>(m_vel_ht.size()) - 2;
-    const int lp1 = lev + 1;
     const amrex::Real* heights = m_vel_ht.data();
+    const amrex::Real* heights_end = m_vel_ht.data() + m_vel_ht.size();
     const amrex::Real* vals = m_vel_vals.data();
 
     amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
         amrex::IntVect iv(i, j, k);
         const amrex::Real ht = problo[idir] + (iv[idir] + 0.5) * dx[idir];
 
-        const int il = amrex::min(k / lp1, nh_max);
-        const int ir = il + 1;
-
         /*const amrex::Real umean =
             vals[3 * il + 0] + ((vals[3 * ir + 0] - vals[3 * il + 0]) /
                                 (heights[ir] - heights[il])) *
                                    (ht - heights[il]);
         */
         const amrex::Real dTdR_z = dTdR * (dTzh - ht) / dTzh;
         const amrex::Real vmean =
-            vals[3 * il + 1] + ((vals[3 * ir + 1] - vals[3 * il + 1]) /
-                                (heights[ir] - heights[il])) *
-                                   (ht - heights[il]);
+            amr_wind::interp::linear(heights, heights_end, vals, ht, 3, 1);
 
         src_term(i, j, k) -= vmean * dTdR_z;
     });