From 5ce79bb568662a1a4b2898b7074869f2fddf98bc Mon Sep 17 00:00:00 2001 From: japham0 <156237070+japham0@users.noreply.github.com> Date: Fri, 21 Jun 2024 16:51:13 -0700 Subject: [PATCH 1/5] Updates to read land/sea mask (#1656) * changes to read land/sea mask * changes to read land/sea mask --- Source/ERF_read_waves.cpp | 51 ++++++++++++++++----------------------- 1 file changed, 21 insertions(+), 30 deletions(-) diff --git a/Source/ERF_read_waves.cpp b/Source/ERF_read_waves.cpp index 9af9410b4..ec9032992 100644 --- a/Source/ERF_read_waves.cpp +++ b/Source/ERF_read_waves.cpp @@ -12,12 +12,10 @@ 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::Print() << " HERE " << bx << std::endl; amrex::Array4 my_H_arr = Hwave_onegrid[lev]->array(mfi); amrex::Array4 my_L_arr = Lwave_onegrid[lev]->array(mfi); @@ -35,17 +33,17 @@ ERF::read_waves (int lev) } - int nx=2147483647; // sanity check + int nx=2147483647; int ny=2147483647; // sanity check - //JUST RECV + //JUST RECEIVED if (amrex::MPMD::MyProc() == this_root) { - if (rank_offset == 0) // the first program + if (rank_offset == 0) // 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 + else // 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); @@ -54,17 +52,9 @@ ERF::read_waves (int lev) ParallelDescriptor::Bcast(&nx, 1); ParallelDescriptor::Bcast(&ny, 1); } + if((nx)*(ny) > 0) { int nsealm = (nx)*ny; - Print()<ParallelCopy(*Hwave_onegrid[lev]); @@ -99,9 +80,19 @@ ERF::read_waves (int lev) Lwave[lev]->ParallelCopy(*Lwave_onegrid[lev]); Lwave[lev]->FillBoundary(geom[lev].periodicity()); amrex::Print() << "HWAVE BOX " << (*Hwave[lev])[0].box() << std::endl; + for (MFIter mfi(*Hwave[lev],TilingIfNotGPU()); mfi.isValid(); ++mfi) { + Box bx = mfi.tilebox(); + const Array4& 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; + } + }); + } - 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 From a0ea7699865abb63bbe44bc504ba444e8d131b76 Mon Sep 17 00:00:00 2001 From: "Dustin (Ting-Hsuan) Ma" Date: Mon, 24 Jun 2024 09:25:09 -0700 Subject: [PATCH 2/5] Turbulent inflow generation technique (#1657) * Box perturbation setup start * perturbation box implementation start * added geometry into CalcTurbPert() * Using index of box to created boxes * Index Perturbation box initialization working * added Perturbation box initialization in init_stuff(), next create class for PerturbationRegion * created 3 functions init_PerturbationRegion(), calc_TurbPert_amplitude(), and calc_TurbPert_updateTime. Now PB boxes initialize at t=0s * cleaned up some portions of the code. Now start implemented t>0s PB * start of struct TurbPert implementation * created TurbPert structure and compiles/runs fines. Now t>=0 uses same function to create turbulent perturbations. * name change to TurbulentPerturbation structure * made small modifications to perturbation funcion call. Start playing with momentum perturbation. * Averaging paraview out put through impermeability enforcement. Added in random number generator in Perturbation amplitude. * Cleaned up some code, and added comments for tomorrow start. Tomorrow: Implement update-time and amplitude * distribution mapping debug w/ Aaron * debugging mf global access * mf now is seen outside of struct * moving on PB implementation to using vectors for per box unique value. * Replaced all perturbation boxes from boxArray and vector to MultiFab. Current version works with np=1, but not np>1. * fixed parallelization problem. Succesfully accessed velocity at initialization. Moving on to calculate update time/think about how to calculate perturbation amplitude on the fly. * Created elapsed time array for mf in tpi. Start writing perturbation amplitude calculations. * mirgrated over to ssh github check * PB u_mean calc per bx using box union (volume averaged). PB items uses vectors. PB frequency triggeres amplitude calculation. PB apply_tpi applied random number perturbation. * Added in Ma&senocak2023 PB amp formulation. Made PB velocity average thread safe. Created some initial constant values for PB. * debugging ubx output (start) * Bug fixed: object.makeSlab() will modify the object, while makeSlab(object,...) will make a copy then edit copy. * Generalized slab average formulation. Can begin some kind of small test to perturb flow. General setup complete * (June 18) TurbulentPerturbation struct now own all computation for PB. Slab average implemented. Start runtime tests. * deleted trailing white space to pass github check * Modified some names, made sure CMake compile. * Changes in attemping to pass workout test * TurbPertStruct.H doesn't like ParallelFor within it? * TurbPertStruct.H ParallelFor loops now using pass-by-value so GPU does not access illegal memory * Update TurbPertStruct.H cuda pass-by- workout test * quick fix for CUDA error * attemping to fix oneAPI SYCL workflow error * Inflow perturbation 1st PR. Cleared all github workflow error * flipped order of solverChoice.pert_type to avoid level error. ready for 1st PR --------- Co-authored-by: Ting-Hsuan (Dustin) Ma Co-authored-by: Aaron M. Lattanzi <103702284+AMLattanzi@users.noreply.github.com> --- CMake/BuildERFExe.cmake | 1 + Exec/ABL/GNUmakefile | 2 +- Source/DataStructs/DataStruct.H | 27 +- Source/DataStructs/Make.package | 1 + Source/DataStructs/TurbPertStruct.H | 328 ++++++++++++++++++ Source/ERF.H | 9 + Source/ERF.cpp | 9 + Source/ERF_make_new_arrays.cpp | 13 + Source/Initialization/ERF_init_TurbPert.cpp | 116 +++++++ Source/Initialization/ERF_init_custom.cpp | 1 - Source/Initialization/Make.package | 1 + Source/SourceTerms/ERF_make_sources.cpp | 15 +- Source/SourceTerms/Src_headers.H | 4 +- Source/TimeIntegration/ERF_Advance.cpp | 10 +- Source/TimeIntegration/ERF_advance_dycore.cpp | 3 + Source/TimeIntegration/TI_slow_rhs_fun.H | 6 +- Source/prob_common.H | 10 + 17 files changed, 547 insertions(+), 9 deletions(-) create mode 100644 Source/DataStructs/TurbPertStruct.H create mode 100644 Source/Initialization/ERF_init_TurbPert.cpp diff --git a/CMake/BuildERFExe.cmake b/CMake/BuildERFExe.cmake index 5f76b5d72..34f7640d4 100644 --- a/CMake/BuildERFExe.cmake +++ b/CMake/BuildERFExe.cmake @@ -140,6 +140,7 @@ function(build_erf_lib erf_lib_name) ${SRC_DIR}/Initialization/ERF_init_from_metgrid.cpp ${SRC_DIR}/Initialization/ERF_init_uniform.cpp ${SRC_DIR}/Initialization/ERF_init1d.cpp + ${SRC_DIR}/Initialization/ERF_init_TurbPert.cpp ${SRC_DIR}/Initialization/ERF_input_sponge.cpp ${SRC_DIR}/IO/Checkpoint.cpp ${SRC_DIR}/IO/ERF_ReadBndryPlanes.cpp diff --git a/Exec/ABL/GNUmakefile b/Exec/ABL/GNUmakefile index c85179014..f26adcb3a 100644 --- a/Exec/ABL/GNUmakefile +++ b/Exec/ABL/GNUmakefile @@ -19,7 +19,7 @@ USE_HIP = FALSE USE_SYCL = FALSE # Debugging -DEBUG = FALSE +DEBUG = TRUE TEST = TRUE USE_ASSERTION = TRUE diff --git a/Source/DataStructs/DataStruct.H b/Source/DataStructs/DataStruct.H index f3ebcb0b4..44a027b15 100644 --- a/Source/DataStructs/DataStruct.H +++ b/Source/DataStructs/DataStruct.H @@ -15,6 +15,7 @@ #include #include #include +#include enum struct ABLDriverType { None, PressureGradient, GeostrophicWind @@ -60,6 +61,11 @@ enum Sponge { ubar_sponge, vbar_sponge, nvars_sponge }; +enum struct PerturbationType { + BPM, CPM, None + // Box perturbation method, Cell perturbation method, None +}; + /** * Container holding many of the algorithmic options and parameters */ @@ -224,12 +230,28 @@ struct SolverChoice { abl_driver_type = ABLDriverType::PressureGradient; } else if (!abl_driver_type_string.compare("GeostrophicWind")) { abl_driver_type = ABLDriverType::GeostrophicWind; - } else if (!abl_driver_type_string.compare("None")){ + } else if (!abl_driver_type_string.compare("None")) { abl_driver_type = ABLDriverType::None; // No ABL driver for simulating classical fluid dynamics problems } else { amrex::Error("Don't know this abl_driver_type"); } + // Which type of inflow turbulent generation + static std::string turb_pert_type_string = "None"; + pp.query("inlet_perturbation_type",turb_pert_type_string); + if (!turb_pert_type_string.compare("BoxPerturbation")) { + pert_type = PerturbationType::BPM; + amrex::Print() << "Box perturbation method selected\n"; + } else if (!turb_pert_type_string.compare("CellPerturbation")) { + pert_type = PerturbationType::CPM; + amrex::Print() << "Cell perturbation method selected\n"; + } else if (!turb_pert_type_string.compare("None")) { + pert_type = PerturbationType::None; + amrex::Print() << "No perturbation method selected\n"; + } else { + amrex::Abort("Dont know this inlet_perturbation_type"); + } + amrex::Vector abl_pressure_grad_in = {0.0, 0.0, 0.0}; pp.queryarr("abl_pressure_grad",abl_pressure_grad_in); for(int i = 0; i < AMREX_SPACEDIM; ++i) abl_pressure_grad[i] = abl_pressure_grad_in[i]; @@ -487,6 +509,9 @@ struct SolverChoice { // User wishes to output time averaged velocity fields bool time_avg_vel = false; + // Type of perturbation + PerturbationType pert_type; + // Numerical diffusion bool use_NumDiff{false}; amrex::Real NumDiffCoeff{0.}; diff --git a/Source/DataStructs/Make.package b/Source/DataStructs/Make.package index b82b7b0ad..c1bc1767f 100644 --- a/Source/DataStructs/Make.package +++ b/Source/DataStructs/Make.package @@ -3,3 +3,4 @@ CEXE_headers += DiffStruct.H CEXE_headers += AdvStruct.H CEXE_headers += SpongeStruct.H CEXE_headers += TurbStruct.H +CEXE_headers += TurbPertStruct.H diff --git a/Source/DataStructs/TurbPertStruct.H b/Source/DataStructs/TurbPertStruct.H new file mode 100644 index 000000000..5d5021248 --- /dev/null +++ b/Source/DataStructs/TurbPertStruct.H @@ -0,0 +1,328 @@ +#ifndef _TURB_PERT_STRUCT_H_ +#define _TURB_PERT_STRUCT_H_ + +#include +#include + +/** + * Container holding quantities related to turbulent perturbation parameters + */ + +struct TurbulentPerturbation { + + public: + + ~TurbulentPerturbation () { + amrex::Print() << "~TurbulentPerturbation destructor\n"; + } + +//#define DIAGONAL_REF +#define LENGTH_REF + +/* void init_tpi_const(const amrex::Real T_0) + { + // Bulk Richardson number + tpi_nonDim = 0.042; + + // Ambient Temperature + tpi_Tinfty = T_0; + }*/ + + // Initializing Perturbation Region + // Note: Box will always be initialized at cell centered + void init_tpi (const int lev, + const amrex::IntVect& nx, + amrex::GpuArray const dx) + + { + // TODO: Use user Param to define this (Currently HARDCODED) + // Perturbation box region setup + int pertBox_offset = 0; + amrex::IntVect regionLo(0+pertBox_offset, 0, 0); + amrex::IntVect regionHi(7+pertBox_offset, nx[1], nx[2]); + amrex::IntVect boxSize(8,8,8); + //amrex::IntVect boxSize(4,4,4); // HACK + + // TODO: This will need to be changed later on + tpi_Lpb = 8*dx[0]; + tpi_Wpb = 8*dx[1]; + tpi_Hpb = 8*dx[2]; + + #ifdef DIAGONAL_REF + tpi_lref = sqrt(boxSize(0)*boxSize(0)*dx[0]*dx[0] + boxSize(1)*boxSize(1)*dx[1]*dx[1]); + #elif defined(LENGTH_REF) + tpi_lref = 8*dx[0]; + #endif + + // Creating structure box array for conserved quantity + amrex::Box tmp_bx(regionLo, regionHi, amrex::IntVect(0,0,0)); + amrex::BoxArray tmp_ba(tmp_bx); + tmp_ba.maxSize(boxSize); + pb_ba.push_back(tmp_ba); + + // Initializing mean magnitude vector + pb_mag.resize(pb_ba[lev].size(), 0.); + + // Set size of vector and initialize + pb_freq.resize(pb_ba[lev].size(), -1.0); + pb_local_etime.resize(pb_ba[lev].size(), 0.0); + pb_amp.resize(pb_ba[lev].size(), 0.0); + + // Function check point message + amrex::Print() << "Turbulent perturbation region initialized with ba " << pb_ba[lev].size() << " boxes\n"; + + // Function method check point message + amrex::Print() << "Using "; + #ifdef DIAGONAL_REF + amrex::Print() << "box diagonal "; + #elif defined(LENGTH_REF) + amrex::Print() << "box length "; + #endif + amrex::Print() << "distance as reference for turbulent inlet perturbation: tpi_lref = " << tpi_lref << "\n"; + } + + // Perturbation update frequency check. + // This function will trigger calc_tpi_meanMag() and calc_tpi_amp(), when + // the local elasped time is greater than the perturbation frequency + void calc_tpi_update (const int lev, + const amrex::Real dt, + amrex::MultiFab& mf_xvel, + amrex::MultiFab& mf_yvel, + amrex::MultiFab& mf_cons) + { + for (int boxIdx = 0; boxIdx < pb_ba[lev].size(); boxIdx++) { + if ( pb_freq[boxIdx] <= pb_local_etime[boxIdx] ) { + // Compute mean velocity magnitude within perturbation box + calc_tpi_meanMag(boxIdx, lev, mf_xvel, mf_yvel, mf_cons); + + // Ma and Senocak (2015) Eq. 15 + if (pb_mag[boxIdx] !=0.) { + pb_freq[boxIdx] = tpi_lref / pb_mag[boxIdx]; + } + + // Trigger amplitude calculation per perturbation box + calc_tpi_amp(boxIdx); + + // Reset local elapsed time + pb_local_etime[boxIdx] = 0.; + } else { + // Increase by timestep + pb_local_etime[boxIdx] += dt; + } // if + } // for + } + +// TODO: Test the difference between these two +//#define USE_SLAB_AVERAGE +#define USE_VOLUME_AVERAGE + + // Perturbation box mean velocity magnitude calculation + // This is pulled into the strucutre to also utilize during runtime + void calc_tpi_meanMag (const int boxIdx, + const int lev, + amrex::MultiFab& mf_xvel, + amrex::MultiFab& mf_yvel, + amrex::MultiFab& mf_cons) + + { + // Creating local copy of PB box array and magnitude + const amrex::BoxArray m_pb_ba = pb_ba[lev]; + amrex::Real* m_pb_mag = get_pb_mag(); + + // Storage of averages per PB + // Index: 0=u (vol/slab_lo), 1=v (vol/slab_lo) + // 2=u (slab_hi), 3=v (slab_hi) + amrex::Gpu::DeviceVector tmp_d(4,0.); + amrex::Real* tmp = tmp_d.data(); + + // Averaging u components + for (amrex::MFIter mfi(mf_xvel, TileNoZ()) ; mfi.isValid(); ++mfi) { + const amrex::Box &vbx = mfi.validbox(); + amrex::Box pbx = convert(m_pb_ba[boxIdx], vbx.ixType()); + amrex::Box ubx = pbx & vbx; + + // Operation over box union + if (ubx.ok()) { + const amrex::Array4 &xvel_arry = mf_xvel.array(mfi); + + #ifdef USE_VOLUME_AVERAGE + int npts = ubx.numPts(); + ParallelFor(amrex::Gpu::KernelInfo().setReduction(true), ubx, [=] + AMREX_GPU_DEVICE(int i, int j, int k, amrex::Gpu::Handler const& handler) noexcept { + amrex::Gpu::deviceReduceSum(&tmp[0], xvel_arry(i,j,k), handler); + }); + tmp[0] /= (amrex::Real) npts; + #endif + + #ifdef USE_SLAB_AVERAGE + amrex::Box ubxSlab_lo = makeSlab(ubx,2,ubx.smallEnd(2)); + amrex::Box ubxSlab_hi = makeSlab(ubx,2,ubx.bigEnd(2)); + int npts_lo = ubxSlab_lo.numPts(); + int npts_hi = ubxSlab_hi.numPts(); + + // Average u in the low slab + ParallelFor(amrex::Gpu::KernelInfo().setReduction(true), ubxSlab_lo, [=] + AMREX_GPU_DEVICE(int i, int j, int k, amrex::Gpu::Handler const& handler) noexcept { + amrex::Gpu::deviceReduceSum(&tmp[0], xvel_arry(i,j,k), handler); + }); + tmp[0] /= (amrex::Real) npts_lo; + + // Average u in the high slab + ParallelFor(amrex::Gpu::KernelInfo().setReduction(true), ubxSlab_hi, [=] + AMREX_GPU_DEVICE(int i, int j, int k, amrex::Gpu::Handler const& handler) noexcept { + amrex::Gpu::deviceReduceSum(&tmp[2], xvel_arry(i,j,k), handler); + }); + tmp[2] /= (amrex::Real) npts_hi; + #endif + } // if + } // MFIter + + // Averaging v components + for (amrex::MFIter mfi(mf_yvel, TileNoZ()); mfi.isValid(); ++mfi) { + amrex::Box const& vbx = mfi.validbox(); + amrex::Box pbx = convert(m_pb_ba[boxIdx], vbx.ixType()); + amrex::Box ubx = pbx & vbx; + + // Operation over box union + if (ubx.ok()) { + const amrex::Array4 &yvel_arry = mf_yvel.array(mfi); + + #ifdef USE_VOLUME_AVERAGE + int npts = ubx.numPts(); + ParallelFor(amrex::Gpu::KernelInfo().setReduction(true), ubx, [=] + AMREX_GPU_DEVICE(int i, int j, int k, amrex::Gpu::Handler const& handler) noexcept { + amrex::Gpu::deviceReduceSum(&tmp[1], yvel_arry(i,j,k), handler); + }); + tmp[1] /= (amrex::Real) npts; + #endif + + #ifdef USE_SLAB_AVERAGE + amrex::Box ubxSlab_lo = makeSlab(ubx,2,ubx.smallEnd(2)); + amrex::Box ubxSlab_hi = makeSlab(ubx,2,ubx.bigEnd(2)); + int npts_lo = ubxSlab_lo.numPts(); + int npts_hi = ubxSlab_hi.numPts(); + + // Average v in the low slab + ParallelFor(amrex::Gpu::KernelInfo().setReduction(true), ubxSlab_lo, [=] + AMREX_GPU_DEVICE(int i, int j, int k, amrex::Gpu::Handler const& handler) noexcept { + amrex::Gpu::deviceReduceSum(&tmp[1], yvel_arry(i,j,k), handler); + }); + tmp[1] /= (amrex::Real) npts_lo; + + // Average v in the high slab + ParallelFor(amrex::Gpu::KernelInfo().setReduction(true), ubxSlab_hi, [=] + AMREX_GPU_DEVICE(int i, int j, int k, amrex::Gpu::Handler const& handler) noexcept { + amrex::Gpu::deviceReduceSum(&tmp[3], yvel_arry(i,j,k), handler); + }); + tmp[3] /= (amrex::Real) npts_hi; + #endif + } // if + } // MFIter + + // Computing the average magnitude within PB + for (amrex::MFIter mfi(mf_cons, TileNoZ()); mfi.isValid(); ++mfi) { + amrex::Box const& vbx = mfi.validbox(); + amrex::Box pbx = convert(m_pb_ba[boxIdx], vbx.ixType()); + amrex::Box ubx = pbx & vbx; + if (ubx.ok()) { + #ifdef USE_SLAB_AVERAGE + m_pb_mag[boxIdx] = 0.5*(sqrt(tmp[0]*tmp[0] + tmp[1]*tmp[1]) + sqrt(tmp[2]*tmp[2] + tmp[3]*tmp[3])); + #endif + + #ifdef USE_VOLUME_AVERAGE + m_pb_mag[boxIdx] = sqrt(tmp[0]*tmp[0] + tmp[1]*tmp[1]); + #endif + } // if + } // MFIter + + } + +//#define ECKERT_FORMULATION +#define BULK_RICHARDSON_FORMULATION + + // Perturbation amplitude calcluation + void calc_tpi_amp (const int boxIdx) + { + amrex::Real g = CONST_GRAV; + amrex::Real beta = 1./tpi_Tinfty; + amrex::Real Um = pb_mag[boxIdx]; + + // Reset the perturbation amplitude + pb_amp[boxIdx] = 0.; + + #ifdef BULK_RICHARDSON_FORMULATION + // Ma and Senocak (2023) Eq. 17 + pb_amp[boxIdx] = tpi_nonDim * Um * Um * Um / (g * beta * tpi_lref * tpi_Hpb); + #endif + + } + +//#define INDEX_PERTURB + + // Applying perturbation amplitude onto source term (Umphrey and Senocak 2016) + // Random perturbation is applied as a white noise on the temperature source term + // it becomes colored through the filtering of temperature transport equation from + // the eddy diffusivity. + void apply_tpi (const int lev, + const amrex::Box& vbx, // box union from upper level + const int comp, // Component to modify + const amrex::IndexType m_ixtype, // IntVect type of src_arr + amrex::Array4 const& src_arr)// Source Array to apply perturbation + { + for (int boxIdx = 0; boxIdx < pb_ba[lev].size(); boxIdx++) { + amrex::Box pbx = convert(pb_ba[lev][boxIdx], m_ixtype); + amrex::Box ubx = pbx & vbx; + if (ubx.ok()) { + amrex::Real amp_copy = pb_amp[boxIdx]; // This exposes a copy so the GPU can capture by value + ParallelForRNG(ubx, [=] AMREX_GPU_DEVICE(int i, int j, int k, const amrex::RandomEngine& engine) noexcept { + // Random number generation between 0 and 1 + amrex::Real rand_double = amrex::Random(engine); + + // Rescale between -1 and 1 and multiple with amplitude + src_arr(i,j,k,comp) = (rand_double*2.0 - 1.0) * amp_copy; + + // For box region debug only + #ifdef INDEX_PERTURB + src_arr(i,j,k,comp) = (amrex::Real) (boxIdx + 1.); + #endif + }); + } + } + } + + + // Quality of life function definitions + amrex::Real* get_pb_mag() { return pb_mag.data(); } + + void debug () + { + for (int i = 0; i < pb_mag.size(); i++) { + amrex::Print() << "[" << i<< "] pb_mag=" << pb_mag[i] + << " | pb_freq=" << pb_freq[i] + << " | pb_local_etime=" << pb_local_etime[i] <<"\n"; + } + } + + // Public constant constants + amrex::Real tpi_Tinfty; // [K] + amrex::Real tpi_nonDim; + + // Public data members + amrex::Vector pb_ba; // PB box array + amrex::Vector pb_mag; // BP mean magnitude [m/s] + + private: + + // Private constant constants + amrex::Real tpi_Hpb; // PB height [m] + amrex::Real tpi_Lpb; // PB length [m] + amrex::Real tpi_Wpb; // PB width [m] + amrex::Real tpi_lref; // PB reference length [m] + + // Storage Arrays + amrex::Vector pb_freq; // PB frequency [1/s] + amrex::Vector pb_local_etime; // PB local elapsed time [s] + amrex::Vector pb_amp; // PB perturbation amplitude Ri:[K] Ec:[K/s] + +}; +#endif diff --git a/Source/ERF.H b/Source/ERF.H index 4cb59e790..926d1d183 100644 --- a/Source/ERF.H +++ b/Source/ERF.H @@ -29,6 +29,7 @@ #include #include +#include #include #include #include @@ -423,6 +424,11 @@ private: amrex::Vector& lev_new, amrex::Vector& lev_old); + // Initialize the Turbulent perturbation + void turbPert_constants (const int lev); + void turbPert_update (const int lev, const amrex::Real dt); + void turbPert_amplitude (const int lev); + // Initialize the integrator object void initialize_integrator (int lev, amrex::MultiFab& cons_mf, amrex::MultiFab& vel_mf); @@ -842,6 +848,9 @@ private: // algorithm choices static SolverChoice solverChoice; + // Turbulent perturbation structure + TurbulentPerturbation turbPert; + #ifdef ERF_USE_PARTICLES // Particle container with all particle species ParticleData particleData; diff --git a/Source/ERF.cpp b/Source/ERF.cpp index 26faab742..7fff9ac0e 100644 --- a/Source/ERF.cpp +++ b/Source/ERF.cpp @@ -1216,6 +1216,15 @@ ERF::init_only (int lev, Real time) input_sponge(lev); } + // Initialize turbulent perturbation + if (solverChoice.pert_type == PerturbationType::BPM) { + if (lev == 0) { + turbPert_constants(lev); + turbPert_update(lev, 0.); + turbPert_amplitude(lev); + } + } + } // read in some parameters from inputs file diff --git a/Source/ERF_make_new_arrays.cpp b/Source/ERF_make_new_arrays.cpp index 7e2e5d561..710abb39e 100644 --- a/Source/ERF_make_new_arrays.cpp +++ b/Source/ERF_make_new_arrays.cpp @@ -304,6 +304,16 @@ ERF::init_stuff (int lev, const BoxArray& ba, const DistributionMapping& dm, } #endif + //********************************************************* + // Turbulent perturbation region initialization + //********************************************************* + // TODO: Test perturbation on multiple levels + if (solverChoice.pert_type == PerturbationType::BPM) { + if (lev == 0) { + turbPert.init_tpi(lev, geom[lev].Domain().bigEnd(), geom[lev].CellSizeArray()); + } + } + // // Define the land mask here and set it to all land // NOTE: the logic below will BREAK if we have any grids not touching the bottom boundary @@ -476,3 +486,6 @@ ERF::initialize_bcs (int lev) (lev, geom[lev], domain_bcs_type, domain_bcs_type_d, m_bc_extdir_vals, m_bc_neumann_vals, use_real_bcs); } + + + diff --git a/Source/Initialization/ERF_init_TurbPert.cpp b/Source/Initialization/ERF_init_TurbPert.cpp new file mode 100644 index 000000000..beb2c924d --- /dev/null +++ b/Source/Initialization/ERF_init_TurbPert.cpp @@ -0,0 +1,116 @@ +/** + * \ERF_init_TurbPert.cpp + */ + +//DUSTIN MA: May 28th, 2024 + +#define DEBUG_PERTBOX_MSG + +#include +#include +#include +#include + +using namespace amrex; + +void +ERF::turbPert_constants(const int lev) +{ + prob->init_turbPert_const(turbPert); +} + +void +ERF::turbPert_update (const int lev, const Real local_dt) +{ + // Grabing data from velocity field + auto& lev_new = vars_new[lev]; + + // Accessing local data + MultiFab cons_data(lev_new[Vars::cons].boxArray(), lev_new[Vars::cons].DistributionMap(), + lev_new[Vars::cons].nComp() , lev_new[Vars::cons].nGrow()); + 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()); + MultiFab::Copy (cons_data, lev_new[Vars::cons], 0, 0, 1, lev_new[Vars::xvel].nGrowVect()); + 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()); + + // Computing perturbation update time + turbPert.calc_tpi_update(lev, local_dt, xvel_data, yvel_data, cons_data); + + #ifdef DEBUG_PERTBOX_MSG + turbPert.debug(); + #endif + + Print() << "Turbulent perturbation update time and amplitude initialized\n"; +} + +// Calculate the perturbation region amplitude. +// This function heavily emmulates the ERF::init_custom () +void +ERF::turbPert_amplitude (int lev) +{ + auto& lev_new = vars_new[lev]; + + MultiFab r_hse(base_state[lev], make_alias, 0, 1); // r_0 is first component + MultiFab p_hse(base_state[lev], make_alias, 1, 1); // p_0 is second component + + MultiFab cons_pert(lev_new[Vars::cons].boxArray(), lev_new[Vars::cons].DistributionMap(), + lev_new[Vars::cons].nComp() , lev_new[Vars::cons].nGrow()); + MultiFab xvel_pert(lev_new[Vars::xvel].boxArray(), lev_new[Vars::xvel].DistributionMap(), 1, lev_new[Vars::xvel].nGrowVect()); + MultiFab yvel_pert(lev_new[Vars::yvel].boxArray(), lev_new[Vars::yvel].DistributionMap(), 1, lev_new[Vars::yvel].nGrowVect()); + MultiFab zvel_pert(lev_new[Vars::zvel].boxArray(), lev_new[Vars::zvel].DistributionMap(), 1, lev_new[Vars::zvel].nGrowVect()); + + // Only storing for conserved quantity for now + auto m_ixtype = cons_pert.boxArray().ixType(); + + // Default perturbations to zero + cons_pert.setVal(0.); + xvel_pert.setVal(0.); + yvel_pert.setVal(0.); + zvel_pert.setVal(0.); + + int fix_random_seed = 0; + ParmParse pp("erf"); pp.query("fix_random_seed", fix_random_seed); + // Note that the value of 1024UL is not significant -- the point here is just to set the + // same seed for all MPI processes for the purpose of regression testing + if (fix_random_seed) { + Print() << "Fixing the random seed" << std::endl; + InitRandom(1024UL); + } + +#ifdef _OPENMP +#pragma omp parallel if (amrex::Gpu::notInLaunchRegion()) +#endif + for (MFIter mfi(lev_new[Vars::cons], TileNoZ()); mfi.isValid(); ++mfi) { + const Box &bx = mfi.validbox(); + const Box &xbx = mfi.tilebox(IntVect(1,0,0)); + const Box &ybx = mfi.tilebox(IntVect(0,1,0)); + const Box &zbx = mfi.tilebox(IntVect(0,0,1)); + + // Perturbation on to different components + const auto &cons_pert_arr = cons_pert.array(mfi); + const auto &xvel_pert_arr = xvel_pert.array(mfi); + const auto &yvel_pert_arr = yvel_pert.array(mfi); + const auto &zvel_pert_arr = zvel_pert.array(mfi); + + Array4 cons_arr = lev_new[Vars::cons].const_array(mfi); + Array4 z_nd_arr = (solverChoice.use_terrain) ? z_phys_nd[lev]->const_array(mfi) : Array4{}; + Array4 z_cc_arr = (solverChoice.use_terrain) ? z_phys_cc[lev]->const_array(mfi) : Array4{}; + + Array4 mf_m = mapfac_m[lev]->array(mfi); + Array4 mf_u = mapfac_m[lev]->array(mfi); + Array4 mf_v = mapfac_m[lev]->array(mfi); + + Array4 r_hse_arr = r_hse.array(mfi); + Array4 p_hse_arr = p_hse.array(mfi); + + turbPert.apply_tpi(lev, bx, RhoTheta_comp, m_ixtype, cons_pert_arr); // Using boxArray + prob->init_custom_pert(bx, xbx, ybx, zbx, cons_arr, cons_pert_arr, + xvel_pert_arr, yvel_pert_arr, zvel_pert_arr, + r_hse_arr, p_hse_arr, z_nd_arr, z_cc_arr, + geom[lev].data(), mf_m, mf_u, mf_v, + solverChoice); + } // mfi + // This initialization adds the initial perturbation direction on the RhoTheta field + MultiFab::Add(lev_new[Vars::cons], cons_pert, RhoTheta_comp, RhoTheta_comp, 1, cons_pert.nGrow()); +} diff --git a/Source/Initialization/ERF_init_custom.cpp b/Source/Initialization/ERF_init_custom.cpp index 98e935d9f..61f3368e8 100644 --- a/Source/Initialization/ERF_init_custom.cpp +++ b/Source/Initialization/ERF_init_custom.cpp @@ -61,7 +61,6 @@ ERF::init_custom (int lev) const Box &ybx = mfi.tilebox(IntVect(0,1,0)); const Box &zbx = mfi.tilebox(IntVect(0,0,1)); - const auto &cons_pert_arr = cons_pert.array(mfi); const auto &xvel_pert_arr = xvel_pert.array(mfi); const auto &yvel_pert_arr = yvel_pert.array(mfi); diff --git a/Source/Initialization/Make.package b/Source/Initialization/Make.package index 685889786..21444ddce 100644 --- a/Source/Initialization/Make.package +++ b/Source/Initialization/Make.package @@ -5,6 +5,7 @@ CEXE_sources += ERF_init_from_input_sounding.cpp CEXE_sources += ERF_input_sponge.cpp CEXE_sources += ERF_init_bcs.cpp CEXE_sources += ERF_init1d.cpp +CEXE_sources += ERF_init_TurbPert.cpp ifeq ($(USE_WINDFARM),TRUE) CEXE_sources += ERF_init_windfarm.cpp diff --git a/Source/SourceTerms/ERF_make_sources.cpp b/Source/SourceTerms/ERF_make_sources.cpp index bb1e38fa6..69fbc5002 100644 --- a/Source/SourceTerms/ERF_make_sources.cpp +++ b/Source/SourceTerms/ERF_make_sources.cpp @@ -44,12 +44,13 @@ void make_sources (int level, const Real* dptr_rhotheta_src, const Real* dptr_rhoqt_src, const Real* dptr_wbar_sub, - const Vector d_rayleigh_ptrs_at_lev) + const Vector d_rayleigh_ptrs_at_lev, + TurbulentPerturbation& turbPert) { BL_PROFILE_REGION("erf_make_sources()"); // ***************************************************************************** - // Initialize source to zero since we re-compute it ever RK stage + // Initialize source to zero since we re-compute it every RK stage // ***************************************************************************** source.setVal(0.0); @@ -303,6 +304,16 @@ void make_sources (int level, if(!(solverChoice.spongeChoice.sponge_type == "input_sponge")){ ApplySpongeZoneBCsForCC(solverChoice.spongeChoice, geom, bx, cell_src, cell_data); } + + // ************************************************************************************* + // Add perturbation + // ************************************************************************************* + if (solverChoice.pert_type == PerturbationType::BPM) { + auto m_ixtype = S_data[IntVars::cons].boxArray().ixType(); + + // Apply stored values onto cell_src + turbPert.apply_tpi(level, bx, RhoTheta_comp, m_ixtype, cell_src); + } } // mfi } // OMP } diff --git a/Source/SourceTerms/Src_headers.H b/Source/SourceTerms/Src_headers.H index da30c721c..5cbb0c7a2 100644 --- a/Source/SourceTerms/Src_headers.H +++ b/Source/SourceTerms/Src_headers.H @@ -3,6 +3,7 @@ #include #include "DataStruct.H" +#include "TurbPertStruct.H" #ifdef ERF_USE_EB #include @@ -36,7 +37,8 @@ void make_sources (int level, int nrk, const amrex::Real* dptr_rhotheta_src, const amrex::Real* dptr_rhoqt_src, const amrex::Real* dptr_wbar_sub, - const amrex::Vector d_rayleigh_ptrs_at_lev); + const amrex::Vector d_rayleigh_ptrs_at_lev, + TurbulentPerturbation& turbPert); void make_mom_sources (int level, int nrk, amrex::Real dt, diff --git a/Source/TimeIntegration/ERF_Advance.cpp b/Source/TimeIntegration/ERF_Advance.cpp index 60fd99139..02448f2bf 100644 --- a/Source/TimeIntegration/ERF_Advance.cpp +++ b/Source/TimeIntegration/ERF_Advance.cpp @@ -35,7 +35,15 @@ ERF::Advance (int lev, Real time, Real dt_lev, int iteration, int /*ncycle*/) MultiFab& V_new = vars_new[lev][Vars::yvel]; MultiFab& W_new = vars_new[lev][Vars::zvel]; - + // TODO: Can test on multiple levels later + // Only apply to level 0 + // DUSTIN MA + if (solverChoice.pert_type == PerturbationType::BPM) { + if (lev == 0) { + turbPert.calc_tpi_update(lev, dt_lev, U_old, V_old, S_old); + turbPert.debug(); + } + } // configure ABLMost params if used MostWall boundary condition if (phys_bc_type[Orientation(Direction::z,Orientation::low)] == ERF_BC::MOST) { if (m_most) { diff --git a/Source/TimeIntegration/ERF_advance_dycore.cpp b/Source/TimeIntegration/ERF_advance_dycore.cpp index d4691bf38..4e59b5c9d 100644 --- a/Source/TimeIntegration/ERF_advance_dycore.cpp +++ b/Source/TimeIntegration/ERF_advance_dycore.cpp @@ -66,6 +66,9 @@ void ERF::advance_dycore(int level, Real* dptr_rhoqt_src = solverChoice.custom_moisture_forcing ? d_rhoqt_src[level].data() : nullptr; Real* dptr_wbar_sub = solverChoice.custom_w_subsidence ? d_w_subsid[level].data() : nullptr; + // Turbulent Perturbation Pointer + //Real* dptr_rhotheta_src = solverChoice.pert_type ? d_rhotheta_src[level].data() : nullptr; + Vector d_rayleigh_ptrs_at_lev; d_rayleigh_ptrs_at_lev.resize(Rayleigh::nvars); bool rayleigh_damp_any = (solverChoice.rayleigh_damp_U ||solverChoice.rayleigh_damp_V || diff --git a/Source/TimeIntegration/TI_slow_rhs_fun.H b/Source/TimeIntegration/TI_slow_rhs_fun.H index 54a709392..0c50a2d22 100644 --- a/Source/TimeIntegration/TI_slow_rhs_fun.H +++ b/Source/TimeIntegration/TI_slow_rhs_fun.H @@ -44,7 +44,8 @@ fine_geom, solverChoice, mapfac_u[level], mapfac_v[level], dptr_rhotheta_src, dptr_rhoqt_src, - dptr_wbar_sub, d_rayleigh_ptrs_at_lev); + dptr_wbar_sub, d_rayleigh_ptrs_at_lev, + turbPert); // Moving terrain if ( solverChoice.use_terrain && (solverChoice.terrain_type == TerrainType::Moving) ) @@ -413,7 +414,8 @@ fine_geom, solverChoice, mapfac_u[level], mapfac_v[level], dptr_rhotheta_src, dptr_rhoqt_src, - dptr_wbar_sub, d_rayleigh_ptrs_at_lev); + dptr_wbar_sub, d_rayleigh_ptrs_at_lev, + turbPert); int n_qstate = micro->Get_Qstate_Size(); make_mom_sources(level, nrk, slow_dt, S_data, S_prim, diff --git a/Source/prob_common.H b/Source/prob_common.H index 0e35a1869..acba961b7 100644 --- a/Source/prob_common.H +++ b/Source/prob_common.H @@ -365,6 +365,16 @@ public: }); } + /** + * Function to set constant values for turbulent pertubation calculation + */ + void + init_turbPert_const(TurbulentPerturbation& turbPert) + { + turbPert.tpi_Tinfty = base_parms.T_0; + turbPert.tpi_nonDim = 0.042; + } + protected: // Struct to store problem parameters ProbParmDefaults base_parms; From e571ffe625c319da2506ba2a7e88cdb5a5ee1a8c Mon Sep 17 00:00:00 2001 From: "Aaron M. Lattanzi" <103702284+AMLattanzi@users.noreply.github.com> Date: Tue, 25 Jun 2024 13:09:50 -0700 Subject: [PATCH 3/5] Terrain From File (#1658) * first pass at terrain initialization from file. * Uncomment time for the new function call. --- Exec/DevTests/MovingTerrain/prob.cpp | 94 +++++----- Exec/RegTests/ParticlesOverWoA/prob.cpp | 106 ++++++----- Exec/RegTests/Terrain2d_Cylinder/prob.cpp | 106 ++++++----- Exec/RegTests/Terrain3d_Hemisphere/prob.cpp | 96 +++++----- Exec/RegTests/WitchOfAgnesi/prob.cpp | 194 ++++++++++---------- Exec/SpongeTest/prob.cpp | 106 ++++++----- Source/prob_common.H | 80 ++++++++ 7 files changed, 457 insertions(+), 325 deletions(-) diff --git a/Exec/DevTests/MovingTerrain/prob.cpp b/Exec/DevTests/MovingTerrain/prob.cpp index c251cda5b..82b7a346b 100644 --- a/Exec/DevTests/MovingTerrain/prob.cpp +++ b/Exec/DevTests/MovingTerrain/prob.cpp @@ -155,47 +155,57 @@ Problem::init_custom_terrain (const Geometry& geom, MultiFab& z_phys_nd, const Real& time) { - // Domain cell size and real bounds - auto dx = geom.CellSizeArray(); - - // Domain valid box (z_nd is nodal) - const amrex::Box& domain = geom.Domain(); - int domlo_x = domain.smallEnd(0); int domhi_x = domain.bigEnd(0) + 1; - int domlo_z = domain.smallEnd(2); - - // Number of ghost cells - int ngrow = z_phys_nd.nGrow(); - - // Populate bottom plane - int k0 = domlo_z; - - Real Ampl = parms.Ampl; - Real wavelength = 100.; - Real kp = 2.0 * PI / wavelength; - Real g = CONST_GRAV; - Real omega = std::sqrt(g * kp); - - for (MFIter mfi(z_phys_nd,amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi) - { - // Grown box with no z range - amrex::Box xybx = mfi.growntilebox(ngrow); - xybx.setRange(2,0); - - amrex::Array4 const& z_arr = z_phys_nd.array(mfi); - - ParallelFor(xybx, [=] AMREX_GPU_DEVICE (int i, int j, int) { - - // Clip indices for ghost-cells - int ii = amrex::min(amrex::max(i,domlo_x),domhi_x); - - // Location of nodes - Real x = ii * dx[0]; - - // Wave height - Real height = Ampl * std::sin(kp * x - omega * time); - - // Populate terrain height - z_arr(i,j,k0) = height; - }); + // Check if a valid csv file exists for the terrain + std::string fname; + amrex::ParmParse pp("erf"); + auto valid_fname = pp.query("terrain_file_name",fname); + if (valid_fname) { + this->read_custom_terrain(fname,geom,z_phys_nd,time); + } else { + + // Domain cell size and real bounds + auto dx = geom.CellSizeArray(); + + // Domain valid box (z_nd is nodal) + const amrex::Box& domain = geom.Domain(); + int domlo_x = domain.smallEnd(0); int domhi_x = domain.bigEnd(0) + 1; + int domlo_z = domain.smallEnd(2); + + // Number of ghost cells + int ngrow = z_phys_nd.nGrow(); + + // Populate bottom plane + int k0 = domlo_z; + + Real Ampl = parms.Ampl; + Real wavelength = 100.; + Real kp = 2.0 * PI / wavelength; + Real g = CONST_GRAV; + Real omega = std::sqrt(g * kp); + + for (MFIter mfi(z_phys_nd,amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi) + { + // Grown box with no z range + amrex::Box xybx = mfi.growntilebox(ngrow); + xybx.setRange(2,0); + + amrex::Array4 const& z_arr = z_phys_nd.array(mfi); + + ParallelFor(xybx, [=] AMREX_GPU_DEVICE (int i, int j, int) + { + + // Clip indices for ghost-cells + int ii = amrex::min(amrex::max(i,domlo_x),domhi_x); + + // Location of nodes + Real x = ii * dx[0]; + + // Wave height + Real height = Ampl * std::sin(kp * x - omega * time); + + // Populate terrain height + z_arr(i,j,k0) = height; + }); + } } } diff --git a/Exec/RegTests/ParticlesOverWoA/prob.cpp b/Exec/RegTests/ParticlesOverWoA/prob.cpp index 422823dde..b65838449 100644 --- a/Exec/RegTests/ParticlesOverWoA/prob.cpp +++ b/Exec/RegTests/ParticlesOverWoA/prob.cpp @@ -126,54 +126,64 @@ void Problem::init_custom_terrain( const Geometry& geom, MultiFab& z_phys_nd, - const Real& /*time*/) + const Real& time) { - // Domain cell size and real bounds - auto dx = geom.CellSizeArray(); - auto ProbLoArr = geom.ProbLoArray(); - auto ProbHiArr = geom.ProbHiArray(); - - // Domain valid box (z_nd is nodal) - const amrex::Box& domain = geom.Domain(); - int domlo_x = domain.smallEnd(0); int domhi_x = domain.bigEnd(0) + 1; - // int domlo_y = domain.smallEnd(1); int domhi_y = domain.bigEnd(1) + 1; - int domlo_z = domain.smallEnd(2); - - // User function parameters - Real a = 0.5; - Real num = 8 * a * a * a; - Real xcen = 0.5 * (ProbLoArr[0] + ProbHiArr[0]); - // Real ycen = 0.5 * (ProbLoArr[1] + ProbHiArr[1]); - - // Number of ghost cells - int ngrow = z_phys_nd.nGrow(); - - // Populate bottom plane - int k0 = domlo_z; - - for ( amrex::MFIter mfi(z_phys_nd,amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi ) - { - // Grown box with no z range - amrex::Box xybx = mfi.growntilebox(ngrow); - xybx.setRange(2,0); - - amrex::Array4 const& z_arr = z_phys_nd.array(mfi); - - ParallelFor(xybx, [=] AMREX_GPU_DEVICE (int i, int j, int) { - - // Clip indices for ghost-cells - int ii = amrex::min(amrex::max(i,domlo_x),domhi_x); - // int jj = amrex::min(amrex::max(j,domlo_y),domhi_y); - - // Location of nodes - Real x = (ii * dx[0] - xcen); - // Real y = (jj * dx[1] - ycen); - - // WoA Hill in x-direction - Real height = num / (x*x + 4 * a * a); - - // Populate terrain height - z_arr(i,j,k0) = height; - }); + // Check if a valid csv file exists for the terrain + std::string fname; + amrex::ParmParse pp("erf"); + auto valid_fname = pp.query("terrain_file_name",fname); + if (valid_fname) { + this->read_custom_terrain(fname,geom,z_phys_nd,time); + } else { + + // Domain cell size and real bounds + auto dx = geom.CellSizeArray(); + auto ProbLoArr = geom.ProbLoArray(); + auto ProbHiArr = geom.ProbHiArray(); + + // Domain valid box (z_nd is nodal) + const amrex::Box& domain = geom.Domain(); + int domlo_x = domain.smallEnd(0); int domhi_x = domain.bigEnd(0) + 1; + // int domlo_y = domain.smallEnd(1); int domhi_y = domain.bigEnd(1) + 1; + int domlo_z = domain.smallEnd(2); + + // User function parameters + Real a = 0.5; + Real num = 8 * a * a * a; + Real xcen = 0.5 * (ProbLoArr[0] + ProbHiArr[0]); + // Real ycen = 0.5 * (ProbLoArr[1] + ProbHiArr[1]); + + // Number of ghost cells + int ngrow = z_phys_nd.nGrow(); + + // Populate bottom plane + int k0 = domlo_z; + + for ( amrex::MFIter mfi(z_phys_nd,amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi ) + { + // Grown box with no z range + amrex::Box xybx = mfi.growntilebox(ngrow); + xybx.setRange(2,0); + + amrex::Array4 const& z_arr = z_phys_nd.array(mfi); + + ParallelFor(xybx, [=] AMREX_GPU_DEVICE (int i, int j, int) + { + + // Clip indices for ghost-cells + int ii = amrex::min(amrex::max(i,domlo_x),domhi_x); + // int jj = amrex::min(amrex::max(j,domlo_y),domhi_y); + + // Location of nodes + Real x = (ii * dx[0] - xcen); + // Real y = (jj * dx[1] - ycen); + + // WoA Hill in x-direction + Real height = num / (x*x + 4 * a * a); + + // Populate terrain height + z_arr(i,j,k0) = height; + }); + } } } diff --git a/Exec/RegTests/Terrain2d_Cylinder/prob.cpp b/Exec/RegTests/Terrain2d_Cylinder/prob.cpp index 948d8c9b2..f99f0e224 100644 --- a/Exec/RegTests/Terrain2d_Cylinder/prob.cpp +++ b/Exec/RegTests/Terrain2d_Cylinder/prob.cpp @@ -87,54 +87,64 @@ void Problem::init_custom_terrain( const Geometry& geom, MultiFab& z_phys_nd, - const Real& /*time*/) + const Real& time) { - // Domain cell size and real bounds - auto dx = geom.CellSizeArray(); - auto ProbLoArr = geom.ProbLoArray(); - auto ProbHiArr = geom.ProbHiArray(); - - // Domain valid box (z_nd is nodal) - const amrex::Box& domain = geom.Domain(); - int domlo_x = domain.smallEnd(0); int domhi_x = domain.bigEnd(0) + 1; - // int domlo_y = domain.smallEnd(1); int domhi_y = domain.bigEnd(1) + 1; - int domlo_z = domain.smallEnd(2); - - // User function parameters - Real a = 0.5; - Real num = 8 * a * a * a; - Real xcen = 0.5 * (ProbLoArr[0] + ProbHiArr[0]); - Real ycen = 0.5 * (ProbLoArr[1] + ProbHiArr[1]); - - // Number of ghost cells - int ngrow = z_phys_nd.nGrow(); - - // Populate bottom plane - int k0 = domlo_z; - - for ( amrex::MFIter mfi(z_phys_nd,amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi ) - { - // Grown box with no z range - amrex::Box xybx = mfi.growntilebox(ngrow); - xybx.setRange(2,0); - - amrex::Array4 const& z_arr = z_phys_nd.array(mfi); - - ParallelFor(xybx, [=] AMREX_GPU_DEVICE (int i, int j, int) { - - // Clip indices for ghost-cells - int ii = amrex::min(amrex::max(i,domlo_x),domhi_x); - - // Location of nodes - Real x = (ii * dx[0] - xcen); - - if(fabs(x)read_custom_terrain(fname,geom,z_phys_nd,time); + } else { + + // Domain cell size and real bounds + auto dx = geom.CellSizeArray(); + auto ProbLoArr = geom.ProbLoArray(); + auto ProbHiArr = geom.ProbHiArray(); + + // Domain valid box (z_nd is nodal) + const amrex::Box& domain = geom.Domain(); + int domlo_x = domain.smallEnd(0); int domhi_x = domain.bigEnd(0) + 1; + // int domlo_y = domain.smallEnd(1); int domhi_y = domain.bigEnd(1) + 1; + int domlo_z = domain.smallEnd(2); + + // User function parameters + Real a = 0.5; + Real num = 8 * a * a * a; + Real xcen = 0.5 * (ProbLoArr[0] + ProbHiArr[0]); + Real ycen = 0.5 * (ProbLoArr[1] + ProbHiArr[1]); + + // Number of ghost cells + int ngrow = z_phys_nd.nGrow(); + + // Populate bottom plane + int k0 = domlo_z; + + for ( amrex::MFIter mfi(z_phys_nd,amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi ) + { + // Grown box with no z range + amrex::Box xybx = mfi.growntilebox(ngrow); + xybx.setRange(2,0); + + amrex::Array4 const& z_arr = z_phys_nd.array(mfi); + + ParallelFor(xybx, [=] AMREX_GPU_DEVICE (int i, int j, int) + { + + // Clip indices for ghost-cells + int ii = amrex::min(amrex::max(i,domlo_x),domhi_x); + + // Location of nodes + Real x = (ii * dx[0] - xcen); + + if(fabs(x)read_custom_terrain(fname,geom,z_phys_nd,time); + } else { + + // Domain cell size and real bounds + auto dx = geom.CellSizeArray(); + auto ProbLoArr = geom.ProbLoArray(); + auto ProbHiArr = geom.ProbHiArray(); + + // Domain valid box (z_nd is nodal) + const amrex::Box& domain = geom.Domain(); + int domlo_x = domain.smallEnd(0); int domhi_x = domain.bigEnd(0) + 1; + int domlo_y = domain.smallEnd(1); int domhi_y = domain.bigEnd(1) + 1; + int domlo_z = domain.smallEnd(2); + + // User function parameters + Real a = 0.5; + Real xcen = 0.5 * (ProbLoArr[0] + ProbHiArr[0]); + Real ycen = 0.5 * (ProbLoArr[1] + ProbHiArr[1]); + + // Number of ghost cells + int ngrow = z_phys_nd.nGrow(); + + // Populate bottom plane + int k0 = domlo_z; + + for ( amrex::MFIter mfi(z_phys_nd,amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi ) + { + // Grown box with no z range + amrex::Box xybx = mfi.growntilebox(ngrow); + xybx.setRange(2,0); - amrex::Array4 const& z_arr = z_phys_nd.array(mfi); + amrex::Array4 const& z_arr = z_phys_nd.array(mfi); - ParallelFor(xybx, [=] AMREX_GPU_DEVICE (int i, int j, int) { + ParallelFor(xybx, [=] AMREX_GPU_DEVICE (int i, int j, int) + { - // Clip indices for ghost-cells - int ii = amrex::min(amrex::max(i,domlo_x),domhi_x); - int jj = amrex::min(amrex::max(j,domlo_y),domhi_y); + // Clip indices for ghost-cells + int ii = amrex::min(amrex::max(i,domlo_x),domhi_x); + int jj = amrex::min(amrex::max(j,domlo_y),domhi_y); - // Location of nodes - Real x = (ii * dx[0] - xcen); - Real y = (jj * dx[1] - ycen); + // Location of nodes + Real x = (ii * dx[0] - xcen); + Real y = (jj * dx[1] - ycen); - if(std::pow(x*x + y*y, 0.5) < a){ - z_arr(i,j,k0) = std::pow(a*a - x*x - y*y , 0.5); - } - else{ - z_arr(i,j,k0) = 0.0; - } + if(std::pow(x*x + y*y, 0.5) < a){ + z_arr(i,j,k0) = std::pow(a*a - x*x - y*y , 0.5); + } + else{ + z_arr(i,j,k0) = 0.0; + } - }); + }); + } } } diff --git a/Exec/RegTests/WitchOfAgnesi/prob.cpp b/Exec/RegTests/WitchOfAgnesi/prob.cpp index 45c64c93d..82e89f3a1 100644 --- a/Exec/RegTests/WitchOfAgnesi/prob.cpp +++ b/Exec/RegTests/WitchOfAgnesi/prob.cpp @@ -4,14 +4,14 @@ using namespace amrex; std::unique_ptr -amrex_probinit( +amrex_probinit ( const amrex_real* /*problo*/, const amrex_real* /*probhi*/) { return std::make_unique(); } -Problem::Problem() +Problem::Problem () { // Parse params amrex::ParmParse pp("prob"); @@ -25,7 +25,7 @@ Problem::Problem() } void -Problem::init_custom_pert( +Problem::init_custom_pert ( const Box& bx, const Box& xbx, const Box& ybx, @@ -35,119 +35,121 @@ Problem::init_custom_pert( Array4 const& x_vel_pert, Array4 const& y_vel_pert, Array4 const& z_vel_pert, - Array4 const& r_hse, - Array4 const& p_hse, + Array4 const& /*r_hse*/, + Array4 const& /*p_hse*/, Array4 const& z_nd, - Array4 const& z_cc, + Array4 const& /*z_cc*/, GeometryData const& geomdata, Array4 const& /*mf_m*/, Array4 const& /*mf_u*/, Array4 const& /*mf_v*/, const SolverChoice& sc) { - const int khi = geomdata.Domain().bigEnd()[2]; + const int khi = geomdata.Domain().bigEnd()[2]; const bool use_moisture = (sc.moisture_type != MoistureType::None); - AMREX_ALWAYS_ASSERT(bx.length()[2] == khi+1); + AMREX_ALWAYS_ASSERT(bx.length()[2] == khi+1); - amrex::ParallelFor(bx, [=, parms=parms] AMREX_GPU_DEVICE(int i, int j, int k) noexcept - { - // Geometry (note we must include these here to get the data on device) - const auto prob_lo = geomdata.ProbLo(); - const auto dx = geomdata.CellSize(); - const Real x = prob_lo[0] + (i + 0.5) * dx[0]; - const Real z = z_cc(i,j,k); + amrex::ParallelFor(bx, [=, parms=parms] AMREX_GPU_DEVICE(int i, int j, int k) noexcept + { + // Set scalar = 0 everywhere + state_pert(i, j, k, RhoScalar_comp) = 0.0; - // Set scalar = 0 everywhere - state_pert(i, j, k, RhoScalar_comp) = 0.0; + if (use_moisture) { + state_pert(i, j, k, RhoQ1_comp) = 0.0; + state_pert(i, j, k, RhoQ2_comp) = 0.0; + } + }); - if (use_moisture) { - state_pert(i, j, k, RhoQ1_comp) = 0.0; - state_pert(i, j, k, RhoQ2_comp) = 0.0; - } - }); - - // Set the x-velocity - amrex::ParallelFor(xbx, [=, parms=parms] AMREX_GPU_DEVICE(int i, int j, int k) noexcept - { - x_vel_pert(i, j, k) = parms.U_0; - }); - - // Set the y-velocity - amrex::ParallelFor(ybx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept - { - y_vel_pert(i, j, k) = 0.0; - }); - - const auto dx = geomdata.CellSize(); - amrex::GpuArray dxInv; - dxInv[0] = 1. / dx[0]; - dxInv[1] = 1. / dx[1]; - dxInv[2] = 1. / dx[2]; - - // Set the z-velocity from impenetrable condition - amrex::ParallelFor(zbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept - { - z_vel_pert(i, j, k) = WFromOmega(i, j, k, 0.0, x_vel_pert, y_vel_pert, z_nd, dxInv); - }); - - amrex::Gpu::streamSynchronize(); + // Set the x-velocity + amrex::ParallelFor(xbx, [=, parms=parms] AMREX_GPU_DEVICE(int i, int j, int k) noexcept + { + x_vel_pert(i, j, k) = parms.U_0; + }); + + // Set the y-velocity + amrex::ParallelFor(ybx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept + { + y_vel_pert(i, j, k) = 0.0; + }); + + const auto dx = geomdata.CellSize(); + amrex::GpuArray dxInv; + dxInv[0] = 1. / dx[0]; + dxInv[1] = 1. / dx[1]; + dxInv[2] = 1. / dx[2]; + + // Set the z-velocity from impenetrable condition + amrex::ParallelFor(zbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept + { + z_vel_pert(i, j, k) = WFromOmega(i, j, k, 0.0, x_vel_pert, y_vel_pert, z_nd, dxInv); + }); + amrex::Gpu::streamSynchronize(); } void -Problem::init_custom_terrain( +Problem::init_custom_terrain ( const Geometry& geom, MultiFab& z_phys_nd, - const Real& /*time*/) + const Real& time) { - // Domain cell size and real bounds - auto dx = geom.CellSizeArray(); - auto ProbLoArr = geom.ProbLoArray(); - auto ProbHiArr = geom.ProbHiArray(); - - // Domain valid box (z_nd is nodal) - const amrex::Box& domain = geom.Domain(); - int domlo_x = domain.smallEnd(0); int domhi_x = domain.bigEnd(0) + 1; - // int domlo_y = domain.smallEnd(1); int domhi_y = domain.bigEnd(1) + 1; - int domlo_z = domain.smallEnd(2); - - // User function parameters - Real a = 0.5; - Real num = 8 * a * a * a; - Real xcen = 0.5 * (ProbLoArr[0] + ProbHiArr[0]); - // Real ycen = 0.5 * (ProbLoArr[1] + ProbHiArr[1]); - - // Number of ghost cells - int ngrow = z_phys_nd.nGrow(); - - // Populate bottom plane - int k0 = domlo_z; - - for ( amrex::MFIter mfi(z_phys_nd,amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi ) - { - // Grown box with no z range - amrex::Box xybx = mfi.growntilebox(ngrow); - xybx.setRange(2,0); - - amrex::Array4 const& z_arr = z_phys_nd.array(mfi); - - ParallelFor(xybx, [=] AMREX_GPU_DEVICE (int i, int j, int) { - - // Clip indices for ghost-cells - int ii = amrex::min(amrex::max(i,domlo_x),domhi_x); - // int jj = amrex::min(amrex::max(j,domlo_y),domhi_y); - - // Location of nodes - Real x = (ii * dx[0] - xcen); - // Real y = (jj * dx[1] - ycen); - - // WoA Hill in x-direction - Real height = num / (x*x + 4 * a * a); - // Populate terrain height - z_arr(i,j,k0) = height; - }); + // Check if a valid csv file exists for the terrain + std::string fname; + amrex::ParmParse pp("erf"); + auto valid_fname = pp.query("terrain_file_name",fname); + if (valid_fname) { + this->read_custom_terrain(fname,geom,z_phys_nd,time); + } else { + // Domain cell size and real bounds + auto dx = geom.CellSizeArray(); + auto ProbLoArr = geom.ProbLoArray(); + auto ProbHiArr = geom.ProbHiArray(); + + // Domain valid box (z_nd is nodal) + const amrex::Box& domain = geom.Domain(); + int domlo_x = domain.smallEnd(0); int domhi_x = domain.bigEnd(0) + 1; + // int domlo_y = domain.smallEnd(1); int domhi_y = domain.bigEnd(1) + 1; + int domlo_z = domain.smallEnd(2); + + // User function parameters + Real a = 0.5; + Real num = 8 * a * a * a; + Real xcen = 0.5 * (ProbLoArr[0] + ProbHiArr[0]); + // Real ycen = 0.5 * (ProbLoArr[1] + ProbHiArr[1]); + + // Number of ghost cells + int ngrow = z_phys_nd.nGrow(); + + // Populate bottom plane + int k0 = domlo_z; + + for ( amrex::MFIter mfi(z_phys_nd,amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi ) + { + // Grown box with no z range + amrex::Box xybx = mfi.growntilebox(ngrow); + xybx.setRange(2,0); + + amrex::Array4 const& z_arr = z_phys_nd.array(mfi); + + ParallelFor(xybx, [=] AMREX_GPU_DEVICE (int i, int j, int) + { + // Clip indices for ghost-cells + int ii = amrex::min(amrex::max(i,domlo_x),domhi_x); + // int jj = amrex::min(amrex::max(j,domlo_y),domhi_y); + + // Location of nodes + Real x = (ii * dx[0] - xcen); + // Real y = (jj * dx[1] - ycen); + + // WoA Hill in x-direction + Real height = num / (x*x + 4 * a * a); + + // Populate terrain height + z_arr(i,j,k0) = height; + }); + } } } diff --git a/Exec/SpongeTest/prob.cpp b/Exec/SpongeTest/prob.cpp index 469ba545a..2cb3b36e2 100644 --- a/Exec/SpongeTest/prob.cpp +++ b/Exec/SpongeTest/prob.cpp @@ -87,54 +87,64 @@ void Problem::init_custom_terrain( const Geometry& geom, MultiFab& z_phys_nd, - const Real& /*time*/) + const Real& time) { - // Domain cell size and real bounds - auto dx = geom.CellSizeArray(); - auto ProbLoArr = geom.ProbLoArray(); - auto ProbHiArr = geom.ProbHiArray(); - - // Domain valid box (z_nd is nodal) - const amrex::Box& domain = geom.Domain(); - int domlo_x = domain.smallEnd(0); int domhi_x = domain.bigEnd(0) + 1; - // int domlo_y = domain.smallEnd(1); int domhi_y = domain.bigEnd(1) + 1; - int domlo_z = domain.smallEnd(2); - - // User function parameters - Real a = 0.5; - Real num = 8 * a * a * a; - Real xcen = 0.5 * (ProbLoArr[0] + ProbHiArr[0]); - Real ycen = 0.5 * (ProbLoArr[1] + ProbHiArr[1]); - - // Number of ghost cells - int ngrow = z_phys_nd.nGrow(); - - // Populate bottom plane - int k0 = domlo_z; - - for ( amrex::MFIter mfi(z_phys_nd,amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi ) - { - // Grown box with no z range - amrex::Box xybx = mfi.growntilebox(ngrow); - xybx.setRange(2,0); - - amrex::Array4 const& z_arr = z_phys_nd.array(mfi); - - ParallelFor(xybx, [=] AMREX_GPU_DEVICE (int i, int j, int) { - - // Clip indices for ghost-cells - int ii = amrex::min(amrex::max(i,domlo_x),domhi_x); - - // Location of nodes - Real x = (ii * dx[0] - xcen); - - if(fabs(x)read_custom_terrain(fname,geom,z_phys_nd,time); + } else { + + // Domain cell size and real bounds + auto dx = geom.CellSizeArray(); + auto ProbLoArr = geom.ProbLoArray(); + auto ProbHiArr = geom.ProbHiArray(); + + // Domain valid box (z_nd is nodal) + const amrex::Box& domain = geom.Domain(); + int domlo_x = domain.smallEnd(0); int domhi_x = domain.bigEnd(0) + 1; + // int domlo_y = domain.smallEnd(1); int domhi_y = domain.bigEnd(1) + 1; + int domlo_z = domain.smallEnd(2); + + // User function parameters + Real a = 0.5; + Real num = 8 * a * a * a; + Real xcen = 0.5 * (ProbLoArr[0] + ProbHiArr[0]); + Real ycen = 0.5 * (ProbLoArr[1] + ProbHiArr[1]); + + // Number of ghost cells + int ngrow = z_phys_nd.nGrow(); + + // Populate bottom plane + int k0 = domlo_z; + + for ( amrex::MFIter mfi(z_phys_nd,amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi ) + { + // Grown box with no z range + amrex::Box xybx = mfi.growntilebox(ngrow); + xybx.setRange(2,0); + + amrex::Array4 const& z_arr = z_phys_nd.array(mfi); + + ParallelFor(xybx, [=] AMREX_GPU_DEVICE (int i, int j, int) + { + + // Clip indices for ghost-cells + int ii = amrex::min(amrex::max(i,domlo_x),domhi_x); + + // Location of nodes + Real x = (ii * dx[0] - xcen); + + if(fabs(x) m_xterrain,m_yterrain,m_zterrain; + amrex::Real value1,value2,value3; + while(file>>value1>>value2>>value3){ + m_xterrain.push_back(value1); + m_yterrain.push_back(value2); + m_zterrain.push_back(value3); + } + file.close(); + + // Copy data to the GPU + int ncell = m_xterrain.size(); + amrex::Gpu::DeviceVector d_xterrain(ncell),d_yterrain(ncell),d_zterrain(ncell); + amrex::Gpu::copy(amrex::Gpu::hostToDevice, m_xterrain.begin(), m_xterrain.end(), d_xterrain.begin()); + amrex::Gpu::copy(amrex::Gpu::hostToDevice, m_yterrain.begin(), m_yterrain.end(), d_yterrain.begin()); + amrex::Gpu::copy(amrex::Gpu::hostToDevice, m_zterrain.begin(), m_zterrain.end(), d_zterrain.begin()); + amrex::Real* d_xt = d_xterrain.data(); + amrex::Real* d_yt = d_yterrain.data(); + amrex::Real* d_zt = d_zterrain.data(); + + // Populate z_phys data + amrex::Real tol = 1.0e-4; + int ngrow = z_phys_nd.nGrow(); + auto dx = geom.CellSizeArray(); + auto ProbLoArr = geom.ProbLoArray(); + int ilo = geom.Domain().smallEnd(0); + int jlo = geom.Domain().smallEnd(1); + int klo = geom.Domain().smallEnd(2); + int ihi = geom.Domain().bigEnd(0); + int jhi = geom.Domain().bigEnd(1); + for (amrex::MFIter mfi(z_phys_nd,amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi ) + { + // Grown box with no z range + amrex::Box xybx = mfi.growntilebox(ngrow); + xybx.setRange(2,0); + + amrex::Array4 const& z_arr = z_phys_nd.array(mfi); + amrex::ParallelFor(xybx, [=] AMREX_GPU_DEVICE (int i, int j, int /*k*/) + { + // Clip indices for ghost-cells + int ii = amrex::min(amrex::max(i,ilo),ihi); + int jj = amrex::min(amrex::max(j,jlo),jhi); + + // Location of nodes + bool found = false; + amrex::Real x = ProbLoArr[0] + ii * dx[0]; + amrex::Real y = ProbLoArr[1] + jj * dx[1]; + amrex::Real zloc = 0.0; + for(int n=0; n Date: Tue, 25 Jun 2024 14:30:18 -0700 Subject: [PATCH 4/5] Limit surface roughness to on land value. (#1659) --- Source/BoundaryConditions/MOSTStress.H | 15 +++++++++------ 1 file changed, 9 insertions(+), 6 deletions(-) diff --git a/Source/BoundaryConditions/MOSTStress.H b/Source/BoundaryConditions/MOSTStress.H index 272213257..b15bd2410 100644 --- a/Source/BoundaryConditions/MOSTStress.H +++ b/Source/BoundaryConditions/MOSTStress.H @@ -256,8 +256,8 @@ struct adiabatic_wave_coupled u_star_arr(i,j,k) = mdata.kappa * umm_arr(i,j,k) / std::log(mdata.zref / z0_arr(i,j,k)); do { ustar = u_star_arr(i,j,k); - z0 = std::max(1200.0 * Hwave_arr(i,j,k) * std::pow( Hwave_arr(i,j,k)/(Lwave_arr(i,j,k)+eps), 4.5 ) - + 0.11 * eta_arr(ie,je,k,EddyDiff::Mom_v) / ustar, z0_eps); + z0 = std::min( std::max(1200.0 * Hwave_arr(i,j,k) * std::pow( Hwave_arr(i,j,k)/(Lwave_arr(i,j,k)+eps), 4.5 ) + + 0.11 * eta_arr(ie,je,k,EddyDiff::Mom_v) / ustar, z0_eps), z0_max ); u_star_arr(i,j,k) = mdata.kappa * umm_arr(i,j,k) / std::log(mdata.zref / z0); ++iter; } while ((std::abs(u_star_arr(i,j,k) - ustar) > tol) && iter <= max_iters); @@ -273,6 +273,7 @@ private: const amrex::Real tol = 1.0e-5; const amrex::Real eps = 1e-15; const amrex::Real z0_eps = 1.0e-6; + const amrex::Real z0_max = 0.1; }; @@ -515,8 +516,8 @@ struct surface_flux_wave_coupled u_star_arr(i,j,k) = mdata.kappa * umm_arr(i,j,k) / std::log(mdata.zref / z0_arr(i,j,k)); do { ustar = u_star_arr(i,j,k); - z0 = std::max(1200.0 * Hwave_arr(i,j,k) * std::pow( Hwave_arr(i,j,k)/(Lwave_arr(i,j,k)+eps), 4.5 ) - + 0.11 * eta_arr(ie,je,k,EddyDiff::Mom_v) / ustar, z0_eps); + z0 = std::min( std::max(1200.0 * Hwave_arr(i,j,k) * std::pow( Hwave_arr(i,j,k)/(Lwave_arr(i,j,k)+eps), 4.5 ) + + 0.11 * eta_arr(ie,je,k,EddyDiff::Mom_v) / ustar, z0_eps), z0_max ); Olen = -ustar * ustar * ustar * tm_arr(i,j,k) / (mdata.kappa * mdata.gravity * mdata.surf_temp_flux); zeta = mdata.zref / Olen; @@ -539,6 +540,7 @@ private: const amrex::Real tol = 1.0e-5; const amrex::Real eps = 1e-15; const amrex::Real z0_eps = 1.0e-6; + const amrex::Real z0_max = 0.1; }; @@ -788,8 +790,8 @@ struct surface_temp_wave_coupled u_star_arr(i,j,k) = mdata.kappa * umm_arr(i,j,k) / std::log(mdata.zref / z0_arr(i,j,k)); do { ustar = u_star_arr(i,j,k); - z0 = std::max(1200.0 * Hwave_arr(i,j,k) * std::pow( Hwave_arr(i,j,k)/(Lwave_arr(i,j,k)+eps), 4.5 ) - + 0.11 * eta_arr(ie,je,k,EddyDiff::Mom_v) / ustar, z0_eps); + z0 = std::min( std::max(1200.0 * Hwave_arr(i,j,k) * std::pow( Hwave_arr(i,j,k)/(Lwave_arr(i,j,k)+eps), 4.5 ) + + 0.11 * eta_arr(ie,je,k,EddyDiff::Mom_v) / ustar, z0_eps), z0_max ); tflux = -(tm_arr(i,j,k) - t_surf_arr(i,j,k)) * ustar * mdata.kappa / (std::log(mdata.zref / z0) - psi_h); Olen = -ustar * ustar * ustar * tm_arr(i,j,k) / @@ -813,6 +815,7 @@ private: const amrex::Real tol = 1.0e-5; const amrex::Real eps = 1e-15; const amrex::Real z0_eps = 1.0e-6; + const amrex::Real z0_max = 0.1; }; From aa98170cc61e8e53eda49bfe1d1c1e339c2222ea Mon Sep 17 00:00:00 2001 From: "Jean M. Sexton" Date: Wed, 26 Jun 2024 19:21:53 -0400 Subject: [PATCH 5/5] Add WW3 submodule to point to mpmd branch (#1660) --- .gitmodules | 3 +++ Submodules/WW3 | 1 + 2 files changed, 4 insertions(+) create mode 160000 Submodules/WW3 diff --git a/.gitmodules b/.gitmodules index 7a95f042f..92e06008c 100644 --- a/.gitmodules +++ b/.gitmodules @@ -7,3 +7,6 @@ [submodule "Submodules/RRTMGP"] path = Submodules/RRTMGP url = https://github.com/E3SM-Project/rte-rrtmgp +[submodule "Submodules/WW3"] + path = Submodules/WW3 + url = https://github.com/erf-model/WW3 diff --git a/Submodules/WW3 b/Submodules/WW3 new file mode 160000 index 000000000..f5fa58228 --- /dev/null +++ b/Submodules/WW3 @@ -0,0 +1 @@ +Subproject commit f5fa5822849e33976cdc935a91b304c60bc2a976