first pass at coupling with WW3 (#1653)

* first pass at coupling with WW3

* Move the ifdef so we don't have to worry about finding mpi.h

---------

Co-authored-by: AMLattanzi <AMLattanzi@lbl.gov>

CMake/BuildERFExe.cmake

@@ -107,6 +107,7 @@ function(build_erf_lib erf_lib_name)
        ${SRC_DIR}/ERF.cpp
        ${SRC_DIR}/ERF_make_new_arrays.cpp
        ${SRC_DIR}/ERF_make_new_level.cpp
+       ${SRC_DIR}/ERF_read_waves.cpp
        ${SRC_DIR}/ERF_Tagging.cpp
        ${SRC_DIR}/Advection/AdvectionSrcForMom.cpp
        ${SRC_DIR}/Advection/AdvectionSrcForState.cpp

Exec/ABL/inputs_mpmd

@@ -0,0 +1,68 @@
+# ------------------  INPUTS TO MAIN PROGRAM  -------------------
+max_step = 5
+
+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.is_periodic = 1 1 0
+
+zlo.type = "NoSlipWall"
+zhi.type = "SlipWall"
+
+# TIME STEP CONTROL
+erf.fixed_dt           = 10  # fixed time step depending on grid resolution
+erf.max_step=1
+# DIAGNOSTICS & VERBOSITY
+erf.sum_interval   = 1       # timesteps between computing mass
+erf.v              = 1       # verbosity in ERF.cpp
+amr.v                = 1       # verbosity in Amr.cpp
+
+# REFINEMENT / REGRIDDING
+amr.max_level       = 0       # maximum level number allowed
+
+# CHECKPOINT FILES
+erf.check_file      = chk        # root name of checkpoint file
+erf.check_int       = 100        # number of timesteps between checkpoints
+
+# PLOTFILES
+erf.plot_file_1     = plt        # prefix of plotfile name
+erf.plot_int_1      = 10         # number of timesteps between plotfiles
+erf.plot_vars_1     = density rhoadv_0 x_velocity y_velocity z_velocity pressure temp theta
+
+# SOLVER CHOICE
+erf.alpha_T = 0.0
+erf.alpha_C = 1.0
+erf.use_gravity = false
+
+erf.molec_diff_type = "None"
+erf.les_type        = "Smagorinsky"
+erf.Cs              = 0.1
+
+erf.init_type = "uniform"
+
+# MOST BOUNDARY (DEFAULT IS ADIABATIC FOR THETA)
+zlo.type      = "Most"
+# erf.most.z0   = 0.1  Don't set z0
+erf.most.zref = 128.0 # set 1/2 of the vertical dz: 1024 / 4 = 256 -> 256/2 = 128
+erf.most.roughness_type_sea = wave_coupled
+
+# PROBLEM PARAMETERS
+prob.rho_0 = 1.0
+prob.A_0 = 1.0
+
+prob.U_0 = 10.0
+prob.V_0 = 0.0
+prob.W_0 = 0.0
+prob.T_0 = 300.0
+
+# Higher values of perturbations lead to instability
+# Instability seems to be coming from BC
+prob.U_0_Pert_Mag = 0.08
+prob.V_0_Pert_Mag = 0.08 #
+prob.W_0_Pert_Mag = 0.0
+

Exec/Make.ERF

@@ -155,6 +155,10 @@ ifeq ($(USE_WINDFARM), TRUE)
 	INCLUDE_LOCATIONS += $(ERF_WINDFARM_EWP_DIR)
 endif
 
+ifeq ($(USE_WW3_COUPLING), TRUE)
+	DEFINES += -DERF_USE_WW3_COUPLING
+endif
+
 ERF_LSM_DIR = $(ERF_SOURCE_DIR)/LandSurfaceModel
 include $(ERF_LSM_DIR)/Make.package
 VPATH_LOCATIONS   += $(ERF_LSM_DIR)

Source/ERF.H

@@ -238,6 +238,10 @@ public:
     // compute dt from CFL considerations
     amrex::Real estTimeStep (int lev, long& dt_fast_ratio) const;
 
+#ifdef ERF_USE_WW3_COUPLING
+    void read_waves (int lev);
+#endif
+
     // Interface for advancing the data at one level by one "slow" timestep
     void advance_dycore (int level,
                          amrex::Vector<amrex::MultiFab>& state_old,
@@ -716,6 +720,9 @@ private:
     // Wave coupling data
     amrex::Vector<std::unique_ptr<amrex::MultiFab>> Hwave;
     amrex::Vector<std::unique_ptr<amrex::MultiFab>> Lwave;
+    amrex::Vector<std::unique_ptr<amrex::MultiFab>> Hwave_onegrid;
+    amrex::Vector<std::unique_ptr<amrex::MultiFab>> Lwave_onegrid;
+    bool finished_wave = false;
 
     // array of flux registers
     amrex::Vector<amrex::YAFluxRegister*> advflux_reg;

Source/ERF.cpp

@@ -276,7 +276,13 @@ ERF::ERF ()
         Hwave[lev] = nullptr;
         Lwave[lev] = nullptr;
     }
-
+    Hwave_onegrid.resize(nlevs_max);
+    Lwave_onegrid.resize(nlevs_max);
+    for (int lev = 0; lev < max_level; ++lev)
+    {
+        Hwave_onegrid[lev] = nullptr;
+        Lwave_onegrid[lev] = nullptr;
+    }
     // Theta prim for MOST
     Theta_prim.resize(nlevs_max);
 
@@ -900,6 +906,11 @@ ERF::InitData ()
         }
     }
 
