diff --git a/CMake/BuildERFExe.cmake b/CMake/BuildERFExe.cmake index 4b0dc4cee..209decc1a 100644 --- a/CMake/BuildERFExe.cmake +++ b/CMake/BuildERFExe.cmake @@ -142,6 +142,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_input_sponge.cpp ${SRC_DIR}/IO/Checkpoint.cpp ${SRC_DIR}/IO/ERF_ReadBndryPlanes.cpp ${SRC_DIR}/IO/ERF_WriteBndryPlanes.cpp @@ -152,6 +153,7 @@ function(build_erf_lib erf_lib_name) ${SRC_DIR}/IO/writeJobInfo.cpp ${SRC_DIR}/IO/console_io.cpp ${SRC_DIR}/SourceTerms/ERF_ApplySpongeZoneBCs.cpp + ${SRC_DIR}/SourceTerms/ERF_ApplySpongeZoneBCs_ReadFromFile.cpp ${SRC_DIR}/SourceTerms/ERF_make_buoyancy.cpp ${SRC_DIR}/SourceTerms/ERF_add_thin_body_sources.cpp ${SRC_DIR}/SourceTerms/ERF_make_mom_sources.cpp @@ -174,7 +176,8 @@ function(build_erf_lib erf_lib_name) ${SRC_DIR}/Utils/MomentumToVelocity.cpp ${SRC_DIR}/Utils/TerrainMetrics.cpp ${SRC_DIR}/Utils/VelocityToMomentum.cpp - ${SRC_DIR}/Utils/InteriorGhostCells.cpp + ${SRC_DIR}/Utils/InteriorGhostCells.cpp + ${SRC_DIR}/Utils/Time_Avg_Vel.cpp ${SRC_DIR}/Microphysics/SAM/Init_SAM.cpp ${SRC_DIR}/Microphysics/SAM/Cloud_SAM.cpp ${SRC_DIR}/Microphysics/SAM/IceFall.cpp diff --git a/Docs/sphinx_doc/BoundaryConditions.rst b/Docs/sphinx_doc/BoundaryConditions.rst index 5073c3a32..1bfdea0ae 100644 --- a/Docs/sphinx_doc/BoundaryConditions.rst +++ b/Docs/sphinx_doc/BoundaryConditions.rst @@ -198,7 +198,8 @@ ERF provides the capability to apply sponge zones at the boundaries to prevent s \frac{dQ}{dt} = \mathrm{RHS} - A\xi^n(Q-Q_\mathrm{target}) -where RHS are the other right-hand side terms. The parameters to be set by the user are -- `A` is the sponge amplitude, `n` is the sponge strength and the :math:`Q_\mathrm{target}` -- the target solution in the sponge. `\xi` is a linear coordinate that is 0 at the beginning of the sponge and 1 at the end. An example of the sponge inputs can be found in ``Exec/RegTests/Terrain2d_Cylinder``. +where RHS are the other right-hand side terms. The parameters to be set by the user are -- `A` is the sponge amplitude, `n` is the sponge strength and the :math:`Q_\mathrm{target}` -- the target solution in the sponge. :math:`\xi` is a linear coordinate that is 0 at the beginning of the sponge and 1 at the end. An example of the sponge inputs can be found in ``Exec/RegTests/Terrain2d_Cylinder`` and is given below. This list of inputs specifies sponge zones in the inlet and outlet of the domain in the x-direction and the outlet of the domain in the z-direction. The `start` and `end` parameters specify the starting and ending of the sponge zones. At the inlet, the sponge starts at :math:`x=0` and at the outlet the sponge ends at :math:`x=L` -- the end of the domain. The sponge amplitude `A` has to be adjust +ed in a problem-specific manner. The density and the :math:`x, y, z` velocities to be used in the sponge zones have to be specified in the inputs list. :: @@ -214,3 +215,14 @@ where RHS are the other right-hand side terms. The parameters to be set by the u erf.sponge_x_velocity = 10.0 erf.sponge_y_velocity = 0.0 erf.sponge_z_velocity = 0.0 + +Another way of specifying sponge zones is by providing the sponge zone data as a text file input. This is currently implemented only for forcing :math:`x` and :math:`y` velocities in the sponge zones. +The sponge data is input as a text file with 3 columns containing :math:`z, u, v` values. An example can be found in ``Exec/SpongeTest`` and a sample inputs list for using this feature is given below. This list specifies a sponge zone in the inlet in the x-direction. The :math:`u` and :math:`v` velocity forcing in the sponge zones will be read in from the text file -- `input_sponge_file.txt`. + +:: + + erf.sponge_type = "input_sponge" + erf.input_sponge_file = "input_sponge_file.txt" + erf.sponge_strength = 1000.0 + erf.use_xlo_sponge_damping = true + erf.xlo_sponge_end = 4.0 diff --git a/Docs/sphinx_doc/Plotfiles.rst b/Docs/sphinx_doc/Plotfiles.rst index 4ff201b9c..a13332a1d 100644 --- a/Docs/sphinx_doc/Plotfiles.rst +++ b/Docs/sphinx_doc/Plotfiles.rst @@ -108,7 +108,9 @@ PlotFile Outputs ================ Plotfiles can include the quantities of several simulation parameters as output. -They are summarized in the list below. +They are summarized in the list below. Note that temporally averaged quantities +(e.g., ``u_t_avg, v_t_avg, w_t_avg, umag_t_avg``) require the user to enable the +storage of the time averaged variables with ``erf.time_avg_vel = true``. Output Options -------------- @@ -201,18 +203,34 @@ Output Options | | | | | | +-----------------------------+------------------+ -| **vorticity_x* | x-component of | +| **vorticity_x** | x-component of | | | vorticity | | | | +-----------------------------+------------------+ -| **vorticity_y* | y-component of | +| **vorticity_y** | y-component of | | | vorticity | | | | +-----------------------------+------------------+ -| **vorticity_z* | z-component of | +| **vorticity_z** | z-component of | | | vorticity | | | | +-----------------------------+------------------+ +| **u_t_avg** | time average of | +| | x-component of | +| | velocity | ++-----------------------------+------------------+ +| **v_t_avg** | time average of | +| | y-component of | +| | velocity | ++-----------------------------+------------------+ +| **w_t_avg** | time average of | +| | z-component of | +| | velocity | ++-----------------------------+------------------+ +| **umag_t_avg** | time average of | +| | velocity mag | +| | | ++-----------------------------+------------------+ | **rhoadv_0** | Conserved scalar | | | | | | | @@ -233,6 +251,14 @@ Output Options | | | | | | +-----------------------------+------------------+ +| **lat_m** | Latitude at mass | +| | points | +| | | ++-----------------------------+------------------+ +| **lon_m** | Longitude at | +| | mass points | +| | | ++-----------------------------+------------------+ | **Kmv** | Vertical | | | Eddy Diffusivity | | | of Momentum | diff --git a/Exec/ABL/inputs_vel_avg b/Exec/ABL/inputs_vel_avg new file mode 100644 index 000000000..9b114b86d --- /dev/null +++ b/Exec/ABL/inputs_vel_avg @@ -0,0 +1,65 @@ +# ------------------ INPUTS TO MAIN PROGRAM ------------------- +max_step = 4000 + +amrex.fpe_trap_invalid = 1 + +fabarray.mfiter_tile_size = 1024 1024 1024 + +# PROBLEM SIZE & GEOMETRY +geometry.prob_extent = 1024 1024 1024 +amr.n_cell = 64 64 64 + +geometry.is_periodic = 1 1 0 + +zlo.type = "NoSlipWall" +zhi.type = "SlipWall" + +# TIME STEP CONTROL +erf.fixed_dt = 0.1 # fixed time step depending on grid resolution + +# 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.time_avg_vel = true +erf.plot_file_1 = plt # prefix of plotfile name +erf.plot_int_1 = 10 # number of timesteps between plotfiles +erf.plot_vars_1 = density rhoKE rhoadv_0 x_velocity y_velocity z_velocity pressure temp theta u_t_avg v_t_avg w_t_avg umag_t_avg + +# SOLVER CHOICE +erf.alpha_T = 0.0 +erf.alpha_C = 1.0 +erf.use_gravity = false + +erf.molec_diff_type = "None" +erf.les_type = "Deardorff" +erf.Ck = 0.1 +erf.sigma_k = 1.0 +erf.Ce = 0.1 + +erf.init_type = "uniform" +erf.KE_0 = 0.1 # for Deardorff + +# 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/DevTests/MetGrid/inputs b/Exec/DevTests/MetGrid/inputs index ccf0ee7dc..6219d7d40 100644 --- a/Exec/DevTests/MetGrid/inputs +++ b/Exec/DevTests/MetGrid/inputs @@ -44,7 +44,7 @@ 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 dens_hse rhoadv_0 x_velocity y_velocity z_velocity pressure temp theta z_phys mapfac pres_hse pert_pres +erf.plot_vars_1 = lat_m lon_m density dens_hse rhoadv_0 x_velocity y_velocity z_velocity pressure temp theta z_phys mapfac pres_hse pert_pres # SOLVER CHOICE erf.alpha_T = 1.0 diff --git a/Exec/EWP/CMakeLists.txt b/Exec/EWP/CMakeLists.txt new file mode 100644 index 000000000..2e6da768c --- /dev/null +++ b/Exec/EWP/CMakeLists.txt @@ -0,0 +1,12 @@ +set(erf_exe_name erf_fitch) + +add_executable(${erf_exe_name} "") +target_sources(${erf_exe_name} + PRIVATE + prob.cpp +) + +target_include_directories(${erf_exe_name} PRIVATE ${CMAKE_CURRENT_SOURCE_DIR}) + +include(${CMAKE_SOURCE_DIR}/CMake/BuildERFExe.cmake) +build_erf_exe(${erf_exe_name}) diff --git a/Exec/EWP/GNUmakefile b/Exec/EWP/GNUmakefile new file mode 100644 index 000000000..991c5d326 --- /dev/null +++ b/Exec/EWP/GNUmakefile @@ -0,0 +1,37 @@ +# AMReX +COMP = gnu +PRECISION = DOUBLE + +# Profiling +PROFILE = FALSE +TINY_PROFILE = TRUE +COMM_PROFILE = FALSE +TRACE_PROFILE = FALSE +MEM_PROFILE = FALSE +USE_GPROF = FALSE + +# Performance +USE_MPI = TRUE +USE_OMP = FALSE + +USE_CUDA = FALSE +USE_HIP = FALSE +USE_SYCL = FALSE + +# Debugging +DEBUG = FALSE + +TEST = TRUE +USE_ASSERTION = TRUE + +USE_WINDFARM = TRUE + + +#USE_POISSON_SOLVE = TRUE + +# GNU Make +Bpack := ./Make.package +Blocs := . +ERF_HOME := ../.. +ERF_PROBLEM_DIR = $(ERF_HOME)/Exec/EWP +include $(ERF_HOME)/Exec/Make.ERF diff --git a/Exec/EWP/Make.package b/Exec/EWP/Make.package new file mode 100644 index 000000000..aeacb72f0 --- /dev/null +++ b/Exec/EWP/Make.package @@ -0,0 +1,2 @@ +CEXE_headers += prob.H +CEXE_sources += prob.cpp diff --git a/Exec/EWP/README b/Exec/EWP/README new file mode 100644 index 000000000..2504f4de7 --- /dev/null +++ b/Exec/EWP/README @@ -0,0 +1,6 @@ +This problem setup is for simulation of the Atmospheric Boundary Layer (ABL) +using one of two turbulence schemes (Smagorinsky or Deardorff) and the bottom +boundary condition possibly specified by Monin Obukhov Similarity Theory (MOST). + +This version of the ABL problem initializes the data using a hydrostatic profile +with random perturbations in velocity and potential temperature. diff --git a/Exec/EWP/inputs_EWP_WindFarm b/Exec/EWP/inputs_EWP_WindFarm new file mode 100644 index 000000000..b7284ee28 --- /dev/null +++ b/Exec/EWP/inputs_EWP_WindFarm @@ -0,0 +1,78 @@ +# ------------------ INPUTS TO MAIN PROGRAM ------------------- +max_step = 100000 + +amrex.fpe_trap_invalid = 1 + +fabarray.mfiter_tile_size = 1024 1024 1024 + +# PROBLEM SIZE & GEOMETRY +erf.latitude_lo = 35.0 +erf.longitude_lo = -100.0 +geometry.prob_extent = 200150 202637 1000 +amr.n_cell = 50 50 40 + +#erf.grid_stretching_ratio = 1.025 +#erf.initial_dz = 16.0 + +geometry.is_periodic = 0 0 0 + +# MOST BOUNDARY (DEFAULT IS ADIABATIC FOR THETA) +#zlo.type = "MOST" +#erf.most.z0 = 0.1 +#erf.most.zref = 8.0 + +zlo.type = "SlipWall" +zhi.type = "SlipWall" +xlo.type = "Outflow" +xhi.type = "Outflow" +ylo.type = "Outflow" +yhi.type = "Outflow" + +# TIME STEP CONTROL +erf.use_native_mri = 1 +erf.fixed_dt = 0.25 # fixed time step depending on grid resolution +#erf.fixed_fast_dt = 0.0025 + +# 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 = 10000 # number of timesteps between checkpoints +#erf.restart = chk02000 + +# PLOTFILES +erf.plot_file_1 = plt # prefix of plotfile name +erf.plot_int_1 = 1000 # number of timesteps between plotfiles +erf.plot_vars_1 = density rhoadv_0 x_velocity y_velocity z_velocity pressure temp theta QKE num_turb vorticity + +# SOLVER CHOICE +erf.alpha_T = 0.0 +erf.alpha_C = 1.0 +erf.use_gravity = false + +erf.molec_diff_type = "ConstantAlpha" +erf.les_type = "None" +erf.Cs = 1.5 +erf.dynamicViscosity = 100.0 + +erf.pbl_type = "None" + +erf.init_type = "uniform" + +erf.windfarm_type = "EWP" + +# PROBLEM PARAMETERS +prob.rho_0 = 1.0 +prob.A_0 = 1.0 + +prob.U_0 = 10.0 +prob.V_0 = 10.0 +prob.W_0 = 0.0 +prob.T_0 = 300.0 + diff --git a/Exec/EWP/plot_profiles.py b/Exec/EWP/plot_profiles.py new file mode 100644 index 000000000..50ab6f76c --- /dev/null +++ b/Exec/EWP/plot_profiles.py @@ -0,0 +1,70 @@ +#!/usr/bin/env python +import os +import glob +import numpy as np +import matplotlib.pyplot as plt +import yt + +pltfiles = [pltfile for pltfile in glob.glob('plt*') + if os.path.isdir(pltfile) and not ('old' in pltfile)] +nsteps = [int(pltfile[3:]) for pltfile in pltfiles] +latestoutput = pltfiles[np.argmax(nsteps)] +print(latestoutput) + +ds = yt.load(latestoutput) + +soln = ds.covering_grid(level=0, left_edge=ds.domain_left_edge, dims=ds.domain_dimensions) +zp = np.linspace(ds.domain_left_edge[2], + ds.domain_right_edge[2], + ds.domain_dimensions[2]+1) +z = (zp[1:] + zp[:-1]).value / 2 + +utot = soln['x_velocity'].value +vtot = soln['y_velocity'].value +wtot = soln['z_velocity'].value + +U = np.mean(utot,axis=(0,1)) +V = np.mean(vtot,axis=(0,1)) +W = np.mean(wtot,axis=(0,1)) +u = utot - U[np.newaxis,np.newaxis,:] +v = vtot - V[np.newaxis,np.newaxis,:] +w = wtot - W[np.newaxis,np.newaxis,:] + +uu = np.var(u, axis=(0,1)) +vv = np.var(v, axis=(0,1)) +ww = np.var(w, axis=(0,1)) + +uw = np.mean(u*w, axis=(0,1)) +vw = np.mean(v*w, axis=(0,1)) +uv = np.mean(u*v, axis=(0,1)) + +fig,ax = plt.subplots(nrows=3,ncols=3,sharey=True,figsize=(8.5,11)) +ax[0,0].plot(U,z) +ax[0,1].plot(V,z) +ax[0,2].plot(W,z) +ax[0,0].set_xlabel(r'$\langle U \rangle$ [m/s]') +ax[0,1].set_xlabel(r'$\langle V \rangle$ [m/s]') +ax[0,2].set_xlabel(r'$\langle W \rangle$ [m/s]') +ax[1,0].plot(uu,z) +ax[1,1].plot(vv,z) +ax[1,2].plot(ww,z) +ax[1,0].set_xlabel(r"$\langle u'u' \rangle$ [m/s]") +ax[1,1].set_xlabel(r"$\langle v'v' \rangle$ [m/s]") +ax[1,2].set_xlabel(r"$\langle w'w' \rangle$ [m/s]") +ax[2,0].plot(uw,z) +ax[2,1].plot(vw,z) +ax[2,2].plot(uv,z) +ax[2,0].set_xlabel(r"$\langle u'w' \rangle$ [m/s]") +ax[2,1].set_xlabel(r"$\langle v'w' \rangle$ [m/s]") +ax[2,2].set_xlabel(r"$\langle u'v' \rangle$ [m/s]") +ax[0,0].set_ylabel('$z$ [m]') +ax[1,0].set_ylabel('$z$ [m]') +ax[2,0].set_ylabel('$z$ [m]') +ax[0,0].set_ylim((ds.domain_left_edge[2], ds.domain_right_edge[2])) +for axi in ax.ravel(): + axi.grid() + +fig.savefig('profiles.png',bbox_inches='tight') + +plt.show() + diff --git a/Exec/EWP/prob.H b/Exec/EWP/prob.H new file mode 100644 index 000000000..aa0c28435 --- /dev/null +++ b/Exec/EWP/prob.H @@ -0,0 +1,79 @@ +#ifndef _PROB_H_ +#define _PROB_H_ + +#include + +#include "AMReX_REAL.H" + +#include "prob_common.H" + +struct ProbParm : ProbParmDefaults { + amrex::Real rho_0 = 0.0; + amrex::Real T_0 = 0.0; + amrex::Real A_0 = 1.0; + amrex::Real QKE_0 = 0.1; + + amrex::Real U_0 = 0.0; + amrex::Real V_0 = 0.0; + amrex::Real W_0 = 0.0; + + // random initial perturbations (legacy code) + amrex::Real U_0_Pert_Mag = 0.0; + amrex::Real V_0_Pert_Mag = 0.0; + amrex::Real W_0_Pert_Mag = 0.0; + amrex::Real T_0_Pert_Mag = 0.0; // perturbation to rho*Theta + + // divergence-free initial perturbations + amrex::Real pert_deltaU = 0.0; + amrex::Real pert_deltaV = 0.0; + amrex::Real pert_periods_U = 5.0; + amrex::Real pert_periods_V = 5.0; + amrex::Real pert_ref_height = 100.0; + + // rayleigh damping + amrex::Real dampcoef = 0.2; // inverse time scale [1/s] + amrex::Real zdamp = 500.0; // damping depth [m] from model top + + // helper vars + amrex::Real aval; + amrex::Real bval; + amrex::Real ufac; + amrex::Real vfac; +}; // namespace ProbParm + +class Problem : public ProblemBase +{ +public: + Problem(const amrex::Real* problo, const amrex::Real* probhi); + +#include "Prob/init_constant_density_hse.H" +#include "Prob/init_rayleigh_damping.H" + + void init_custom_pert ( + const amrex::Box& bx, + const amrex::Box& xbx, + const amrex::Box& ybx, + const amrex::Box& zbx, + amrex::Array4 const& state, + amrex::Array4 const& state_pert, + amrex::Array4 const& x_vel_pert, + amrex::Array4 const& y_vel_pert, + amrex::Array4 const& z_vel_pert, + amrex::Array4 const& r_hse, + amrex::Array4 const& p_hse, + amrex::Array4 const& z_nd, + amrex::Array4 const& z_cc, + amrex::GeometryData const& geomdata, + amrex::Array4 const& mf_m, + amrex::Array4 const& mf_u, + amrex::Array4 const& mf_v, + const SolverChoice& sc) override; + +protected: + std::string name() override { return "ABL"; } + +private: + ProbParm parms; +}; + +#endif diff --git a/Exec/EWP/prob.cpp b/Exec/EWP/prob.cpp new file mode 100644 index 000000000..736211998 --- /dev/null +++ b/Exec/EWP/prob.cpp @@ -0,0 +1,167 @@ +#include "prob.H" +#include "AMReX_Random.H" + +using namespace amrex; + +std::unique_ptr +amrex_probinit(const amrex_real* problo, const amrex_real* probhi) +{ + return std::make_unique(problo, probhi); +} + +Problem::Problem(const amrex::Real* problo, const amrex::Real* probhi) +{ + // Parse params + ParmParse pp("prob"); + pp.query("rho_0", parms.rho_0); + pp.query("T_0", parms.T_0); + pp.query("A_0", parms.A_0); + pp.query("QKE_0", parms.QKE_0); + + pp.query("U_0", parms.U_0); + pp.query("V_0", parms.V_0); + pp.query("W_0", parms.W_0); + pp.query("U_0_Pert_Mag", parms.U_0_Pert_Mag); + pp.query("V_0_Pert_Mag", parms.V_0_Pert_Mag); + pp.query("W_0_Pert_Mag", parms.W_0_Pert_Mag); + pp.query("T_0_Pert_Mag", parms.T_0_Pert_Mag); + + pp.query("pert_deltaU", parms.pert_deltaU); + pp.query("pert_deltaV", parms.pert_deltaV); + pp.query("pert_periods_U", parms.pert_periods_U); + pp.query("pert_periods_V", parms.pert_periods_V); + pp.query("pert_ref_height", parms.pert_ref_height); + parms.aval = parms.pert_periods_U * 2.0 * PI / (probhi[1] - problo[1]); + parms.bval = parms.pert_periods_V * 2.0 * PI / (probhi[0] - problo[0]); + parms.ufac = parms.pert_deltaU * std::exp(0.5) / parms.pert_ref_height; + parms.vfac = parms.pert_deltaV * std::exp(0.5) / parms.pert_ref_height; + + pp.query("dampcoef", parms.dampcoef); + pp.query("zdamp", parms.zdamp); + + init_base_parms(parms.rho_0, parms.T_0); +} + +void +Problem::init_custom_pert( + const amrex::Box& bx, + const amrex::Box& xbx, + const amrex::Box& ybx, + const amrex::Box& zbx, + amrex::Array4 const& /*state*/, + amrex::Array4 const& state_pert, + amrex::Array4 const& x_vel_pert, + amrex::Array4 const& y_vel_pert, + amrex::Array4 const& z_vel_pert, + amrex::Array4 const& /*r_hse*/, + amrex::Array4 const& /*p_hse*/, + amrex::Array4 const& /*z_nd*/, + amrex::Array4 const& /*z_cc*/, + amrex::GeometryData const& geomdata, + amrex::Array4 const& /*mf_m*/, + amrex::Array4 const& /*mf_u*/, + amrex::Array4 const& /*mf_v*/, + const SolverChoice& sc) +{ + const bool use_moisture = (sc.moisture_type != MoistureType::None); + + ParallelForRNG(bx, [=, parms=parms] AMREX_GPU_DEVICE(int i, int j, int k, const amrex::RandomEngine& engine) noexcept { + // Geometry + const Real* prob_lo = geomdata.ProbLo(); + const Real* prob_hi = geomdata.ProbHi(); + const Real* dx = geomdata.CellSize(); + const Real x = prob_lo[0] + (i + 0.5) * dx[0]; + const Real y = prob_lo[1] + (j + 0.5) * dx[1]; + const Real z = prob_lo[2] + (k + 0.5) * dx[2]; + + // Define a point (xc,yc,zc) at the center of the domain + const Real xc = 0.5 * (prob_lo[0] + prob_hi[0]); + const Real yc = 0.5 * (prob_lo[1] + prob_hi[1]); + const Real zc = 0.5 * (prob_lo[2] + prob_hi[2]); + + const Real r = std::sqrt((x-xc)*(x-xc) + (y-yc)*(y-yc) + (z-zc)*(z-zc)); + + // Add temperature perturbations + if ((z <= parms.pert_ref_height) && (parms.T_0_Pert_Mag != 0.0)) { + Real rand_double = amrex::Random(engine); // Between 0.0 and 1.0 + state_pert(i, j, k, RhoTheta_comp) = (rand_double*2.0 - 1.0)*parms.T_0_Pert_Mag; + } + + // Set scalar = A_0*exp(-10r^2), where r is distance from center of domain + state_pert(i, j, k, RhoScalar_comp) = parms.A_0 * exp(-10.*r*r); + + // Set an initial value for QKE + state_pert(i, j, k, RhoQKE_comp) = parms.QKE_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 + ParallelForRNG(xbx, [=, parms=parms] AMREX_GPU_DEVICE(int i, int j, int k, const amrex::RandomEngine& engine) noexcept { + const Real* prob_lo = geomdata.ProbLo(); + const Real* dx = geomdata.CellSize(); + const Real y = prob_lo[1] + (j + 0.5) * dx[1]; + const Real z = prob_lo[2] + (k + 0.5) * dx[2]; + + // Set the x-velocity + x_vel_pert(i, j, k) = parms.U_0; + if ((z <= parms.pert_ref_height) && (parms.U_0_Pert_Mag != 0.0)) + { + Real rand_double = amrex::Random(engine); // Between 0.0 and 1.0 + Real x_vel_prime = (rand_double*2.0 - 1.0)*parms.U_0_Pert_Mag; + x_vel_pert(i, j, k) += x_vel_prime; + } + if (parms.pert_deltaU != 0.0) + { + const amrex::Real yl = y - prob_lo[1]; + const amrex::Real zl = z / parms.pert_ref_height; + const amrex::Real damp = std::exp(-0.5 * zl * zl); + x_vel_pert(i, j, k) += parms.ufac * damp * z * std::cos(parms.aval * yl); + } + }); + + // Set the y-velocity + ParallelForRNG(ybx, [=, parms=parms] AMREX_GPU_DEVICE(int i, int j, int k, const amrex::RandomEngine& engine) noexcept { + const Real* prob_lo = geomdata.ProbLo(); + const Real* dx = geomdata.CellSize(); + const Real x = prob_lo[0] + (i + 0.5) * dx[0]; + const Real z = prob_lo[2] + (k + 0.5) * dx[2]; + + // Set the y-velocity + y_vel_pert(i, j, k) = parms.V_0; + if ((z <= parms.pert_ref_height) && (parms.V_0_Pert_Mag != 0.0)) + { + Real rand_double = amrex::Random(engine); // Between 0.0 and 1.0 + Real y_vel_prime = (rand_double*2.0 - 1.0)*parms.V_0_Pert_Mag; + y_vel_pert(i, j, k) += y_vel_prime; + } + if (parms.pert_deltaV != 0.0) + { + const amrex::Real xl = x - prob_lo[0]; + const amrex::Real zl = z / parms.pert_ref_height; + const amrex::Real damp = std::exp(-0.5 * zl * zl); + y_vel_pert(i, j, k) += parms.vfac * damp * z * std::cos(parms.bval * xl); + } + }); + + // Set the z-velocity + ParallelForRNG(zbx, [=, parms=parms] AMREX_GPU_DEVICE(int i, int j, int k, const amrex::RandomEngine& engine) noexcept { + const int dom_lo_z = geomdata.Domain().smallEnd()[2]; + const int dom_hi_z = geomdata.Domain().bigEnd()[2]; + + // Set the z-velocity + if (k == dom_lo_z || k == dom_hi_z+1) + { + z_vel_pert(i, j, k) = 0.0; + } + else if (parms.W_0_Pert_Mag != 0.0) + { + Real rand_double = amrex::Random(engine); // Between 0.0 and 1.0 + Real z_vel_prime = (rand_double*2.0 - 1.0)*parms.W_0_Pert_Mag; + z_vel_pert(i, j, k) = parms.W_0 + z_vel_prime; + } + }); +} diff --git a/Exec/EWP/wind-turbine-1.tbl b/Exec/EWP/wind-turbine-1.tbl new file mode 100644 index 000000000..6d678c3d1 --- /dev/null +++ b/Exec/EWP/wind-turbine-1.tbl @@ -0,0 +1,24 @@ +22 +75. 85. 0.130 2.0 +4. 0.805 50.0 +5. 0.805 150.0 +6. 0.805 280.0 +7. 0.805 460.0 +8. 0.805 700.0 +9. 0.805 990.0 +10. 0.790 1300.0 +11. 0.740 1600.0 +12. 0.700 1850.0 +13. 0.400 1950.0 +14. 0.300 1990.0 +15. 0.250 1995.0 +16. 0.200 2000.0 +17. 0.160 2000.0 +18. 0.140 2000.0 +19. 0.120 2000.0 +20. 0.100 2000.0 +21. 0.080 2000.0 +22. 0.070 2000.0 +23. 0.060 2000.0 +24. 0.055 2000.0 +25. 0.050 2000.0 diff --git a/Exec/EWP/windturbines.txt b/Exec/EWP/windturbines.txt new file mode 100644 index 000000000..4c68601f6 --- /dev/null +++ b/Exec/EWP/windturbines.txt @@ -0,0 +1,97 @@ +35.7828828829 -99.0168 1 +35.8219219219 -99.0168 1 +35.860960961 -99.0168 1 +35.9 -99.0168 1 +35.939039039 -99.0168 1 +35.9780780781 -99.0168 1 +36.0171171171 -99.0168 1 +35.7828828829 -98.9705333333 1 +35.8219219219 -98.9705333333 1 +35.860960961 -98.9705333333 1 +35.9 -98.9705333333 1 +35.939039039 -98.9705333333 1 +35.9780780781 -98.9705333333 1 +36.0171171171 -98.9705333333 1 +35.7828828829 -98.9242666667 1 +35.8219219219 -98.9242666667 1 +35.860960961 -98.9242666667 1 +35.9 -98.9242666667 1 +35.939039039 -98.9242666667 1 +35.9780780781 -98.9242666667 1 +36.0171171171 -98.9242666667 1 +35.7828828829 -98.878 1 +35.8219219219 -98.878 1 +35.860960961 -98.878 1 +35.9 -98.878 1 +35.939039039 -98.878 1 +35.9780780781 -98.878 1 +36.0171171171 -98.878 1 +35.7828828829 -98.8317333333 1 +35.8219219219 -98.8317333333 1 +35.860960961 -98.8317333333 1 +35.9 -98.8317333333 1 +35.939039039 -98.8317333333 1 +35.9780780781 -98.8317333333 1 +36.0171171171 -98.8317333333 1 +35.7828828829 -98.7854666667 1 +35.8219219219 -98.7854666667 1 +35.860960961 -98.7854666667 1 +35.9 -98.7854666667 1 +35.939039039 -98.7854666667 1 +35.9780780781 -98.7854666667 1 +36.0171171171 -98.7854666667 1 +35.7828828829 -98.7392 1 +35.8219219219 -98.7392 1 +35.860960961 -98.7392 1 +35.9 -98.7392 1 +35.939039039 -98.7392 1 +35.9780780781 -98.7392 1 +36.0171171171 -98.7392 1 +35.463963964 -98.2228 1 +35.487987988 -98.2228 1 +35.512012012 -98.2228 1 +35.536036036 -98.2228 1 +35.463963964 -98.1929333333 1 +35.487987988 -98.1929333333 1 +35.512012012 -98.1929333333 1 +35.536036036 -98.1929333333 1 +35.463963964 -98.1630666667 1 +35.487987988 -98.1630666667 1 +35.512012012 -98.1630666667 1 +35.536036036 -98.1630666667 1 +35.463963964 -98.1332 1 +35.487987988 -98.1332 1 +35.512012012 -98.1332 1 +35.536036036 -98.1332 1 +36.363963964 -99.4228 1 +36.387987988 -99.4228 1 +36.412012012 -99.4228 1 +36.436036036 -99.4228 1 +36.363963964 -99.3929333333 1 +36.387987988 -99.3929333333 1 +36.412012012 -99.3929333333 1 +36.436036036 -99.3929333333 1 +36.363963964 -99.3630666667 1 +36.387987988 -99.3630666667 1 +36.412012012 -99.3630666667 1 +36.436036036 -99.3630666667 1 +36.363963964 -99.3332 1 +36.387987988 -99.3332 1 +36.412012012 -99.3332 1 +36.436036036 -99.3332 1 +35.263963964 -99.5228 1 +35.287987988 -99.5228 1 +35.312012012 -99.5228 1 +35.336036036 -99.5228 1 +35.263963964 -99.4929333333 1 +35.287987988 -99.4929333333 1 +35.312012012 -99.4929333333 1 +35.336036036 -99.4929333333 1 +35.263963964 -99.4630666667 1 +35.287987988 -99.4630666667 1 +35.312012012 -99.4630666667 1 +35.336036036 -99.4630666667 1 +35.263963964 -99.4332 1 +35.287987988 -99.4332 1 +35.312012012 -99.4332 1 +35.336036036 -99.4332 1 diff --git a/Exec/EWP/windturbines.txt_LargeDomain b/Exec/EWP/windturbines.txt_LargeDomain new file mode 100644 index 000000000..586e94b81 --- /dev/null +++ b/Exec/EWP/windturbines.txt_LargeDomain @@ -0,0 +1,30 @@ +36.732 -97.128 1 +36.824 -97.367 1 +36.967 -96.754 1 +36.607 -96.846 1 +36.719 -97.249 1 +36.889 -97.031 1 +36.621 -96.932 1 +36.752 -96.788 1 +36.855 -97.421 1 +36.694 -96.656 1 +36.579 -97.223 1 +36.912 -97.379 1 +36.721 -96.815 1 +36.685 -97.437 1 +36.781 -97.049 1 +36.936 -97.203 1 +36.541 -96.791 1 +36.655 -97.284 1 +36.798 -96.695 1 +36.887 -96.956 1 +36.633 -96.716 1 +36.912 -96.813 1 +36.725 -96.579 1 +36.594 -97.184 1 +36.902 -97.354 1 +36.628 -96.871 1 +36.762 -97.428 1 +36.811 -97.285 1 +36.549 -97.492 1 +36.642 -96.673 1 diff --git a/Exec/EWP/windturbines_WindFarm.txt b/Exec/EWP/windturbines_WindFarm.txt new file mode 100644 index 000000000..4c68601f6 --- /dev/null +++ b/Exec/EWP/windturbines_WindFarm.txt @@ -0,0 +1,97 @@ +35.7828828829 -99.0168 1 +35.8219219219 -99.0168 1 +35.860960961 -99.0168 1 +35.9 -99.0168 1 +35.939039039 -99.0168 1 +35.9780780781 -99.0168 1 +36.0171171171 -99.0168 1 +35.7828828829 -98.9705333333 1 +35.8219219219 -98.9705333333 1 +35.860960961 -98.9705333333 1 +35.9 -98.9705333333 1 +35.939039039 -98.9705333333 1 +35.9780780781 -98.9705333333 1 +36.0171171171 -98.9705333333 1 +35.7828828829 -98.9242666667 1 +35.8219219219 -98.9242666667 1 +35.860960961 -98.9242666667 1 +35.9 -98.9242666667 1 +35.939039039 -98.9242666667 1 +35.9780780781 -98.9242666667 1 +36.0171171171 -98.9242666667 1 +35.7828828829 -98.878 1 +35.8219219219 -98.878 1 +35.860960961 -98.878 1 +35.9 -98.878 1 +35.939039039 -98.878 1 +35.9780780781 -98.878 1 +36.0171171171 -98.878 1 +35.7828828829 -98.8317333333 1 +35.8219219219 -98.8317333333 1 +35.860960961 -98.8317333333 1 +35.9 -98.8317333333 1 +35.939039039 -98.8317333333 1 +35.9780780781 -98.8317333333 1 +36.0171171171 -98.8317333333 1 +35.7828828829 -98.7854666667 1 +35.8219219219 -98.7854666667 1 +35.860960961 -98.7854666667 1 +35.9 -98.7854666667 1 +35.939039039 -98.7854666667 1 +35.9780780781 -98.7854666667 1 +36.0171171171 -98.7854666667 1 +35.7828828829 -98.7392 1 +35.8219219219 -98.7392 1 +35.860960961 -98.7392 1 +35.9 -98.7392 1 +35.939039039 -98.7392 1 +35.9780780781 -98.7392 1 +36.0171171171 -98.7392 1 +35.463963964 -98.2228 1 +35.487987988 -98.2228 1 +35.512012012 -98.2228 1 +35.536036036 -98.2228 1 +35.463963964 -98.1929333333 1 +35.487987988 -98.1929333333 1 +35.512012012 -98.1929333333 1 +35.536036036 -98.1929333333 1 +35.463963964 -98.1630666667 1 +35.487987988 -98.1630666667 1 +35.512012012 -98.1630666667 1 +35.536036036 -98.1630666667 1 +35.463963964 -98.1332 1 +35.487987988 -98.1332 1 +35.512012012 -98.1332 1 +35.536036036 -98.1332 1 +36.363963964 -99.4228 1 +36.387987988 -99.4228 1 +36.412012012 -99.4228 1 +36.436036036 -99.4228 1 +36.363963964 -99.3929333333 1 +36.387987988 -99.3929333333 1 +36.412012012 -99.3929333333 1 +36.436036036 -99.3929333333 1 +36.363963964 -99.3630666667 1 +36.387987988 -99.3630666667 1 +36.412012012 -99.3630666667 1 +36.436036036 -99.3630666667 1 +36.363963964 -99.3332 1 +36.387987988 -99.3332 1 +36.412012012 -99.3332 1 +36.436036036 -99.3332 1 +35.263963964 -99.5228 1 +35.287987988 -99.5228 1 +35.312012012 -99.5228 1 +35.336036036 -99.5228 1 +35.263963964 -99.4929333333 1 +35.287987988 -99.4929333333 1 +35.312012012 -99.4929333333 1 +35.336036036 -99.4929333333 1 +35.263963964 -99.4630666667 1 +35.287987988 -99.4630666667 1 +35.312012012 -99.4630666667 1 +35.336036036 -99.4630666667 1 +35.263963964 -99.4332 1 +35.287987988 -99.4332 1 +35.312012012 -99.4332 1 +35.336036036 -99.4332 1 diff --git a/Exec/Fitch/inputs_Fitch_WindFarm b/Exec/Fitch/inputs_Fitch_WindFarm new file mode 100644 index 000000000..e901933cb --- /dev/null +++ b/Exec/Fitch/inputs_Fitch_WindFarm @@ -0,0 +1,78 @@ +# ------------------ INPUTS TO MAIN PROGRAM ------------------- +max_step = 100000 + +amrex.fpe_trap_invalid = 1 + +fabarray.mfiter_tile_size = 1024 1024 1024 + +# PROBLEM SIZE & GEOMETRY +erf.latitude_lo = 35.0 +erf.longitude_lo = -100.0 +geometry.prob_extent = 200150 202637 1000 +amr.n_cell = 50 50 40 + +#erf.grid_stretching_ratio = 1.025 +#erf.initial_dz = 16.0 + +geometry.is_periodic = 0 0 0 + +# MOST BOUNDARY (DEFAULT IS ADIABATIC FOR THETA) +#zlo.type = "MOST" +#erf.most.z0 = 0.1 +#erf.most.zref = 8.0 + +zlo.type = "SlipWall" +zhi.type = "SlipWall" +xlo.type = "Outflow" +xhi.type = "Outflow" +ylo.type = "Outflow" +yhi.type = "Outflow" + +# TIME STEP CONTROL +erf.use_native_mri = 1 +erf.fixed_dt = 0.25 # fixed time step depending on grid resolution +#erf.fixed_fast_dt = 0.0025 + +# 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 = 10000 # number of timesteps between checkpoints +#erf.restart = chk02000 + +# PLOTFILES +erf.plot_file_1 = plt # prefix of plotfile name +erf.plot_int_1 = 1000 # number of timesteps between plotfiles +erf.plot_vars_1 = density rhoadv_0 x_velocity y_velocity z_velocity pressure temp theta QKE num_turb vorticity + +# SOLVER CHOICE +erf.alpha_T = 0.0 +erf.alpha_C = 1.0 +erf.use_gravity = false + +erf.molec_diff_type = "ConstantAlpha" +erf.les_type = "None" +erf.Cs = 1.5 +erf.dynamicViscosity = 100.0 + +erf.pbl_type = "None" + +erf.init_type = "uniform" + +erf.windfarm_type = "Fitch" + +# PROBLEM PARAMETERS +prob.rho_0 = 1.0 +prob.A_0 = 1.0 + +prob.U_0 = 10.0 +prob.V_0 = 10.0 +prob.W_0 = 0.0 +prob.T_0 = 300.0 + diff --git a/Exec/Fitch/windturbines_WindFarm.txt b/Exec/Fitch/windturbines_WindFarm.txt new file mode 100644 index 000000000..4c68601f6 --- /dev/null +++ b/Exec/Fitch/windturbines_WindFarm.txt @@ -0,0 +1,97 @@ +35.7828828829 -99.0168 1 +35.8219219219 -99.0168 1 +35.860960961 -99.0168 1 +35.9 -99.0168 1 +35.939039039 -99.0168 1 +35.9780780781 -99.0168 1 +36.0171171171 -99.0168 1 +35.7828828829 -98.9705333333 1 +35.8219219219 -98.9705333333 1 +35.860960961 -98.9705333333 1 +35.9 -98.9705333333 1 +35.939039039 -98.9705333333 1 +35.9780780781 -98.9705333333 1 +36.0171171171 -98.9705333333 1 +35.7828828829 -98.9242666667 1 +35.8219219219 -98.9242666667 1 +35.860960961 -98.9242666667 1 +35.9 -98.9242666667 1 +35.939039039 -98.9242666667 1 +35.9780780781 -98.9242666667 1 +36.0171171171 -98.9242666667 1 +35.7828828829 -98.878 1 +35.8219219219 -98.878 1 +35.860960961 -98.878 1 +35.9 -98.878 1 +35.939039039 -98.878 1 +35.9780780781 -98.878 1 +36.0171171171 -98.878 1 +35.7828828829 -98.8317333333 1 +35.8219219219 -98.8317333333 1 +35.860960961 -98.8317333333 1 +35.9 -98.8317333333 1 +35.939039039 -98.8317333333 1 +35.9780780781 -98.8317333333 1 +36.0171171171 -98.8317333333 1 +35.7828828829 -98.7854666667 1 +35.8219219219 -98.7854666667 1 +35.860960961 -98.7854666667 1 +35.9 -98.7854666667 1 +35.939039039 -98.7854666667 1 +35.9780780781 -98.7854666667 1 +36.0171171171 -98.7854666667 1 +35.7828828829 -98.7392 1 +35.8219219219 -98.7392 1 +35.860960961 -98.7392 1 +35.9 -98.7392 1 +35.939039039 -98.7392 1 +35.9780780781 -98.7392 1 +36.0171171171 -98.7392 1 +35.463963964 -98.2228 1 +35.487987988 -98.2228 1 +35.512012012 -98.2228 1 +35.536036036 -98.2228 1 +35.463963964 -98.1929333333 1 +35.487987988 -98.1929333333 1 +35.512012012 -98.1929333333 1 +35.536036036 -98.1929333333 1 +35.463963964 -98.1630666667 1 +35.487987988 -98.1630666667 1 +35.512012012 -98.1630666667 1 +35.536036036 -98.1630666667 1 +35.463963964 -98.1332 1 +35.487987988 -98.1332 1 +35.512012012 -98.1332 1 +35.536036036 -98.1332 1 +36.363963964 -99.4228 1 +36.387987988 -99.4228 1 +36.412012012 -99.4228 1 +36.436036036 -99.4228 1 +36.363963964 -99.3929333333 1 +36.387987988 -99.3929333333 1 +36.412012012 -99.3929333333 1 +36.436036036 -99.3929333333 1 +36.363963964 -99.3630666667 1 +36.387987988 -99.3630666667 1 +36.412012012 -99.3630666667 1 +36.436036036 -99.3630666667 1 +36.363963964 -99.3332 1 +36.387987988 -99.3332 1 +36.412012012 -99.3332 1 +36.436036036 -99.3332 1 +35.263963964 -99.5228 1 +35.287987988 -99.5228 1 +35.312012012 -99.5228 1 +35.336036036 -99.5228 1 +35.263963964 -99.4929333333 1 +35.287987988 -99.4929333333 1 +35.312012012 -99.4929333333 1 +35.336036036 -99.4929333333 1 +35.263963964 -99.4630666667 1 +35.287987988 -99.4630666667 1 +35.312012012 -99.4630666667 1 +35.336036036 -99.4630666667 1 +35.263963964 -99.4332 1 +35.287987988 -99.4332 1 +35.312012012 -99.4332 1 +35.336036036 -99.4332 1 diff --git a/Exec/Make.ERF b/Exec/Make.ERF index 6fd4eca97..328613d09 100644 --- a/Exec/Make.ERF +++ b/Exec/Make.ERF @@ -142,9 +142,13 @@ INCLUDE_LOCATIONS += $(ERF_MOISTURE_KESSLER_DIR) ifeq ($(USE_WINDFARM), TRUE) DEFINES += -DERF_USE_WINDFARM ERF_WINDFARM_FITCH_DIR = $(ERF_SOURCE_DIR)/WindFarmParametrization/Fitch + ERF_WINDFARM_EWP_DIR = $(ERF_SOURCE_DIR)/WindFarmParametrization/EWP include $(ERF_WINDFARM_FITCH_DIR)/Make.package + include $(ERF_WINDFARM_EWP_DIR)/Make.package VPATH_LOCATIONS += $(ERF_WINDFARM_FITCH_DIR) INCLUDE_LOCATIONS += $(ERF_WINDFARM_FITCH_DIR) + VPATH_LOCATIONS += $(ERF_WINDFARM_EWP_DIR) + INCLUDE_LOCATIONS += $(ERF_WINDFARM_EWP_DIR) endif ERF_LSM_DIR = $(ERF_SOURCE_DIR)/LandSurfaceModel diff --git a/Exec/RegTests/WPS_Test/inputs_real_ChisholmView b/Exec/RegTests/WPS_Test/inputs_real_ChisholmView index dfd5efb4f..de034cab9 100644 --- a/Exec/RegTests/WPS_Test/inputs_real_ChisholmView +++ b/Exec/RegTests/WPS_Test/inputs_real_ChisholmView @@ -39,7 +39,7 @@ erf.check_int = -100 # number of timesteps between checkpoints # PLOTFILES erf.plot_file_1 = plt # prefix of plotfile name erf.plot_int_1 = 1 # number of timesteps between plotfiles -erf.plot_vars_1 = density rhoQ1 x_velocity y_velocity z_velocity pressure temp theta pres_hse z_phys qt qp qv qc +erf.plot_vars_1 = lat_m lon_m density rhoQ1 x_velocity y_velocity z_velocity pressure temp theta pres_hse z_phys qt qp qv qc # SOLVER CHOICE erf.alpha_T = 0.0 diff --git a/Exec/SpongeTest/CMakeLists.txt b/Exec/SpongeTest/CMakeLists.txt new file mode 100644 index 000000000..54678203c --- /dev/null +++ b/Exec/SpongeTest/CMakeLists.txt @@ -0,0 +1,15 @@ +set(erf_exe_name erf_sponge_test) + +add_executable(${erf_exe_name} "") +target_sources(${erf_exe_name} + PRIVATE + prob.cpp +) + +target_include_directories(${erf_exe_name} PRIVATE ${CMAKE_CURRENT_SOURCE_DIR}) + +include(${CMAKE_SOURCE_DIR}/CMake/BuildERFExe.cmake) +build_erf_exe(${erf_exe_name}) + +#find_package(AMReX) +#target_link_libraries( ${_target} AMReX::amrex) diff --git a/Exec/SpongeTest/GNUmakefile b/Exec/SpongeTest/GNUmakefile new file mode 100644 index 000000000..aa571732d --- /dev/null +++ b/Exec/SpongeTest/GNUmakefile @@ -0,0 +1,32 @@ +# AMReX +COMP = gnu +PRECISION = DOUBLE + +# Profiling +PROFILE = FALSE +TINY_PROFILE = FALSE +COMM_PROFILE = FALSE +TRACE_PROFILE = FALSE +MEM_PROFILE = FALSE +USE_GPROF = FALSE + +# Performance +USE_MPI = TRUE +USE_OMP = FALSE + +USE_CUDA = FALSE +USE_HIP = FALSE +USE_SYCL = FALSE + +# Debugging +DEBUG = FALSE + +TEST = TRUE +USE_ASSERTION = TRUE + +# GNU Make +Bpack := ./Make.package +Blocs := . +ERF_HOME := ../.. +ERF_PROBLEM_DIR = $(ERF_HOME)/Exec/SpongeTest +include $(ERF_HOME)/Exec/Make.ERF diff --git a/Exec/SpongeTest/Make.package b/Exec/SpongeTest/Make.package new file mode 100644 index 000000000..aeacb72f0 --- /dev/null +++ b/Exec/SpongeTest/Make.package @@ -0,0 +1,2 @@ +CEXE_headers += prob.H +CEXE_sources += prob.cpp diff --git a/Exec/SpongeTest/input_sponge_file.txt b/Exec/SpongeTest/input_sponge_file.txt new file mode 100644 index 000000000..0a32f830c --- /dev/null +++ b/Exec/SpongeTest/input_sponge_file.txt @@ -0,0 +1,11 @@ +0.0 10.0 0.0 +1.0 10.0 0.0 +2.0 10.0 0.0 +3.0 10.0 0.0 +4.0 10.0 0.0 +5.0 10.0 0.0 +6.0 10.0 0.0 +7.0 10.0 0.0 +8.0 10.0 0.0 +9.0 10.0 0.0 +10.0 10.0 0.0 diff --git a/Exec/SpongeTest/inputs_sponge_test b/Exec/SpongeTest/inputs_sponge_test new file mode 100644 index 000000000..8d72f7cfb --- /dev/null +++ b/Exec/SpongeTest/inputs_sponge_test @@ -0,0 +1,99 @@ +# ------------------ INPUTS TO MAIN PROGRAM ------------------- +max_step = 6000 + +amrex.fpe_trap_invalid = 1 + +fabarray.mfiter_tile_size = 1024 1024 1024 + +# PROBLEM SIZE & GEOMETRY +geometry.prob_lo = 0 0. 0. +geometry.prob_hi = 30.0 1. 10.0 + +amr.n_cell = 256 4 128 # dx=dy=dz=100 m, Straka et al 1993 / Xue et al 2000 + +geometry.is_periodic = 0 1 0 + +xlo.type = "Outflow" +xhi.type = "Outflow" +xlo.velocity = 1. 0. 0. +xlo.density = 1.16 +xlo.theta = 300. +xlo.scalar = 0. + +zlo.type = "SlipWall" +zhi.type = "Outflow" + +# SPONGE ZONE BCS + +erf.sponge_type = "input_sponge" +erf.input_sponge_file = "input_sponge_file.txt" +erf.sponge_strength = 1000.0 +erf.use_xlo_sponge_damping = true +erf.xlo_sponge_end = 4.0 +#erf.use_xhi_sponge_damping = true +#erf.xhi_sponge_start = 26.0 +#erf.use_zhi_sponge_damping = true +#erf.zhi_sponge_start = 8.0 + + +# TIME STEP CONTROL +erf.no_substepping = 1 +#erf.fixed_dt = 1E-5 +erf.cfl = 0.3 + +# 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 checkpoint +#erf.restart = chk06000 + +# PLOTFILES +erf.plot_file_1 = plt # prefix of plotfile name +erf.plot_int_1 = 100 # number of timesteps between plotfiles +erf.plot_vars_1 = density x_velocity y_velocity z_velocity pressure theta pres_hse dens_hse pert_pres pert_dens z_phys detJ dpdx dpdy pres_hse_x pres_hse_y + +# SOLVER CHOICE +erf.use_gravity = true +erf.use_coriolis = false +erf.les_type = "None" + +# MULTILEVEL +amr.max_level = 0 +amr.ref_ratio_vect = 2 2 1 + +erf.refinement_indicators = box1 +erf.box1.max_level = 1 +erf.box1.in_box_lo = 2. 0.25 +erf.box1.in_box_hi = 8. 0.75 + +# TERRRAIN GRID TYPE +erf.use_terrain = true +erf.terrain_smoothing = 0 + +#erf.terrain_z_levels = 0 0.0236057 0.0486738 0.0752909 0.103548 0.133541 0.165372 0.199145 0.234973 0.272972 0.313264 0.355978 0.401247 0.449212 0.500017 0.553815 0.610764 0.671027 0.734774 0.802181 0.873428 0.948703 1.0282 1.1121 1.20063 1.29397 1.39234 1.49595 1.605 1.71971 1.84029 1.96694 2.09987 2.23928 2.38537 2.53831 2.69828 2.86545 3.03996 3.22194 3.41151 3.60875 3.81373 4.02649 4.24706 4.47539 4.71146 4.95516 5.20639 5.46498 5.73072 6.00339 6.28271 6.56834 6.85995 7.15712 7.45942 7.76639 8.07751 8.39226 8.71007 9.03036 9.35252 9.67594 10 + +# USE WENO FOR ADVECTION +#erf.all_use_WENO = true +#erf.spatial_order_WENO = 5 + +# Diffusion coefficient from Straka, K = 75 m^2/s +erf.molec_diff_type = "ConstantAlpha" +erf.rho0_trans = 1.0 # [kg/m^3], used to convert input diffusivities +erf.dynamicViscosity = 0.0 # [kg/(m-s)] ==> nu = 75.0 m^2/s +erf.alpha_T = 0.0 # [m^2/s] + +#erf.abl_driver_type = "PressureGradient" +#erf.abl_pressure_grad = -0.2 0. 0. + +# PROBLEM PARAMETERS (optional) +prob.T_0 = 300.0 +prob.U_0 = 0.0 +prob.rho_0 = 1.16 diff --git a/Exec/SpongeTest/prob.H b/Exec/SpongeTest/prob.H new file mode 100644 index 000000000..3de42008b --- /dev/null +++ b/Exec/SpongeTest/prob.H @@ -0,0 +1,56 @@ +#ifndef _PROB_H_ +#define _PROB_H_ + +#include + +#include "AMReX_REAL.H" + +#include "prob_common.H" + +struct ProbParm : ProbParmDefaults { + amrex::Real rho_0 = 1.2; + amrex::Real T_0 = 300.0; // surface temperature == mean potential temperature + amrex::Real U_0 = 0.0; + amrex::Real V_0 = 0.0; +}; // namespace ProbParm + +class Problem : public ProblemBase +{ +public: + Problem(); + +#include "Prob/init_density_hse_dry.H" + + void init_custom_pert ( + const amrex::Box& bx, + const amrex::Box& xbx, + const amrex::Box& ybx, + const amrex::Box& zbx, + amrex::Array4 const& state, + amrex::Array4 const& state_pert, + amrex::Array4 const& x_vel_pert, + amrex::Array4 const& y_vel_pert, + amrex::Array4 const& z_vel_pert, + amrex::Array4 const& r_hse, + amrex::Array4 const& p_hse, + amrex::Array4 const& z_nd, + amrex::Array4 const& z_cc, + amrex::GeometryData const& geomdata, + amrex::Array4 const& mf_m, + amrex::Array4 const& mf_u, + amrex::Array4 const& mf_v, + const SolverChoice& sc) override; + + void init_custom_terrain ( + const amrex::Geometry& geom, + amrex::MultiFab& z_phys_nd, + const amrex::Real& time) override; + +protected: + std::string name() override { return "2D Terrain - Cylinder"; } + +private: + ProbParm parms; +}; + +#endif diff --git a/Exec/SpongeTest/prob.cpp b/Exec/SpongeTest/prob.cpp new file mode 100644 index 000000000..469ba545a --- /dev/null +++ b/Exec/SpongeTest/prob.cpp @@ -0,0 +1,140 @@ +#include "prob.H" +#include "EOS.H" +#include "TerrainMetrics.H" + +using namespace amrex; + +std::unique_ptr +amrex_probinit( + const amrex_real* /*problo*/, + const amrex_real* /*probhi*/) +{ + return std::make_unique(); +} + +Problem::Problem() +{ + // Parse params + ParmParse pp("prob"); + pp.query("rho_0", parms.rho_0); + pp.query("U_0", parms.U_0); + + init_base_parms(parms.rho_0, parms.T_0); +} + +void +Problem::init_custom_pert( + const Box& bx, + const Box& xbx, + const Box& ybx, + const Box& zbx, + Array4 const& /*state*/, + Array4 const& state_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& z_nd, + Array4 const& z_cc, + GeometryData const& geomdata, + Array4 const& /*mf_m*/, + Array4 const& /*mf_u*/, + Array4 const& /*mf_v*/, + const SolverChoice& sc) +{ + const bool use_moisture = (sc.moisture_type != MoistureType::None); + + // Geometry (note we must include these here to get the data on device) + ParallelFor(bx, [=, parms=parms] AMREX_GPU_DEVICE(int i, int j, int k) noexcept + { + 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); + + // 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; + } + }); + + // Set the x-velocity + ParallelFor(xbx, [=, parms=parms] AMREX_GPU_DEVICE(int i, int j, int k) noexcept + { + x_vel_pert(i, j, k) = 0.0; + }); + + // Set the y-velocity + ParallelFor(ybx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept + { + y_vel_pert(i, j, k) = 0.0; + }); + + // Set the z-velocity from impenetrable condition + ParallelFor(zbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept + { + z_vel_pert(i, j, k) = 0.0; + }); + + amrex::Gpu::streamSynchronize(); +} + +void +Problem::init_custom_terrain( + const Geometry& geom, + MultiFab& z_phys_nd, + 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) #include #include +#include #include #include #include @@ -307,6 +308,8 @@ public: void init_from_input_sounding (int lev); + void input_sponge (int lev); + void init_from_hse (int lev); #ifdef ERF_USE_MULTIBLOCK @@ -336,6 +339,11 @@ public: MultiBlockContainer *m_mbc = nullptr; amrex::Vector > vars_new; amrex::Vector > vars_old; + + // Velocity time averaged field + amrex::Vector> vel_t_avg; + amrex::Vector t_avg_cnt; + #endif std::string pp_prefix {"erf"}; @@ -454,9 +462,15 @@ private: //! Initialize Rayleigh damping profiles void initRayleigh (); + //! Initialize sponge profiles + void initSponge (); + //! Set Rayleigh mean profiles from input sounding void setRayleighRefFromSounding (bool restarting); + //! Set sponge mean profiles from input sounding + void setSpongeRefFromSounding (bool restarting); + // a wrapper for estTimeStep() void ComputeDt (); @@ -523,6 +537,7 @@ private: amrex::Vector> bdy_data_yhi; amrex::Real bdy_time_interval; + amrex::Vector> lat_m, lon_m; amrex::Real Latitude; amrex::Real Longitude; #endif // ERF_USE_NETCDF @@ -530,6 +545,9 @@ private: // Struct for working with the sounding data we take as an input InputSoundingData input_sounding_data; + // Struct for working with the sponge data we take as an input + InputSpongeData input_sponge_data; + // write checkpoint file to disk void WriteCheckpointFile () const; @@ -578,6 +596,10 @@ private: #ifndef ERF_USE_MULTIBLOCK amrex::Vector > vars_new; amrex::Vector > vars_old; + + // Velocity time averaged field + amrex::Vector> vel_t_avg; + amrex::Vector t_avg_cnt; #endif amrex::Vector > > > mri_integrator_mem; @@ -606,6 +628,7 @@ private: //#if defined(ERF_USE_WINDFARM) amrex::Vector Nturb; amrex::Vector vars_fitch; // Vabs, Vabsdt, dudt, dvdt, dTKEdt + amrex::Vector vars_ewp; // dudt, dvdt, dTKEdt amrex::Real hub_height, rotor_dia, thrust_coeff_standing, nominal_power; amrex::Vector wind_speed, thrust_coeff, power; @@ -616,7 +639,7 @@ private: #if defined(ERF_USE_RRTMGP) Radiation rad; - amrex::Vector qheating_rates; // radiation heating rate source terms + amrex::Vector> qheating_rates; // radiation heating rate source terms #endif // Fillpatcher classes for coarse-fine boundaries @@ -770,7 +793,9 @@ private: "magvel", "divU", "pres_hse", "dens_hse", "pressure", "pert_pres", "pert_dens", "eq_pot_temp", "num_turb", "dpdx", "dpdy", "pres_hse_x", "pres_hse_y", - "z_phys", "detJ" , "mapfac", + "z_phys", "detJ" , "mapfac", "lat_m", "lon_m", + // Time averaged velocity + "u_t_avg", "v_t_avg", "w_t_avg", "umag_t_avg", // eddy viscosity "Kmv","Kmh", // eddy diffusivity of heat @@ -821,6 +846,9 @@ private: // init_type: "ideal", "real", "input_sounding", "metgrid" or "" static std::string init_type; + // sponge_type: "input_sponge" + static std::string sponge_type; + // use_real_bcs: only true if 1) ( (init_type == real) or (init_type == metgrid) ) // AND 2) we want to use the bc's from the WRF bdy file static bool use_real_bcs; @@ -836,6 +864,9 @@ private: // Text input_sounding file static std::string input_sounding_file; + // Text input_sponge file + static std::string input_sponge_file; + // Flag to trigger initialization from input_sounding like WRF's ideal.exe // used with init_type == "input_sounding" static bool init_sounding_ideal; @@ -879,9 +910,15 @@ private: // This is a vector over levels of vectors across quantities of Vectors amrex::Vector > > h_rayleigh_ptrs; + // This is a vector over levels of vectors across quantities of Vectors + amrex::Vector > > h_sponge_ptrs; + // This is a vector over levels of vectors across quantities of DeviceVectors amrex::Vector > > d_rayleigh_ptrs; + // This is a vector over levels of vectors across quantities of DeviceVectors + amrex::Vector > > d_sponge_ptrs; + amrex::Vector h_havg_density; amrex::Vector h_havg_temperature; amrex::Vector h_havg_pressure; diff --git a/Source/ERF.cpp b/Source/ERF.cpp index 51ab507bc..aec7f0ac3 100644 --- a/Source/ERF.cpp +++ b/Source/ERF.cpp @@ -62,6 +62,9 @@ std::string ERF::nc_bdy_file; // Must provide via input // Text input_sounding file std::string ERF::input_sounding_file = "input_sounding"; +// Text input_sponge file +std::string ERF::input_sponge_file = "input_sponge_file.txt"; + // Flag to trigger initialization from input_sounding like WRF's ideal.exe bool ERF::init_sounding_ideal = false; @@ -112,10 +115,9 @@ ERF::ERF () int nlevs_max = max_level + 1; #ifdef ERF_USE_WINDFARM - if(solverChoice.windfarm_type == WindFarmType::Fitch){ - Nturb.resize(nlevs_max); - vars_fitch.resize(nlevs_max); - } + Nturb.resize(nlevs_max); + vars_fitch.resize(nlevs_max); + vars_ewp.resize(nlevs_max); #endif #if defined(ERF_USE_RRTMGP) @@ -263,6 +265,18 @@ ERF::ERF () // Qv prim for MOST Qv_prim.resize(nlevs_max); + // Time averaged velocity field + vel_t_avg.resize(nlevs_max); + t_avg_cnt.resize(nlevs_max); + +#ifdef ERF_USE_NETCDF + // Longitude and latitude (only filled for use_real_bcs==True) + if (use_real_bcs) { + lat_m.resize(nlevs_max); + lon_m.resize(nlevs_max); + } +#endif + // Initialize tagging criteria for mesh refinement refinement_criteria_setup(); @@ -321,7 +335,7 @@ ERF::Evolve () cur_time += dt[0]; Print() << "Coarse STEP " << step+1 << " ends." << " TIME = " << cur_time - << " DT = " << dt[0] << std::endl; + << " DT = " << dt[0] << std::endl; post_timestep(step, cur_time, dt[0]); @@ -761,7 +775,14 @@ ERF::InitData () bool restarting = (!restart_chkfile.empty()); setRayleighRefFromSounding(restarting); } + } + // Read in sponge data from input file + if(solverChoice.spongeChoice.sponge_type == "input_sponge") + { + initSponge(); + bool restarting = (!restart_chkfile.empty()); + setSpongeRefFromSounding(restarting); } if (is_it_time_for_action(istep[0], t_new[0], dt[0], sum_interval, sum_per)) { @@ -902,6 +923,16 @@ ERF::InitData () for (int lev = 0; lev <= finest_level; ++lev) micro->Update_Micro_Vars_Lev(lev, vars_new[lev][Vars::cons]); } + // Fill time averaged velocities before first plot file + if (solverChoice.time_avg_vel) { + for (int lev = 0; lev <= finest_level; ++lev) { + Time_Avg_Vel_atCC(dt[lev], t_avg_cnt[lev], vel_t_avg[lev].get(), + vars_new[lev][Vars::xvel], + vars_new[lev][Vars::yvel], + vars_new[lev][Vars::zvel]); + } + } + // check for additional plotting variables that are available after particle containers // are setup. const std::string& pv1 = "plot_vars_1"; appendPlotVariables(pv1,plot_var_names_1); @@ -1149,6 +1180,11 @@ ERF::init_only (int lev, Real time) #ifdef ERF_USE_WINDFARM init_windfarm(lev); #endif + + if(solverChoice.spongeChoice.sponge_type == "input_sponge"){ + input_sponge(lev); + } + } // read in some parameters from inputs file @@ -1295,6 +1331,9 @@ ERF::ReadParameters () // Text input_sounding file pp.query("input_sounding_file", input_sounding_file); + // Text input_sounding file + pp.query("input_sponge_file", input_sponge_file); + // Flag to trigger initialization from input_sounding like WRF's ideal.exe pp.query("init_sounding_ideal", init_sounding_ideal); @@ -1741,10 +1780,9 @@ ERF::ERF (const RealBox& rb, int max_level_in, int nlevs_max = max_level + 1; #ifdef ERF_USE_WINDFARM - if(solverChoice.windfarm_type == WindFarmType::Fitch){ - Nturb.resize(nlevs_max); - vars_fitch.resize(nlevs_max); - } + Nturb.resize(nlevs_max); + vars_fitch.resize(nlevs_max); + vars_ewp.resize(nlevs_max); #endif #if defined(ERF_USE_RRTMGP) @@ -1847,6 +1885,10 @@ ERF::ERF (const RealBox& rb, int max_level_in, // Theta prim for MOST Theta_prim.resize(nlevs_max); + // Time averaged velocity field + vel_t_avg.resize(nlevs_max); + t_avg_cnt.resize(nlevs_max); + // Initialize tagging criteria for mesh refinement refinement_criteria_setup(); diff --git a/Source/ERF_make_new_arrays.cpp b/Source/ERF_make_new_arrays.cpp index 316217410..3051e3d48 100644 --- a/Source/ERF_make_new_arrays.cpp +++ b/Source/ERF_make_new_arrays.cpp @@ -158,6 +158,21 @@ ERF::init_stuff (int lev, const BoxArray& ba, const DistributionMapping& dm, rV_new[lev].setVal(3.4e22); rW_new[lev].setVal(5.6e23); + // ******************************************************************************************** + // These are just time averaged fields for diagnostics + // ******************************************************************************************** + + // NOTE: We are not completing a fillpach call on the time averaged data; + // which would copy on intersection and interpolate from coarse. + // Therefore, we are restarting the averaging when the ba changes, + // this may give poor statistics for dynamic mesh refinment. + vel_t_avg[lev] = nullptr; + if (solverChoice.time_avg_vel) { + vel_t_avg[lev] = std::make_unique(ba, dm, 4, 0); // Each vel comp and the mag + vel_t_avg[lev]->setVal(0.0); + t_avg_cnt[lev] = 0.0; + } + // ******************************************************************************************** // Initialize flux registers whenever we create/re-create a level // ******************************************************************************************** @@ -169,7 +184,7 @@ ERF::init_stuff (int lev, const BoxArray& ba, const DistributionMapping& dm, advflux_reg[lev] = new YAFluxRegister(ba , grids[lev-1], dm , dmap[lev-1], geom[lev], geom[lev-1], - ref_ratio[lev-1], lev, ncomp_reflux); + ref_ratio[lev-1], lev, ncomp_reflux); } } @@ -177,15 +192,15 @@ ERF::init_stuff (int lev, const BoxArray& ba, const DistributionMapping& dm, // Define Theta_prim storage if using MOST BC // ******************************************************************************************** if (phys_bc_type[Orientation(Direction::z,Orientation::low)] == ERF_BC::MOST) { - Theta_prim[lev] = std::make_unique(ba,dm,1,IntVect(ngrow_state,ngrow_state,0)); - if (solverChoice.moisture_type != MoistureType::None) { - Qv_prim[lev] = std::make_unique(ba,dm,1,IntVect(ngrow_state,ngrow_state,0)); - } else { - Qv_prim[lev] = nullptr; - } + Theta_prim[lev] = std::make_unique(ba,dm,1,IntVect(ngrow_state,ngrow_state,0)); + if (solverChoice.moisture_type != MoistureType::None) { + Qv_prim[lev] = std::make_unique(ba,dm,1,IntVect(ngrow_state,ngrow_state,0)); + } else { + Qv_prim[lev] = nullptr; + } } else { - Theta_prim[lev] = nullptr; - Qv_prim[lev] = nullptr; + Theta_prim[lev] = nullptr; + Qv_prim[lev] = nullptr; } // ******************************************************************************************** @@ -220,14 +235,19 @@ ERF::init_stuff (int lev, const BoxArray& ba, const DistributionMapping& dm, vars_fitch[lev].define(ba, dm, 5, ngrow_state); // V, dVabsdt, dudt, dvdt, dTKEdt Nturb[lev].define(ba, dm, 1, ngrow_state); // Number of turbines in a cell } + if (solverChoice.windfarm_type == WindFarmType::EWP){ + int ngrow_state = ComputeGhostCells(solverChoice.advChoice, solverChoice.use_NumDiff) + 1; + vars_ewp[lev].define(ba, dm, 3, ngrow_state); // dudt, dvdt, dTKEdt + Nturb[lev].define(ba, dm, 1, ngrow_state); // Number of turbines in a cell + } #endif #if defined(ERF_USE_RRTMGP) //********************************************************* // Radiation heating source terms //********************************************************* - qheating_rates[lev].define(ba, dm, 2, ngrow_state); - qheating_rates[lev].setVal(0.); + qheating_rates[lev] = std::make_unique(ba, dm, 2, ngrow_state); + qheating_rates[lev]->setVal(0.); #endif } diff --git a/Source/ERF_make_new_level.cpp b/Source/ERF_make_new_level.cpp index 23be8b4fa..43f8c7dff 100644 --- a/Source/ERF_make_new_level.cpp +++ b/Source/ERF_make_new_level.cpp @@ -381,7 +381,7 @@ ERF::RemakeLevel (int lev, Real time, const BoxArray& ba, const DistributionMapp //******************************************************************************************** // This allocates all kinds of things, including but not limited to: solution arrays, - // terrain arrays and metrics,a nd base state. + // terrain arrays and metrics, and base state. // ******************************************************************************************* init_stuff(lev, ba, dm, temp_lev_new, temp_lev_old); diff --git a/Source/IO/Checkpoint.cpp b/Source/IO/Checkpoint.cpp index 9bc2d8a12..ee4d22564 100644 --- a/Source/IO/Checkpoint.cpp +++ b/Source/IO/Checkpoint.cpp @@ -172,7 +172,7 @@ ERF::WriteCheckpointFile () const #if defined(ERF_USE_WINDFARM) - if(solverChoice.windfarm_type == WindFarmType::Fitch){ + if(solverChoice.windfarm_type == WindFarmType::Fitch or solverChoice.windfarm_type == WindFarmType::EWP){ ng = Nturb[lev].nGrowVect(); MultiFab mf_Nturb(grids[lev],dmap[lev],1,ng); MultiFab::Copy(mf_Nturb,Nturb[lev],0,0,1,ng); @@ -424,7 +424,7 @@ ERF::ReadCheckpointFile () } #if defined(ERF_USE_WINDFARM) - if(solverChoice.windfarm_type == WindFarmType::Fitch){ + if(solverChoice.windfarm_type == WindFarmType::Fitch or solverChoice.windfarm_type == WindFarmType::EWP){ ng = Nturb[lev].nGrowVect(); MultiFab mf_Nturb(grids[lev],dmap[lev],1,ng); VisMF::Read(mf_Nturb, amrex::MultiFabFileFullPrefix(lev, restart_chkfile, "Level_", "NumTurb")); diff --git a/Source/IO/ERF_WriteScalarProfiles.cpp b/Source/IO/ERF_WriteScalarProfiles.cpp index a5885161e..78a61f9ec 100644 --- a/Source/IO/ERF_WriteScalarProfiles.cpp +++ b/Source/IO/ERF_WriteScalarProfiles.cpp @@ -30,9 +30,9 @@ ERF::sum_integrated_quantities (Real time) Real scal_ml = 0.0; #if 1 - mass_sl = volWgtSumMF(0,vars_new[0][Vars::cons],Rho_comp,*mapfac_m[0],false,false); + mass_sl = volWgtSumMF(0,vars_new[0][Vars::cons],Rho_comp,*mapfac_m[0],true,false); for (int lev = 0; lev <= finest_level; lev++) { - mass_ml += volWgtSumMF(lev,vars_new[lev][Vars::cons],Rho_comp,*mapfac_m[lev],false,true); + mass_ml += volWgtSumMF(lev,vars_new[lev][Vars::cons],Rho_comp,*mapfac_m[lev],true,true); } #else for (int lev = 0; lev <= finest_level; lev++) { @@ -51,110 +51,107 @@ ERF::sum_integrated_quantities (Real time) }); } if (lev == 0) { - mass_sl = volWgtSumMF(0,pert_dens,0,*mapfac_m[0],false,false); + mass_sl = volWgtSumMF(0,pert_dens,0,*mapfac_m[0],true,false); } - mass_ml += volWgtSumMF(lev,pert_dens,0,*mapfac_m[lev],false,true); + mass_ml += volWgtSumMF(lev,pert_dens,0,*mapfac_m[lev],true,true); } // lev #endif - Real rhth_sl = volWgtSumMF(0,vars_new[0][Vars::cons], RhoTheta_comp,*mapfac_m[0],false,false); - Real scal_sl = volWgtSumMF(0,vars_new[0][Vars::cons],RhoScalar_comp,*mapfac_m[0],false,false); + Real rhth_sl = volWgtSumMF(0,vars_new[0][Vars::cons], RhoTheta_comp,*mapfac_m[0],true,false); + Real scal_sl = volWgtSumMF(0,vars_new[0][Vars::cons],RhoScalar_comp,*mapfac_m[0],true,false); for (int lev = 0; lev <= finest_level; lev++) { - rhth_ml += volWgtSumMF(lev,vars_new[lev][Vars::cons], RhoTheta_comp,*mapfac_m[lev],false,true); - scal_ml += volWgtSumMF(lev,vars_new[lev][Vars::cons],RhoScalar_comp,*mapfac_m[lev],false,true); + rhth_ml += volWgtSumMF(lev,vars_new[lev][Vars::cons], RhoTheta_comp,*mapfac_m[lev],true,true); + scal_ml += volWgtSumMF(lev,vars_new[lev][Vars::cons],RhoScalar_comp,*mapfac_m[lev],true,true); } - if (verbose > 0) { - - Gpu::HostVector h_avg_ustar; h_avg_ustar.resize(1); - Gpu::HostVector h_avg_tstar; h_avg_tstar.resize(1); - Gpu::HostVector h_avg_olen; h_avg_olen.resize(1); - if ((m_most != nullptr) && (NumDataLogs() > 0)) { - Box domain = geom[0].Domain(); - int zdir = 2; - h_avg_ustar = sumToLine(*m_most->get_u_star(0),0,1,domain,zdir); - h_avg_tstar = sumToLine(*m_most->get_t_star(0),0,1,domain,zdir); - h_avg_olen = sumToLine(*m_most->get_olen(0),0,1,domain,zdir); - - // Divide by the total number of cells we are averaging over - Real area_z = static_cast(domain.length(0)*domain.length(1)); - h_avg_ustar[0] /= area_z; - h_avg_tstar[0] /= area_z; - h_avg_olen[0] /= area_z; - - } else { - h_avg_ustar[0] = 0.; - h_avg_tstar[0] = 0.; - h_avg_olen[0] = 0.; - } + Gpu::HostVector h_avg_ustar; h_avg_ustar.resize(1); + Gpu::HostVector h_avg_tstar; h_avg_tstar.resize(1); + Gpu::HostVector h_avg_olen; h_avg_olen.resize(1); + if ((m_most != nullptr) && (NumDataLogs() > 0)) { + Box domain = geom[0].Domain(); + int zdir = 2; + h_avg_ustar = sumToLine(*m_most->get_u_star(0),0,1,domain,zdir); + h_avg_tstar = sumToLine(*m_most->get_t_star(0),0,1,domain,zdir); + h_avg_olen = sumToLine(*m_most->get_olen(0),0,1,domain,zdir); + + // Divide by the total number of cells we are averaging over + Real area_z = static_cast(domain.length(0)*domain.length(1)); + h_avg_ustar[0] /= area_z; + h_avg_tstar[0] /= area_z; + h_avg_olen[0] /= area_z; + + } else { + h_avg_ustar[0] = 0.; + h_avg_tstar[0] = 0.; + h_avg_olen[0] = 0.; + } - const int nfoo = 6; - Real foo[nfoo] = {mass_sl,rhth_sl,scal_sl,mass_ml,rhth_ml,scal_ml}; + const int nfoo = 6; + Real foo[nfoo] = {mass_sl,rhth_sl,scal_sl,mass_ml,rhth_ml,scal_ml}; #ifdef AMREX_LAZY - Lazy::QueueReduction([=]() mutable { + Lazy::QueueReduction([=]() mutable { #endif - ParallelDescriptor::ReduceRealSum( - foo, nfoo, ParallelDescriptor::IOProcessorNumber()); - - if (ParallelDescriptor::IOProcessor()) { - int i = 0; - mass_sl = foo[i++]; - rhth_sl = foo[i++]; - scal_sl = foo[i++]; - mass_ml = foo[i++]; - rhth_ml = foo[i++]; - scal_ml = foo[i++]; - - Print() << '\n'; - if (finest_level == 0) { + ParallelDescriptor::ReduceRealSum( + foo, nfoo, ParallelDescriptor::IOProcessorNumber()); + + if (ParallelDescriptor::IOProcessor()) { + int i = 0; + mass_sl = foo[i++]; + rhth_sl = foo[i++]; + scal_sl = foo[i++]; + mass_ml = foo[i++]; + rhth_ml = foo[i++]; + scal_ml = foo[i++]; + + Print() << '\n'; + if (finest_level == 0) { #if 1 - Print() << "TIME= " << time << " MASS = " << mass_sl << '\n'; + Print() << "TIME= " << time << " MASS = " << mass_sl << '\n'; #else - Print() << "TIME= " << time << " PERT MASS = " << mass_sl << '\n'; + Print() << "TIME= " << time << " PERT MASS = " << mass_sl << '\n'; #endif - Print() << "TIME= " << time << " RHO THETA = " << rhth_sl << '\n'; - Print() << "TIME= " << time << " RHO SCALAR = " << scal_sl << '\n'; - } else { + Print() << "TIME= " << time << " RHO THETA = " << rhth_sl << '\n'; + Print() << "TIME= " << time << " RHO SCALAR = " << scal_sl << '\n'; + } else { #if 1 - Print() << "TIME= " << time << " MASS SL/ML = " << mass_sl << " " << mass_ml << '\n'; + Print() << "TIME= " << time << " MASS SL/ML = " << mass_sl << " " << mass_ml << '\n'; #else - Print() << "TIME= " << time << " PERT MASS SL/ML = " << mass_sl << " " << mass_ml << '\n'; + Print() << "TIME= " << time << " PERT MASS SL/ML = " << mass_sl << " " << mass_ml << '\n'; #endif - Print() << "TIME= " << time << " RHO THETA SL/ML = " << rhth_sl << " " << rhth_ml << '\n'; - Print() << "TIME= " << time << " RHO SCALAR SL/ML = " << scal_sl << " " << scal_ml << '\n'; - } - - // The first data log only holds scalars - if (NumDataLogs() > 0) - { - int nd = 0; - std::ostream& data_log1 = DataLog(nd); - if (data_log1.good()) { - if (time == 0.0) { - data_log1 << std::setw(datwidth) << " time"; - data_log1 << std::setw(datwidth) << " u_star"; - data_log1 << std::setw(datwidth) << " t_star"; - data_log1 << std::setw(datwidth) << " olen"; - data_log1 << std::endl; - } // time = 0 - - // Write the quantities at this time - data_log1 << std::setw(datwidth) << time; - data_log1 << std::setw(datwidth) << std::setprecision(datprecision) - << h_avg_ustar[0]; - data_log1 << std::setw(datwidth) << std::setprecision(datprecision) - << h_avg_tstar[0]; - data_log1 << std::setw(datwidth) << std::setprecision(datprecision) - << h_avg_olen[0]; - data_log1 << std::endl; - } // if good - } // loop over i - } // if IOProcessor + Print() << "TIME= " << time << " RHO THETA SL/ML = " << rhth_sl << " " << rhth_ml << '\n'; + Print() << "TIME= " << time << " RHO SCALAR SL/ML = " << scal_sl << " " << scal_ml << '\n'; + } + + // The first data log only holds scalars + if (NumDataLogs() > 0) + { + int nd = 0; + std::ostream& data_log1 = DataLog(nd); + if (data_log1.good()) { + if (time == 0.0) { + data_log1 << std::setw(datwidth) << " time"; + data_log1 << std::setw(datwidth) << " u_star"; + data_log1 << std::setw(datwidth) << " t_star"; + data_log1 << std::setw(datwidth) << " olen"; + data_log1 << std::endl; + } // time = 0 + + // Write the quantities at this time + data_log1 << std::setw(datwidth) << time; + data_log1 << std::setw(datwidth) << std::setprecision(datprecision) + << h_avg_ustar[0]; + data_log1 << std::setw(datwidth) << std::setprecision(datprecision) + << h_avg_tstar[0]; + data_log1 << std::setw(datwidth) << std::setprecision(datprecision) + << h_avg_olen[0]; + data_log1 << std::endl; + } // if good + } // loop over i + } // if IOProcessor #ifdef AMREX_LAZY - }); + }); #endif - } // if verbose // This is just an alias for convenience int lev = 0; diff --git a/Source/IO/Plotfile.cpp b/Source/IO/Plotfile.cpp index 253d353ce..e486d5314 100644 --- a/Source/IO/Plotfile.cpp +++ b/Source/IO/Plotfile.cpp @@ -87,6 +87,9 @@ ERF::setPlotVariables (const std::string& pp_plot_var_names, Vector #endif plot_var_names = tmp_plot_names; + + for (int i=0; i plot_var_names) mf_comp ++; } +#ifdef ERF_USE_NETCDF + if (use_real_bcs) { + if (containerHasElement(plot_var_names, "lat_m")) { +#ifdef _OPENMP +#pragma omp parallel if (amrex::Gpu::notInLaunchRegion()) +#endif + for ( MFIter mfi(mf[lev],TilingIfNotGPU()); mfi.isValid(); ++mfi) + { + const Box& bx = mfi.tilebox(); + const Array4& derdat = mf[lev].array(mfi); + const Array4& data = lat_m[lev]->array(mfi); + ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + derdat(i, j, k, mf_comp) = data(i,j,0); + }); + } + mf_comp ++; + } // lat_m + if (containerHasElement(plot_var_names, "lon_m")) { +#ifdef _OPENMP +#pragma omp parallel if (amrex::Gpu::notInLaunchRegion()) +#endif + for ( MFIter mfi(mf[lev],TilingIfNotGPU()); mfi.isValid(); ++mfi) + { + const Box& bx = mfi.tilebox(); + const Array4& derdat = mf[lev].array(mfi); + const Array4& data = lon_m[lev]->array(mfi); + ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + derdat(i, j, k, mf_comp) = data(i,j,0); + }); + } + mf_comp ++; + } // lon_m + } // use_real_bcs +#endif + + + if (solverChoice.time_avg_vel) { + if (containerHasElement(plot_var_names, "u_t_avg")) { +#ifdef _OPENMP +#pragma omp parallel if (amrex::Gpu::notInLaunchRegion()) +#endif + for ( MFIter mfi(mf[lev],TilingIfNotGPU()); mfi.isValid(); ++mfi) + { + const Box& bx = mfi.tilebox(); + const Array4& derdat = mf[lev].array(mfi); + const Array4& data = vel_t_avg[lev]->array(mfi); + const Real norm = t_avg_cnt[lev]; + ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept + { + derdat(i ,j ,k, mf_comp) = data(i,j,k,0) / norm; + }); + } + mf_comp ++; + } + + if (containerHasElement(plot_var_names, "v_t_avg")) { +#ifdef _OPENMP +#pragma omp parallel if (amrex::Gpu::notInLaunchRegion()) +#endif + for ( MFIter mfi(mf[lev],TilingIfNotGPU()); mfi.isValid(); ++mfi) + { + const Box& bx = mfi.tilebox(); + const Array4& derdat = mf[lev].array(mfi); + const Array4& data = vel_t_avg[lev]->array(mfi); + const Real norm = t_avg_cnt[lev]; + ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept + { + derdat(i ,j ,k, mf_comp) = data(i,j,k,1) / norm; + }); + } + mf_comp ++; + } + + if (containerHasElement(plot_var_names, "w_t_avg")) { +#ifdef _OPENMP +#pragma omp parallel if (amrex::Gpu::notInLaunchRegion()) +#endif + for ( MFIter mfi(mf[lev],TilingIfNotGPU()); mfi.isValid(); ++mfi) + { + const Box& bx = mfi.tilebox(); + const Array4& derdat = mf[lev].array(mfi); + const Array4& data = vel_t_avg[lev]->array(mfi); + const Real norm = t_avg_cnt[lev]; + ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept + { + derdat(i ,j ,k, mf_comp) = data(i,j,k,2) / norm; + }); + } + mf_comp ++; + } + + if (containerHasElement(plot_var_names, "umag_t_avg")) { +#ifdef _OPENMP +#pragma omp parallel if (amrex::Gpu::notInLaunchRegion()) +#endif + for ( MFIter mfi(mf[lev],TilingIfNotGPU()); mfi.isValid(); ++mfi) + { + const Box& bx = mfi.tilebox(); + const Array4& derdat = mf[lev].array(mfi); + const Array4& data = vel_t_avg[lev]->array(mfi); + const Real norm = t_avg_cnt[lev]; + ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept + { + derdat(i ,j ,k, mf_comp) = data(i,j,k,3) / norm; + }); + } + mf_comp ++; + } + } + if (containerHasElement(plot_var_names, "Kmv")) { MultiFab::Copy(mf[lev],*eddyDiffs_lev[lev],EddyDiff::Mom_v,mf_comp,1,0); mf_comp ++; diff --git a/Source/Initialization/ERF_init1d.cpp b/Source/Initialization/ERF_init1d.cpp index 59bdf973b..c41ab360a 100644 --- a/Source/Initialization/ERF_init1d.cpp +++ b/Source/Initialization/ERF_init1d.cpp @@ -114,6 +114,88 @@ ERF::setRayleighRefFromSounding (bool restarting) } } +/** + * Initialization function for host and device vectors + * used to store the effects of sponge Damping. + */ +void +ERF::initSponge () +{ + h_sponge_ptrs.resize(max_level+1); + d_sponge_ptrs.resize(max_level+1); + + for (int lev = 0; lev <= finest_level; lev++) + { + // These have 2 components: ubar, vbar + h_sponge_ptrs[lev].resize(Sponge::nvars_sponge); + d_sponge_ptrs[lev].resize(Sponge::nvars_sponge); + + const int zlen_sponge = geom[lev].Domain().length(2); + + // Allocate space for these 1D vectors + for (int n = 0; n < Sponge::nvars_sponge; n++) { + h_sponge_ptrs[lev][n].resize(zlen_sponge, 0.0_rt); + d_sponge_ptrs[lev][n].resize(zlen_sponge, 0.0_rt); + } + + } +} + +/** + * Sets the sponge damping averaged quantities from an + * externally supplied input sponge data file. + * + * @param[in] restarting Boolean parameter that indicates whether + we are currently restarting from a checkpoint file. + */ +void +ERF::setSpongeRefFromSounding (bool restarting) +{ + + // If we are restarting then we haven't read the input_sponge file yet + // so we need to read it here + // TODO: should we store this information in the checkpoint file instead? + if (restarting) { + input_sponge_data.read_from_file(input_sponge_file, geom[0], zlevels_stag); + } + + const Real* z_inp_sponge = input_sponge_data.z_inp_sponge.dataPtr(); + const Real* U_inp_sponge = input_sponge_data.U_inp_sponge.dataPtr(); + const Real* V_inp_sponge = input_sponge_data.V_inp_sponge.dataPtr(); + const int inp_sponge_size = input_sponge_data.size(); + + for (int lev = 0; lev <= finest_level; lev++) + { + const int khi = geom[lev].Domain().bigEnd()[2]; + Vector zcc(khi+1); + + if (z_phys_cc[lev]) { + // use_terrain=1 + // calculate the damping strength based on the max height at each k + reduce_to_max_per_level(zcc, z_phys_cc[lev]); + } else { + const auto *const prob_lo = geom[lev].ProbLo(); + const auto *const dx = geom[lev].CellSize(); + for (int k = 0; k <= khi; k++) + { + zcc[k] = prob_lo[2] + (k+0.5) * dx[2]; + } + } + + for (int k = 0; k <= khi; k++) + { + h_sponge_ptrs[lev][Sponge::ubar_sponge][k] = interpolate_1d(z_inp_sponge, U_inp_sponge, zcc[k], inp_sponge_size); + h_sponge_ptrs[lev][Sponge::vbar_sponge][k] = interpolate_1d(z_inp_sponge, V_inp_sponge, zcc[k], inp_sponge_size); + } + + // Copy from host version to device version + Gpu::copy(Gpu::hostToDevice, h_sponge_ptrs[lev][Sponge::ubar_sponge].begin(), h_sponge_ptrs[lev][Sponge::ubar_sponge].end(), + d_sponge_ptrs[lev][Sponge::ubar_sponge].begin()); + Gpu::copy(Gpu::hostToDevice, h_sponge_ptrs[lev][Sponge::vbar_sponge].begin(), h_sponge_ptrs[lev][Sponge::vbar_sponge].end(), + d_sponge_ptrs[lev][Sponge::vbar_sponge].begin()); + } +} + /** * Initialize density and pressure base state in * hydrostatic equilibrium. diff --git a/Source/Initialization/ERF_init_from_metgrid.cpp b/Source/Initialization/ERF_init_from_metgrid.cpp index 2a0a6a505..71dcc1748 100644 --- a/Source/Initialization/ERF_init_from_metgrid.cpp +++ b/Source/Initialization/ERF_init_from_metgrid.cpp @@ -140,7 +140,7 @@ ERF::init_from_metgrid (int lev) // This defines all the z(i,j,k) values given z(i,j,0) from above. init_terrain_grid(lev, geom[lev], *z_phys, zlevels_stag); - // Copy SST and LANDMASK data into MF and iMF data structures + // Copy LATITUDE, LONGITUDE, SST and LANDMASK data into MF and iMF data structures auto& ba = lev_new[Vars::cons].boxArray(); auto& dm = lev_new[Vars::cons].DistributionMap(); auto ngv = lev_new[Vars::cons].nGrowVect(); ngv[2] = 0; @@ -193,6 +193,34 @@ ERF::init_from_metgrid (int lev) } else { for (int it = 0; it < ntimes; ++it) lmask_lev[lev][it] = nullptr; } + lat_m[lev] = std::make_unique(ba2d,dm,1,ngv); + for ( MFIter mfi(*(lat_m[lev]), TilingIfNotGPU()); mfi.isValid(); ++mfi ) { + Box gtbx = mfi.growntilebox(); + FArrayBox& dst = (*(lat_m[lev]))[mfi]; + FArrayBox& src = NC_LAT_fab[0]; + const Array4< Real>& dst_arr = dst.array(); + const Array4& src_arr = src.const_array(); + ParallelFor(gtbx, [=] AMREX_GPU_DEVICE (int i, int j, int) noexcept + { + int li = amrex::min(amrex::max(i, i_lo), i_hi); + int lj = amrex::min(amrex::max(j, j_lo), j_hi); + dst_arr(i,j,0) = src_arr(li,lj,0); + }); + } + lon_m[lev] = std::make_unique(ba2d,dm,1,ngv); + for ( MFIter mfi(*(lon_m[lev]), TilingIfNotGPU()); mfi.isValid(); ++mfi ) { + Box gtbx = mfi.growntilebox(); + FArrayBox& dst = (*(lon_m[lev]))[mfi]; + FArrayBox& src = NC_LON_fab[0]; + const Array4< Real>& dst_arr = dst.array(); + const Array4& src_arr = src.const_array(); + ParallelFor(gtbx, [=] AMREX_GPU_DEVICE (int i, int j, int) noexcept + { + int li = amrex::min(amrex::max(i, i_lo), i_hi); + int lj = amrex::min(amrex::max(j, j_lo), j_hi); + dst_arr(i,j,0) = src_arr(li,lj,0); + }); + } for (int it = 0; it < ntimes; it++) { // Verify that the grid size and resolution from met_em file matches that in geom (from ERF inputs file). diff --git a/Source/Initialization/ERF_init_from_wrfinput.cpp b/Source/Initialization/ERF_init_from_wrfinput.cpp index 9e59f2adf..fb35b79b1 100644 --- a/Source/Initialization/ERF_init_from_wrfinput.cpp +++ b/Source/Initialization/ERF_init_from_wrfinput.cpp @@ -185,6 +185,45 @@ ERF::init_from_wrfinput (int lev) solverChoice.moisture_type, n_qstate); } // mf + auto& ba = lev_new[Vars::cons].boxArray(); + auto& dm = lev_new[Vars::cons].DistributionMap(); + auto ngv = lev_new[Vars::cons].nGrowVect(); ngv[2] = 0; + BoxList bl2d = ba.boxList(); + for (auto& b : bl2d) { + b.setRange(2,0); + } + BoxArray ba2d(std::move(bl2d)); + int i_lo = geom[lev].Domain().smallEnd(0); int i_hi = geom[lev].Domain().bigEnd(0); + int j_lo = geom[lev].Domain().smallEnd(1); int j_hi = geom[lev].Domain().bigEnd(1); + lat_m[lev] = std::make_unique(ba2d,dm,1,ngv); + for ( MFIter mfi(*(lat_m[lev]), TilingIfNotGPU()); mfi.isValid(); ++mfi ) { + Box gtbx = mfi.growntilebox(); + FArrayBox& dst = (*(lat_m[lev]))[mfi]; + FArrayBox& src = NC_LAT_fab[0]; + const Array4< Real>& dst_arr = dst.array(); + const Array4& src_arr = src.const_array(); + ParallelFor(gtbx, [=] AMREX_GPU_DEVICE (int i, int j, int) noexcept + { + int li = amrex::min(amrex::max(i, i_lo), i_hi); + int lj = amrex::min(amrex::max(j, j_lo), j_hi); + dst_arr(i,j,0) = src_arr(li,lj,0); + }); + } + lon_m[lev] = std::make_unique(ba2d,dm,1,ngv); + for ( MFIter mfi(*(lon_m[lev]), TilingIfNotGPU()); mfi.isValid(); ++mfi ) { + Box gtbx = mfi.growntilebox(); + FArrayBox& dst = (*(lon_m[lev]))[mfi]; + FArrayBox& src = NC_LON_fab[0]; + const Array4< Real>& dst_arr = dst.array(); + const Array4& src_arr = src.const_array(); + ParallelFor(gtbx, [=] AMREX_GPU_DEVICE (int i, int j, int) noexcept + { + int li = amrex::min(amrex::max(i, i_lo), i_hi); + int lj = amrex::min(amrex::max(j, j_lo), j_hi); + dst_arr(i,j,0) = src_arr(li,lj,0); + }); + } + #ifdef _OPENMP #pragma omp parallel if (amrex::Gpu::notInLaunchRegion()) #endif diff --git a/Source/Initialization/ERF_input_sponge.cpp b/Source/Initialization/ERF_input_sponge.cpp new file mode 100644 index 000000000..aee3c8bab --- /dev/null +++ b/Source/Initialization/ERF_input_sponge.cpp @@ -0,0 +1,32 @@ +/** + * \file ERF_input_sponge.cpp + */ + +#include +#include +#include +#include +#include + +using namespace amrex; + +/** + * High level wrapper for sponge x and y velocities + * level data from input sponge data. + * + * @param lev Integer specifying the current level + */ +void +ERF::input_sponge (int lev) +{ + // We only want to read the file once + if (lev == 0) { + if (input_sponge_file.empty()) + Error("input_sounding file name must be provided via input"); + + // this will interpolate the input profiles to the nominal height levels + // (ranging from 0 to the domain top) + input_sponge_data.read_from_file(input_sponge_file, geom[lev], zlevels_stag); + + } +} diff --git a/Source/Initialization/InputSpongeData.H b/Source/Initialization/InputSpongeData.H new file mode 100644 index 000000000..73fbdd44f --- /dev/null +++ b/Source/Initialization/InputSpongeData.H @@ -0,0 +1,109 @@ +#ifndef _INPUT_SPONGE_DATA_H_ +#define _INPUT_SPONGE_DATA_H_ + +#include +#include + +#include +#include +#include +#include + +#include +#include + +/** + * Data structure storing input sponge data. Also + * handles reading the input file for sponge data + */ +struct InputSpongeData { +public: + InputSpongeData () {} + + void read_from_file (const std::string input_sponge_file, + const amrex::Geometry &geom, + const amrex::Vector& zlevels_stag) + { + const int klo = 0; + const int khi = geom.Domain().bigEnd()[AMREX_SPACEDIM-1]; + const int Nz = geom.Domain().size()[AMREX_SPACEDIM-1]; + const amrex::Real dz = geom.CellSize()[AMREX_SPACEDIM-1]; + + const bool use_terrain = (zlevels_stag.size() > 0); + const amrex::Real zbot = (use_terrain) ? zlevels_stag[klo] : geom.ProbLo(AMREX_SPACEDIM-1); + const amrex::Real ztop = (use_terrain) ? zlevels_stag[khi+1] : geom.ProbHi(AMREX_SPACEDIM-1); + + z_inp_sponge.resize(Nz+2); + U_inp_sponge.resize(Nz+2); + V_inp_sponge.resize(Nz+2); + + // Read the input_sponge file + amrex::Print() << "input_sponge file location : " << input_sponge_file << std::endl; + std::ifstream input_sponge_reader(input_sponge_file); + if(!input_sponge_reader.is_open()) { + amrex::Error("Error opening the input_sponge file.\n"); + } + else { + // Read the contents of the input_sponge file + amrex::Print() << "Successfully opened the input_sponge file. Now reading... " << std::endl; + std::string line; + + // First, read the input data into temp vectors; then, interpolate vectors to the + // domain lo/hi and cell centers (from level 0) + amrex::Vector z_inp_sponge_tmp, U_inp_sponge_tmp, V_inp_sponge_tmp; + + // Add surface + z_inp_sponge_tmp.push_back(zbot); // height above sea level [m] + U_inp_sponge_tmp.push_back(0); + V_inp_sponge_tmp.push_back(0); + + // Read the vertical profile at each given height + amrex::Real z, U, V; + while(std::getline(input_sponge_reader, line)) { + std::istringstream iss_z(line); + iss_z >> z >> U >> V; + if (z == zbot) { + U_inp_sponge_tmp[0] = U; + V_inp_sponge_tmp[0] = V; + } else { + AMREX_ALWAYS_ASSERT(z > z_inp_sponge_tmp[z_inp_sponge_tmp.size()-1]); // sounding is increasing in height + z_inp_sponge_tmp.push_back(z); + U_inp_sponge_tmp.push_back(U); + V_inp_sponge_tmp.push_back(V); + if (z >= ztop) break; + } + } + + // At this point, we have an input_sponge from zbot up to + // z_inp_sponge_tmp[N-1] >= ztop. Now, interpolate to grid level 0 heights + const int Ninp = z_inp_sponge_tmp.size(); + z_inp_sponge[0] = zbot; + U_inp_sponge[0] = U_inp_sponge_tmp[0]; + V_inp_sponge[0] = V_inp_sponge_tmp[0]; + for (int k=0; k < Nz; ++k) { + z_inp_sponge[k+1] = (use_terrain) ? 0.5 * (zlevels_stag[k] + zlevels_stag[k+1]) + : zbot + (k + 0.5) * dz; + U_inp_sponge[k+1] = interpolate_1d(z_inp_sponge_tmp.dataPtr(), U_inp_sponge_tmp.dataPtr(), z_inp_sponge[k+1], Ninp); + V_inp_sponge[k+1] = interpolate_1d(z_inp_sponge_tmp.dataPtr(), V_inp_sponge_tmp.dataPtr(), z_inp_sponge[k+1], Ninp); + } + z_inp_sponge[Nz+1] = ztop; + U_inp_sponge[Nz+1] = interpolate_1d(z_inp_sponge_tmp.dataPtr(), U_inp_sponge_tmp.dataPtr(), ztop, Ninp); + V_inp_sponge[Nz+1] = interpolate_1d(z_inp_sponge_tmp.dataPtr(), V_inp_sponge_tmp.dataPtr(), ztop, Ninp); + } + + amrex::Print() << "Successfully read the input_sponge file..." << std::endl; + input_sponge_reader.close(); + } + + int size () const + { + AMREX_ALWAYS_ASSERT(z_inp_sponge.size() == U_inp_sponge.size()); + AMREX_ALWAYS_ASSERT(z_inp_sponge.size() == V_inp_sponge.size()); + return z_inp_sponge.size(); + } + + // Members + // - read from file + amrex::Vector z_inp_sponge, U_inp_sponge, V_inp_sponge; +}; +#endif diff --git a/Source/Initialization/Make.package b/Source/Initialization/Make.package index f84d76f20..685889786 100644 --- a/Source/Initialization/Make.package +++ b/Source/Initialization/Make.package @@ -2,6 +2,7 @@ CEXE_sources += ERF_init_uniform.cpp CEXE_sources += ERF_init_from_hse.cpp CEXE_sources += ERF_init_custom.cpp 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 @@ -10,6 +11,7 @@ CEXE_sources += ERF_init_windfarm.cpp endif CEXE_headers += InputSoundingData.H +CEXE_headers += InputSpongeData.H ifeq ($(USE_NETCDF),TRUE) CEXE_headers += Metgrid_utils.H diff --git a/Source/Radiation/Modal_aero_wateruptake.H b/Source/Radiation/Modal_aero_wateruptake.H index 19533b05c..2900550a9 100644 --- a/Source/Radiation/Modal_aero_wateruptake.H +++ b/Source/Radiation/Modal_aero_wateruptake.H @@ -311,13 +311,12 @@ class ModalAeroWateruptake { r=rdry*(1.+p*third/(1.-slog*rdry/a)); } else { - // AML TODO: FIX WHATEVER THIS IS TRYING TO DO... amrex::GpuComplex cx4[4]; makoh_quartic(cx4,p43,p42,p41,p40); //find smallest real(r8) solution r = 1000.*rdry; - auto nsol = 0; - for(auto n=0; n<4; ++n) { + int nsol = -1; + for(int n=0; n<4; ++n) { auto xr=cx4[n].real(); auto xi=cx4[n].imag(); if(abs(xi) > abs(xr)*eps) continue; @@ -327,7 +326,7 @@ class ModalAeroWateruptake { r=xr; nsol=n; } - if(nsol != 0) { + if(nsol == -1) { amrex::Print() << "kohlerc: no real solution found(quartic), nsol= " << nsol << "\n"; r=rdry; } @@ -346,7 +345,7 @@ class ModalAeroWateruptake { makoh_cubic(cx3,p32,p31,p30); //find smallest real(r8) solution r=1000.*rdry; - auto nsol = 0; + int nsol = -1; for(auto n=0; n<3; ++n) { auto xr = cx3[n].real(); auto xi = cx3[n].imag(); @@ -357,8 +356,8 @@ class ModalAeroWateruptake { r=xr; nsol=n; } - if(nsol == 0) { - amrex::Print() << "kohlerc: no real solution found (cubic)\n"; + if(nsol == -1) { + amrex::Print() << "kohlerc: no real solution found (cubic), nsol= " << nsol << "\n"; r=rdry; } } @@ -374,10 +373,10 @@ class ModalAeroWateruptake { // solves x**3 + p2 x**2 + p1 x + p0 = 0 // where p0, p1, p2 are real static - void makoh_cubic(amrex::GpuComplex cx[], - const real& p2, - const real& p1, - const real& p0) + void makoh_cubic (amrex::GpuComplex cx[], + const real& p2, + const real& p1, + const real& p0) { const real eps = 1.0e-20; const real third=1./3.; @@ -413,11 +412,11 @@ class ModalAeroWateruptake { // solves x**4 + p3 x**3 + p2 x**2 + p1 x + p0 = 0 // where p0, p1, p2, p3 are real static - void makoh_quartic(amrex::GpuComplex cx[], - const real& p3, - const real& p2, - const real& p1, - const real& p0) + void makoh_quartic (amrex::GpuComplex cx[], + const real& p3, + const real& p2, + const real& p1, + const real& p0) { const real third=1./3.; auto q = -p2*p2/36.+(p3*p1-4*p0)/12.; diff --git a/Source/Radiation/Phys_prop.H b/Source/Radiation/Phys_prop.H index e5de2832a..54c30ba34 100644 --- a/Source/Radiation/Phys_prop.H +++ b/Source/Radiation/Phys_prop.H @@ -699,6 +699,8 @@ class PhysProp { // Read optics data for modal aerosols void modal_optics_init (physprop_t& phys_prop, yakl::SimpleNetCDF& prop) { + // NOTE: Definitions for real arrays come from rrtmgp_const.h + // and they default to styleFortran ordering. realHost5d extpsw, abspsw, asmpsw, absplw; auto nlwbnd = prop.getDimSize( "lw_band" ); auto nswbnd = prop.getDimSize( "sw_band" ); @@ -714,23 +716,24 @@ class PhysProp { prop.read(asmpsw, "asmpsw" ); prop.read(absplw, "absplw" ); - phys_prop.extpsw = real4d("extpsw", nswbnd, prefi, prefr, ncoef); - phys_prop.abspsw = real4d("abspsw", nswbnd, prefi, prefr, ncoef); - phys_prop.asmpsw = real4d("asmpsw", nswbnd, prefi, prefr, ncoef); - phys_prop.absplw = real4d("absplw", nlwbnd, prefi, prefr, ncoef); + // styleFortran ordering to be consistent with realHost5d definition + phys_prop.extpsw = real4d("extpsw", ncoef, prefr, prefi, nswbnd); + phys_prop.abspsw = real4d("abspsw", ncoef, prefr, prefi, nswbnd); + phys_prop.asmpsw = real4d("asmpsw", ncoef, prefr, prefi, nswbnd); + phys_prop.absplw = real4d("absplw", ncoef, prefr, prefi, nswbnd); - parallel_for (SimpleBounds<4>(nswbnd, prefi, prefr, ncoef), + parallel_for (SimpleBounds<4>(nswbnd, prefr, prefi, ncoef), YAKL_LAMBDA (int i, int j, int k, int l) { - phys_prop.extpsw(i,j,k,l) = extpsw(i,1,j,k,l); - phys_prop.abspsw(i,j,k,l) = abspsw(i,1,j,k,l); - phys_prop.asmpsw(i,j,k,l) = asmpsw(i,1,j,k,l); + phys_prop.extpsw(i,j,k,l) = extpsw(i,j,k,1,l); + phys_prop.abspsw(i,j,k,l) = abspsw(i,j,k,1,l); + phys_prop.asmpsw(i,j,k,l) = asmpsw(i,j,k,1,l); }); - parallel_for (SimpleBounds<4>(nlwbnd, prefi, prefr, ncoef), + parallel_for (SimpleBounds<4>(nlwbnd, prefr, prefi, ncoef), YAKL_LAMBDA (int i, int j, int k, int l) { - phys_prop.absplw(i,j,k,l) = absplw(i,1,j,k,l); + phys_prop.absplw(i,j,k,l) = absplw(i,j,k,1,l); }); prop.read(phys_prop.refrtabsw, "refindex_real_sw" ); diff --git a/Source/Radiation/Radiation.H b/Source/Radiation/Radiation.H index 6f86be73a..cbc731a37 100644 --- a/Source/Radiation/Radiation.H +++ b/Source/Radiation/Radiation.H @@ -44,6 +44,7 @@ class Radiation { // init void initialize (const amrex::MultiFab& cons_in, + amrex::MultiFab* qheating_rates, amrex::Vector qmoist, const amrex::BoxArray& grids, const amrex::Geometry& geom, @@ -93,6 +94,9 @@ class Radiation { // valid boxes on which to evolve the solution amrex::BoxArray m_box; + // Pointer to the radiation source terms + amrex::MultiFab* qrad_src; + // number of vertical levels int nlev, zlo, zhi; @@ -176,9 +180,11 @@ class Radiation { // k-distribution coefficients files to read from. These are set via namelist // variables. - std::string rrtmgp_file_path; - std::string rrtmgp_coefficients_file_sw = "rrtmgp_coefficients_sw_20181204.nc"; - std::string rrtmgp_coefficients_file_lw = "rrtmgp_coefficients_lw_20181204.nc"; + std::string rrtmgp_data_path; + std::string rrtmgp_coefficients_file_sw; + std::string rrtmgp_coefficients_file_lw; + std::string rrtmgp_coefficients_file_name_sw = "rrtmgp_coefficients_sw_20181204.nc"; + std::string rrtmgp_coefficients_file_name_lw = "rrtmgp_coefficients_lw_20181204.nc"; // Band midpoints; these need to be module variables because of how cam_history works; // add_hist_coord sets up pointers to these, so they need to persist. diff --git a/Source/Radiation/Radiation.cpp b/Source/Radiation/Radiation.cpp index cc5ffddc1..aa2ac16f1 100644 --- a/Source/Radiation/Radiation.cpp +++ b/Source/Radiation/Radiation.cpp @@ -36,13 +36,14 @@ using yakl::fortran::SimpleBounds; namespace internal { void initial_fluxes (int nz, int nlay, int nbands, FluxesByband& fluxes) { - fluxes.flux_up = real2d("flux_up", nz, nlay+1); - fluxes.flux_dn = real2d("flux_dn", nz, nlay+1); - fluxes.flux_net = real2d("flux_net", nz, nlay+1); + fluxes.flux_up = real2d("flux_up" , nz, nlay+1); + fluxes.flux_dn = real2d("flux_dn" , nz, nlay+1); + fluxes.flux_net = real2d("flux_net" , nz, nlay+1); fluxes.flux_dn_dir = real2d("flux_dn_dir", nz, nlay+1); - fluxes.bnd_flux_up = real3d("flux_up", nz, nlay+1, nbands); - fluxes.bnd_flux_dn = real3d("flux_dn", nz, nlay+1, nbands); - fluxes.bnd_flux_net = real3d("flux_net", nz, nlay+1, nbands); + + fluxes.bnd_flux_up = real3d("flux_up" , nz, nlay+1, nbands); + fluxes.bnd_flux_dn = real3d("flux_dn" , nz, nlay+1, nbands); + fluxes.bnd_flux_net = real3d("flux_net" , nz, nlay+1, nbands); fluxes.bnd_flux_dn_dir = real3d("flux_dn_dir", nz, nlay+1, nbands); } @@ -54,32 +55,32 @@ namespace internal { auto nlev = size(daytime_fluxes.bnd_flux_up, 2); auto nbnds = size(daytime_fluxes.bnd_flux_up, 3); - int1d nday_1d("nday_1d", 1), - nday_host("nday_host",1); + int1d nday_1d("nday_1d", 1),nday_host("nday_host",1); yakl::memset(nday_1d, 0); parallel_for(SimpleBounds<1>(ncol), YAKL_LAMBDA (int icol) { if (day_indices(icol) > 0) nday_1d(1)++; - printf("daynight indices(check): %d, %d, %d\n",icol,day_indices(icol),nday_1d(1)); + //printf("daynight indices(check): %d, %d, %d\n",icol,day_indices(icol),nday_1d(1)); }); nday_1d.deep_copy_to(nday_host); auto nday = nday_host(1); + AMREX_ASSERT_WITH_MESSAGE((nday>0) && (nday<=ncol), "RADIATION: Invalid number of days!"); parallel_for(SimpleBounds<3>(nday, nlev, nbnds), YAKL_LAMBDA (int iday, int ilev, int ibnd) { // Map daytime index to proper column index - // auto icol = day_indices(iday); - auto icol = iday; + auto icol = day_indices(iday); + //auto icol = iday; // Expand broadband fluxes - expanded_fluxes.flux_up(icol,ilev) = daytime_fluxes.flux_up(iday,ilev); - expanded_fluxes.flux_dn(icol,ilev) = daytime_fluxes.flux_dn(iday,ilev); - expanded_fluxes.flux_net(icol,ilev) = daytime_fluxes.flux_net(iday,ilev); + expanded_fluxes.flux_up(icol,ilev) = daytime_fluxes.flux_up(iday,ilev); + expanded_fluxes.flux_dn(icol,ilev) = daytime_fluxes.flux_dn(iday,ilev); + expanded_fluxes.flux_net(icol,ilev) = daytime_fluxes.flux_net(iday,ilev); expanded_fluxes.flux_dn_dir(icol,ilev) = daytime_fluxes.flux_dn_dir(iday,ilev); // Expand band-by-band fluxes - expanded_fluxes.bnd_flux_up(icol,ilev,ibnd) = daytime_fluxes.bnd_flux_up(iday,ilev,ibnd); - expanded_fluxes.bnd_flux_dn(icol,ilev,ibnd) = daytime_fluxes.bnd_flux_dn(iday,ilev,ibnd); - expanded_fluxes.bnd_flux_net(icol,ilev,ibnd) = daytime_fluxes.bnd_flux_net(iday,ilev,ibnd); + expanded_fluxes.bnd_flux_up(icol,ilev,ibnd) = daytime_fluxes.bnd_flux_up(iday,ilev,ibnd); + expanded_fluxes.bnd_flux_dn(icol,ilev,ibnd) = daytime_fluxes.bnd_flux_dn(iday,ilev,ibnd); + expanded_fluxes.bnd_flux_net(icol,ilev,ibnd) = daytime_fluxes.bnd_flux_net(iday,ilev,ibnd); expanded_fluxes.bnd_flux_dn_dir(icol,ilev,ibnd) = daytime_fluxes.bnd_flux_dn_dir(iday,ilev,ibnd); }); } @@ -97,6 +98,7 @@ namespace internal { // init void Radiation::initialize (const MultiFab& cons_in, + MultiFab* qheating_rates, Vector qmoist, const BoxArray& grids, const Geometry& geom, @@ -110,6 +112,8 @@ void Radiation::initialize (const MultiFab& cons_in, m_geom = geom; m_box = grids; + qrad_src = qheating_rates; + auto dz = m_geom.CellSize(2); auto lowz = m_geom.ProbLo(2); @@ -121,9 +125,9 @@ void Radiation::initialize (const MultiFab& cons_in, do_snow_optics = do_snow_opt; is_cmip6_volc = is_cmip6_volcano; - rrtmgp_file_path = getRadiationDataDir() + "/"; - rrtmgp_coefficients_file_sw.insert(0,rrtmgp_file_path); - rrtmgp_coefficients_file_lw.insert(0,rrtmgp_file_path); + rrtmgp_data_path = getRadiationDataDir() + "/"; + rrtmgp_coefficients_file_sw = rrtmgp_data_path + rrtmgp_coefficients_file_name_sw; + rrtmgp_coefficients_file_lw = rrtmgp_data_path + rrtmgp_coefficients_file_name_lw; for ( MFIter mfi(cons_in, TilingIfNotGPU()); mfi.isValid(); ++mfi) { @@ -173,6 +177,7 @@ void Radiation::initialize (const MultiFab& cons_in, for ( MFIter mfi(cons_in, false); mfi.isValid(); ++mfi) { auto states_array = cons_in.array(mfi); auto qt_array = (qmoist[0]) ? qmoist[0]->array(mfi) : Array4 {}; + auto qv_array = (qmoist[1]) ? qmoist[1]->array(mfi) : Array4 {}; auto qc_array = (qmoist[2]) ? qmoist[2]->array(mfi) : Array4 {}; auto qi_array = (qmoist.size()>=8) ? qmoist[3]->array(mfi) : Array4 {}; @@ -184,26 +189,32 @@ void Radiation::initialize (const MultiFab& cons_in, { auto icol = j*nx+i+1; auto ilev = k+1; + Real qv = (qv_array) ? qv_array(i,j,k): 0.0; qt(icol,ilev) = (qt_array) ? qt_array(i,j,k): 0.0; qc(icol,ilev) = (qc_array) ? qc_array(i,j,k): 0.0; qi(icol,ilev) = (qi_array) ? qi_array(i,j,k): 0.0; qn(icol,ilev) = qc(icol,ilev) + qi(icol,ilev); - tmid(icol,ilev) = getTgivenRandRTh(states_array(i,j,k,Rho_comp),states_array(i,j,k,RhoTheta_comp)); - pmid(icol,ilev) = getPgivenRTh(states_array(i,j,k,RhoTheta_comp))/1.0e3; + tmid(icol,ilev) = getTgivenRandRTh(states_array(i,j,k,Rho_comp),states_array(i,j,k,RhoTheta_comp),qv); + // NOTE: RRTMGP code expects pressure in mb so we convert it here + pmid(icol,ilev) = getPgivenRTh(states_array(i,j,k,RhoTheta_comp),qv) * 1.0e-2; }); } parallel_for(SimpleBounds<2>(ncol, nlev+1), YAKL_LAMBDA (int icol, int ilev) { if (ilev == 1) { - pint(icol, 1) = 2.*pmid(icol, 2) - pmid(icol, 1); - tint(icol, 1) = 2.*tmid(icol, 2) - tmid(icol, 1); + //pint(icol, 1) = 2.*pmid(icol, 2) - pmid(icol, 1); + //tint(icol, 1) = 2.*tmid(icol, 2) - tmid(icol, 1); + pint(icol, 1) = -0.5*pmid(icol, 2) + 1.5*pmid(icol, 1); + tint(icol, 1) = -0.5*tmid(icol, 2) + 1.5*tmid(icol, 1); } else if (ilev <= nlev) { pint(icol, ilev) = 0.5*(pmid(icol, ilev-1) + pmid(icol, ilev)); tint(icol, ilev) = 0.5*(tmid(icol, ilev-1) + tmid(icol, ilev)); } else { - pint(icol, nlev+1) = 2.*pmid(icol, nlev-1) - pmid(icol, nlev); - tint(icol, nlev+1) = 2.*tmid(icol, nlev-1) - tmid(icol, nlev); + //pint(icol, nlev+1) = 2.*pmid(icol, nlev-1) - pmid(icol, nlev); + //tint(icol, nlev+1) = 2.*tmid(icol, nlev-1) - tmid(icol, nlev); + pint(icol, nlev+1) = -0.5*pmid(icol, nlev-1) + 1.5*pmid(icol, nlev); + tint(icol, nlev+1) = -0.5*tmid(icol, nlev-1) + 1.5*tmid(icol, nlev); } }); @@ -248,42 +259,42 @@ void Radiation::initialize (const MultiFab& cons_in, // run radiation model void Radiation::run () { - // Temporary variable for heating rate output - real2d hr("hr", ncol, nlev); - // Cosine solar zenith angle for all columns in chunk real1d coszrs("coszrs", ncol); // Pointers to fields on the physics buffer - real2d cld("cld", ncol, nlev), cldfsnow("cldfsnow", ncol, nlev), - iclwp("iclwp", ncol, nlev), iciwp("iciwp", ncol, nlev), - icswp("icswp", ncol, nlev), dei("dei", ncol, nlev), - des("des", ncol, nlev), lambdac("lambdac", ncol, nlev), - mu("mu", ncol, nlev), rei("rei", ncol, nlev), rel("rel", ncol, nlev); + real2d cld("cld" , ncol, nlev), cldfsnow("cldfsnow", ncol, nlev), + iclwp("iclwp", ncol, nlev), iciwp("iciwp" , ncol, nlev), + icswp("icswp", ncol, nlev), dei("dei" , ncol, nlev), + des("des" , ncol, nlev), lambdac("lambdac" , ncol, nlev), + mu("mu" , ncol, nlev), rei("rei" , ncol, nlev), + rel("rel" , ncol, nlev); // Cloud, snow, and aerosol optical properties real3d cld_tau_gpt_sw("cld_tau_gpt_sw", ncol, nlev, nswgpts), - cld_ssa_gpt_sw("cld_ssa_gpt_sw", ncol, nlev, nswgpts), - cld_asm_gpt_sw("cld_asm_gpt_sw", ncol, nlev, nswgpts); + cld_ssa_gpt_sw("cld_ssa_gpt_sw", ncol, nlev, nswgpts), + cld_asm_gpt_sw("cld_asm_gpt_sw", ncol, nlev, nswgpts); + real3d cld_tau_bnd_sw("cld_tau_bnd_sw", ncol, nlev, nswbands), - cld_ssa_bnd_sw("cld_ssa_bnd_sw", ncol, nlev, nswbands), - cld_asm_bnd_sw("cld_asm_bnd_sw", ncol, nlev, nswbands); + cld_ssa_bnd_sw("cld_ssa_bnd_sw", ncol, nlev, nswbands), + cld_asm_bnd_sw("cld_asm_bnd_sw", ncol, nlev, nswbands); real3d aer_tau_bnd_sw("aer_tau_bnd_sw", ncol, nlev, nswbands), - aer_ssa_bnd_sw("aer_ssa_bnd_sw", ncol, nlev, nswbands), - aer_asm_bnd_sw("aer_asm_bnd_sw", ncol, nlev, nswbands); + aer_ssa_bnd_sw("aer_ssa_bnd_sw", ncol, nlev, nswbands), + aer_asm_bnd_sw("aer_asm_bnd_sw", ncol, nlev, nswbands); real3d cld_tau_bnd_lw("cld_tau_bnd_lw", ncol, nlev, nlwbands), - aer_tau_bnd_lw("aer_tau_bnd_lw", ncol, nlev, nlwbands); + aer_tau_bnd_lw("aer_tau_bnd_lw", ncol, nlev, nlwbands); + real3d cld_tau_gpt_lw("cld_tau_gpt_lw", ncol, nlev, nlwgpts); // NOTE: these are diagnostic only real3d liq_tau_bnd_sw("liq_tau_bnd_sw", ncol, nlev, nswbands), - ice_tau_bnd_sw("ice_tau_bnd_sw", ncol, nlev, nswbands), - snw_tau_bnd_sw("snw_tau_bnd_sw", ncol, nlev, nswbands); + ice_tau_bnd_sw("ice_tau_bnd_sw", ncol, nlev, nswbands), + snw_tau_bnd_sw("snw_tau_bnd_sw", ncol, nlev, nswbands); real3d liq_tau_bnd_lw("liq_tau_bnd_lw", ncol, nlev, nlwbands), - ice_tau_bnd_lw("ice_tau_bnd_lw", ncol, nlev, nlwbands), - snw_tau_bnd_lw("snw_tau_bnd_lw", ncol, nlev, nlwbands); + ice_tau_bnd_lw("ice_tau_bnd_lw", ncol, nlev, nlwbands), + snw_tau_bnd_lw("snw_tau_bnd_lw", ncol, nlev, nlwbands); // Gas volume mixing ratios real3d gas_vmr("gas_vmr", ngas, ncol, nlev); @@ -309,8 +320,8 @@ void Radiation::run () if (do_short_wave_rad) { // Radiative fluxes FluxesByband fluxes_allsky, fluxes_clrsky; - internal::initial_fluxes(ncol, nlev, nlwbands, fluxes_allsky); - internal::initial_fluxes(ncol, nlev, nlwbands, fluxes_clrsky); + internal::initial_fluxes(ncol, nlev+1, nswbands, fluxes_allsky); + internal::initial_fluxes(ncol, nlev+1, nswbands, fluxes_clrsky); // Get cosine solar zenith angle for current time step. ( still NOT YET implemented here) // set_cosine_solar_zenith_angle(state, dt_avg, coszrs(1:ncol)) @@ -508,18 +519,20 @@ void Radiation::run () } } // dolw - // Compute heating rate for dtheta/dt - parallel_for(SimpleBounds<2>(ncol, nlev), YAKL_LAMBDA (int icol, int ilay) - { - hr(icol,ilay) = (qrs(icol,ilay) + qrl(icol,ilay)) / Cp_d * (1.e5 / std::pow(pmid(icol,ilay), R_d/Cp_d)); - }); - - // convert radiative heating rates to Q*dp for energy conservation - if (conserve_energy) { - parallel_for(SimpleBounds<2>(ncol, nlev), YAKL_LAMBDA (int icol, int ilev) + // Populate source term for theta dycore variable + for (MFIter mfi(*(qrad_src)); mfi.isValid(); ++mfi) { + auto qrad_src_array = qrad_src->array(mfi); + const auto& box3d = mfi.tilebox(); + auto nx = box3d.length(0); + amrex::ParallelFor(box3d, [=] AMREX_GPU_DEVICE (int i, int j, int k) { - qrs(icol,ilev) = qrs(icol,ilev)*pdel(icol,ilev); - qrl(icol,ilev) = qrl(icol,ilev)*pdel(icol,ilev); + // Map (col,lev) to (i,j,k) + auto icol = j*nx+i+1; + auto ilev = k+1; + + // SW and LW sources + qrad_src_array(i,j,k,0) = qrs(icol,ilev); + qrad_src_array(i,j,k,1) = qrl(icol,ilev); }); } } @@ -696,14 +709,12 @@ void Radiation::radiation_driver_lw (int ncol, int nlev, FluxesByband& fluxes_allsky, const real2d& qrl, const real2d& qrlc) { real3d cld_tau_gpt_rad("cld_tau_gpt_rad", ncol, nlev+1, nlwgpts); - real3d aer_tau_bnd_rad("aer_tau_bnd-rad", ncol, nlev+1, nlwgpts); + real3d aer_tau_bnd_rad("aer_tau_bnd_rad", ncol, nlev+1, nlwgpts); // Surface emissivity needed for longwave real2d surface_emissivity("surface_emissivity", nlwbands, ncol); // Temporary heating rates on radiation vertical grid - real2d qrl_rad("qrl_rad", ncol, nlev); - real2d qrlc_rad("qrlc_rad", ncol, nlev); real3d gas_vmr_rad("gas_vmr_rad", ngas, ncol, nlev); // Set surface emissivity to 1 here. There is a note in the RRTMG @@ -718,7 +729,7 @@ void Radiation::radiation_driver_lw (int ncol, int nlev, yakl::memset(cld_tau_gpt_rad, 0.); yakl::memset(aer_tau_bnd_rad, 0.); - parallel_for(SimpleBounds<3>(ncol, nlev, nswgpts), YAKL_LAMBDA (int icol, int ilev, int igpt) + parallel_for(SimpleBounds<3>(ncol, nlev, nlwgpts), YAKL_LAMBDA (int icol, int ilev, int igpt) { cld_tau_gpt_rad(icol,ilev,igpt) = cld_tau_gpt(icol,ilev,igpt); aer_tau_bnd_rad(icol,ilev,igpt) = aer_tau_bnd(icol,ilev,igpt); @@ -737,18 +748,11 @@ void Radiation::radiation_driver_lw (int ncol, int nlev, // Calculate heating rates calculate_heating_rate(fluxes_allsky.flux_up, fluxes_allsky.flux_dn, - pint, qrl_rad); + pint, qrl); calculate_heating_rate(fluxes_allsky.flux_up, fluxes_allsky.flux_dn, - pint,qrlc_rad); - - // Map heating rates to CAM columns and levels - parallel_for(SimpleBounds<2>(ncol, nlev), YAKL_LAMBDA (int icol, int ilev) - { - qrl(icol,ilev) = qrl_rad(icol,ilev); - qrlc(icol,ilev) = qrlc_rad(icol,ilev); - }); + pint, qrlc); } // Initialize array of daytime indices to be all zero. If any zeros exist when @@ -885,10 +889,29 @@ void Radiation::calculate_heating_rate (const real2d& flux_up, const real2d& pint, const real2d& heating_rate) { + // NOTE: The pressure is in [mb] for RRTMGP to use. + // The fluxes are in [W/m^2] and gravity is [m/s^2]. + // We need to convert pressure from [mb] -> [Pa] and divide by Cp [J/kg*K] + // The heating rate is {dF/dP * g / Cp} with units [K/s] + real1d heatfac("heatfac",1); + yakl::memset(heatfac, 1.0e-2/Cp_d); parallel_for(SimpleBounds<2>(ncol, nlev), YAKL_LAMBDA (int icol, int ilev) { - heating_rate(icol,ilev) = (flux_up(icol,ilev+1) - flux_up(icol,ilev)- flux_dn(icol,ilev+1)+flux_dn(icol,ilev)) - *CONST_GRAV/(pint(icol,ilev+1)-pint(icol,ilev)); + heating_rate(icol,ilev) = heatfac(1) * ( (flux_up(icol,ilev+1) - flux_dn(icol,ilev+1)) + - (flux_up(icol,ilev ) - flux_dn(icol,ilev )) ) + * CONST_GRAV/(pint(icol,ilev+1)-pint(icol,ilev)); + /* + if (icol==1) { + amrex::Print() << "HR: " << ilev << ' ' + << heating_rate(icol,ilev) << ' ' + << (flux_up(icol,ilev+1) - flux_dn(icol,ilev+1)) << ' ' + << (flux_up(icol,ilev ) - flux_dn(icol,ilev )) << ' ' + << (flux_up(icol,ilev+1) - flux_dn(icol,ilev+1)) + - (flux_up(icol,ilev ) - flux_dn(icol,ilev )) << ' ' + << (pint(icol,ilev+1)-pint(icol,ilev)) << ' ' + << CONST_GRAV << "\n"; + } + */ }); } diff --git a/Source/Radiation/Run_longwave_rrtmgp.cpp b/Source/Radiation/Run_longwave_rrtmgp.cpp index d66a6719f..ddab671ff 100644 --- a/Source/Radiation/Run_longwave_rrtmgp.cpp +++ b/Source/Radiation/Run_longwave_rrtmgp.cpp @@ -98,12 +98,13 @@ void Rrtmgp::run_longwave_rrtmgp (int ngas, int ncol, int nlay, // Do the clearsky calculation before adding in clouds FluxesByband fluxes_clrsky; - fluxes_clrsky.flux_up = real2d("clrsky_flux_up", ncol, nlay+1); // clrsky_flux_up; - fluxes_clrsky.flux_dn = real2d("clrsky_flux_up", ncol, nlay+1); //clrsky_flux_dn; - fluxes_clrsky.flux_net = real2d("clrsky_flux_up", ncol, nlay+1); //clrsky_flux_net; - fluxes_clrsky.bnd_flux_up = real3d("clrsky_flux_up", ncol, nlay+1, nlwbands); //clrsky_bnd_flux_up; - fluxes_clrsky.bnd_flux_dn = real3d("clrsky_flux_up", ncol, nlay+1, nlwbands); //clrsky_bnd_flux_dn; - fluxes_clrsky.bnd_flux_net = real3d("clrsky_flux_up", ncol, nlay+1, nlwbands); //clrsky_bnd_flux_net; + fluxes_clrsky.flux_up = real2d("clrsky_flux_up" , ncol, nlay+1); // clrsky_flux_up; + fluxes_clrsky.flux_dn = real2d("clrsky_flux_dn" , ncol, nlay+1); //clrsky_flux_dn; + fluxes_clrsky.flux_net = real2d("clrsky_flux_net", ncol, nlay+1); //clrsky_flux_net; + fluxes_clrsky.bnd_flux_up = real3d("clrsky_bnd_flux_up" , ncol, nlay+1, nlwbands); //clrsky_bnd_flux_up; + fluxes_clrsky.bnd_flux_dn = real3d("clrsky_bnd_flux_dn" , ncol, nlay+1, nlwbands); //clrsky_bnd_flux_dn; + fluxes_clrsky.bnd_flux_net = real3d("clrsky_bnd_flux_net", ncol, nlay+1, nlwbands); //clrsky_bnd_flux_net; + rte_lw(max_gauss_pts, gauss_Ds, gauss_wts, combined_optics, top_at_1, lw_sources, emis_sfc, fluxes_clrsky); // Copy fluxes back out of FluxesByband object @@ -125,12 +126,13 @@ void Rrtmgp::run_longwave_rrtmgp (int ngas, int ncol, int nlay, // Call LW flux driver FluxesByband fluxes_allsky; - fluxes_allsky.flux_up = real2d("flux_up", ncol, nlay+1); //allsky_flux_up; - fluxes_allsky.flux_dn = real2d("flux_dn", ncol, nlay+1); //allsky_flux_dn; - fluxes_allsky.flux_net = real2d("flux_net", ncol, nlay+1); //allsky_flux_net; - fluxes_allsky.bnd_flux_up = real3d("flux_up", ncol, nlay+1, nlwbands); //allsky_bnd_flux_up; - fluxes_allsky.bnd_flux_dn = real3d("flux_dn", ncol, nlay+1, nlwbands); //allsky_bnd_flux_dn; - fluxes_allsky.bnd_flux_net = real3d("flux_net", ncol, nlay+1, nlwbands); //allsky_bnd_flux_net; + fluxes_allsky.flux_up = real2d("allsky_flux_up" , ncol, nlay+1); //allsky_flux_up; + fluxes_allsky.flux_dn = real2d("allsky_flux_dn" , ncol, nlay+1); //allsky_flux_dn; + fluxes_allsky.flux_net = real2d("allsky_flux_net", ncol, nlay+1); //allsky_flux_net; + fluxes_allsky.bnd_flux_up = real3d("allsky_bnd_flux_up" , ncol, nlay+1, nlwbands); //allsky_bnd_flux_up; + fluxes_allsky.bnd_flux_dn = real3d("allsky_bnd_flux_dn" , ncol, nlay+1, nlwbands); //allsky_bnd_flux_dn; + fluxes_allsky.bnd_flux_net = real3d("allsky_bnd_flux_net", ncol, nlay+1, nlwbands); //allsky_bnd_flux_net; + rte_lw(max_gauss_pts, gauss_Ds, gauss_wts, combined_optics, top_at_1, lw_sources, emis_sfc, fluxes_allsky); // Copy fluxes back out of FluxesByband object diff --git a/Source/Radiation/Run_shortwave_rrtmgp.cpp b/Source/Radiation/Run_shortwave_rrtmgp.cpp index 9d7fcaa6e..086a90688 100644 --- a/Source/Radiation/Run_shortwave_rrtmgp.cpp +++ b/Source/Radiation/Run_shortwave_rrtmgp.cpp @@ -34,16 +34,18 @@ void Rrtmgp::run_shortwave_rrtmgp (int ngas, int ncol, int nlay, bool1d top_at_1_g("top_at_1_g",1); boolHost1d top_at_1_h("top_at_1_h",1); bool top_at_1; + parallel_for(SimpleBounds<1>(1), YAKL_LAMBDA (int icol) { // HACK: Single loop kernel is not efficient + top_at_1_g(1) = pmid(1, 1) < pmid (1, 2); + }); + top_at_1_g.deep_copy_to(top_at_1_h); + top_at_1 = top_at_1_h(1); real2d toa_flux("toa_flux", ncol, nswgpts); - k_dist_sw.gas_optics(ncol, nlay, &top_at_1, pmid, pint, tmid, gas_concs, combined_optics, toa_flux); // TODO: top_at_1 is not a valid address and this doesn't match longwave call + k_dist_sw.gas_optics(ncol, nlay, top_at_1, pmid, pint, tmid, gas_concs, combined_optics, toa_flux); // Apply TOA flux scaling parallel_for(SimpleBounds<2>(nswgpts,ncol), YAKL_LAMBDA (int igpt, int icol) { toa_flux(icol, igpt) = tsi_scaling * toa_flux(icol,igpt); - top_at_1_g(1) = pmid(1, 1) < pmid (1, 2); }); - top_at_1_g.deep_copy_to(top_at_1_h); - top_at_1 = top_at_1_h(1); // Add in aerosol // TODO: should we avoid allocating here? @@ -74,25 +76,26 @@ void Rrtmgp::run_shortwave_rrtmgp (int ngas, int ncol, int nlay, // Do the clearsky calculation before adding in clouds FluxesByband fluxes_clrsky; - fluxes_clrsky.flux_up = real2d("clrsky_flux_up", ncol, nlay+1); // clrsky_flux_up; - fluxes_clrsky.flux_dn = real2d("clrsky_flux_up", ncol, nlay+1); //clrsky_flux_dn; - fluxes_clrsky.flux_dn_dir = real2d("clrsky_flux_up", ncol, nlay+1); //clrsky_flux_dn_dir; - fluxes_clrsky.flux_net = real2d("clrsky_flux_up", ncol, nlay+1); //clrsky_flux_net; - fluxes_clrsky.bnd_flux_up = real3d("clrsky_flux_up", ncol, nlay+1, nswbands); //clrsky_bnd_flux_up; - fluxes_clrsky.bnd_flux_dn = real3d("clrsky_flux_up", ncol, nlay+1, nswbands); //clrsky_bnd_flux_dn; - fluxes_clrsky.bnd_flux_dn_dir = real3d("clrsky_flux_up", ncol, nlay+1, nswbands); //clrsky_bnd_flux_dn_dir; - fluxes_clrsky.bnd_flux_net = real3d("clrsky_flux_up", ncol, nlay+1, nswbands); //clrsky_bnd_flux_net; + fluxes_clrsky.flux_up = real2d("clrsky_flux_up" , ncol, nlay+1); // clrsky_flux_up; + fluxes_clrsky.flux_dn = real2d("clrsky_flux_nd" , ncol, nlay+1); //clrsky_flux_dn; + fluxes_clrsky.flux_net = real2d("clrsky_flux_net", ncol, nlay+1); //clrsky_flux_net; + fluxes_clrsky.flux_dn_dir = real2d("clrsky_flux_dn_dir", ncol, nlay+1); //clrsky_flux_dn_dir; + fluxes_clrsky.bnd_flux_up = real3d("clrsky_bnd_flux_up" , ncol, nlay+1, nswbands); //clrsky_bnd_flux_up; + fluxes_clrsky.bnd_flux_dn = real3d("clrsky_bnd_flux_dn" , ncol, nlay+1, nswbands); //clrsky_bnd_flux_dn; + fluxes_clrsky.bnd_flux_net = real3d("clrsky_bnd_flux_net", ncol, nlay+1, nswbands); //clrsky_bnd_flux_net; + fluxes_clrsky.bnd_flux_dn_dir = real3d("clrsky_bnd_flux_dn_dir", ncol, nlay+1, nswbands); //clrsky_bnd_flux_dn_dir; + rte_sw(combined_optics, top_at_1, coszrs, toa_flux, albedo_dir, albedo_dif, fluxes_clrsky); // Copy fluxes back out of FluxesByband object - fluxes_clrsky.flux_up.deep_copy_to(clrsky_flux_up); - fluxes_clrsky.flux_dn.deep_copy_to(clrsky_flux_dn); - fluxes_clrsky.flux_dn_dir.deep_copy_to(clrsky_flux_dn_dir); + fluxes_clrsky.flux_up.deep_copy_to (clrsky_flux_up); + fluxes_clrsky.flux_dn.deep_copy_to (clrsky_flux_dn); fluxes_clrsky.flux_net.deep_copy_to(clrsky_flux_net); - fluxes_clrsky.bnd_flux_up.deep_copy_to(clrsky_bnd_flux_up); - fluxes_clrsky.bnd_flux_dn.deep_copy_to(clrsky_bnd_flux_dn); - fluxes_clrsky.bnd_flux_dn_dir.deep_copy_to(clrsky_bnd_flux_dn_dir); + fluxes_clrsky.flux_dn_dir.deep_copy_to(clrsky_flux_dn_dir); + fluxes_clrsky.bnd_flux_up.deep_copy_to (clrsky_bnd_flux_up); + fluxes_clrsky.bnd_flux_dn.deep_copy_to (clrsky_bnd_flux_dn); fluxes_clrsky.bnd_flux_net.deep_copy_to(clrsky_bnd_flux_net); + fluxes_clrsky.bnd_flux_dn_dir.deep_copy_to(clrsky_bnd_flux_dn_dir); // Add in clouds OpticalProps2str cloud_optics; @@ -110,25 +113,26 @@ void Rrtmgp::run_shortwave_rrtmgp (int ngas, int ncol, int nlay, // Call SW flux driver FluxesByband fluxes_allsky; - fluxes_allsky.flux_up = real2d("allsky_flux_up", ncol, nlay+1); // allsky_flux_up; - fluxes_allsky.flux_dn = real2d("allsky_flux_up", ncol, nlay+1); //allsky_flux_dn; - fluxes_allsky.flux_dn_dir = real2d("allsky_flux_up", ncol, nlay+1); //allsky_flux_dn_dir; - fluxes_allsky.flux_net = real2d("allsky_flux_up", ncol, nlay+1); //allsky_flux_net; - fluxes_allsky.bnd_flux_up = real3d("allsky_flux_up", ncol, nlay+1, nswbands); //allsky_bnd_flux_up; - fluxes_allsky.bnd_flux_dn = real3d("allsky_flux_up", ncol, nlay+1, nswbands); //allsky_bnd_flux_dn; - fluxes_allsky.bnd_flux_dn_dir = real3d("allsky_flux_up", ncol, nlay+1, nswbands); //allsky_bnd_flux_dn_dir; - fluxes_allsky.bnd_flux_net = real3d("allsky_flux_up", ncol, nlay+1, nswbands); //allsky_bnd_flux_net; + fluxes_allsky.flux_up = real2d("allsky_flux_up" , ncol, nlay+1); //allsky_flux_up; + fluxes_allsky.flux_dn = real2d("allsky_flux_dn" , ncol, nlay+1); //allsky_flux_dn; + fluxes_allsky.flux_net = real2d("allsky_flux_net", ncol, nlay+1); //allsky_flux_net; + fluxes_allsky.flux_dn_dir = real2d("allsky_flux_dn_dir", ncol, nlay+1); //allsky_flux_dn_dir; + fluxes_allsky.bnd_flux_up = real3d("allsky_bnd_flux_up" , ncol, nlay+1, nswbands); //allsky_bnd_flux_up; + fluxes_allsky.bnd_flux_dn = real3d("allsky_bnd_flux_dn" , ncol, nlay+1, nswbands); //allsky_bnd_flux_dn; + fluxes_allsky.bnd_flux_net = real3d("allsky_bnd_flux_net", ncol, nlay+1, nswbands); //allsky_bnd_flux_net; + fluxes_allsky.bnd_flux_dn_dir = real3d("allsky_bnd_flux_dn_dir", ncol, nlay+1, nswbands); //allsky_bnd_flux_dn_dir; + rte_sw(combined_optics, top_at_1, coszrs, toa_flux, albedo_dir, albedo_dif, fluxes_allsky); // Copy fluxes back out of FluxesByband object fluxes_allsky.flux_up.deep_copy_to(allsky_flux_up); fluxes_allsky.flux_dn.deep_copy_to(allsky_flux_dn); - fluxes_allsky.flux_dn_dir.deep_copy_to(allsky_flux_dn_dir); fluxes_allsky.flux_net.deep_copy_to(allsky_flux_net); + fluxes_allsky.flux_dn_dir.deep_copy_to(allsky_flux_dn_dir); fluxes_allsky.bnd_flux_up.deep_copy_to(allsky_bnd_flux_up); fluxes_allsky.bnd_flux_dn.deep_copy_to(allsky_bnd_flux_dn); - fluxes_allsky.bnd_flux_dn_dir.deep_copy_to(allsky_bnd_flux_dn_dir); fluxes_allsky.bnd_flux_net.deep_copy_to(allsky_bnd_flux_net); + fluxes_allsky.bnd_flux_dn_dir.deep_copy_to(allsky_bnd_flux_dn_dir); yakl::fence(); } diff --git a/Source/SourceTerms/ERF_ApplySpongeZoneBCs.cpp b/Source/SourceTerms/ERF_ApplySpongeZoneBCs.cpp index d11788907..c8c391446 100644 --- a/Source/SourceTerms/ERF_ApplySpongeZoneBCs.cpp +++ b/Source/SourceTerms/ERF_ApplySpongeZoneBCs.cpp @@ -59,9 +59,9 @@ ApplySpongeZoneBCsForCC ( int jj = amrex::min(amrex::max(j, domlo_y), domhi_y); int kk = amrex::min(amrex::max(k, domlo_z), domhi_z); - Real x = (ii+0.5) * dx[0]; - Real y = (jj+0.5) * dx[1]; - Real z = (kk+0.5) * dx[2]; + Real x = ProbLoArr[0] + (ii+0.5) * dx[0]; + Real y = ProbLoArr[1] + (jj+0.5) * dx[1]; + Real z = ProbLoArr[2] + (kk+0.5) * dx[2]; // x left sponge if(use_xlo_sponge_damping){ @@ -171,9 +171,9 @@ ApplySpongeZoneBCsForMom ( int jj = amrex::min(amrex::max(j, domlo_y), domhi_y); int kk = amrex::min(amrex::max(k, domlo_z), domhi_z); - Real x = ii * dx[0]; - Real y = (jj+0.5) * dx[1]; - Real z = (kk+0.5) * dx[2]; + Real x = ProbLoArr[0] + ii * dx[0]; + Real y = ProbLoArr[1] + (jj+0.5) * dx[1]; + Real z = ProbLoArr[2] + (kk+0.5) * dx[2]; // x lo sponge if(use_xlo_sponge_damping){ @@ -230,9 +230,9 @@ ApplySpongeZoneBCsForMom ( int jj = amrex::min(amrex::max(j, domlo_y), domhi_y); int kk = amrex::min(amrex::max(k, domlo_z), domhi_z); - Real x = (ii+0.5) * dx[0]; - Real y = jj * dx[1]; - Real z = (kk+0.5) * dx[2]; + Real x = ProbLoArr[0] + (ii+0.5) * dx[0]; + Real y = ProbLoArr[1] + jj * dx[1]; + Real z = ProbLoArr[2] + (kk+0.5) * dx[2]; // x lo sponge if(use_xlo_sponge_damping){ @@ -289,9 +289,9 @@ ApplySpongeZoneBCsForMom ( int jj = amrex::min(amrex::max(j, domlo_y), domhi_y); int kk = amrex::min(amrex::max(k, domlo_z), domhi_z); - Real x = (ii+0.5) * dx[0]; - Real y = (jj+0.5) * dx[1]; - Real z = kk * dx[2]; + Real x = ProbLoArr[0] + (ii+0.5) * dx[0]; + Real y = ProbLoArr[1] + (jj+0.5) * dx[1]; + Real z = ProbLoArr[2] + kk * dx[2]; // x left sponge if(use_xlo_sponge_damping){ diff --git a/Source/SourceTerms/ERF_ApplySpongeZoneBCs_ReadFromFile.cpp b/Source/SourceTerms/ERF_ApplySpongeZoneBCs_ReadFromFile.cpp new file mode 100644 index 000000000..ea350f1d9 --- /dev/null +++ b/Source/SourceTerms/ERF_ApplySpongeZoneBCs_ReadFromFile.cpp @@ -0,0 +1,179 @@ +#include +#include +#include + +//#include +//#include +//#include + +using namespace amrex; + +void +ApplySpongeZoneBCsForMom_ReadFromFile ( + const SpongeChoice& spongeChoice, + const Geometry geom, + const Box& tbx, + const Box& tby, + const Array4& cell_data, + const Array4& rho_u_rhs, + const Array4& rho_v_rhs, + const Array4& rho_u, + const Array4& rho_v, + const Vector d_sponge_ptrs_at_lev) +{ + // Domain cell size and real bounds + auto dx = geom.CellSizeArray(); + auto ProbHiArr = geom.ProbHiArray(); + auto ProbLoArr = geom.ProbLoArray(); + + const Real sponge_strength = spongeChoice.sponge_strength; + const int use_xlo_sponge_damping = spongeChoice.use_xlo_sponge_damping; + const int use_xhi_sponge_damping = spongeChoice.use_xhi_sponge_damping; + const int use_ylo_sponge_damping = spongeChoice.use_ylo_sponge_damping; + const int use_yhi_sponge_damping = spongeChoice.use_yhi_sponge_damping; + const int use_zlo_sponge_damping = spongeChoice.use_zlo_sponge_damping; + const int use_zhi_sponge_damping = spongeChoice.use_zhi_sponge_damping; + + const Real xlo_sponge_end = spongeChoice.xlo_sponge_end; + const Real xhi_sponge_start = spongeChoice.xhi_sponge_start; + const Real ylo_sponge_end = spongeChoice.ylo_sponge_end; + const Real yhi_sponge_start = spongeChoice.yhi_sponge_start; + const Real zlo_sponge_end = spongeChoice.zlo_sponge_end; + const Real zhi_sponge_start = spongeChoice.zhi_sponge_start; + + // Domain valid box + const 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); + int domhi_z = domain.bigEnd(2) + 1; + + Real* ubar_sponge = d_sponge_ptrs_at_lev[Sponge::ubar_sponge]; + Real* vbar_sponge = d_sponge_ptrs_at_lev[Sponge::vbar_sponge]; + + if(use_xlo_sponge_damping)AMREX_ALWAYS_ASSERT(xlo_sponge_end > ProbLoArr[0]); + if(use_xhi_sponge_damping)AMREX_ALWAYS_ASSERT(xhi_sponge_start < ProbHiArr[0]); + if(use_ylo_sponge_damping)AMREX_ALWAYS_ASSERT(ylo_sponge_end > ProbLoArr[1]); + if(use_yhi_sponge_damping)AMREX_ALWAYS_ASSERT(yhi_sponge_start < ProbHiArr[1]); + if(use_zlo_sponge_damping)AMREX_ALWAYS_ASSERT(zlo_sponge_end > ProbLoArr[2]); + if(use_zhi_sponge_damping)AMREX_ALWAYS_ASSERT(zhi_sponge_start < ProbHiArr[2]); + + ParallelFor(tbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) + { + int ii = amrex::min(amrex::max(i, domlo_x), domhi_x); + int jj = amrex::min(amrex::max(j, domlo_y), domhi_y); + int kk = amrex::min(amrex::max(k, domlo_z), domhi_z); + + Real x = ProbLoArr[0] + ii * dx[0]; + Real y = ProbLoArr[1] + (jj+0.5) * dx[1]; + Real z = ProbLoArr[2] + (kk+0.5) * dx[2]; + + // x lo sponge + if(use_xlo_sponge_damping){ + if (x < xlo_sponge_end) { + Real xi = (xlo_sponge_end - x) / (xlo_sponge_end - ProbLoArr[0]); + rho_u_rhs(i, j, k) -= sponge_strength * xi * xi * (rho_u(i, j, k) - cell_data(i,j,k,0)*ubar_sponge[k]); + } + } + // x hi sponge + if(use_xhi_sponge_damping){ + if (x > xhi_sponge_start) { + Real xi = (x - xhi_sponge_start) / (ProbHiArr[0] - xhi_sponge_start); + rho_u_rhs(i, j, k) -= sponge_strength * xi * xi * (rho_u(i, j, k) - cell_data(i,j,k,0)*ubar_sponge[k]); + } + } + + // y lo sponge + if(use_ylo_sponge_damping){ + if (y < ylo_sponge_end) { + Real xi = (ylo_sponge_end - y) / (ylo_sponge_end - ProbLoArr[1]); + rho_u_rhs(i, j, k) -= sponge_strength * xi * xi * (rho_u(i, j, k) - cell_data(i,j,k,0)*ubar_sponge[k]); + } + } + // x right sponge + if(use_yhi_sponge_damping){ + if (y > yhi_sponge_start) { + Real xi = (y - yhi_sponge_start) / (ProbHiArr[1] - yhi_sponge_start); + rho_u_rhs(i, j, k) -= sponge_strength * xi * xi * (rho_u(i, j, k) - cell_data(i,j,k,0)*ubar_sponge[k]); + } + } + + // z lo sponge + if(use_zlo_sponge_damping){ + if (z < zlo_sponge_end) { + Real xi = (zlo_sponge_end - z) / (zlo_sponge_end - ProbLoArr[2]); + rho_u_rhs(i, j, k) -= sponge_strength * xi * xi * (rho_u(i, j, k) - cell_data(i,j,k,0)*ubar_sponge[k]); + } + } + + + // z hi sponge + if(use_zhi_sponge_damping){ + if (z > zhi_sponge_start) { + Real xi = (z - zhi_sponge_start) / (ProbHiArr[2] - zhi_sponge_start); + rho_u_rhs(i, j, k) -= sponge_strength * xi * xi * (rho_u(i, j, k) - cell_data(i,j,k,0)*ubar_sponge[k]); + } + } + }); + + + ParallelFor(tby, [=] AMREX_GPU_DEVICE(int i, int j, int k) + { + int ii = amrex::min(amrex::max(i, domlo_x), domhi_x); + int jj = amrex::min(amrex::max(j, domlo_y), domhi_y); + int kk = amrex::min(amrex::max(k, domlo_z), domhi_z); + + Real x = ProbLoArr[0] + (ii+0.5) * dx[0]; + Real y = ProbLoArr[1] + jj * dx[1]; + Real z = ProbLoArr[2] + (kk+0.5) * dx[2]; + + // x lo sponge + if(use_xlo_sponge_damping){ + if (x < xlo_sponge_end) { + Real xi = (xlo_sponge_end - x) / (xlo_sponge_end - ProbLoArr[0]); + rho_v_rhs(i, j, k) -= sponge_strength * xi * xi * (rho_v(i, j, k) - cell_data(i,j,k,0)*vbar_sponge[k]); + } + } + // x hi sponge + if(use_xhi_sponge_damping){ + if (x > xhi_sponge_start) { + Real xi = (x - xhi_sponge_start) / (ProbHiArr[0] - xhi_sponge_start); + rho_v_rhs(i, j, k) -= sponge_strength * xi * xi * (rho_v(i, j, k) - cell_data(i,j,k,0)*vbar_sponge[k]); + } + } + + // y lo sponge + if(use_ylo_sponge_damping){ + if (y < ylo_sponge_end) { + Real xi = (ylo_sponge_end - y) / (ylo_sponge_end - ProbLoArr[1]); + rho_v_rhs(i, j, k) -= sponge_strength * xi * xi * (rho_v(i, j, k) - cell_data(i,j,k,0)*vbar_sponge[k]); + } + } + // x right sponge + if(use_yhi_sponge_damping){ + if (y > yhi_sponge_start) { + Real xi = (y - yhi_sponge_start) / (ProbHiArr[1] - yhi_sponge_start); + rho_v_rhs(i, j, k) -= sponge_strength * xi * xi * (rho_v(i, j, k) - cell_data(i,j,k,0)*vbar_sponge[k]); + } + } + + // z lo sponge + if(use_zlo_sponge_damping){ + if (z < zlo_sponge_end) { + Real xi = (zlo_sponge_end - z) / (zlo_sponge_end - ProbLoArr[2]); + rho_v_rhs(i, j, k) -= sponge_strength * xi * xi * (rho_v(i, j, k) - cell_data(i,j,k,0)*vbar_sponge[k]); + } + } + + + // z hi sponge + if(use_zhi_sponge_damping){ + if (z > zhi_sponge_start) { + Real xi = (z - zhi_sponge_start) / (ProbHiArr[2] - zhi_sponge_start); + rho_v_rhs(i, j, k) -= sponge_strength * xi * xi * (rho_v(i, j, k) - cell_data(i,j,k,0)*vbar_sponge[k]); + } + } + }); +} diff --git a/Source/SourceTerms/ERF_make_mom_sources.cpp b/Source/SourceTerms/ERF_make_mom_sources.cpp index 602feafb8..cbfcecbaf 100644 --- a/Source/SourceTerms/ERF_make_mom_sources.cpp +++ b/Source/SourceTerms/ERF_make_mom_sources.cpp @@ -56,6 +56,7 @@ void make_mom_sources (int /*level*/, const Real* dptr_v_geos, const Real* dptr_wbar_sub, const Vector d_rayleigh_ptrs_at_lev, + const Vector d_sponge_ptrs_at_lev, int n_qstate) { BL_PROFILE_REGION("erf_make_mom_sources()"); @@ -304,8 +305,16 @@ void make_mom_sources (int /*level*/, // ***************************************************************************** // Add SPONGING // ***************************************************************************** - ApplySpongeZoneBCsForMom(solverChoice.spongeChoice, geom, tbx, tby, tbz, + if(solverChoice.spongeChoice.sponge_type == "input_sponge") + { + ApplySpongeZoneBCsForMom_ReadFromFile(solverChoice.spongeChoice, geom, tbx, tby, cell_data, + xmom_src_arr, ymom_src_arr, rho_u, rho_v, d_sponge_ptrs_at_lev); + } + else + { + ApplySpongeZoneBCsForMom(solverChoice.spongeChoice, geom, tbx, tby, tbz, xmom_src_arr, ymom_src_arr, zmom_src_arr, rho_u, rho_v, rho_w); + } } // mfi } diff --git a/Source/SourceTerms/ERF_make_sources.cpp b/Source/SourceTerms/ERF_make_sources.cpp index cb87201a6..c4ea1db24 100644 --- a/Source/SourceTerms/ERF_make_sources.cpp +++ b/Source/SourceTerms/ERF_make_sources.cpp @@ -35,7 +35,7 @@ void make_sources (int level, const MultiFab & S_prim, MultiFab & source, #ifdef ERF_USE_RRTMGP - const MultiFab& qheating_rates, + const MultiFab* qheating_rates, #endif const Geometry geom, const SolverChoice& solverChoice, @@ -150,9 +150,10 @@ void make_sources (int level, // Add radiation source terms to (rho theta) // ************************************************************************************* { - auto const& qheating_arr = qheating_rates.const_array(mfi); + auto const& qheating_arr = qheating_rates->const_array(mfi); ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { + // Short-wavelength and long-wavelength radiation source terms cell_src(i,j,k,RhoTheta_comp) += qheating_arr(i,j,k,0) + qheating_arr(i,j,k,1); }); } @@ -272,8 +273,9 @@ void make_sources (int level, // ************************************************************************************* // Add sponging // ************************************************************************************* - ApplySpongeZoneBCsForCC(solverChoice.spongeChoice, geom, bx, cell_src, cell_data); - + if(!(solverChoice.spongeChoice.sponge_type == "input_sponge")){ + ApplySpongeZoneBCsForCC(solverChoice.spongeChoice, geom, bx, cell_src, cell_data); + } } // mfi } // OMP } diff --git a/Source/SourceTerms/Make.package b/Source/SourceTerms/Make.package index e115d5074..0f1df167e 100644 --- a/Source/SourceTerms/Make.package +++ b/Source/SourceTerms/Make.package @@ -3,6 +3,7 @@ CEXE_sources += ERF_make_buoyancy.cpp CEXE_sources += ERF_make_mom_sources.cpp CEXE_sources += ERF_make_sources.cpp CEXE_sources += ERF_ApplySpongeZoneBCs.cpp +CEXE_sources += ERF_ApplySpongeZoneBCs_ReadFromFile.cpp CEXE_sources += NumericalDiffusion.cpp ifeq ($(USE_NETCDF),TRUE) diff --git a/Source/SourceTerms/Src_headers.H b/Source/SourceTerms/Src_headers.H index 0e1eab1d3..da30c721c 100644 --- a/Source/SourceTerms/Src_headers.H +++ b/Source/SourceTerms/Src_headers.H @@ -27,7 +27,7 @@ void make_sources (int level, int nrk, const amrex::MultiFab& S_prim, amrex::MultiFab& cc_source, #if defined(ERF_USE_RRTMGP) - const amrex::MultiFab& qheating_rates, + const amrex::MultiFab* qheating_rates, #endif const amrex::Geometry geom, const SolverChoice& solverChoice, @@ -57,6 +57,7 @@ void make_mom_sources (int level, int nrk, const amrex::Real* dptr_rhoqt_src, const amrex::Real* dptr_wbar_sub, const amrex::Vector d_rayleigh_ptrs_at_lev, + const amrex::Vector d_sponge_ptrs_at_lev, const int n_qstate); void add_thin_body_sources (amrex::MultiFab& xmom_source, @@ -105,4 +106,15 @@ void ApplySpongeZoneBCsForMom (const SpongeChoice& spongeChoice, const amrex::Array4& rho_u, const amrex::Array4& rho_v, const amrex::Array4& rho_w); + +void ApplySpongeZoneBCsForMom_ReadFromFile (const SpongeChoice& spongeChoice, + const amrex::Geometry geom, + const amrex::Box& tbx, + const amrex::Box& tby, + const amrex::Array4& cell_data, + const amrex::Array4& rho_u_rhs, + const amrex::Array4& rho_v_rhs, + const amrex::Array4& rho_u, + const amrex::Array4& rho_v, + const amrex::Vector d_sponge_ptrs_at_lev); #endif diff --git a/Source/TimeIntegration/ERF_Advance.cpp b/Source/TimeIntegration/ERF_Advance.cpp index d8b901e03..eff2ea0a8 100644 --- a/Source/TimeIntegration/ERF_Advance.cpp +++ b/Source/TimeIntegration/ERF_Advance.cpp @@ -3,6 +3,7 @@ #ifdef ERF_USE_WINDFARM #include +#include #endif using namespace amrex; @@ -74,6 +75,11 @@ ERF::Advance (int lev, Real time, Real dt_lev, int iteration, int /*ncycle*/) fitch_advance(lev, Geom(lev), dt_lev, S_old, U_old, V_old, W_old, vars_fitch[lev], Nturb[lev]); } + else if (solverChoice.windfarm_type == WindFarmType::EWP) { + ewp_advance(lev, Geom(lev), dt_lev, S_old, + U_old, V_old, W_old, vars_ewp[lev], Nturb[lev]); + } + #endif const BoxArray& ba = S_old.boxArray(); @@ -179,4 +185,11 @@ ERF::Advance (int lev, Real time, Real dt_lev, int iteration, int /*ncycle*/) {time, time + dt_lev}); } } + + // *********************************************************************************************** + // Update the time averaged velocities if they are requested + // *********************************************************************************************** + if (solverChoice.time_avg_vel) { + Time_Avg_Vel_atCC(dt[lev], t_avg_cnt[lev], vel_t_avg[lev].get(), U_new, V_new, W_new); + } } diff --git a/Source/TimeIntegration/ERF_advance_dycore.cpp b/Source/TimeIntegration/ERF_advance_dycore.cpp index 8c4736e77..5f434b02d 100644 --- a/Source/TimeIntegration/ERF_advance_dycore.cpp +++ b/Source/TimeIntegration/ERF_advance_dycore.cpp @@ -49,8 +49,9 @@ void ERF::advance_dycore(int level, const Box& domain = fine_geom.Domain(); - DiffChoice dc = solverChoice.diffChoice; - TurbChoice tc = solverChoice.turbChoice[level]; + DiffChoice dc = solverChoice.diffChoice; + TurbChoice tc = solverChoice.turbChoice[level]; + SpongeChoice sc = solverChoice.spongeChoice; MultiFab r_hse (base_state[level], make_alias, 0, 1); // r_0 is first component MultiFab p_hse (base_state[level], make_alias, 1, 1); // p_0 is second component @@ -75,6 +76,14 @@ void ERF::advance_dycore(int level, d_rayleigh_ptrs_at_lev[Rayleigh::wbar] = solverChoice.rayleigh_damp_W ? d_rayleigh_ptrs[level][Rayleigh::wbar].data() : nullptr; d_rayleigh_ptrs_at_lev[Rayleigh::thetabar] = solverChoice.rayleigh_damp_T ? d_rayleigh_ptrs[level][Rayleigh::thetabar].data() : nullptr; + Vector d_sponge_ptrs_at_lev; + if(sc.sponge_type=="input_sponge") + { + d_sponge_ptrs_at_lev.resize(Sponge::nvars_sponge); + d_sponge_ptrs_at_lev[Sponge::ubar_sponge] = d_sponge_ptrs[level][Sponge::ubar_sponge].data(); + d_sponge_ptrs_at_lev[Sponge::vbar_sponge] = d_sponge_ptrs[level][Sponge::vbar_sponge].data(); + } + bool l_use_terrain = solverChoice.use_terrain; bool l_use_diff = ( (dc.molec_diff_type != MolecDiffType::None) || (tc.les_type != LESType::None) || diff --git a/Source/TimeIntegration/ERF_advance_radiation.cpp b/Source/TimeIntegration/ERF_advance_radiation.cpp index 38d071a10..159451e3a 100644 --- a/Source/TimeIntegration/ERF_advance_radiation.cpp +++ b/Source/TimeIntegration/ERF_advance_radiation.cpp @@ -1,4 +1,3 @@ - #include using namespace amrex; @@ -14,7 +13,9 @@ void ERF::advance_radiation (int lev, bool do_snow_opt {true}; bool is_cmip6_volcano {false}; - rad.initialize(cons, qmoist[lev], + rad.initialize(cons, + qheating_rates[lev].get(), + qmoist[lev], grids[lev], Geom(lev), dt_advance, diff --git a/Source/TimeIntegration/TI_slow_rhs_fun.H b/Source/TimeIntegration/TI_slow_rhs_fun.H index b07d4c237..98e730866 100644 --- a/Source/TimeIntegration/TI_slow_rhs_fun.H +++ b/Source/TimeIntegration/TI_slow_rhs_fun.H @@ -43,7 +43,7 @@ // Construct the source terms for the cell-centered (conserved) variables make_sources(level, nrk, slow_dt, S_data, S_prim, cc_src, #if defined(ERF_USE_RRTMGP) - qheating_rates[level], + qheating_rates[level].get(), #endif fine_geom, solverChoice, mapfac_u[level], mapfac_v[level], @@ -117,7 +117,7 @@ r0_new, fine_geom, solverChoice, mapfac_m[level], mapfac_u[level], mapfac_v[level], dptr_u_geos, dptr_v_geos, dptr_wbar_sub, - d_rayleigh_ptrs_at_lev, n_qstate); + d_rayleigh_ptrs_at_lev, d_sponge_ptrs_at_lev, n_qstate); erf_slow_rhs_pre(level, finest_level, nrk, slow_dt, S_rhs, S_data, S_prim, S_scratch, xvel_new, yvel_new, zvel_new, @@ -217,7 +217,7 @@ r0, fine_geom, solverChoice, mapfac_m[level], mapfac_u[level], mapfac_v[level], dptr_u_geos, dptr_v_geos, dptr_wbar_sub, - d_rayleigh_ptrs_at_lev, n_qstate); + d_rayleigh_ptrs_at_lev, d_sponge_ptrs_at_lev, n_qstate); erf_slow_rhs_pre(level, finest_level, nrk, slow_dt, S_rhs, S_data, S_prim, S_scratch, xvel_new, yvel_new, zvel_new, diff --git a/Source/Utils/Make.package b/Source/Utils/Make.package index f0b87ba00..2100badf8 100644 --- a/Source/Utils/Make.package +++ b/Source/Utils/Make.package @@ -1,7 +1,3 @@ -CEXE_sources += MomentumToVelocity.cpp -CEXE_sources += VelocityToMomentum.cpp -CEXE_sources += InteriorGhostCells.cpp - CEXE_headers += TerrainMetrics.H CEXE_headers += Microphysics_Utils.H CEXE_headers += TileNoZ.H @@ -11,7 +7,6 @@ CEXE_headers += Interpolation_WENO.H CEXE_headers += Interpolation_WENO_Z.H CEXE_headers += Interpolation.H CEXE_headers += Interpolation_1D.H -CEXE_sources += TerrainMetrics.cpp CEXE_headers += ParFunctions.H @@ -19,6 +14,12 @@ CEXE_headers += Sat_methods.H CEXE_headers += Water_vapor_saturation.H CEXE_headers += DirectionSelector.H +CEXE_sources += MomentumToVelocity.cpp +CEXE_sources += VelocityToMomentum.cpp +CEXE_sources += InteriorGhostCells.cpp +CEXE_sources += TerrainMetrics.cpp +CEXE_sources += Time_Avg_Vel.cpp + ifeq ($(USE_POISSON_SOLVE),TRUE) CEXE_sources += ERF_PoissonSolve.cpp CEXE_sources += ERF_PoissonSolve_tb.cpp diff --git a/Source/Utils/MomentumToVelocity.cpp b/Source/Utils/MomentumToVelocity.cpp index 281c9567b..9ca10654c 100644 --- a/Source/Utils/MomentumToVelocity.cpp +++ b/Source/Utils/MomentumToVelocity.cpp @@ -86,8 +86,8 @@ MomentumToVelocity (MultiFab& xvel, MultiFab& yvel, MultiFab& zvel, vely(i,j,k) = momy(i,j,k) / dens_arr(i,j-1,k,Rho_comp); }); } - if ( (bx.bigEnd(0) == domain.bigEnd(0)) && - (bc_ptr_h[BCVars::cons_bc].hi(0) == ERFBCType::ext_dir) ) { + if ( (bx.bigEnd(1) == domain.bigEnd(1)) && + (bc_ptr_h[BCVars::cons_bc].hi(1) == ERFBCType::ext_dir) ) { ParallelFor(makeSlab(tby,1,domain.bigEnd(1)+1), [=] AMREX_GPU_DEVICE (int i, int j, int k) { vely(i,j,k) = momy(i,j,k) / dens_arr(i,j,k,Rho_comp); }); diff --git a/Source/Utils/Time_Avg_Vel.cpp b/Source/Utils/Time_Avg_Vel.cpp new file mode 100644 index 000000000..5121ede60 --- /dev/null +++ b/Source/Utils/Time_Avg_Vel.cpp @@ -0,0 +1,47 @@ +#include + +using namespace amrex; + +/* + * Accumulate time averaged velocity fields + */ +void +Time_Avg_Vel_atCC (const Real& dt, + Real& t_avg_cnt, + MultiFab* vel_t_avg, + MultiFab& xvel, + MultiFab& yvel, + MultiFab& zvel) +{ + // Augment the counter + t_avg_cnt += dt; + +#ifdef _OPENMP +#pragma omp parallel if (Gpu::notInLaunchRegion()) +#endif + for ( MFIter mfi(*(vel_t_avg),TilingIfNotGPU()); mfi.isValid(); ++mfi) + { + // CC tilebox + Box tbx = mfi.tilebox(); + + // Velocity on faces + const Array4& velx = xvel.array(mfi); + const Array4& vely = yvel.array(mfi); + const Array4& velz = zvel.array(mfi); + + // Time average at CC + Array4 vel_t_avg_arr = vel_t_avg->array(mfi); + + ParallelFor(tbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) + { + Real u_cc = 0.5 * ( velx(i,j,k) + velx(i+1,j ,k ) ); + Real v_cc = 0.5 * ( vely(i,j,k) + vely(i ,j+1,k ) ); + Real w_cc = 0.5 * ( velz(i,j,k) + velz(i ,j ,k+1) ); + Real umag_cc = std::sqrt(u_cc*u_cc + v_cc*v_cc + w_cc*w_cc); + vel_t_avg_arr(i,j,k,0) += u_cc; + vel_t_avg_arr(i,j,k,1) += v_cc; + vel_t_avg_arr(i,j,k,2) += w_cc; + vel_t_avg_arr(i,j,k,3) += umag_cc; + }); + } +} diff --git a/Source/Utils/Utils.H b/Source/Utils/Utils.H index 4f73531a1..42c25ad1b 100644 --- a/Source/Utils/Utils.H +++ b/Source/Utils/Utils.H @@ -109,6 +109,17 @@ fine_compute_interior_ghost_rhs (const amrex::Real& time, amrex::Vector& S_rhs_f, amrex::Vector& S_data_f); +/* + * Accumulate time averaged velocity fields + */ +void +Time_Avg_Vel_atCC (const amrex::Real& dt, + amrex::Real& t_avg_cnt, + amrex::MultiFab* vel_t_avg, + amrex::MultiFab& xvel, + amrex::MultiFab& yvel, + amrex::MultiFab& zvel); + /** * Zero RHS in the set region * diff --git a/Source/WindFarmParametrization/EWP/AdvanceEWP.cpp b/Source/WindFarmParametrization/EWP/AdvanceEWP.cpp new file mode 100644 index 000000000..13c360363 --- /dev/null +++ b/Source/WindFarmParametrization/EWP/AdvanceEWP.cpp @@ -0,0 +1,117 @@ +#include +#include +using namespace amrex; + +Real R_ewp = 30.0; +Real z_c_ewp = 100.0; +Real C_T_ewp = 0.8, C_TKE_ewp = 0.0; + +void ewp_advance (int lev, + const Geometry& geom, + const Real& dt_advance, + MultiFab& cons_in, + MultiFab& U_old, MultiFab& V_old, MultiFab& W_old, + MultiFab& mf_vars_ewp, const amrex::MultiFab& mf_Nturb) +{ + ewp_source_terms_cellcentered(geom, cons_in, U_old, V_old, W_old, mf_vars_ewp, mf_Nturb); + ewp_update(dt_advance, cons_in, U_old, V_old, mf_vars_ewp); +} + + +void ewp_update (const Real& dt_advance, + MultiFab& cons_in, + MultiFab& U_old, MultiFab& V_old, + const MultiFab& mf_vars_ewp) +{ + + for ( MFIter mfi(cons_in,TilingIfNotGPU()); mfi.isValid(); ++mfi) { + + Box bx = mfi.tilebox(); + Box tbx = mfi.nodaltilebox(0); + Box tby = mfi.nodaltilebox(1); + + auto cons_array = cons_in.array(mfi); + auto ewp_array = mf_vars_ewp.array(mfi); + auto u_vel = U_old.array(mfi); + auto v_vel = V_old.array(mfi); + + ParallelFor(tbx, tby, bx, + [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept + { + u_vel(i,j,k) = u_vel(i,j,k) + (ewp_array(i-1,j,k,0) + ewp_array(i,j,k,0))/2.0*dt_advance; + }, + [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept + { + v_vel(i,j,k) = v_vel(i,j,k) + (ewp_array(i,j-1,k,1) + ewp_array(i,j,k,1))/2.0*dt_advance; + }, + [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept + { + cons_array(i,j,k,RhoQKE_comp) = cons_array(i,j,k,RhoQKE_comp) + ewp_array(i,j,k,2)*dt_advance; + }); + } +} + +void ewp_source_terms_cellcentered (const Geometry& geom, + const MultiFab& cons_in, + const MultiFab& U_old, const MultiFab& V_old, const MultiFab& W_old, + MultiFab& mf_vars_ewp, const amrex::MultiFab& mf_Nturb) +{ + + auto dx = geom.CellSizeArray(); + auto ProbHiArr = geom.ProbHiArray(); + auto ProbLoArr = geom.ProbLoArray(); + + // Domain valid box + 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); + int domhi_z = domain.bigEnd(2) + 1; + + // The order of variables are - Vabs dVabsdt, dudt, dvdt, dTKEdt + mf_vars_ewp.setVal(0.0); + + for ( MFIter mfi(cons_in,TilingIfNotGPU()); mfi.isValid(); ++mfi) { + + const Box& bx = mfi.tilebox(); + const Box& gbx = mfi.growntilebox(1); + auto cons_array = cons_in.array(mfi); + auto ewp_array = mf_vars_ewp.array(mfi); + auto Nturb_array = mf_Nturb.array(mfi); + auto u_vel = U_old.array(mfi); + auto v_vel = V_old.array(mfi); + auto w_vel = W_old.array(mfi); + + amrex::IntVect lo = bx.smallEnd(); + + ParallelFor(gbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + int ii = amrex::min(amrex::max(i, domlo_x), domhi_x); + int jj = amrex::min(amrex::max(j, domlo_y), domhi_y); + int kk = amrex::min(amrex::max(k, domlo_z), domhi_z); + + + Real x = ProbLoArr[0] + (ii+0.5) * dx[0]; + Real y = ProbLoArr[1] + (jj+0.5) * dx[1]; + Real z = ProbLoArr[2] + (kk+0.5) * dx[2]; + + // Compute Fitch source terms + + Real Vabs = std::pow(u_vel(i,j,k)*u_vel(i,j,k) + + v_vel(i,j,k)*v_vel(i,j,k) + + w_vel(i,j,kk)*w_vel(i,j,kk), 0.5); + + Real L_wake = std::pow(dx[0]*dx[1],0.5); + Real K_turb = 100.0; + Real sigma_0 = 60.0; + Real sigma_e = Vabs/(3.0*K_turb*L_wake)*(std::pow(2.0*K_turb*L_wake/Vabs + std::pow(sigma_0,2),3.0/2.0) - std::pow(sigma_0,3)); + + Real phi = std::atan2(v_vel(i,j,k),u_vel(i,j,k)); // Wind direction w.r.t the x-dreiction + Real fac = -std::pow(M_PI/8.0,0.5)*C_T_ewp*std::pow(R_ewp,2)*std::pow(Vabs,2)/(dx[0]*dx[1]*sigma_e)*std::exp(-0.5*std::pow((z - z_c_ewp)/sigma_e,2)); + ewp_array(i,j,k,0) = fac*std::cos(phi)*Nturb_array(i,j,k); + ewp_array(i,j,k,1) = fac*std::sin(phi)*Nturb_array(i,j,k); + ewp_array(i,j,k,2) = 0.0; + }); + } +} diff --git a/Source/WindFarmParametrization/EWP/EWP.H b/Source/WindFarmParametrization/EWP/EWP.H new file mode 100644 index 000000000..7bd5db3cf --- /dev/null +++ b/Source/WindFarmParametrization/EWP/EWP.H @@ -0,0 +1,25 @@ +#ifndef EWP_H +#define EWP_H + +#include +#include + +void ewp_advance (int lev, + const amrex::Geometry& geom, + const amrex::Real& dt_advance, + amrex::MultiFab& cons_in, + amrex::MultiFab& U_old, amrex::MultiFab& V_old, amrex::MultiFab& W_old, + amrex::MultiFab& mf_vars_ewp, const amrex::MultiFab& mf_Nturb); + +void ewp_source_terms_cellcentered (const amrex::Geometry& geom, + const amrex::MultiFab& cons_in, + const amrex::MultiFab& U_old, const amrex::MultiFab& V_old, const amrex::MultiFab& W_old, + amrex::MultiFab& mf_vars_ewp, const amrex::MultiFab& mf_Nturb); + +void ewp_update (const amrex::Real& dt_advance, + amrex::MultiFab& cons_in, + amrex::MultiFab& U_old, amrex::MultiFab& V_old, + const amrex::MultiFab& mf_vars_ewp); + +#endif + diff --git a/Source/WindFarmParametrization/EWP/Make.package b/Source/WindFarmParametrization/EWP/Make.package new file mode 100644 index 000000000..56789f3d1 --- /dev/null +++ b/Source/WindFarmParametrization/EWP/Make.package @@ -0,0 +1,2 @@ +CEXE_sources += AdvanceEWP.cpp +CEXE_headers += EWP.H