From 2851437de54787174878d0923b2c2c9f33d44fcf Mon Sep 17 00:00:00 2001 From: japham0 <156237070+japham0@users.noreply.github.com> Date: Thu, 20 Jun 2024 19:14:24 -0700 Subject: [PATCH] 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 --- CMake/BuildERFExe.cmake | 1 + Exec/ABL/inputs_mpmd | 68 +++++++++++++++ Exec/Make.ERF | 4 + Source/ERF.H | 7 ++ Source/ERF.cpp | 13 ++- Source/ERF_make_new_arrays.cpp | 26 ++++++ Source/ERF_read_waves.cpp | 107 ++++++++++++++++++++++++ Source/Make.package | 4 + Source/TimeIntegration/ERF_TimeStep.cpp | 4 + Source/main.cpp | 19 ++++- 10 files changed, 251 insertions(+), 2 deletions(-) create mode 100644 Exec/ABL/inputs_mpmd create mode 100644 Source/ERF_read_waves.cpp diff --git a/CMake/BuildERFExe.cmake b/CMake/BuildERFExe.cmake index 412c1d926..5f76b5d72 100644 --- a/CMake/BuildERFExe.cmake +++ b/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 diff --git a/Exec/ABL/inputs_mpmd b/Exec/ABL/inputs_mpmd new file mode 100644 index 000000000..d03df316a --- /dev/null +++ b/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 + diff --git a/Exec/Make.ERF b/Exec/Make.ERF index 7efb8b1fa..8f81e8092 100644 --- a/Exec/Make.ERF +++ b/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) diff --git a/Source/ERF.H b/Source/ERF.H index d92977148..4cb59e790 100644 --- a/Source/ERF.H +++ b/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& state_old, @@ -716,6 +720,9 @@ private: // Wave coupling data amrex::Vector> Hwave; amrex::Vector> Lwave; + amrex::Vector> Hwave_onegrid; + amrex::Vector> Lwave_onegrid; + bool finished_wave = false; // array of flux registers amrex::Vector advflux_reg; diff --git a/Source/ERF.cpp b/Source/ERF.cpp index a884317f8..26faab742 100644 --- a/Source/ERF.cpp +++ b/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). diff --git a/Source/ERF_make_new_arrays.cpp b/Source/ERF_make_new_arrays.cpp index 08ccf435a..7e2e5d561 100644 --- a/Source/ERF_make_new_arrays.cpp +++ b/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 pmap; + pmap.resize(1); + pmap[0]=0; + DistributionMapping dm_onegrid(ba2d_onegrid); + dm_onegrid.define(pmap); + + Hwave_onegrid[lev] = std::make_unique(ba2d_onegrid,dm_onegrid,1,IntVect(1,1,0)); + Lwave_onegrid[lev] = std::make_unique(ba2d_onegrid,dm_onegrid,1,IntVect(1,1,0)); + Hwave[lev] = std::make_unique(ba2d,dm,1,IntVect(3,3,0)); + Lwave[lev] = std::make_unique(ba2d,dm,1,IntVect(3,3,0)); + std::cout< +#include +#include +#include + +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 my_H_arr = Hwave_onegrid[lev]->array(mfi); + amrex::Array4 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()<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 + diff --git a/Source/Make.package b/Source/Make.package index 119eb03f0..b62909334 100644 --- a/Source/Make.package +++ b/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 diff --git a/Source/TimeIntegration/ERF_TimeStep.cpp b/Source/TimeIntegration/ERF_TimeStep.cpp index 2212ffcbf..4084498b7 100644 --- a/Source/TimeIntegration/ERF_TimeStep.cpp +++ b/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]); diff --git a/Source/main.cpp b/Source/main.cpp index a79f97c75..2f473b85f 100644 --- a/Source/main.cpp +++ b/Source/main.cpp @@ -11,6 +11,11 @@ #include #endif +#ifdef ERF_USE_WW3_COUPLING +#include +#include +#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 }