+#ifdef ERF_USE_WW3_COUPLING
+    int lev = 0;
+    read_waves(lev);
+#endif
+
     // Configure ABLMost params if used MostWall boundary condition
     // NOTE: we must set up the MOST routine after calling FillPatch
     //       in order to have lateral ghost cells filled (MOST + terrain interp).

Source/ERF_make_new_arrays.cpp

@@ -245,6 +245,32 @@ ERF::init_stuff (int lev, const BoxArray& ba, const DistributionMapping& dm,
     }
 #endif
 
+
+#ifdef ERF_USE_WW3_COUPLING
+    // create a new BoxArray and DistributionMapping for a MultiFab with 1 box
+    BoxArray ba_onegrid(geom[lev].Domain());
+    BoxList bl2d_onegrid = ba_onegrid.boxList();
+    for (auto& b : bl2d_onegrid) {
+        b.setRange(2,0);
+    }
+    BoxArray ba2d_onegrid(std::move(bl2d_onegrid));
+    Vector<int> pmap;
+    pmap.resize(1);
+    pmap[0]=0;
+    DistributionMapping dm_onegrid(ba2d_onegrid);
+    dm_onegrid.define(pmap);
+
+    Hwave_onegrid[lev] = std::make_unique<MultiFab>(ba2d_onegrid,dm_onegrid,1,IntVect(1,1,0));
+    Lwave_onegrid[lev] = std::make_unique<MultiFab>(ba2d_onegrid,dm_onegrid,1,IntVect(1,1,0));
+    Hwave[lev] = std::make_unique<MultiFab>(ba2d,dm,1,IntVect(3,3,0));
+    Lwave[lev] = std::make_unique<MultiFab>(ba2d,dm,1,IntVect(3,3,0));
+    std::cout<<ba_onegrid<<std::endl;
+    std::cout<<ba2d_onegrid<<std::endl;
+    std::cout<<dm_onegrid<<std::endl;
+    std::cout<<dm_onegrid<<std::endl;
+#endif
+
+
 #if defined(ERF_USE_RRTMGP)
     //*********************************************************
     // Radiation heating source terms

Source/ERF_read_waves.cpp

