Merge branch 'development' into dg/erf_particles

Docs/sphinx_doc/CouplingToWW3.rst

@@ -0,0 +1,48 @@
+
+ .. role:: cpp(code)
+    :language: c++
+
+ .. _CouplingToWW3:
+
+Coupling To WW3
+===============
+
+Coupling with WaveWatch III is currently a work in progress.
+Currently, we have a one-way coupling between ERF and WaveWatch III (WW3), where WW3 sends ERF Hwave (significant wave height) and Lwave (mean wavelength) over a grid.
+
+One-way coupling WW3 to ERF
+---------------------------
+
+The values are used to compute the surface roughness :math:`\overline{z_{0}}` through a fixed-point iteration:
+
+.. math::
+  z_{0} = 1200.0 Hwave (\frac{Hwave}{Lwave})^{4.5} + \frac{0.11 \mu}{u_*}
+
+To run the coupled model:
+
+.. code-block:: bash
+
+    git clone --recursive [email protected]:erf-model/ERF
+    cd ERF/Exec/ABL
+    make -j4 USE_WW3_COUPLING=TRUE
+    cd ../../Submodules/WW3
+    ./model/bin/w3_setup model -c gnu -s Ifremer1
+    cd regtests
+    ./bin/run_cmake_test -C MPMD -n 2 -p mpirun -f -s PR1_MPI ../model ww3_tp2.2
+
+Modifications to the problem size and geometry, as well as other parameters can be done in the `inputs_mpmd` file. The `plt` files as well as relevant outputs can be viewed in the `regtests/ww3_tp2.2/work` directory.
+
+Two-way coupling:
+-----------------
+
+Disclaimer: Two-way coupling is currently a work in progress. Two-way coupling involves sending the wind velocity and direction to WW3. We convert the x and y velocities from ERF to a wind speed (:math:`U_{wind}`) and wind direction (:math:`\theta`) from the reference height.
+
+.. math::
+
+  U_{wind} = \sqrt{u^{2} + v^2}
+
+.. math::
+
+  \theta = \mathrm{arctan}{\frac{v}{u}}
+
+Both :math:`U_{wind}` and :math:`\theta` are then used in the wind source term :math:`S_{in}` in the ST6 subroutine in WW3

Docs/sphinx_doc/index.rst

@@ -74,10 +74,11 @@ In addition to this documentation, there is API documentation for ERF generated
    Visualization.rst
 
 .. toctree::
-   :caption: COUPLING TO AMR-WIND
+   :caption: EXTERNAL COUPLING
    :maxdepth: 1
 
    CouplingToAMRWind.rst
+   CouplingToWW3.rst
 
 .. toctree::
    :caption: ERF vs WRF

Exec/ABL/inputs_mpmd

@@ -6,8 +6,8 @@ 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.prob_extent =  191  91    1024
+amr.n_cell           =  191  91       4
 
 geometry.is_periodic = 1 1 0
 

Source/ERF.H

@@ -241,6 +241,7 @@ public:
 
 #ifdef ERF_USE_WW3_COUPLING
     void read_waves (int lev);
+    void send_waves (int lev);
 #endif
 
     // Interface for advancing the data at one level by one "slow" timestep

Source/ERF.cpp

@@ -909,6 +909,7 @@ ERF::InitData ()
 #ifdef ERF_USE_WW3_COUPLING
     int lev = 0;
     read_waves(lev);
+    send_waves(lev);
 #endif
 
     // Configure ABLMost params if used MostWall boundary condition

Source/ERF_read_waves.cpp

@@ -4,6 +4,8 @@
 #include <Utils.H>
 #include <mpi.h>
 #include <AMReX_MPMD.H>
+#include <AMReX_MultiFab.H>
+#include <AMReX_Geometry.H>
 
 using namespace amrex;
 
@@ -69,6 +71,9 @@ ERF::read_waves (int lev)
                  }
              }
 
+             // amrex::AllPrintToFile("output_HS_cpp.txt")<<FArrayBox(my_H_arr)<<std::endl;
+             // amrex::AllPrintToFile("output_L_cpp.txt")<<FArrayBox(my_L_arr)<<std::endl;
+
          }
     }
 
@@ -93,6 +98,149 @@ ERF::read_waves (int lev)
         });
     }
 
