Skip to content

Commit

Permalink
Merge branch 'development' into metgrid_interpolator
Browse files Browse the repository at this point in the history
  • Loading branch information
AMLattanzi authored Dec 2, 2024
2 parents 8d8b5af + fc0f5a9 commit 3a9f526
Show file tree
Hide file tree
Showing 27 changed files with 609 additions and 213 deletions.
Original file line number Diff line number Diff line change
@@ -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}
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -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
Original file line number Diff line number Diff line change
@@ -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})
Original file line number Diff line number Diff line change
@@ -0,0 +1,80 @@
#ifndef ERF_PROB_H_
#define ERF_PROB_H_

#include <string>

#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<amrex::Real const> const& state,
amrex::Array4<amrex::Real > const& state_pert,
amrex::Array4<amrex::Real > const& x_vel_pert,
amrex::Array4<amrex::Real > const& y_vel_pert,
amrex::Array4<amrex::Real > const& z_vel_pert,
amrex::Array4<amrex::Real > const& r_hse,
amrex::Array4<amrex::Real > const& p_hse,
amrex::Array4<amrex::Real const> const& z_nd,
amrex::Array4<amrex::Real const> const& z_cc,
amrex::GeometryData const& geomdata,
amrex::Array4<amrex::Real const> const& mf_m,
amrex::Array4<amrex::Real const> const& mf_u,
amrex::Array4<amrex::Real const> 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
Original file line number Diff line number Diff line change
@@ -0,0 +1,238 @@
#include "ERF_prob.H"
#include "AMReX_Random.H"

using namespace amrex;

std::unique_ptr<ProblemBase>
amrex_probinit(const amrex_real* problo, const amrex_real* probhi)
{
return std::make_unique<Problem>(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<amrex::Real const> const& /*state*/,
amrex::Array4<amrex::Real > const& state_pert,
amrex::Array4<amrex::Real > const& x_vel_pert,
amrex::Array4<amrex::Real > const& y_vel_pert,
amrex::Array4<amrex::Real > const& z_vel_pert,
amrex::Array4<amrex::Real > const& /*r_hse*/,
amrex::Array4<amrex::Real > const& /*p_hse*/,
amrex::Array4<amrex::Real const> const& /*z_nd*/,
amrex::Array4<amrex::Real const> const& /*z_cc*/,
amrex::GeometryData const& geomdata,
amrex::Array4<amrex::Real const> const& /*mf_m*/,
amrex::Array4<amrex::Real const> const& /*mf_u*/,
amrex::Array4<amrex::Real const> 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<Real> 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);
});
}
}
}


Original file line number Diff line number Diff line change
@@ -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
Original file line number Diff line number Diff line change
@@ -0,0 +1,2 @@
CEXE_headers += ERF_prob.H
CEXE_sources += ERF_prob.cpp
Original file line number Diff line number Diff line change
@@ -0,0 +1,2 @@
This folder contains an example of the Simple actuator disk wind farm parametetrization
for a single turbine case.
Loading

0 comments on commit 3a9f526

Please sign in to comment.