diff --git a/Exec/WindFarmTests/SingleTurbine/SimpleActuatorDisk/CMakeLists.txt b/Exec/WindFarmTests/SingleTurbine/SimpleActuatorDisk/NoTerrain/CMakeLists.txt similarity index 87% rename from Exec/WindFarmTests/SingleTurbine/SimpleActuatorDisk/CMakeLists.txt rename to Exec/WindFarmTests/SingleTurbine/SimpleActuatorDisk/NoTerrain/CMakeLists.txt index d6967632f..5531ae3c8 100644 --- a/Exec/WindFarmTests/SingleTurbine/SimpleActuatorDisk/CMakeLists.txt +++ b/Exec/WindFarmTests/SingleTurbine/SimpleActuatorDisk/NoTerrain/CMakeLists.txt @@ -1,4 +1,4 @@ -set(erf_exe_name erf_fitch) +set(erf_exe_name erf_sad_no_terrain) add_executable(${erf_exe_name} "") target_sources(${erf_exe_name} diff --git a/Exec/WindFarmTests/SingleTurbine/SimpleActuatorDisk/ERF_prob.H b/Exec/WindFarmTests/SingleTurbine/SimpleActuatorDisk/NoTerrain/ERF_prob.H similarity index 100% rename from Exec/WindFarmTests/SingleTurbine/SimpleActuatorDisk/ERF_prob.H rename to Exec/WindFarmTests/SingleTurbine/SimpleActuatorDisk/NoTerrain/ERF_prob.H diff --git a/Exec/WindFarmTests/SingleTurbine/SimpleActuatorDisk/ERF_prob.cpp b/Exec/WindFarmTests/SingleTurbine/SimpleActuatorDisk/NoTerrain/ERF_prob.cpp similarity index 100% rename from Exec/WindFarmTests/SingleTurbine/SimpleActuatorDisk/ERF_prob.cpp rename to Exec/WindFarmTests/SingleTurbine/SimpleActuatorDisk/NoTerrain/ERF_prob.cpp diff --git a/Exec/WindFarmTests/SingleTurbine/SimpleActuatorDisk/GNUmakefile b/Exec/WindFarmTests/SingleTurbine/SimpleActuatorDisk/NoTerrain/GNUmakefile similarity index 90% rename from Exec/WindFarmTests/SingleTurbine/SimpleActuatorDisk/GNUmakefile rename to Exec/WindFarmTests/SingleTurbine/SimpleActuatorDisk/NoTerrain/GNUmakefile index 0e9dba57b..ac6ef05f2 100644 --- a/Exec/WindFarmTests/SingleTurbine/SimpleActuatorDisk/GNUmakefile +++ b/Exec/WindFarmTests/SingleTurbine/SimpleActuatorDisk/NoTerrain/GNUmakefile @@ -29,6 +29,6 @@ USE_WINDFARM = TRUE # GNU Make Bpack := ./Make.package Blocs := . -ERF_HOME := ../../../.. -ERF_PROBLEM_DIR = $(ERF_HOME)/Exec/WindFarmTests/SingleTurbine/SimpleActuatorDisk +ERF_HOME := ../../../../.. +ERF_PROBLEM_DIR = $(ERF_HOME)/Exec/WindFarmTests/SingleTurbine/SimpleActuatorDisk/NoTerrain include $(ERF_HOME)/Exec/Make.ERF diff --git a/Exec/WindFarmTests/SingleTurbine/SimpleActuatorDisk/Make.package b/Exec/WindFarmTests/SingleTurbine/SimpleActuatorDisk/NoTerrain/Make.package similarity index 100% rename from Exec/WindFarmTests/SingleTurbine/SimpleActuatorDisk/Make.package rename to Exec/WindFarmTests/SingleTurbine/SimpleActuatorDisk/NoTerrain/Make.package diff --git a/Exec/WindFarmTests/SingleTurbine/SimpleActuatorDisk/README b/Exec/WindFarmTests/SingleTurbine/SimpleActuatorDisk/NoTerrain/README similarity index 100% rename from Exec/WindFarmTests/SingleTurbine/SimpleActuatorDisk/README rename to Exec/WindFarmTests/SingleTurbine/SimpleActuatorDisk/NoTerrain/README diff --git a/Exec/WindFarmTests/SingleTurbine/SimpleActuatorDisk/inputs_1WT_x_y b/Exec/WindFarmTests/SingleTurbine/SimpleActuatorDisk/NoTerrain/inputs_1WT_x_y similarity index 100% rename from Exec/WindFarmTests/SingleTurbine/SimpleActuatorDisk/inputs_1WT_x_y rename to Exec/WindFarmTests/SingleTurbine/SimpleActuatorDisk/NoTerrain/inputs_1WT_x_y diff --git a/Exec/WindFarmTests/SingleTurbine/SimpleActuatorDisk/windturbines_loc_x_y_1WT.txt b/Exec/WindFarmTests/SingleTurbine/SimpleActuatorDisk/NoTerrain/windturbines_loc_x_y_1WT.txt similarity index 100% rename from Exec/WindFarmTests/SingleTurbine/SimpleActuatorDisk/windturbines_loc_x_y_1WT.txt rename to Exec/WindFarmTests/SingleTurbine/SimpleActuatorDisk/NoTerrain/windturbines_loc_x_y_1WT.txt diff --git a/Exec/WindFarmTests/SingleTurbine/SimpleActuatorDisk/windturbines_spec_1WT.tbl b/Exec/WindFarmTests/SingleTurbine/SimpleActuatorDisk/NoTerrain/windturbines_spec_1WT.tbl similarity index 100% rename from Exec/WindFarmTests/SingleTurbine/SimpleActuatorDisk/windturbines_spec_1WT.tbl rename to Exec/WindFarmTests/SingleTurbine/SimpleActuatorDisk/NoTerrain/windturbines_spec_1WT.tbl diff --git a/Exec/WindFarmTests/SingleTurbine/SimpleActuatorDisk/WithTerrain/CMakeLists.txt b/Exec/WindFarmTests/SingleTurbine/SimpleActuatorDisk/WithTerrain/CMakeLists.txt new file mode 100644 index 000000000..5dba1aa31 --- /dev/null +++ b/Exec/WindFarmTests/SingleTurbine/SimpleActuatorDisk/WithTerrain/CMakeLists.txt @@ -0,0 +1,12 @@ +set(erf_exe_name erf_sad_with_terrain) + +add_executable(${erf_exe_name} "") +target_sources(${erf_exe_name} + PRIVATE + ERF_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/WindFarmTests/SingleTurbine/SimpleActuatorDisk/WithTerrain/ERF_prob.H b/Exec/WindFarmTests/SingleTurbine/SimpleActuatorDisk/WithTerrain/ERF_prob.H new file mode 100644 index 000000000..b20ea6b88 --- /dev/null +++ b/Exec/WindFarmTests/SingleTurbine/SimpleActuatorDisk/WithTerrain/ERF_prob.H @@ -0,0 +1,80 @@ +#ifndef ERF_PROB_H_ +#define ERF_PROB_H_ + +#include + +#include "AMReX_REAL.H" + +#include "ERF_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 KE_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; + + // 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/ERF_init_constant_density_hse.H" +#include "Prob/ERF_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; + + void init_custom_terrain ( + const amrex::Geometry& geom, + amrex::MultiFab& z_phys_nd, + const amrex::Real& time) override; + +protected: + std::string name() override { return "ABL"; } + +private: + ProbParm parms; +}; + +#endif diff --git a/Exec/WindFarmTests/SingleTurbine/SimpleActuatorDisk/WithTerrain/ERF_prob.cpp b/Exec/WindFarmTests/SingleTurbine/SimpleActuatorDisk/WithTerrain/ERF_prob.cpp new file mode 100644 index 000000000..4df40c464 --- /dev/null +++ b/Exec/WindFarmTests/SingleTurbine/SimpleActuatorDisk/WithTerrain/ERF_prob.cpp @@ -0,0 +1,238 @@ +#include "ERF_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("KE_0", parms.KE_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; + + 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_d=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_d.pert_ref_height) && (parms_d.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_d.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_d.A_0 * exp(-10.*r*r); + + // Set an initial value for KE + state_pert(i, j, k, RhoKE_comp) = parms_d.KE_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_d=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_d.U_0; + if ((z <= parms_d.pert_ref_height) && (parms_d.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_d.U_0_Pert_Mag; + x_vel_pert(i, j, k) += x_vel_prime; + } + if (parms_d.pert_deltaU != 0.0) + { + const amrex::Real yl = y - prob_lo[1]; + const amrex::Real zl = z / parms_d.pert_ref_height; + const amrex::Real damp = std::exp(-0.5 * zl * zl); + x_vel_pert(i, j, k) += parms_d.ufac * damp * z * std::cos(parms_d.aval * yl); + } + }); + + // Set the y-velocity + ParallelForRNG(ybx, [=, parms_d=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_d.V_0; + if ((z <= parms_d.pert_ref_height) && (parms_d.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_d.V_0_Pert_Mag; + y_vel_pert(i, j, k) += y_vel_prime; + } + if (parms_d.pert_deltaV != 0.0) + { + const amrex::Real xl = x - prob_lo[0]; + const amrex::Real zl = z / parms_d.pert_ref_height; + const amrex::Real damp = std::exp(-0.5 * zl * zl); + y_vel_pert(i, j, k) += parms_d.vfac * damp * z * std::cos(parms_d.bval * xl); + } + }); + + // Set the z-velocity + ParallelForRNG(zbx, [=, parms_d=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_d.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_d.W_0_Pert_Mag; + z_vel_pert(i, j, k) = parms_d.W_0 + z_vel_prime; + } + }); +} + +void +Problem::init_custom_terrain ( + const Geometry& geom, + MultiFab& z_phys_nd, + const Real& time) +{ + + // Check if a valid text file exists for the terrain + std::string fname; + ParmParse pp("erf"); + auto valid_fname = pp.query("terrain_file_name",fname); + if (valid_fname) { + this->read_custom_terrain(fname,geom,z_phys_nd,time); + } else { + // Domain cell size and real bounds + auto dx = geom.CellSizeArray(); + auto ProbLoArr = geom.ProbLoArray(); + auto ProbHiArr = geom.ProbHiArray(); + + // Domain valid box (z_nd is nodal) + const amrex::Box& domain = geom.Domain(); + int domlo_x = domain.smallEnd(0); int domhi_x = domain.bigEnd(0) + 1; + int domlo_y = domain.smallEnd(1); int domhi_y = domain.bigEnd(1) + 1; + int domlo_z = domain.smallEnd(2); + + // User function parameters + Real a = 0.5; + Real num = 8 * a * a * a; + Real xcen = 500.0; + Real ycen = 500.0; + + // if hm is nonzero, then use alternate hill definition + + //Real hm = parms.hmax; + //Real L = parms.L; + + // 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 ) + { + amrex::Box zbx = mfi.nodaltilebox(2); + if (zbx.smallEnd(2) > k0) continue; + + // Grown box with no z range + amrex::Box xybx = mfi.growntilebox(ngrow); + xybx.setRange(2,0); + + amrex::Array4 const& z_arr = z_phys_nd.array(mfi); + + ParallelFor(xybx, [=] AMREX_GPU_DEVICE (int i, int j, int) + { + // Clip indices for ghost-cells + int ii = amrex::min(amrex::max(i,domlo_x),domhi_x); + int jj = amrex::min(amrex::max(j,domlo_y),domhi_y); + + // Location of nodes + Real x = ProbLoArr[0] + ii * dx[0] - xcen; + Real y = ProbLoArr[1] + jj * dx[1] - ycen; + // Real y = (jj * dx[1] - ycen); + + Real x_L = x/100.0; + Real y_L = y/100.0; + z_arr(i,j,k0) = 100.0 / (1.0 + x_L*x_L + y_L*y_L); + }); + } + } +} + + diff --git a/Exec/WindFarmTests/SingleTurbine/SimpleActuatorDisk/WithTerrain/GNUmakefile b/Exec/WindFarmTests/SingleTurbine/SimpleActuatorDisk/WithTerrain/GNUmakefile new file mode 100644 index 000000000..5536699e0 --- /dev/null +++ b/Exec/WindFarmTests/SingleTurbine/SimpleActuatorDisk/WithTerrain/GNUmakefile @@ -0,0 +1,34 @@ +# 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 + +# GNU Make +Bpack := ./Make.package +Blocs := . +ERF_HOME := ../../../../.. +ERF_PROBLEM_DIR = $(ERF_HOME)/Exec/WindFarmTests/SingleTurbine/SimpleActuatorDisk/WithTerrain +include $(ERF_HOME)/Exec/Make.ERF diff --git a/Exec/WindFarmTests/SingleTurbine/SimpleActuatorDisk/WithTerrain/Make.package b/Exec/WindFarmTests/SingleTurbine/SimpleActuatorDisk/WithTerrain/Make.package new file mode 100644 index 000000000..5fc21f61c --- /dev/null +++ b/Exec/WindFarmTests/SingleTurbine/SimpleActuatorDisk/WithTerrain/Make.package @@ -0,0 +1,2 @@ +CEXE_headers += ERF_prob.H +CEXE_sources += ERF_prob.cpp diff --git a/Exec/WindFarmTests/SingleTurbine/SimpleActuatorDisk/WithTerrain/README b/Exec/WindFarmTests/SingleTurbine/SimpleActuatorDisk/WithTerrain/README new file mode 100644 index 000000000..b8e994ac1 --- /dev/null +++ b/Exec/WindFarmTests/SingleTurbine/SimpleActuatorDisk/WithTerrain/README @@ -0,0 +1,2 @@ +This folder contains an example of the Simple actuator disk wind farm parametetrization +for a single turbine case. diff --git a/Exec/WindFarmTests/SingleTurbine/SimpleActuatorDisk/inputs_1WT_lat_lon b/Exec/WindFarmTests/SingleTurbine/SimpleActuatorDisk/WithTerrain/inputs_1WT_x_y similarity index 58% rename from Exec/WindFarmTests/SingleTurbine/SimpleActuatorDisk/inputs_1WT_lat_lon rename to Exec/WindFarmTests/SingleTurbine/SimpleActuatorDisk/WithTerrain/inputs_1WT_x_y index 4781b7b8a..2e95a59a4 100644 --- a/Exec/WindFarmTests/SingleTurbine/SimpleActuatorDisk/inputs_1WT_lat_lon +++ b/Exec/WindFarmTests/SingleTurbine/SimpleActuatorDisk/WithTerrain/inputs_1WT_x_y @@ -1,20 +1,26 @@ # ------------------ INPUTS TO MAIN PROGRAM ------------------- -max_step = 3000 +max_step = 1000 amrex.fpe_trap_invalid = 1 fabarray.mfiter_tile_size = 1024 1024 1024 # PROBLEM SIZE & GEOMETRY -geometry.prob_extent = 100000.0 100000.0 1000 -amr.n_cell = 50 50 80 +geometry.prob_extent = 1000.0 1000.0 500.0 +amr.n_cell = 100 100 50 # WINDFARM PARAMETRIZATION PARAMETERS -erf.windfarm_type = "SimpleActuatorDisk" -erf.windfarm_loc_type = "lat_lon" -erf.latitude_lo = 35.0 -erf.longitude_lo = -100.0 -erf.windfarm_loc_table = "windturbines_loc_lat_lon_1WT.txt" +erf.use_terrain = 1 +erf.terrain_type = Static +erf.terrain_z_levels = 0.000000 5.000000 11.040962 17.202744 23.487761 29.898478 36.437410 43.107120 49.910224 56.849391 63.927341 71.146850 78.510749 86.021926 93.683327 101.497956 109.468877 117.599217 125.892163 134.350968 142.978950 151.779491 160.756043 169.912126 179.251331 188.777319 198.493828 208.404667 218.513722 228.824959 239.342420 250.070231 261.012597 272.173811 283.558250 295.170377 307.014747 319.096004 331.418886 343.988226 356.808952 369.886094 383.224778 396.830235 410.707802 424.862920 439.301141 454.028126 469.049651 484.371606 500.000000 + +erf.windfarm_type = "SimpleAD" +erf.windfarm_loc_type = "x_y" +erf.sampling_distance_by_D = 0.5 +erf.turb_disk_angle_from_x = 90.0 +erf.windfarm_x_shift = 200.0 +erf.windfarm_y_shift = 200.0 +erf.windfarm_loc_table = "windturbines_loc_x_y_1WT.txt" erf.windfarm_spec_table = "windturbines_spec_1WT.tbl" #erf.grid_stretching_ratio = 1.025 @@ -51,7 +57,8 @@ xlo.theta = 300. # TIME STEP CONTROL -erf.fixed_dt = 3.0 # fixed time step depending on grid resolution +erf.use_native_mri = 1 +erf.fixed_dt = 0.1 # fixed time step depending on grid resolution #erf.fixed_fast_dt = 0.0025 # DIAGNOSTICS & VERBOSITY @@ -65,20 +72,20 @@ amr.max_level = 0 # maximum level number allowed # CHECKPOINT FILES erf.check_file = chk # root name of checkpoint file erf.check_int = 1000 # number of timesteps between checkpoints -#erf.restart = chk02000 +#erf.restart = chk01000 # 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 rhoadv_0 x_velocity y_velocity z_velocity pressure temp theta KE num_turb vorticity +erf.plot_vars_1 = density rhoadv_0 x_velocity y_velocity z_velocity pressure temp theta QKE num_turb SMark0 SMark1 vorticity_x vorticity_y vorticity_z # ADVECTION SCHEMES -erf.dycore_horiz_adv_type = "Centered_2nd" -erf.dycore_vert_adv_type = "Centered_2nd" -erf.dryscal_horiz_adv_type = "Centered_2nd" -erf.dryscal_vert_adv_type = "Centered_2nd" -erf.moistscal_horiz_adv_type = "Centered_2nd" -erf.moistscal_vert_adv_type = "Centered_2nd" +#erf.dycore_horiz_adv_type = "Centered_2nd" +#erf.dycore_vert_adv_type = "Centered_2nd" +#erf.dryscal_horiz_adv_type = "Centered_2nd" +#erf.dryscal_vert_adv_type = "Centered_2nd" +#erf.moistscal_horiz_adv_type = "Centered_2nd" +#erf.moistscal_vert_adv_type = "Centered_2nd" # SOLVER CHOICE erf.alpha_T = 0.0 diff --git a/Exec/WindFarmTests/SingleTurbine/SimpleActuatorDisk/WithTerrain/windturbines_loc_x_y_1WT.txt b/Exec/WindFarmTests/SingleTurbine/SimpleActuatorDisk/WithTerrain/windturbines_loc_x_y_1WT.txt new file mode 100644 index 000000000..792d5b75b --- /dev/null +++ b/Exec/WindFarmTests/SingleTurbine/SimpleActuatorDisk/WithTerrain/windturbines_loc_x_y_1WT.txt @@ -0,0 +1,3 @@ +300.0 500.0 +500.0 500.0 +700.0 500.0 diff --git a/Exec/WindFarmTests/SingleTurbine/SimpleActuatorDisk/WithTerrain/windturbines_spec_1WT.tbl b/Exec/WindFarmTests/SingleTurbine/SimpleActuatorDisk/WithTerrain/windturbines_spec_1WT.tbl new file mode 100644 index 000000000..e1abbf4cf --- /dev/null +++ b/Exec/WindFarmTests/SingleTurbine/SimpleActuatorDisk/WithTerrain/windturbines_spec_1WT.tbl @@ -0,0 +1,6 @@ +4 +119.0 178.0 0.130 2.0 +9 0.805 50.0 +10 0.805 50.0 +11 0.805 50.0 +12 0.805 50.0 diff --git a/Exec/WindFarmTests/SingleTurbine/SimpleActuatorDisk/inputs_WindFarm_lat_lon b/Exec/WindFarmTests/SingleTurbine/SimpleActuatorDisk/inputs_WindFarm_lat_lon deleted file mode 100644 index 8c6e9c77d..000000000 --- a/Exec/WindFarmTests/SingleTurbine/SimpleActuatorDisk/inputs_WindFarm_lat_lon +++ /dev/null @@ -1,84 +0,0 @@ -# ------------------ INPUTS TO MAIN PROGRAM ------------------- -max_step = 100000 - -amrex.fpe_trap_invalid = 1 - -fabarray.mfiter_tile_size = 1024 1024 1024 - -# PROBLEM SIZE & GEOMETRY -geometry.prob_extent = 200150 202637 1000 -amr.n_cell = 50 50 40 - -# WINDFARM PARAMETRIZATION PARAMETERS -erf.windfarm_type = "Fitch" -erf.windfarm_loc_type = "lat_lon" -erf.latitude_lo = 35.0 -erf.longitude_lo = -100.0 -erf.windfarm_loc_table = "windturbines_loc_lat_lon_WindFarm.txt" -erf.windfarm_spec_table = "windturbines_spec_WindFarm.tbl" - - -#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.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 KE 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/WindFarmTests/SingleTurbine/SimpleActuatorDisk/inputs_WindFarm_x_y b/Exec/WindFarmTests/SingleTurbine/SimpleActuatorDisk/inputs_WindFarm_x_y deleted file mode 100644 index 1733cfe54..000000000 --- a/Exec/WindFarmTests/SingleTurbine/SimpleActuatorDisk/inputs_WindFarm_x_y +++ /dev/null @@ -1,82 +0,0 @@ -# ------------------ INPUTS TO MAIN PROGRAM ------------------- -max_step = 100000 - -amrex.fpe_trap_invalid = 1 - -fabarray.mfiter_tile_size = 1024 1024 1024 - -# PROBLEM SIZE & GEOMETRY -geometry.prob_extent = 200150 202637 1000 -amr.n_cell = 50 50 40 - -# WINDFARM PARAMETRIZATION PARAMETERS -erf.windfarm_type = "Fitch" -erf.windfarm_loc_type = "x_y" -erf.windfarm_loc_table = "windturbines_loc_x_y_WindFarm.txt" -erf.windfarm_spec_table = "windturbines_spec_WindFarm.tbl" - - -#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.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 KE 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/WindFarmTests/SingleTurbine/SimpleActuatorDisk/windturbines_loc_lat_lon_1WT.txt b/Exec/WindFarmTests/SingleTurbine/SimpleActuatorDisk/windturbines_loc_lat_lon_1WT.txt deleted file mode 100644 index b78b36903..000000000 --- a/Exec/WindFarmTests/SingleTurbine/SimpleActuatorDisk/windturbines_loc_lat_lon_1WT.txt +++ /dev/null @@ -1,2 +0,0 @@ -35.45 -99.5 1 - diff --git a/Source/ERF_make_new_arrays.cpp b/Source/ERF_make_new_arrays.cpp index 6a9a249c3..9ff78921e 100644 --- a/Source/ERF_make_new_arrays.cpp +++ b/Source/ERF_make_new_arrays.cpp @@ -251,7 +251,7 @@ ERF::init_stuff (int lev, const BoxArray& ba, const DistributionMapping& dm, vars_windfarm[lev].define(ba, dm, 3, ngrow_state);// dudt, dvdt, dwdt } Nturb[lev].define(ba, dm, 1, ngrow_state); // Number of turbines in a cell - SMark[lev].define(ba, dm, 2, ngrow_state); // Free stream velocity/source term + SMark[lev].define(ba, dm, 2, 1); // Free stream velocity/source term // sampling marker in a cell - 2 components #endif @@ -355,7 +355,7 @@ ERF::init_stuff (int lev, const BoxArray& ba, const DistributionMapping& dm, // write out the vtk files for wind turbine location and/or // actuator disks #ifdef ERF_USE_WINDFARM - init_windfarm(lev); + //init_windfarm(lev); #endif } diff --git a/Source/ERF_make_new_level.cpp b/Source/ERF_make_new_level.cpp index c7c529e3c..5ce7a74eb 100644 --- a/Source/ERF_make_new_level.cpp +++ b/Source/ERF_make_new_level.cpp @@ -131,6 +131,13 @@ void ERF::MakeNewLevelFromScratch (int lev, Real time, const BoxArray& ba_in, } } + // Read in tables needed for windfarm simulations + // fill in Nturb multifab - number of turbines in each mesh cell + // write out the vtk files for wind turbine location and/or + // actuator disks + #ifdef ERF_USE_WINDFARM + init_windfarm(lev); + #endif // ******************************************************************************************** // Build the data structures for canopy model (depends upon z_phys) // ******************************************************************************************** diff --git a/Source/Initialization/ERF_init_windfarm.cpp b/Source/Initialization/ERF_init_windfarm.cpp index b87ededd4..bdaae7bdf 100644 --- a/Source/Initialization/ERF_init_windfarm.cpp +++ b/Source/Initialization/ERF_init_windfarm.cpp @@ -29,8 +29,7 @@ ERF::init_windfarm (int lev) true, false); } - - windfarm->fill_Nturb_multifab(geom[lev], Nturb[lev]); + windfarm->fill_Nturb_multifab(geom[lev], Nturb[lev], z_phys_cc[lev]); windfarm->write_turbine_locations_vtk(); @@ -38,8 +37,10 @@ ERF::init_windfarm (int lev) solverChoice.windfarm_type == WindFarmType::GeneralAD) { windfarm->fill_SMark_multifab(geom[lev], SMark[lev], solverChoice.sampling_distance_by_D, - solverChoice.turb_disk_angle); - windfarm->write_actuator_disks_vtk(geom[lev]); + solverChoice.turb_disk_angle, + z_phys_cc[lev]); + windfarm->write_actuator_disks_vtk(geom[lev], + solverChoice.sampling_distance_by_D); } if(solverChoice.windfarm_type == WindFarmType::GeneralAD) { diff --git a/Source/WindFarmParametrization/ERF_InitWindFarm.cpp b/Source/WindFarmParametrization/ERF_InitWindFarm.cpp index 933e1ada8..482a39379 100644 --- a/Source/WindFarmParametrization/ERF_InitWindFarm.cpp +++ b/Source/WindFarmParametrization/ERF_InitWindFarm.cpp @@ -5,6 +5,7 @@ #include #include #include // For POSIX directory handling +#include // For std::sort using namespace amrex; @@ -349,18 +350,104 @@ WindFarm::read_windfarm_airfoil_tables (const std::string windfarm_airfoil_table set_blade_airfoil_spec(bld_airfoil_aoa, bld_airfoil_Cl, bld_airfoil_Cd); } +void +WindFarm::gatherKeyValuePairs(const std::vector>& localData, + std::vector>& globalData) +{ + int myRank = amrex::ParallelDescriptor::MyProc(); + int nProcs = amrex::ParallelDescriptor::NProcs(); + + int localSize = localData.size(); + + // Prepare separate vectors for keys and values + std::vector localKeys(localSize); + std::vector localValues(localSize); + + for (size_t i = 0; i < localData.size(); ++i) { + localKeys[i] = localData[i].first; + localValues[i] = localData[i].second; + } + + // Gather sizes of the arrays from all processes + std::vector allSizes(nProcs); + amrex::ParallelDescriptor::Gather(&localSize, 1, allSizes.data(), 1, 0); + + // Compute displacements for global data + std::vector displs(nProcs, 0); + int totalSize = 0; + if (amrex::ParallelDescriptor::MyProc() == 0) { + for (int i = 0; i < nProcs; ++i) { + displs[i] = totalSize; + totalSize += allSizes[i]; + } + globalData.resize(totalSize); + } + + // Gather keys + std::vector globalKeys(totalSize); + amrex::ParallelDescriptor::Gatherv(localKeys.data(), localSize, + globalKeys.data(), allSizes, + displs, 0); + + // Gather values + std::vector globalValues(totalSize); + amrex::ParallelDescriptor::Gatherv(localValues.data(), localSize, + globalValues.data(), allSizes, + displs, 0); + + // Rank 0 combines keys and values into globalData + if (myRank == 0) { + globalData.clear(); + for (int i = 0; i < totalSize; ++i) { + globalData.emplace_back(globalKeys[i], globalValues[i]); + } + + // Sort global data by keys + std::sort(globalData.begin(), globalData.end()); + } + + // Broadcast the global data to all processes + int globalCount = globalData.size(); + amrex::ParallelDescriptor::Bcast(&globalCount, 1, 0); + + if (myRank != 0) { + globalData.resize(globalCount); + } + + // Broadcast the actual pairs + for (int i = 0; i < globalCount; ++i) { + amrex::ParallelDescriptor::Bcast(&globalData[i].first, 1, 0); + amrex::ParallelDescriptor::Bcast(&globalData[i].second, 1, 0); + } + + for (const auto& kv : globalData) { + std::cout << "Rank " << myRank << "Key: " << kv.first << ", Value: " << kv.second << std::endl; + } +} + + void WindFarm::fill_Nturb_multifab (const Geometry& geom, - MultiFab& mf_Nturb) + MultiFab& mf_Nturb, + std::unique_ptr& z_phys_cc) { + zloc.resize(xloc.size(),-1.0); + turb_index.resize(xloc.size(),-1); + amrex::Gpu::DeviceVector d_xloc(xloc.size()); amrex::Gpu::DeviceVector d_yloc(yloc.size()); + amrex::Gpu::DeviceVector d_zloc(xloc.size()); + amrex::Gpu::DeviceVector d_turb_index(xloc.size()); amrex::Gpu::copyAsync(amrex::Gpu::hostToDevice, xloc.begin(), xloc.end(), d_xloc.begin()); amrex::Gpu::copyAsync(amrex::Gpu::hostToDevice, yloc.begin(), yloc.end(), d_yloc.begin()); + amrex::Gpu::copyAsync(amrex::Gpu::hostToDevice, zloc.begin(), zloc.end(), d_zloc.begin()); + amrex::Gpu::copyAsync(amrex::Gpu::hostToDevice, turb_index.begin(), turb_index.end(), d_turb_index.begin()); - Real* d_xloc_ptr = d_xloc.data(); - Real* d_yloc_ptr = d_yloc.data(); + Real* d_xloc_ptr = d_xloc.data(); + Real* d_yloc_ptr = d_yloc.data(); + Real* d_zloc_ptr = d_zloc.data(); + int* d_turb_index_ptr = d_turb_index.data(); mf_Nturb.setVal(0); @@ -374,10 +461,14 @@ WindFarm::fill_Nturb_multifab (const Geometry& geom, auto ProbLoArr = geom.ProbLoArray(); int num_turb = xloc.size(); + bool is_terrain = z_phys_cc ? true: false; + // Initialize wind farm for ( MFIter mfi(mf_Nturb,TilingIfNotGPU()); mfi.isValid(); ++mfi) { const Box& bx = mfi.tilebox(); auto Nturb_array = mf_Nturb.array(mfi); + const Array4& z_cc_arr = (z_phys_cc) ? z_phys_cc->const_array(mfi) : Array4{}; + int k0 = bx.smallEnd()[2]; ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { int li = amrex::min(amrex::max(i, i_lo), i_hi); int lj = amrex::min(amrex::max(j, j_lo), j_hi); @@ -391,22 +482,56 @@ WindFarm::fill_Nturb_multifab (const Geometry& geom, if( d_xloc_ptr[it]+1e-3 > x1 and d_xloc_ptr[it]+1e-3 < x2 and d_yloc_ptr[it]+1e-3 > y1 and d_yloc_ptr[it]+1e-3 < y2){ Nturb_array(i,j,k,0) = Nturb_array(i,j,k,0) + 1; + if(is_terrain) { + d_zloc_ptr[it] = z_cc_arr(i,j,k0); + d_turb_index_ptr[it] = it; + } + else { + d_zloc_ptr[it] = 0.0; + d_turb_index_ptr[it] = it; + } } } }); } + + Gpu::copy(Gpu::deviceToHost, d_zloc.begin(), d_zloc.end(), zloc.begin()); + Gpu::copy(Gpu::deviceToHost, d_turb_index.begin(), d_turb_index.end(), turb_index.begin()); + + if(is_terrain) { + + std::vector> turb_index_zloc; + for(int it=0;it> turb_index_zloc_glob; + + gatherKeyValuePairs(turb_index_zloc, turb_index_zloc_glob); + + // Each process now has the global array + for (const auto& kv : turb_index_zloc_glob) { + std::cout << "Rank " << amrex::ParallelDescriptor::MyProc() << "Global data" << kv.first << " " << kv.second << "\n"; + zloc[kv.first] = kv.second; + } + } } void WindFarm::fill_SMark_multifab (const Geometry& geom, MultiFab& mf_SMark, const Real& sampling_distance_by_D, - const Real& turb_disk_angle) + const Real& turb_disk_angle, + std::unique_ptr& z_phys_cc) { amrex::Gpu::DeviceVector d_xloc(xloc.size()); amrex::Gpu::DeviceVector d_yloc(yloc.size()); + amrex::Gpu::DeviceVector d_zloc(yloc.size()); amrex::Gpu::copyAsync(amrex::Gpu::hostToDevice, xloc.begin(), xloc.end(), d_xloc.begin()); amrex::Gpu::copyAsync(amrex::Gpu::hostToDevice, yloc.begin(), yloc.end(), d_yloc.begin()); + amrex::Gpu::copyAsync(amrex::Gpu::hostToDevice, zloc.begin(), zloc.end(), d_zloc.begin()); Real d_rotor_rad = rotor_rad; Real d_hub_height = hub_height; @@ -414,6 +539,7 @@ WindFarm::fill_SMark_multifab (const Geometry& geom, Real* d_xloc_ptr = d_xloc.data(); Real* d_yloc_ptr = d_yloc.data(); + Real* d_zloc_ptr = d_zloc.data(); mf_SMark.setVal(-1.0); @@ -434,9 +560,13 @@ WindFarm::fill_SMark_multifab (const Geometry& geom, // Initialize wind farm for ( MFIter mfi(mf_SMark,TilingIfNotGPU()); mfi.isValid(); ++mfi) { - const Box& bx = mfi.tilebox(); + //const Box& bx = mfi.tilebox(); + const Box& gbx = mfi.growntilebox(1); auto SMark_array = mf_SMark.array(mfi); - ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + + const Array4& z_cc_arr = (z_phys_cc) ? z_phys_cc->const_array(mfi) : Array4{}; + + ParallelFor(gbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { int ii = amrex::min(amrex::max(i, i_lo), i_hi); int jj = amrex::min(amrex::max(j, j_lo), j_hi); int kk = amrex::min(amrex::max(k, k_lo), k_hi); @@ -446,7 +576,8 @@ WindFarm::fill_SMark_multifab (const Geometry& geom, Real y1 = ProbLoArr[1] + jj*dx[1]; Real y2 = ProbLoArr[1] + (jj+1)*dx[1]; - Real z = ProbLoArr[2] + (kk+0.5) * dx[2]; + //Real z = ProbLoArr[2] + (kk+0.5) * dx[2]; + Real z = (z_cc_arr) ? z_cc_arr(ii,jj,kk) : ProbLoArr[2] + (kk+0.5) * dx[2]; int turb_indices_overlap[2]; int check_int = 0; @@ -454,17 +585,21 @@ WindFarm::fill_SMark_multifab (const Geometry& geom, Real x0 = d_xloc_ptr[it] + d_sampling_distance*nx; Real y0 = d_yloc_ptr[it] + d_sampling_distance*ny; + Real z0 = 0.0; + if(z_cc_arr) { + z0 = d_zloc_ptr[it]; + } + bool is_cell_marked = find_if_marked(x1, x2, y1, y2, x0, y0, - nx, ny, d_hub_height, d_rotor_rad, z); + nx, ny, d_hub_height+z0, d_rotor_rad, z); if(is_cell_marked) { SMark_array(i,j,k,0) = it; } x0 = d_xloc_ptr[it]; y0 = d_yloc_ptr[it]; - //printf("Values are %d, %0.15g, %0.15g\n", it, x0, y0); is_cell_marked = find_if_marked(x1, x2, y1, y2, x0, y0, - nx, ny, d_hub_height, d_rotor_rad, z); + nx, ny, d_hub_height+z0, d_rotor_rad, z); if(is_cell_marked) { SMark_array(i,j,k,1) = it; turb_indices_overlap[check_int] = it; @@ -501,11 +636,14 @@ WindFarm::write_turbine_locations_vtk () void -WindFarm::write_actuator_disks_vtk (const Geometry& geom) +WindFarm::write_actuator_disks_vtk (const Geometry& geom, + const Real& sampling_distance_by_D) { + Real sampling_distance = sampling_distance_by_D*2.0*rotor_rad; + if (ParallelDescriptor::IOProcessor()){ - FILE *file_actuator_disks_all, *file_actuator_disks_in_dom; + FILE *file_actuator_disks_all, *file_actuator_disks_in_dom, *file_averaging_disks_in_dom; file_actuator_disks_all = fopen("actuator_disks_all.vtk","w"); fprintf(file_actuator_disks_all, "%s\n","# vtk DataFile Version 3.0"); fprintf(file_actuator_disks_all, "%s\n","Actuator Disks"); @@ -518,6 +656,13 @@ WindFarm::write_actuator_disks_vtk (const Geometry& geom) fprintf(file_actuator_disks_in_dom, "%s\n","ASCII"); fprintf(file_actuator_disks_in_dom, "%s\n","DATASET POLYDATA"); + file_averaging_disks_in_dom = fopen("averaging_disks_in_dom.vtk","w"); + fprintf(file_averaging_disks_in_dom, "%s\n","# vtk DataFile Version 3.0"); + fprintf(file_averaging_disks_in_dom, "%s\n","Actuator Disks"); + fprintf(file_averaging_disks_in_dom, "%s\n","ASCII"); + fprintf(file_averaging_disks_in_dom, "%s\n","DATASET POLYDATA"); + + int npts = 100; fprintf(file_actuator_disks_all, "%s %ld %s\n", "POINTS", xloc.size()*npts, "float"); auto ProbLoArr = geom.ProbLoArray(); @@ -534,25 +679,35 @@ WindFarm::write_actuator_disks_vtk (const Geometry& geom) } } fprintf(file_actuator_disks_in_dom, "%s %ld %s\n", "POINTS", static_cast(num_turb_in_dom*npts), "float"); + fprintf(file_averaging_disks_in_dom, "%s %ld %s\n", "POINTS", static_cast(num_turb_in_dom*npts), "float"); Real nx = std::cos(my_turb_disk_angle+0.5*M_PI); Real ny = std::sin(my_turb_disk_angle+0.5*M_PI); + Real nx1 = -std::cos(my_turb_disk_angle); + Real ny1 = -std::sin(my_turb_disk_angle); + for(int it=0; it ProbLoArr[0] and xloc[it] < ProbHiArr[0] and yloc[it] > ProbLoArr[1] and yloc[it] < ProbHiArr[1]) { fprintf(file_actuator_disks_in_dom, "%0.15g %0.15g %0.15g\n", x, y, z); + fprintf(file_averaging_disks_in_dom, "%0.15g %0.15g %0.15g\n", xavg, yavg, z); } } } fprintf(file_actuator_disks_all, "%s %ld %ld\n", "LINES", xloc.size()*(npts-1), static_cast(xloc.size()*(npts-1)*3)); fprintf(file_actuator_disks_in_dom, "%s %ld %ld\n", "LINES", static_cast(num_turb_in_dom*(npts-1)), static_cast(num_turb_in_dom*(npts-1)*3)); + fprintf(file_averaging_disks_in_dom, "%s %ld %ld\n", "LINES", static_cast(num_turb_in_dom*(npts-1)), static_cast(num_turb_in_dom*(npts-1)*3)); for(int it=0; it(2), + static_cast(it*npts+pt), + static_cast(it*npts+pt+1)); + } + } + fclose(file_actuator_disks_all); fclose(file_actuator_disks_in_dom); + fclose(file_averaging_disks_in_dom); } } diff --git a/Source/WindFarmParametrization/ERF_WindFarm.H b/Source/WindFarmParametrization/ERF_WindFarm.H index a94778fdf..68265a6db 100644 --- a/Source/WindFarmParametrization/ERF_WindFarm.H +++ b/Source/WindFarmParametrization/ERF_WindFarm.H @@ -71,16 +71,22 @@ public: void read_windfarm_spec_table_extra (const std::string windfarm_spec_table_extra); void fill_Nturb_multifab (const amrex::Geometry& geom, - amrex::MultiFab& mf_Nturb); + amrex::MultiFab& mf_Nturb, + std::unique_ptr& z_phys_cc); void fill_SMark_multifab (const amrex::Geometry& geom, amrex::MultiFab& mf_SMark, const amrex::Real& sampling_distance_by_D, - const amrex::Real& turb_disk_angle); + const amrex::Real& turb_disk_angle, + std::unique_ptr& z_phys_cc); void write_turbine_locations_vtk (); - void write_actuator_disks_vtk (const amrex::Geometry& geom); + void write_actuator_disks_vtk (const amrex::Geometry& geom, + const amrex::Real& sampling_distance_by_D); + + void gatherKeyValuePairs(const std::vector>& localData, + std::vector>& globalData); void advance (const amrex::Geometry& a_geom, const amrex::Real& dt_advance, @@ -142,7 +148,8 @@ public: protected: - amrex::Vector xloc, yloc; + amrex::Vector xloc, yloc, zloc; + amrex::Vector turb_index; amrex::Real my_turb_disk_angle; amrex::Real hub_height, rotor_rad, thrust_coeff_standing, nominal_power; amrex::Vector wind_speed, thrust_coeff, power; diff --git a/Submodules/AMReX b/Submodules/AMReX index 1c8b54530..c35aad273 160000 --- a/Submodules/AMReX +++ b/Submodules/AMReX @@ -1 +1 @@ -Subproject commit 1c8b5453023acb66bf65580b56e2f2f1dfd40c8e +Subproject commit c35aad273352c2526801f7c987fa4b3c551207f6