Skip to content

Commit

Permalink
WIP: Generalize wind farm parametrization (#1511)
Browse files Browse the repository at this point in the history
* Generalizing Fitch

* WIP: Generalizing Fitch

* Minor changes

* Reverting changes for inputs

* Remove unwanted function

* Adding ifdef to plotfile

* Some minor changes to comments

* Adding params to inputs

* Adding params to inputs

---------

Co-authored-by: Mahesh Natarajan <[email protected]>
Co-authored-by: Aaron M. Lattanzi <[email protected]>
  • Loading branch information
3 people authored Mar 20, 2024
1 parent 5258b98 commit c8677ef
Show file tree
Hide file tree
Showing 9 changed files with 241 additions and 2 deletions.
76 changes: 76 additions & 0 deletions Exec/Fitch/inputs_Fitch_LargeDomain
Original file line number Diff line number Diff line change
@@ -0,0 +1,76 @@
# ------------------ INPUTS TO MAIN PROGRAM -------------------
max_step = 10

amrex.fpe_trap_invalid = 1

fabarray.mfiter_tile_size = 1024 1024 1024

# PROBLEM SIZE & GEOMETRY
erf.windfarm_type = "Fitch"
erf.latitude_lo = 35.0
erf.longitude_lo = -100.0
geometry.prob_extent = 534000 389000 20000
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 num_turb

# 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"


# 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

30 changes: 30 additions & 0 deletions Exec/Fitch/windturbines.txt
Original file line number Diff line number Diff line change
@@ -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
5 changes: 5 additions & 0 deletions Source/DataStructs/DataStruct.H
Original file line number Diff line number Diff line change
Expand Up @@ -255,7 +255,10 @@ struct SolverChoice {

// Which type of windfarm model
static std::string windfarm_type_string = "Fitch";

pp.query("windfarm_type", windfarm_type_string);
pp.query("latitude_lo", latitude_lo);
pp.query("longitude_lo", longitude_lo);
if (windfarm_type_string == "Fitch") {
windfarm_type = WindFarmType::Fitch;
}
Expand Down Expand Up @@ -421,5 +424,7 @@ struct SolverChoice {
bool do_cloud {true};
bool do_precip {true};
bool use_moist_background {false};

amrex::Real latitude_lo, longitude_lo;
};
#endif
8 changes: 6 additions & 2 deletions Source/ERF.H
Original file line number Diff line number Diff line change
Expand Up @@ -342,6 +342,10 @@ public:
void init_from_metgrid (int lev);
#endif // ERF_USE_NETCDF

#ifdef ERF_USE_WINDFARM
void init_windfarm(int lev);
#endif

// more flexible version of AverageDown() that lets you average down across multiple levels
void AverageDownTo (int crse_lev, int scomp, int ncomp); // NOLINT

Expand Down Expand Up @@ -712,7 +716,7 @@ private:

// Note that the order of variable names here must match the order in Derive.cpp
const amrex::Vector<std::string> derived_names {"soundspeed", "temp", "theta", "KE", "QKE", "scalar",
"pres_hse", "dens_hse", "pressure", "pert_pres", "pert_dens", "eq_pot_temp",
"pres_hse", "dens_hse", "pressure", "pert_pres", "pert_dens", "eq_pot_temp", "num_turb",
"dpdx", "dpdy", "pres_hse_x", "pres_hse_y",
"z_phys", "detJ" , "mapfac",
// eddy viscosity
Expand All @@ -722,7 +726,7 @@ private:
// mynn pbl lengthscale
"Lpbl",
// moisture vars
"qt", "qv", "qc", "qi", "qp", "qrain", "qsnow", "qgraup", "rain_accum"
"qt", "qv", "qc", "qi", "qp", "qrain", "qsnow", "qgraup", "rain_accum",
#ifdef ERF_COMPUTE_ERROR
,"xvel_err", "yvel_err", "zvel_err", "pp_err"
#endif
Expand Down
6 changes: 6 additions & 0 deletions Source/ERF.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -1091,6 +1091,12 @@ ERF::init_only (int lev, Real time)
lev_new[Vars::xvel].OverrideSync(geom[lev].periodicity());
lev_new[Vars::yvel].OverrideSync(geom[lev].periodicity());
lev_new[Vars::zvel].OverrideSync(geom[lev].periodicity());

// Initialize wind farm

#ifdef ERF_USE_WINDFARM
init_windfarm(lev);
#endif
}

// read in some parameters from inputs file
Expand Down
1 change: 1 addition & 0 deletions Source/ERF_make_new_level.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -699,6 +699,7 @@ void ERF::init_stuff (int lev, const BoxArray& ba, const DistributionMapping& dm
if (solverChoice.windfarm_type == WindFarmType::Fitch){
int ngrow_state = ComputeGhostCells(solverChoice.advChoice, solverChoice.use_NumDiff) + 1;
vars_fitch[lev].define(ba, dm, 5, ngrow_state); // V, dVabsdt, dudt, dvdt, dTKEdt
Nturb[lev].define(ba, dm, 1, 0); // Number of turbines in a cell
}
#endif
}
Expand Down
20 changes: 20 additions & 0 deletions Source/IO/Plotfile.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -336,6 +336,26 @@ ERF::WritePlotFile (int which, Vector<std::string> plot_var_names)
mf_comp ++;
}

#ifdef ERF_USE_WINDFARM
if (containerHasElement(plot_var_names, "num_turb"))
{
std::cout << "Plotting num_turb" << "\n";
#ifdef _OPENMP
#pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
#endif
for ( MFIter mfi(mf[lev],TilingIfNotGPU()); mfi.isValid(); ++mfi)
{
const Box& bx = mfi.tilebox();
const Array4<Real>& derdat = mf[lev].array(mfi);
const Array4<Real const>& Nturb_array = Nturb[lev].const_array(mfi);
ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
derdat(i, j, k, mf_comp) = Nturb_array(i,j,k,0);
});
}
mf_comp ++;
}
#endif

