diff --git a/Docs/sphinx_doc/CouplingToWW3.rst b/Docs/sphinx_doc/CouplingToWW3.rst new file mode 100644 index 000000000..5a9934d0c --- /dev/null +++ b/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 git@github.com: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 diff --git a/Docs/sphinx_doc/index.rst b/Docs/sphinx_doc/index.rst index 178d0b7fe..7f73dec2a 100644 --- a/Docs/sphinx_doc/index.rst +++ b/Docs/sphinx_doc/index.rst @@ -79,6 +79,12 @@ In addition to this documentation, there is API documentation for ERF generated CouplingToAMRWind.rst +.. toctree:: + :caption: ERF vs WRF + :maxdepth: 1 + + CouplingToWW3.rst + .. toctree:: :caption: ERF vs WRF :maxdepth: 1 diff --git a/Exec/ABL/inputs_mpmd b/Exec/ABL/inputs_mpmd index d03df316a..a4ffbd1dc 100644 --- a/Exec/ABL/inputs_mpmd +++ b/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 diff --git a/Source/ERF.H b/Source/ERF.H index 926d1d183..4ca77a03a 100644 --- a/Source/ERF.H +++ b/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 diff --git a/Source/ERF.cpp b/Source/ERF.cpp index 7fff9ac0e..fa9eb9ca3 100644 --- a/Source/ERF.cpp +++ b/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 diff --git a/Source/ERF_read_waves.cpp b/Source/ERF_read_waves.cpp index ec9032992..d180a3e7b 100644 --- a/Source/ERF_read_waves.cpp +++ b/Source/ERF_read_waves.cpp @@ -4,6 +4,8 @@ #include #include #include +#include +#include using namespace amrex; @@ -69,6 +71,9 @@ ERF::read_waves (int lev) } } + // amrex::AllPrintToFile("output_HS_cpp.txt")<& 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; + } + }); + } + +} + +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& u_vel = x_avg.array(mfi); + const Array4& 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) << " [" <& v_vel = y_avg.array(mfi); + const Array4& 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) << " ["<& magnitude = u_mag.array(mfi); + const Array4& u = x_avg.array(mfi); + const Array4& v = y_avg.array(mfi); + const Array4& 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) <& magnitude = u_mag.array(mfi); + const Array4& 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 diff --git a/Source/TimeIntegration/ERF_TimeStep.cpp b/Source/TimeIntegration/ERF_TimeStep.cpp index 4084498b7..f8c6baa59 100644 --- a/Source/TimeIntegration/ERF_TimeStep.cpp +++ b/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