@@ -0,0 +1,107 @@
+#ifdef ERF_USE_WW3_COUPLING
+
+#include <ERF.H>
+#include <Utils.H>
+#include <mpi.h>
+#include <AMReX_MPMD.H>
+
+using namespace amrex;
+
+void
+ERF::read_waves (int lev)
+{
+    for ( MFIter mfi(*Hwave_onegrid[lev],false); mfi.isValid(); ++mfi)
+    {
+         //auto my_H_ptr = Hwave[lev]->array(mfi);
+         //auto my_L_ptr = Lwave[lev]->array(mfi);
+         const auto & bx = mfi.validbox();
+
+        amrex::Print() <<  " HERE " << bx << std::endl;
+         // How to declare my_H_ptr directly?
+         amrex::Array4<Real> my_H_arr = Hwave_onegrid[lev]->array(mfi);
+         amrex::Array4<Real> my_L_arr = Lwave_onegrid[lev]->array(mfi);
+
+         Real* my_H_ptr = my_H_arr.dataPtr();
+         Real* my_L_ptr = my_L_arr.dataPtr();
+
+         int rank_offset = amrex::MPMD::MyProc() - amrex::ParallelDescriptor::MyProc();
+         int this_root, other_root;
+         if (rank_offset == 0) { // First program
+             this_root = 0;
+             other_root = amrex::ParallelDescriptor::NProcs();
+         } else {
+             this_root = rank_offset;
+             other_root = 0;
+         }
+
+
+         int nx=2147483647; // sanity check
+         int ny=2147483647; // sanity check
+
+         //JUST RECV
+         if (amrex::MPMD::MyProc() == this_root) {
+             if (rank_offset == 0) // the first program
+             {
+                 MPI_Recv(&nx, 1, MPI_INT, other_root, 1, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+                 MPI_Recv(&ny, 1, MPI_INT, other_root, 7, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+             }
+             else // the second program
+             {
+                 MPI_Recv(&nx, 1, MPI_INT, other_root, 0, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+                 MPI_Recv(&ny, 1, MPI_INT, other_root, 6, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+             }
+             //This may not be necessary
+             ParallelDescriptor::Bcast(&nx, 1);
+             ParallelDescriptor::Bcast(&ny, 1);
+         }
+         if((nx)*(ny) > 0) {
+             int nsealm = (nx)*ny;
+             Print()<<nsealm<<std::endl;
+             Print()<<" NX = " << nx<<std::endl;
+             Print()<<" Ny = " << ny<<std::endl;
+             Print()<<" BX  NPTS = " << bx.numPts() <<std::endl;
+             Print()<<" BX       = " << bx <<std::endl;
+             // Print()<<" FAB NPTS = " << Hwave_onegrid[0]->array(mfi).size() <<std::endl;
+             Print()<<" FAB NPTS = " << (*Hwave_onegrid[0])[0].box().numPts()<<std::endl;
+             Print()<<" FAB      = " << (*Hwave_onegrid[0])[0].box() <<std::endl;
+             // AMREX_ALWAYS_ASSERT_WITH_MESSAGE((nx)*ny <= bx.numPts() || bx.numPts() == 0, "total number of points being filled exceeds the size of the current box\n");
+
+             if (amrex::MPMD::MyProc() == this_root) {
+                 if (rank_offset == 0) // the first program
+                 {
+                     MPI_Recv(my_H_ptr, nsealm, MPI_DOUBLE, other_root, 3, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+                     MPI_Recv(my_L_ptr, nsealm, MPI_DOUBLE, other_root, 5, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+                 }
+                 else // the second program
+                 {
+                     MPI_Recv(my_H_ptr, nsealm, MPI_DOUBLE, other_root, 2, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+                     MPI_Recv(my_L_ptr, nsealm, MPI_DOUBLE, other_root, 4, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+                 }
+             }
+             amrex::Print()<<"Just recieved "<<nsealm<<"as a double*"<<std::endl;
+
+             if(bx.contains(IntVect(192,2,0))) {
+                 std::cout<<my_H_arr(192,92,0)<<std::endl;
+                 std::cout<<my_L_arr(192,92,0)<<std::endl;
+                 std::cout<<my_H_ptr[192-0+(92-0)*193]<<std::endl;
+                 std::cout<<my_L_ptr[192-0+(92-0)*193]<<std::endl;
+             }
+             amrex::AllPrintToFile("output_HS_cpp.txt")<<FArrayBox(my_H_arr)<<std::endl;
+             amrex::AllPrintToFile("output_L_cpp.txt")<<FArrayBox(my_L_arr)<<std::endl;
+         } else {
+             finished_wave = true;
+         }
+    }
+    //May need to be Redistribute
+    //    ParallelCopy(Hwave[lev],Hwave_onegrid[lev],0,0,1,0);
+    Hwave[lev]->ParallelCopy(*Hwave_onegrid[lev]);
+    Hwave[lev]->FillBoundary(geom[lev].periodicity());
+    Lwave[lev]->ParallelCopy(*Lwave_onegrid[lev]);
+    Lwave[lev]->FillBoundary(geom[lev].periodicity());
+    amrex::Print() << "HWAVE BOX " << (*Hwave[lev])[0].box() << std::endl;
+
+    print_state(*Hwave[lev],IntVect(103,-3,0),0,IntVect(3,3,0));
+    print_state(*Hwave[lev],IntVect(103,88,0),0,IntVect(3,3,0));
+}
+#endif
+

Source/Make.package

@@ -14,3 +14,7 @@ CEXE_sources += ERF_make_new_level.cpp
 CEXE_sources += ERF_make_new_arrays.cpp
 CEXE_sources += Derive.cpp
 CEXE_headers += Derive.H
+
+ifeq ($(USE_WW3_COUPLING),TRUE)
+CEXE_sources += ERF_read_waves.cpp
+endif

Source/TimeIntegration/ERF_TimeStep.cpp

@@ -63,6 +63,10 @@ ERF::timeStep (int lev, Real time, int /*iteration*/)
                        << " with dt = " << dt[lev] << std::endl;
     }
 
+#ifdef ERF_USE_WW3_COUPLING
+    read_waves(lev);
+#endif
+
     // Advance a single level for a single time step
     Advance(lev, time, dt[lev], istep[lev], nsubsteps[lev]);
 

Source/main.cpp

@@ -11,6 +11,11 @@
 #include <MultiBlockContainer.H>
 #endif
 
+#ifdef ERF_USE_WW3_COUPLING
+#include <mpi.h>
+#include <AMReX_MPMD.H>
+#endif
+
 std::string inputs_name;
 
 using namespace amrex;
@@ -55,6 +60,7 @@ void add_par () {
 */
 int main (int argc, char* argv[])
 {
+
 #ifdef AMREX_USE_MPI
     MPI_Init(&argc, &argv);
 #endif
@@ -96,7 +102,12 @@ int main (int argc, char* argv[])
             }
         }
     }
+#ifdef ERF_USE_WW3_COUPLING
+    MPI_Comm comm = amrex::MPMD::Initialize(argc, argv);
+    amrex::Initialize(argc,argv,true,comm,add_par);
+#else
     amrex::Initialize(argc,argv,true,MPI_COMM_WORLD,add_par);
+#endif
 
     // Save the inputs file name for later.
     if (!strchr(argv[1], '=')) {
@@ -228,9 +239,15 @@ int main (int argc, char* argv[])
 
     // destroy timer for profiling
     BL_PROFILE_VAR_STOP(pmain);
-
+#ifdef ERF_USE_WW3_COUPLING
+    MPI_Barrier(MPI_COMM_WORLD);
+#endif
     amrex::Finalize();
 #ifdef AMREX_USE_MPI
+#ifdef ERF_USE_WW3_COUPLING
+    amrex::MPMD::Finalize();
+#else
     MPI_Finalize();
 #endif
+#endif
 }