int klo = geom[lev].Domain().smallEnd(2);
int khi = geom[lev].Domain().bigEnd(2);

Expand Down
96 changes: 96 additions & 0 deletions Source/Initialization/ERF_init_windfarm.cpp
Original file line number Diff line number Diff line change
@@ -0,0 +1,96 @@
/**
* \file ERF_init_windfarm.cpp
*/

#include <ERF.H>

using namespace amrex;

/**
* Read in the turbine locations in latitude-longitude from windturbines.txt
* and convert it into x and y coordinates in metres
*
* @param lev Integer specifying the current level
*/
void
ERF::init_windfarm (int lev)
{
// Read turbine locations from windturbines.txt
FILE *file_windturbines;
file_windturbines = fopen("windturbines.txt","r");
std::ifstream file("windturbines.txt");
if (!file.is_open()) {
amrex::Error("Wind turbines location file windturbines.txt not found");
}
// Vector of vectors to store the matrix
std::vector<Real> lat, lon, xloc, yloc;
Real value1, value2, value3;
while (file >> value1 >> value2 >> value3) {
lat.push_back(value1);
lon.push_back(value2);
}
file.close();

Real rad_earth = 6371.0e3; // Radius of the earth
Real lat_lo = solverChoice.latitude_lo*M_PI/180.0;
Real lon_lo = solverChoice.longitude_lo*M_PI/180.0;

// (lat_lo, lon_lo) is mapped to (0,0)

for(int it=0;it<lat.size();it++){
lat[it] = lat[it]*M_PI/180.0;
lon[it] = lon[it]*M_PI/180.0;
Real delta_lat = (lat[it] - lat_lo);
Real delta_lon = (lon[it] - lon_lo);

Real term1 = std::pow(sin(delta_lat/2.0),2);
Real term2 = cos(lat[it])*cos(lat_lo)*std::pow(sin(delta_lon/2.0),2);
Real dist = 2.0*rad_earth*std::asin(std::sqrt(term1 + term2));
Real dy_turb = (lat[it] - lat_lo) * 111000.0 * 180.0/M_PI ;
yloc.push_back(dy_turb);
Real dx_turb = std::sqrt(std::pow(dist,2) - std::pow(dy_turb,2));
xloc.push_back(dx_turb);
}

// Write out a vtk file for turbine locations
if (ParallelDescriptor::IOProcessor()){
FILE* file_turbloc_vtk;
file_turbloc_vtk = fopen("turbine_locations.vtk","w");
fprintf(file_turbloc_vtk, "%s\n","# vtk DataFile Version 3.0");
fprintf(file_turbloc_vtk, "%s\n","Wind turbine locations");
fprintf(file_turbloc_vtk, "%s\n","ASCII");
fprintf(file_turbloc_vtk, "%s\n","DATASET POLYDATA");
fprintf(file_turbloc_vtk, "%s %ld %s\n", "POINTS", xloc.size(), "float");
for(int it=0; it<xloc.size(); it++){
fprintf(file_turbloc_vtk, "%0.15g %0.15g %0.15g\n", xloc[it], yloc[it], 1e-12);
}
fclose(file_turbloc_vtk);
}

Nturb[lev].setVal(0);

int i_lo = geom[lev].Domain().smallEnd(0); int i_hi = geom[lev].Domain().bigEnd(0);
int j_lo = geom[lev].Domain().smallEnd(1); int j_hi = geom[lev].Domain().bigEnd(1);

// Initialize wind farm
for ( MFIter mfi(Nturb[lev],TilingIfNotGPU()); mfi.isValid(); ++mfi) {
const Box& bx = mfi.tilebox();
auto Nturb_array = Nturb[lev].array(mfi);
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);

auto dx = geom[lev].CellSizeArray();
Real x1 = li*dx[0], x2 = (li+1)*dx[0];
Real y1 = lj*dx[1], y2 = (lj+1)*dx[1];

for(int it=0; it<xloc.size(); it++){
if( xloc[it]+1e-12 > x1 and xloc[it]+1e-12 < x2 and
yloc[it]+1e-12 > y1 and yloc[it]+1e-12 < y2){
Nturb_array(i,j,k,0) = Nturb_array(i,j,k,0) + 1;
}
}
});
}
}

1 change: 1 addition & 0 deletions Source/Initialization/Make.package
Original file line number Diff line number Diff line change
Expand Up @@ -2,6 +2,7 @@ CEXE_sources += ERF_init_uniform.cpp
CEXE_sources += ERF_init_from_hse.cpp
CEXE_sources += ERF_init_custom.cpp
CEXE_sources += ERF_init_from_input_sounding.cpp
CEXE_sources += ERF_init_windfarm.cpp

CEXE_headers += InputSoundingData.H
CEXE_sources += ERF_init_bcs.cpp
Expand Down

0 comments on commit c8677ef

Please sign in to comment.