From f44544627d59fc46c82472c32889d00738eae22f Mon Sep 17 00:00:00 2001 From: Mahesh Natarajan Date: Wed, 18 Dec 2024 09:31:55 -0800 Subject: [PATCH] Generalizing terrain capability (#2029) * Correcting filenames * Generalizing terrain capability * Removing unused variables --------- Co-authored-by: Mahesh Natarajan Co-authored-by: Mahesh Natarajan Co-authored-by: Aaron M. Lattanzi <103702284+AMLattanzi@users.noreply.github.com> --- .../SimpleActuatorDisk/ReadTerrainUSGS.py | 208 ++++++++++++++++++ .../SimpleActuatorDisk/domain_bounds.txt | 5 + Source/ERF_ProbCommon.H | 187 +++++++++++++++- Source/Initialization/Make.package | 2 +- ...{ERF_InitWindFarm.cpp => ERF_WindFarm.cpp} | 73 ++++-- Source/WindFarmParametrization/Make.package | 2 +- 6 files changed, 452 insertions(+), 25 deletions(-) create mode 100644 Exec/WindFarmTests/WindFarm/SimpleActuatorDisk/ReadTerrainUSGS.py create mode 100644 Exec/WindFarmTests/WindFarm/SimpleActuatorDisk/domain_bounds.txt rename Source/WindFarmParametrization/{ERF_InitWindFarm.cpp => ERF_WindFarm.cpp} (93%) diff --git a/Exec/WindFarmTests/WindFarm/SimpleActuatorDisk/ReadTerrainUSGS.py b/Exec/WindFarmTests/WindFarm/SimpleActuatorDisk/ReadTerrainUSGS.py new file mode 100644 index 000000000..379727aea --- /dev/null +++ b/Exec/WindFarmTests/WindFarm/SimpleActuatorDisk/ReadTerrainUSGS.py @@ -0,0 +1,208 @@ +import math +import numpy as np +import sys +import os +import rasterio +from rasterio.transform import xy +from rasterio.warp import transform_bounds +from rasterio.transform import rowcol + +# Check if two filename arguments are provided +if len(sys.argv) != 3: + print("Usage: python3 ") + sys.exit(1) + +# Get the filenames from the command-line arguments +filename1 = sys.argv[1] +filename2 = sys.argv[2] + +# Check if both files exist, throw an error if not +if not os.path.isfile(filename1): + raise FileNotFoundError(f"Error: The file '{filename1}' does not exist.") +if not os.path.isfile(filename2): + raise FileNotFoundError(f"Error: The file '{filename2}' does not exist.") + + +# Path to the GeoTIFF file +geotiff_file = filename1 + +domain_bounds = np.loadtxt(filename2) + +# Define domain bounds +domain_lon_min = domain_bounds[0] +domain_lon_max = domain_bounds[1] +domain_lat_min = domain_bounds[2] +domain_lat_max = domain_bounds[3] + +# Constants +rad_earth = 6371000.0 # Radius of Earth in meters +M_PI = math.pi +dist_per_deg_lat = 6371000.0*2.0*M_PI/(2.0*180.0) + +row = -1 +col = -1 + +domain_row_min = 0 +domain_row_max = 0 +domain_col_min = 0 +domain_col_max = 0 + +# Open the GeoTIFF file +with rasterio.open(geotiff_file) as src: + # Get the affine transform of the raster + transform = src.transform + + # Calculate the row corresponding to the given latitude + domain_row_min = int((src.bounds.top - domain_lat_max) / abs(transform[4])) # transform[4] is the pixel height (negative) + domain_row_max = int((src.bounds.top - domain_lat_min) / abs(transform[4])) # transform[4] is the pixel height (negative) + domain_col_min = int((domain_lon_min - transform[2]) / transform[0]) # transform[2] is the x-coordinate of the top-left corner + domain_col_max = int((domain_lon_max - transform[2]) / transform[0]) # transform[2] is the x-coordinate of the top-left corner + + print("src.bounds_top are %0.15g %0.15g\n"%(src.bounds.top, transform[4])); + print("The row and cols are %d, %d, %d, %d"%(domain_row_min, domain_row_max, domain_col_min, domain_col_max)); + + # Convert row and column back to latitude and longitude + lon_min = transform[2] + domain_col_min * transform[0] + lon_max = transform[2] + domain_col_max * transform[0] + + lat_min = src.bounds.top + domain_row_max * transform[4] + lat_max = src.bounds.top + domain_row_min * transform[4] + + print("The lon min and max are %0.15g, %0.15g"%(lon_min, lon_max)); + print("The lat min and max are %0.15g, %0.15g"%(lat_min, lat_max)); + + + print("Values are %.15g, %0.15g, %0.15g, %0.15g\n"%(src.bounds.top,domain_row_max,transform[4], lat_min)); +# Open the GeoTIFF file +with rasterio.open(geotiff_file) as src: + + # Read the raster shape (rows, columns) + rows, cols = src.shape + transform = src.transform + + # Print CRS and file bounds for context + print(f"CRS: {src.crs}") + print(f"Bounds: {src.bounds}") + print(f"Raster Shape (rows, cols): {rows}, {cols}\n") + + if(domain_bounds[0] <= src.bounds.left or domain_bounds[0] >= src.bounds.right): + print("The longitude min (the first entry) specified in %s is %0.15g," + " which is not within the data longitude bounds (%0.15g, %0.15g). Exiting....\n"% + (filename2, domain_bounds[0], src.bounds.left,src.bounds.right)); + sys.exit(); + + if(domain_bounds[1] <= src.bounds.left or domain_bounds[1] >= src.bounds.right): + print("The longitude max (the second entry) specified in %s is %0.15g," + " which is not within the data longitude bounds (%0.15g, %0.15g). Exiting....\n"% + (filename2, domain_bounds[1], src.bounds.left,src.bounds.right)); + sys.exit(); + + if(domain_bounds[2] <= src.bounds.bottom or domain_bounds[2] >= src.bounds.top): + print("The latitude min (the third entry) specified in %s is %0.15g," + " which is not within the data longitude bounds (%0.15g, %0.15g). Exiting....\n"% + (filename2, domain_bounds[2], src.bounds.bottom,src.bounds.top)); + sys.exit(); + + if(domain_bounds[3] <= src.bounds.bottom or domain_bounds[3] >= src.bounds.top): + print("The latitude max (the fourth entry) specified in %s is %0.15g," + " which is not within the data longitude bounds (%0.15g, %0.15g). Exiting....\n"% + (filename2, domain_bounds[3], src.bounds.bottom,src.bounds.top)); + sys.exit(); + + # Get the bounds in the source CRS + bounds = src.bounds + + # Transform bounds to latitude and longitude (EPSG:4326) + lat_lon_bounds = transform_bounds(src.crs, "EPSG:4326", bounds.left, bounds.bottom, bounds.right, bounds.top) + + print("Data lat min %0.15g\n"%bounds.bottom); + print("Data lat max %0.15g\n"%bounds.top); + print("Data lon min %0.15g\n"%bounds.left); + print("Data lon max %0.15g\n"%bounds.right); + + lat_min = lat_min*M_PI / 180.0 + lon_min = lon_min*M_PI / 180.0 + + # Read the first band of elevation data + elevation_data = src.read(1) + + xloc = 0; + yloc = [] # To store y-coordinates + + vtk_file = open("terrain_mesh.vtk",'w') + + vtk_file.write("# vtk DataFile Version 2.0\n") + vtk_file.write("Structured Grid Example\n") + vtk_file.write("ASCII\n\n") + + # Write the structured grid header + vtk_file.write(f"DATASET STRUCTURED_GRID\n") + nskip = 1 + rows_in_output = int((domain_row_max-domain_row_min)/nskip); + cols_in_output = int((domain_col_max-domain_col_min)/nskip); + vtk_file.write(f"DIMENSIONS {cols_in_output} {rows_in_output} 1\n") # Adjust for structured grid dimensions + vtk_file.write(f"POINTS {rows_in_output*cols_in_output} float\n") # Adjust for structured grid dimensions + + file_for_erf = open("ERF_terrain_file.txt",'w') + file_for_erf.write("%0.15g %0.15g\n"%(lon_min*180/M_PI,lat_min*180/M_PI)) + file_for_erf.write("%d %d\n"%(domain_col_max-domain_col_min,domain_row_max-domain_row_min)); + + for row in range(domain_row_max, domain_row_min, -1): # Sample every 10% of rows + print("Doing %d of %d"%(row-domain_row_min, domain_row_max-domain_row_min)) + for col in range(domain_col_min, domain_col_max, nskip): # Sample every 10% of columns + + # Read the value at this pixel location + value = elevation_data[row, col] + + # Convert row and column back to latitude and longitude + lon = transform[2] + col * transform[0] # lon at the given column + lat = src.bounds.top + row * transform[4] # lat at the given row + + # Get the elevation value + elevation = elevation_data[row, col] + + if(elevation==0.0): + print("Elevation is zero...cannot be\n"); + sys.exit() + + # Conversion and calculations + lat = lat * M_PI / 180.0 + lon = lon * M_PI / 180.0 + delta_lat = lat - lat_min + delta_lon = lon - lon_min + + term1 = math.pow(math.sin(delta_lat / 2.0), 2) + term2 = math.cos(lat) * math.cos(lat_min) * math.pow(math.sin(delta_lon / 2.0), 2) + dist = 2.0 * rad_earth * math.asin(math.sqrt(term1 + term2)) + dy = (lat - lat_min) * dist_per_deg_lat * 180.0 / M_PI + fac = math.pow(dist, 2) - math.pow(dy, 2) + if(dist < dy): + if(abs(dist-dy) > 1e-10): + print("Error in calculation dist cannot be less than dy %0.15g, %0.15g"%(dist, dy)) + sys.exit() + + if(abs(dist-dy) < 1e-10): + fac = 0.0 + + dx = math.sqrt(fac) + + xloc = dx + file_for_erf.write("%0.15g %0.15g %0.15g\n"%(xloc,dy,elevation)); + + vtk_file.write(f"{xloc} {dy} {elevation}\n") + if(col==domain_col_min): + yloc.append(dy) + + file_for_erf.close() + +# Write the elevation as a scalar field + vtk_file.write(f"\nPOINT_DATA {rows_in_output * cols_in_output}\n") + vtk_file.write("SCALARS Elevation float 1\n") + vtk_file.write("LOOKUP_TABLE default\n") + for row in range(domain_row_max, domain_row_min, -1): + for col in range(domain_col_min, domain_col_max, nskip): + elevation_value = elevation_data[row, col] # Elevation for the current point + vtk_file.write(f"{elevation_value}\n") + +file_for_erf.close() + diff --git a/Exec/WindFarmTests/WindFarm/SimpleActuatorDisk/domain_bounds.txt b/Exec/WindFarmTests/WindFarm/SimpleActuatorDisk/domain_bounds.txt new file mode 100644 index 000000000..2c2f2954e --- /dev/null +++ b/Exec/WindFarmTests/WindFarm/SimpleActuatorDisk/domain_bounds.txt @@ -0,0 +1,5 @@ +-97.57 +-97.34 +36.35 +36.45 + diff --git a/Source/ERF_ProbCommon.H b/Source/ERF_ProbCommon.H index 75c4464a4..6211fe02c 100644 --- a/Source/ERF_ProbCommon.H +++ b/Source/ERF_ProbCommon.H @@ -274,11 +274,14 @@ public: const amrex::Real& time) { // Check if a valid text file exists for the terrain - std::string fname; + std::string fname, fname_usgs, ftype; amrex::ParmParse pp("erf"); auto valid_fname = pp.query("terrain_file_name",fname); + auto valid_fname_USGS = pp.query("terrain_file_name_USGS",fname_usgs); if (valid_fname) { this->read_custom_terrain(fname,geom,z_phys_nd,time); + } else if(valid_fname_USGS) { + this->read_custom_terrain_USGS(fname_usgs,geom,z_phys_nd,time); } else { // Note that this only sets the terrain value at the ground IF k=0 is in the box amrex::Print() << "Initializing flat terrain at z=0" << std::endl; @@ -308,6 +311,7 @@ public: } } + /** * Function to perform custom initialization of terrain * @@ -423,6 +427,187 @@ public: } } + virtual void + read_custom_terrain_USGS (const std::string& fname_usgs, + const amrex::Geometry& geom, + amrex::MultiFab& z_phys_nd, + const amrex::Real& /*time*/) + { + // Read terrain file on the host + amrex::Print()<<"Reading terrain file: "<< fname_usgs << std::endl; + std::ifstream file(fname_usgs); + + if (!file.is_open()) { + amrex::Abort("Error: Could not open the file " + fname_usgs + "\n"); + } + + // Check if file is empty + if (file.peek() == std::ifstream::traits_type::eof()) { + amrex::Abort("Error: The file " + fname_usgs + " is empty.\n"); + } + + amrex::Gpu::HostVector m_xterrain,m_yterrain,m_zterrain; + amrex::Real value1,value2,value3; + amrex::Real lat_min, lon_min; + int nlons, nlats; + + file >> lon_min >> lat_min; + if(std::fabs(lon_min) > 180.0) { + amrex::Error("The value of longitude for entry in the first line in " + fname_usgs + + " should not exceed 180. It is " + std::to_string(lon_min)); + } + if(std::fabs(lat_min) > 90.0) { + amrex::Error("The value of latitude for entry in the first line in " + fname_usgs + + " should not exceed 90. It is " + std::to_string(lat_min)); + } + + file >> nlons >> nlats; + + int counter = 0; + while(file >> value1 >> value2 >> value3){ + m_xterrain.push_back(value1); + if(counter%nlons==0) { + m_yterrain.push_back(value2); + } + m_zterrain.push_back(value3); + counter += 1; + } + file.close(); + + AMREX_ASSERT(m_xterrain.size() == static_cast(nlons*nlats)); + AMREX_ASSERT(m_yterrain.size() == static_cast(nlats)); + AMREX_ASSERT(m_zterrain.size() == static_cast(nlons*nlats)); + + // Shift the terrain so that there is some flat region for the + // inflow to come in + + amrex::ParmParse pp("erf"); + amrex::Real x_shift, y_shift; + pp.query("windfarm_x_shift",x_shift); + pp.query("windfarm_x_shift",y_shift); + + + for(auto& v : m_xterrain){ + v += x_shift; + } + for(auto& v : m_yterrain){ + v += y_shift; + } + + if (m_xterrain.empty()) { + amrex::Abort("Error: m_xterrain is empty!\n"); + } + + auto index = std::max_element(m_xterrain.begin(), m_xterrain.end()); + amrex::Real xmax_terrain = *index; + index = std::max_element(m_yterrain.begin(), m_yterrain.end()); + amrex::Real ymax_terrain = *index; + + index = std::min_element(m_zterrain.begin(), m_zterrain.end()); + amrex::Real zmin_terrain = *index; + + for(auto& v : m_zterrain){ + v -= zmin_terrain; + } + + // Copy data to the GPU + int nnode = nlons*nlats; + amrex::Gpu::DeviceVector d_xterrain(nnode),d_yterrain(nlats),d_zterrain(nnode); + amrex::Gpu::copy(amrex::Gpu::hostToDevice, m_xterrain.begin(), m_xterrain.end(), d_xterrain.begin()); + amrex::Gpu::copy(amrex::Gpu::hostToDevice, m_yterrain.begin(), m_yterrain.end(), d_yterrain.begin()); + amrex::Gpu::copy(amrex::Gpu::hostToDevice, m_zterrain.begin(), m_zterrain.end(), d_zterrain.begin()); + amrex::Real* d_xt = d_xterrain.data(); + amrex::Real* d_yt = d_yterrain.data(); + amrex::Real* d_zt = d_zterrain.data(); + + // Populate z_phys data + int ngrow = z_phys_nd.nGrow(); + auto dx = geom.CellSizeArray(); + auto ProbLoArr = geom.ProbLoArray(); + int ilo = geom.Domain().smallEnd(0); + int jlo = geom.Domain().smallEnd(1); + int klo = geom.Domain().smallEnd(2); + int ihi = geom.Domain().bigEnd(0) + 1; + int jhi = geom.Domain().bigEnd(1) + 1; + + for (amrex::MFIter mfi(z_phys_nd,amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi ) + { + // 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); + amrex::ParallelFor(xybx, [=] AMREX_GPU_DEVICE (int i, int j, int /*k*/) + { + // Clip indices for ghost-cells + int ii = amrex::min(amrex::max(i,ilo),ihi); + int jj = amrex::min(amrex::max(j,jlo),jhi); + + // Location of nodes + amrex::Real x = ProbLoArr[0] + ii * dx[0] + 1e-3; + amrex::Real y = ProbLoArr[1] + jj * dx[1] + 1e-3; + + if(x < x_shift or x > xmax_terrain or + y < y_shift or y > ymax_terrain) { + z_arr(i,j,klo) = 0.0; + } else { + int jindex_terrain = -1; + for (int it = 0; it < nlats && jindex_terrain == -1; it++) { + if (d_yt[it] > y) { + jindex_terrain = it-1; + } + } + //if(jindex_terrain == -1){ + // printf("Could not find the jindex for a y value of %0.15g\n",y); + //} + // Now the index to check for x goes from jindex_terrain*nlons to + // (jindex_terrain+1)*nlons-1 + int iindex_terrain=-1; + int gstart = jindex_terrain*nlons; + int gend = (jindex_terrain+1)*nlons-1; + for (int it = gstart; it <= gend && iindex_terrain == -1; it++) { + if (d_xt[it] > x) { + iindex_terrain = it-gstart; + } + } + //if(iindex_terrain == -1){ + // printf("Could not find the iindex for a x value of %0.15g for a y index of %d and y value of %0.15g\n",x, jindex_terrain,y); + //} + + // Now do an averaging + int ind11, ind12, ind21, ind22; + ind11 = jindex_terrain*nlons + iindex_terrain; + ind12 = ind11+nlons; + ind21 = ind11+1; + ind22 = ind12+1; + amrex::Real x1 = d_xt[ind11]; + amrex::Real x2 = d_xt[ind21]; + amrex::Real y1 = d_yt[jindex_terrain]; + amrex::Real y2 = d_yt[jindex_terrain+1]; + // Do bilinear interpolation + amrex::Real denom = (x2-x1)*(y2-y1); + amrex::Real w_11 = (x2-x)*(y2-y)/denom; + amrex::Real w_12 = (x2-x)*(y-y1)/denom; + amrex::Real w_21 = (x-x1)*(y2-y)/denom; + amrex::Real w_22 = (x-x1)*(y-y1)/denom; + + z_arr(i,j,klo) = w_11*d_zt[ind11] + w_12*d_zt[ind12] + w_21*d_zt[ind21] + w_22*d_zt[ind22]; + + //z_arr(i,j,klo) = 0.25*(d_zt[ind1] + d_zt[ind2] + d_zt[ind3] + d_zt[ind4]); + // Smoothen out the edges using Gaussian + /*amrex::Real sigma = 600.0; + amrex::Real fac = 2.0*sigma*sigma; + z_arr(i,j,klo) = z_arr(i,j,klo)*(1.0-exp(-(x-x_shift)*(x-x_shift)/fac))* + (1.0-exp(-(x-xmax_terrain)*(x-xmax_terrain)/fac))* + (1.0-exp(-(y-y_shift)*(y-y_shift)/fac))* + (1.0-exp(-(y-ymax_terrain)*(y-ymax_terrain)/fac)); */ + //std::cout << "Values of z is " << d_zt[ind1] << " " << d_zt[ind2] << " " << d_zt[ind3] << " " << d_zt[ind4] << "\n"; + } + }); + } + } + + #ifdef ERF_USE_TERRAIN_VELOCITY virtual amrex::Real compute_terrain_velocity(const amrex::Real /*time*/) { diff --git a/Source/Initialization/Make.package b/Source/Initialization/Make.package index 114b59052..3cc1d623d 100644 --- a/Source/Initialization/Make.package +++ b/Source/Initialization/Make.package @@ -10,7 +10,7 @@ CEXE_sources += ERF_Init1D.cpp CEXE_sources += ERF_InitTurbPert.cpp ifeq ($(USE_WINDFARM),TRUE) -CEXE_sources += ERFInitWindFarm.cpp +CEXE_sources += ERF_InitWindFarm.cpp endif ifeq ($(USE_NETCDF),TRUE) diff --git a/Source/WindFarmParametrization/ERF_InitWindFarm.cpp b/Source/WindFarmParametrization/ERF_WindFarm.cpp similarity index 93% rename from Source/WindFarmParametrization/ERF_InitWindFarm.cpp rename to Source/WindFarmParametrization/ERF_WindFarm.cpp index df1fe76d9..948478d54 100644 --- a/Source/WindFarmParametrization/ERF_InitWindFarm.cpp +++ b/Source/WindFarmParametrization/ERF_WindFarm.cpp @@ -84,32 +84,57 @@ WindFarm::init_windfarm_lat_lon (const std::string windfarm_loc_table, file.close(); Real rad_earth = 6371.0e3; // Radius of the earth + Real m_per_deg_lat = rad_earth*2.0*M_PI/(2.0*180.0); // Find the coordinates of average of min and max of the farm // Rotate about that point + ParmParse pp("erf"); + std::string fname_usgs; + auto valid_fname_USGS = pp.query("terrain_file_name_USGS",fname_usgs); + Real lon_ref, lat_ref; + + if (valid_fname_USGS) { + std::ifstream file_usgs(fname_usgs); + file_usgs >> lon_ref >> lat_ref; + file_usgs.close(); + lon_ref = lon_ref*M_PI/180.0; + lat_ref = lat_ref*M_PI/180.0; + } else { + Real lat_min = *std::min_element(lat.begin(), lat.end()); + Real lon_min = *std::min_element(lon.begin(), lon.end()); - Real lat_min = *std::min_element(lat.begin(), lat.end()); - Real lat_max = *std::max_element(lat.begin(), lat.end()); - Real lon_min = *std::min_element(lon.begin(), lon.end()); - Real lon_max = *std::max_element(lon.begin(), lon.end()); - - Real lat_cen = 0.5*(lat_min+lat_max)*M_PI/180.0; - Real lon_cen = 0.5*(lon_min+lon_max)*M_PI/180.0; - - // (lat_lo, lon_lo) is mapped to (0,0) + lon_ref = lon_min*M_PI/180.0; + lat_ref = lat_min*M_PI/180.0; + } for(int it=0;it= 0.0) { xloc.push_back(dx_turb); } @@ -119,12 +144,9 @@ WindFarm::init_windfarm_lat_lon (const std::string windfarm_loc_table, yloc.push_back(dy_turb); } - Real xloc_min = *std::min_element(xloc.begin(),xloc.end()); - Real yloc_min = *std::min_element(yloc.begin(),yloc.end()); - for(int it = 0;it ProbLoArr[0] and + xloc[it] < ProbHiArr[0] and + yloc[it] > ProbLoArr[1] and + yloc[it] < ProbHiArr[1] ) { + if(is_counted[it] != 1) { + Abort("Wind turbine " + std::to_string(it) + "has been counted " + std::to_string(is_counted[it]) + " times" + + " It should have been counted only once. Aborting...."); + } } } diff --git a/Source/WindFarmParametrization/Make.package b/Source/WindFarmParametrization/Make.package index 10f6ffb07..78c13af5d 100644 --- a/Source/WindFarmParametrization/Make.package +++ b/Source/WindFarmParametrization/Make.package @@ -1,2 +1,2 @@ CEXE_headers += ERF_WindFarm.H -CEXE_sources += ERF_InitWindFarm.cpp +CEXE_sources += ERF_WindFarm.cpp