WIP: Fitch wind farm parametrization

Exec/Fitch/CMakeLists.txt

@@ -0,0 +1,12 @@
+set(erf_exe_name erf_abl)
+
+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})

Exec/Fitch/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/Fitch
+include $(ERF_HOME)/Exec/Make.ERF

Exec/Fitch/Make.package

@@ -0,0 +1,2 @@
+CEXE_headers += prob.H
+CEXE_sources += prob.cpp

Exec/Fitch/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.

Exec/Fitch/inputs_Fitch

@@ -0,0 +1,74 @@
+# ------------------  INPUTS TO MAIN PROGRAM  -------------------
+max_step = 1800000
+
+amrex.fpe_trap_invalid = 1
+
+fabarray.mfiter_tile_size = 1024 1024 1024
+
+# PROBLEM SIZE & GEOMETRY
+geometry.prob_extent =  18000    18000  1000
+amr.n_cell           =   32     32      64
+
+#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.005  # fixed time step depending on grid resolution
+
+# DIAGNOSTICS & VERBOSITY
+erf.sum_interval   = 1       # timesteps between computing mass
+erf.v              = 1       # verbosity in ERF.cpp
+amr.v                = 1       # verbosity in Amr.cpp
+
+# REFINEMENT / REGRIDDING
+amr.max_level       = 0       # maximum level number allowed
+
+# CHECKPOINT FILES
+erf.check_file      = chk       # root name of checkpoint file
+erf.check_int       = -1        # number of timesteps between checkpoints
+
+# 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
+
+# SOLVER CHOICE
+erf.alpha_T = 0.0
+erf.alpha_C = 1.0
+erf.use_gravity = false
+
+erf.molec_diff_type = "None"
+erf.les_type        = "None"
+erf.Cs              = 1.5
+erf.dynamicViscosity = 0.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
+

Exec/Fitch/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()
+

Exec/Fitch/prob.H

@@ -0,0 +1,78 @@
+#ifndef _PROB_H_
+#define _PROB_H_
+
+#include <string>
+
+#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<amrex::Real      > const& state,
+        amrex::Array4<amrex::Real      > const& x_vel,
+        amrex::Array4<amrex::Real      > const& y_vel,
+        amrex::Array4<amrex::Real      > const& z_vel,
+        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;
+
+protected:
+    std::string name() override { return "ABL"; }
+
+private:
+    ProbParm parms;
+};
+
+#endif