+    for (MFIter mfi(*Hwave[lev],TilingIfNotGPU()); mfi.isValid(); ++mfi) {
+        Box bx = mfi.tilebox();
+        const Array4<Real const>& Hwave_arr = Hwave[lev]->const_array(mfi);
+        const Array4<int>& 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;
+            }
+        });
+    }
+
+}
+
+void
+ERF::send_waves (int lev)
+{
+    int ncomp = 1; // number components
+    auto& lev_new = vars_new[lev];
+    const double PI = 3.1415926535897932384626433832795028841971693993751058209;
+
+    // Access xvel, yvel from ABL
+    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());
+
+    // Make local copy of xvel, yvel
+    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());
+
+
+    MultiFab x_avg(lev_new[Vars::cons].boxArray(), lev_new[Vars::cons].DistributionMap(),
+                       ncomp, lev_new[Vars::cons].nGrow());
+    MultiFab y_avg(lev_new[Vars::cons].boxArray(), lev_new[Vars::cons].DistributionMap(),
+                       ncomp, lev_new[Vars::cons].nGrow());
+
+    MultiFab u_mag(lev_new[Vars::cons].boxArray(), lev_new[Vars::cons].DistributionMap(),
+                       ncomp, lev_new[Vars::cons].nGrow());
+
+    MultiFab u_dir(lev_new[Vars::cons].boxArray(), lev_new[Vars::cons].DistributionMap(),
+                       ncomp, lev_new[Vars::cons].nGrow());
+    x_avg.setVal(0.);
+    y_avg.setVal(0.);
+    u_mag.setVal(0.);
+    u_dir.setVal(0.);
+
+    for (MFIter mfi(x_avg,TilingIfNotGPU()); mfi.isValid(); ++mfi) {
+
+        Box bx = mfi.tilebox();
+        const Array4<Real>& u_vel = x_avg.array(mfi);
+        const Array4<const Real>& velx_arr = xvel_data.array(mfi);
+
+        ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k){
+            u_vel(i,j,k)  = 0.5 *( velx_arr(i,j,k) + velx_arr(i+1,j,k) );
+
+            // amrex::AllPrintToFile("uvel.txt") << amrex::IntVect(i,j,k) << " [" <<velx_arr(i,j,k) << "| avg:  " << u_vel(i,j,k)<< " | " << velx_arr(i+1,j,k) << "]" <<std::endl;
+        });
+    }
+
+    for (MFIter mfi(y_avg,TilingIfNotGPU()); mfi.isValid(); ++mfi) {
+
+        Box bx = mfi.tilebox();
+        const Array4<Real>& v_vel = y_avg.array(mfi);
+        const Array4<const Real>& vely_arr = yvel_data.array(mfi);
+
+        ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k){
+
+            v_vel(i,j,k)  = 0.5 *( vely_arr(i,j,k) + vely_arr(i,j+1,k) );
+
+            // amrex::AllPrintToFile("vvel.txt") << amrex::IntVect(i,j,k) << " ["<<vely_arr(i,j,k)<<"| avg: " << v_vel(i,j,k)<< " | " << vely_arr(i,j+1,k) << "]" <<std::endl;
+        });
+    }
+
+    for (MFIter mfi(u_mag, TilingIfNotGPU()); mfi.isValid(); ++mfi) {
+
+        Box bx = mfi.tilebox();
+        const Array4<Real>& magnitude = u_mag.array(mfi);
+        const Array4<const Real>& u = x_avg.array(mfi);
+        const Array4<const Real>& v = y_avg.array(mfi);
+        const Array4<Real>& theta = u_dir.array(mfi);
+
+        ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k){
+
+            magnitude(i,j,k)  = std::sqrt( pow(u(i,j,k), 2) + pow(v(i,j,k), 2) );
+
+            double u_val = u(i, j, k);
+            double v_val = v(i, j, k);
+            if ( u_val == 0 ) {
+                u_val = std::max( u_val, 1e-15 );  // Ensure u_val is non-zero
+            }
+
+
+            if ( u_val < 0 && v_val > 0 || u_val < 0 && v_val < 0 ) {
+
+                theta(i,j,k) = PI + ( atan( v_val / u_val ) );
+
+            } else {
+
+                theta(i,j,k) = atan ( v_val / u_val );
+            }
+
+
+            // amrex::AllPrintToFile("mag_theta.txt") << amrex::IntVect(i,j,k) <<  " Magnitude: " << magnitude(i,j,k) << " Theta: " << theta(i,j,k) <<std::endl;
+        });
+
+
+            // MPI_Send to WW3
+    // Calculate the number of elements in the current box
+    int n_elements = bx.numPts();
+
+    // Get the data pointers for the arrays
+    Real* magnitude_ptr = &magnitude(bx.smallEnd());
+    Real* theta_ptr = &theta(bx.smallEnd());
+
+    // Initialize other_root as needed
+    int other_root = 0; // Example initialization, replace with appropriate logic
+
+    // Send magnitude array
+//    MPI_Send(magnitude_ptr, n_elements, MPI_DOUBLE, other_root, 0, MPI_COMM_WORLD);
+    // Send theta array
+//    MPI_Send(theta_ptr, n_elements, MPI_DOUBLE, other_root, 1, MPI_COMM_WORLD);
+// amrex::AllPrintToFile("debug_send.txt") << n_elements << std::endl;
+
+}
+/*
+            for (MFIter mfi(u_mag); mfi.isValid(); ++mfi) {
+
+                const Array4<Real>& magnitude = u_mag.array(mfi);
+                const Array4<Real>& theta = u_dir.array(mfi);
+
+                Real* magnitude_ptr = magnitude.dataPtr();
+                Real* theta_ptr = theta.dataPtr();
+
+                // Get number of elements in arrays
+                int n_elements = mfi.validbox().numPts();
+                int this_root = 0;
+                int other_root = 1;
+
+                // Send magnitude and theta
+                MPI_Send(magnitude_ptr. n_elements, MPI_DOUBLE, this_root, 0, MPI_COMM_WORLD)
+                MPI_Send(theta_ptr, n_elements, MPI_DOUBLE, other_root, 1, MPI_COMM_WORLD);
+            }
+*/
 }
 #endif
 

Source/TimeIntegration/ERF_TimeStep.cpp

@@ -65,6 +65,7 @@ ERF::timeStep (int lev, Real time, int /*iteration*/)
 
 #ifdef ERF_USE_WW3_COUPLING
     read_waves(lev);
+    send_waves(lev);
 #endif
 
     // Advance a single level for a single time step