From 0062ec9c32e78c99004d0614f5a7511540c6c6af Mon Sep 17 00:00:00 2001 From: Mahesh Natarajan Date: Mon, 22 Apr 2024 17:25:59 -0700 Subject: [PATCH] WIP: Explicit Wake Parametrization model for wind farm parametrization (#1583) * WIP:First draft of EWP * WIP: Updates to EWP model * WIP: Updates to EWP model * WIP: Removing unwanted file * Updating Make.ERF * Correcting Fitch input file * Removing unwanted file --------- Co-authored-by: Mahesh Natarajan Co-authored-by: Aaron M. Lattanzi <103702284+AMLattanzi@users.noreply.github.com> --- Exec/EWP/CMakeLists.txt | 12 ++ Exec/EWP/GNUmakefile | 37 ++++ Exec/EWP/Make.package | 2 + Exec/EWP/README | 6 + Exec/EWP/inputs_EWP_WindFarm | 78 ++++++++ Exec/EWP/plot_profiles.py | 70 ++++++++ Exec/EWP/prob.H | 79 +++++++++ Exec/EWP/prob.cpp | 167 ++++++++++++++++++ Exec/EWP/wind-turbine-1.tbl | 24 +++ Exec/EWP/windturbines.txt | 97 ++++++++++ Exec/EWP/windturbines.txt_LargeDomain | 30 ++++ Exec/EWP/windturbines_WindFarm.txt | 97 ++++++++++ Exec/Fitch/inputs_Fitch_WindFarm | 78 ++++++++ Exec/Fitch/windturbines_WindFarm.txt | 97 ++++++++++ Exec/Make.ERF | 4 + Source/DataStructs/DataStruct.H | 7 +- Source/ERF.H | 1 + Source/ERF.cpp | 14 +- Source/ERF_make_new_arrays.cpp | 5 + Source/IO/Checkpoint.cpp | 4 +- Source/TimeIntegration/ERF_Advance.cpp | 6 + .../EWP/AdvanceEWP.cpp | 117 ++++++++++++ Source/WindFarmParametrization/EWP/EWP.H | 25 +++ .../WindFarmParametrization/EWP/Make.package | 2 + 24 files changed, 1047 insertions(+), 12 deletions(-) create mode 100644 Exec/EWP/CMakeLists.txt create mode 100644 Exec/EWP/GNUmakefile create mode 100644 Exec/EWP/Make.package create mode 100644 Exec/EWP/README create mode 100644 Exec/EWP/inputs_EWP_WindFarm create mode 100644 Exec/EWP/plot_profiles.py create mode 100644 Exec/EWP/prob.H create mode 100644 Exec/EWP/prob.cpp create mode 100644 Exec/EWP/wind-turbine-1.tbl create mode 100644 Exec/EWP/windturbines.txt create mode 100644 Exec/EWP/windturbines.txt_LargeDomain create mode 100644 Exec/EWP/windturbines_WindFarm.txt create mode 100644 Exec/Fitch/inputs_Fitch_WindFarm create mode 100644 Exec/Fitch/windturbines_WindFarm.txt create mode 100644 Source/WindFarmParametrization/EWP/AdvanceEWP.cpp create mode 100644 Source/WindFarmParametrization/EWP/EWP.H create mode 100644 Source/WindFarmParametrization/EWP/Make.package 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/Source/DataStructs/DataStruct.H b/Source/DataStructs/DataStruct.H index df38dbb62..daa3d8724 100644 --- a/Source/DataStructs/DataStruct.H +++ b/Source/DataStructs/DataStruct.H @@ -37,7 +37,7 @@ enum struct MoistureType { }; enum struct WindFarmType { - Fitch + Fitch, EWP }; enum struct LandSurfaceType { @@ -275,8 +275,11 @@ struct SolverChoice { if (windfarm_type_string == "Fitch") { windfarm_type = WindFarmType::Fitch; } + else if(windfarm_type_string == "EWP"){ + windfarm_type = WindFarmType::EWP; + } else { - amrex::Abort("Dont know this windfarm_type. Currently the only option is Fitch."); + amrex::Abort("Dont know this windfarm_type. Currently only Fitch and EWP models are supported."); } // Test if time averaged data is to be output diff --git a/Source/ERF.H b/Source/ERF.H index e080ef666..a0effbcd1 100644 --- a/Source/ERF.H +++ b/Source/ERF.H @@ -615,6 +615,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; diff --git a/Source/ERF.cpp b/Source/ERF.cpp index e31ec67df..b7c88af7b 100644 --- a/Source/ERF.cpp +++ b/Source/ERF.cpp @@ -112,10 +112,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) @@ -1755,10 +1754,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) diff --git a/Source/ERF_make_new_arrays.cpp b/Source/ERF_make_new_arrays.cpp index 4b3f1dca7..1f2036bd7 100644 --- a/Source/ERF_make_new_arrays.cpp +++ b/Source/ERF_make_new_arrays.cpp @@ -235,6 +235,11 @@ 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) 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/TimeIntegration/ERF_Advance.cpp b/Source/TimeIntegration/ERF_Advance.cpp index cc2f4694d..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(); 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