Skip to content
This repository has been archived by the owner on Nov 27, 2023. It is now read-only.

Commit

Permalink
IO: TurbSimFile: adding fromAMRWind_PD and AMRWind class
Browse files Browse the repository at this point in the history
  • Loading branch information
ebranlard committed Jul 31, 2023
1 parent 1bb67f9 commit cd94e85
Show file tree
Hide file tree
Showing 2 changed files with 272 additions and 43 deletions.
125 changes: 125 additions & 0 deletions pyFAST/input_output/amrwind_file.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,125 @@
"""Read AMR-Wind NETCDF file
"""
import xarray as xr
import numpy as np

class AMRWindFile(dict):
"""
Read a AMR-Wind output file (.nc)
"""

@staticmethod
def defaultExtensions():
""" List of file extensions expected for this fileformat"""
return ['.nc']

@staticmethod
def formatName():
""" Short string (~100 char) identifying the file format"""
return 'NetCDF plane sampling file from AMRWind'

@staticmethod
def priority(): return 60 # Priority in weio.read fileformat list between 0=high and 100:low

def __init__(self, filename=None, timestep=None, output_frequency=None, **kwargs):
self.filename = filename
self.amrwind_dt = timestep
self.output_dt = timestep * output_frequency

if filename:
self.read(**kwargs)

def read(self, group_name):
"""
Parameters
----------
group_name : str,
group name inside netcdf file that you want to read, e.g. p_slice
TODO: see if group_name can be avoided, and add a read_group function
"""
# --- Standard tests and exceptions (generic code)
if filename:
self.filename = filename
if not self.filename:
raise Exception('No filename provided')
if not os.path.isfile(self.filename):
raise OSError(2,'File not found:',self.filename)
if os.stat(self.filename).st_size == 0:
raise Exception('File is empty:',self.filename)


ds = xr.open_dataset(self.filename,group=group_name)

coordinates = {"x":(0,"axial"), "y":(1,"lateral"),"z":(2,"vertical")}
c = {}
for coordinate,(i,desc) in coordinates.items():
c[coordinate] = xr.IndexVariable(
dims=[coordinate],
data=np.sort(np.unique(ds['coordinates'].isel(ndim=i))),
attrs={"description":"{0} coordinate".format(desc),"units":"m"}
)
c["t"] = xr.IndexVariable(
dims=["t"],
data=ds.num_time_steps*self.output_dt,
attrs={"description":"time from start of simulation","units":"s"}
)

self.nt = len(c["t"])
self.nx = len(c["x"])
self.ny = len(c["y"])
self.nz = len(c["z"])

coordinates = {"x":(0,"axial","u"), "y":(1,"lateral","v"),"z":(2,"vertical","w")}
v = {}
for coordinate,(i,desc,u) in coordinates.items():
v[u] = xr.DataArray(np.reshape(getattr(ds,"velocity{0}".format(coordinate)).values,(self.nt,self.nx,self.ny,self.nz)),
coords=c,
dims=["t","x","y","z"],
name="{0} velocity".format(desc),
attrs={"description":"velocity along {0}".format(coordinate),"units":"m/s"})

ds = xr.Dataset(data_vars=v, coords=v[u].coords)
ds.attrs = {"original file":self.filename}

self.data = ds


def write(self, filename=None):
""" Rewrite object to file, or write object to `filename` if provided """
if filename:
self.filename = filename
if not self.filename:
raise Exception('No filename provided')
raise NotImplementedError()

def toDataFrame(self):
""" Returns object into one DataFrame, or a dictionary of DataFrames"""
# --- Example (returning one DataFrame):
# return pd.DataFrame(data=np.zeros((10,2)),columns=['Col1','Col2'])
# --- Example (returning dict of DataFrames):
#dfs={}
#cols=['Alpha_[deg]','Cl_[-]','Cd_[-]','Cm_[-]']
#dfs['Polar1'] = pd.DataFrame(data=..., columns=cols)
#dfs['Polar1'] = pd.DataFrame(data=..., columns=cols)
# return dfs
raise NotImplementedError()

# --- Optional functions
def __repr__(self):
""" String that is written to screen when the user calls `print()` on the object.
Provide short and relevant information to save time for the user.
"""
s='<{} object>:\n'.format(type(self).__name__)
s+='|Main attributes:\n'
s+='| - filename: {}\n'.format(self.filename)
# --- Example printing some relevant information for user
#s+='|Main keys:\n'
#s+='| - ID: {}\n'.format(self['ID'])
#s+='| - data : shape {}\n'.format(self['data'].shape)
s+='|Main methods:\n'
s+='| - read, write, toDataFrame, keys'
return s

190 changes: 147 additions & 43 deletions pyFAST/input_output/turbsim_file.py
Original file line number Diff line number Diff line change
Expand Up @@ -58,7 +58,7 @@ def __init__(self, filename=None, **kwargs):
if filename:
self.read(filename, **kwargs)

def read(self, filename=None, header_only=False):
def read(self, filename=None, header_only=False, tdecimals=8):
""" read BTS file, with field:
u (3 x nt x ny x nz)
uTwr (3 x nt x nTwr)
Expand Down Expand Up @@ -98,11 +98,11 @@ def read(self, filename=None, header_only=False):
self['uTwr'] = uTwr
self['info'] = info
self['ID'] = ID
self['dt'] = np.round(dt,3) # dt is stored in single precision in the TurbSim output
self['dt'] = np.round(dt, tdecimals) # dt is stored in single precision in the TurbSim output
self['y'] = np.arange(ny)*dy
self['y'] -= np.mean(self['y']) # y always centered on 0
self['z'] = np.arange(nz)*dz +zBottom
self['t'] = np.round(np.arange(nt)*dt, 3)
self['t'] = np.round(np.arange(nt)*dt, tdecimals)
self['zTwr'] =-np.arange(nTwr)*dz + zBottom
self['zRef'] = zHub
self['uRef'] = uHub
Expand Down Expand Up @@ -592,45 +592,13 @@ def __repr__(self):
s+=' ux: min: {}, max: {}, mean: {} \n'.format(np.min(ux), np.max(ux), np.mean(ux))
s+=' uy: min: {}, max: {}, mean: {} \n'.format(np.min(uy), np.max(uy), np.mean(uy))
s+=' uz: min: {}, max: {}, mean: {} \n'.format(np.min(uz), np.max(uz), np.mean(uz))

s += ' Useful methods:\n'
s += ' - read, write, toDataFrame, keys\n'
s += ' - valuesAt, vertProfile, horizontalPlane, verticalPlane, closestPoint\n'
s += ' - fitPowerLaw\n'
s += ' - makePeriodic, checkPeriodic\n'
return s

def toDataSet(self, datetime=False):
import xarray as xr

if datetime:
timearray = pd.to_datetime(self['t'], unit='s', origin=pd.to_datetime('2000-01-01 00:00:00'))
timestr = 'datetime'
else:
timearray = self['t']
timestr = 'time'

ds = xr.Dataset(
data_vars=dict(
u=([timestr,'y','z'], self['u'][0,:,:,:]),
v=([timestr,'y','z'], self['u'][1,:,:,:]),
w=([timestr,'y','z'], self['u'][2,:,:,:]),
),
coords={
timestr : timearray,
'y' : self['y'],
'z' : self['z'],
},
)

# Add mean computations
ds['up'] = ds['u'] - ds['u'].mean(dim=timestr)
ds['vp'] = ds['v'] - ds['v'].mean(dim=timestr)
ds['wp'] = ds['w'] - ds['w'].mean(dim=timestr)

if datetime:
# Add time (in s) to the variable list
ds['time'] = (('datetime'), self['t'])

return ds



def toDataFrame(self):
dfs={}

Expand Down Expand Up @@ -713,31 +681,167 @@ def toDataFrame(self):
# pass
return dfs

def toDataset(self):
"""
Convert the data that was read in into a xarray Dataset
# TODO SORT OUT THE DIFFERENCE WITH toDataSet
"""
from xarray import IndexVariable, DataArray, Dataset

print('[TODO] pyFAST.input_output.turbsim_file.toDataset: merge with function toDataSet')

y = IndexVariable("y", self.y, attrs={"description":"lateral coordinate","units":"m"})
zround = np.asarray([np.round(zz,6) for zz in self.z]) #the open function here returns something like *.0000000001 which is annoying
z = IndexVariable("z", zround, attrs={"description":"vertical coordinate","units":"m"})
time = IndexVariable("time", self.t, attrs={"description":"time since start of simulation","units":"s"})

da = {}
for component,direction,velname in zip([0,1,2],["x","y","z"],["u","v","w"]):
# the dataset produced here has y/z axes swapped relative to data stored in original object
velocity = np.swapaxes(self["u"][component,...],1,2)
da[velname] = DataArray(velocity,
coords={"time":time,"y":y,"z":z},
dims=["time","y","z"],
name="velocity",
attrs={"description":"velocity along {0}".format(direction),"units":"m/s"})

return Dataset(data_vars=da, coords={"time":time,"y":y,"z":z})

def toDataSet(self, datetime=False):
"""
Convert the data that was read in into a xarray Dataset
# TODO SORT OUT THE DIFFERENCE WITH toDataset
"""
import xarray as xr

print('[TODO] pyFAST.input_output.turbsim_file.toDataSet: should be discontinued')
print('[TODO] pyFAST.input_output.turbsim_file.toDataSet: merge with function toDataset')

if datetime:
timearray = pd.to_datetime(self['t'], unit='s', origin=pd.to_datetime('2000-01-01 00:00:00'))
timestr = 'datetime'
else:
timearray = self['t']
timestr = 'time'

ds = xr.Dataset(
data_vars=dict(
u=([timestr,'y','z'], self['u'][0,:,:,:]),
v=([timestr,'y','z'], self['u'][1,:,:,:]),
w=([timestr,'y','z'], self['u'][2,:,:,:]),
),
coords={
timestr : timearray,
'y' : self['y'],
'z' : self['z'],
},
)

# Add mean computations
ds['up'] = ds['u'] - ds['u'].mean(dim=timestr)
ds['vp'] = ds['v'] - ds['v'].mean(dim=timestr)
ds['wp'] = ds['w'] - ds['w'].mean(dim=timestr)

if datetime:
# Add time (in s) to the variable list
ds['time'] = (('datetime'), self['t'])

return ds

# Useful converters
def fromAMRWind(self, filename, dt, nt):
def fromAMRWind(self, filename, timestep, output_frequency, sampling_identifier, verbose=1, fileout=None, zref=None, xloc=None):
"""
Reads a AMRWind netcdf file, grabs a group of sampling planes (e.g. p_slice),
return an instance of TurbSimFile, optionally write turbsim file to disk
Parameters
----------
filename : str,
full path to netcdf file generated by amrwind
timestep : float,
amr-wind code timestep (time.fixed_dt)
output_frequency : int,
frequency chosen for sampling output in amrwind input file (sampling.output_frequency)
sampling_identifier : str,
identifier of the sampling being requested (an entry of sampling.labels in amrwind input file)
zref : float,
height to be written to turbsim as the reference height. if none is given, it is taken as the vertical centerpoint of the slice
"""
try:
from pyFAST.input_output.amrwind_file import AMRWindFile
except:
try:
from .amrwind_file import AMRWindFile
except:
from amrwind_file import AMRWindFile

obj = AMRWindFile(filename,timestep,output_frequency, group_name=sampling_identifier)

self["u"] = np.ndarray((3,obj.nt,obj.ny,obj.nz))

xloc = float(obj.data.x[0]) if xloc is None else xloc
if verbose:
print("Grabbing the slice at x={0} m".format(xloc))
self['u'][0,:,:,:] = np.swapaxes(obj.data.u.sel(x=xloc).values,1,2)
self['u'][1,:,:,:] = np.swapaxes(obj.data.v.sel(x=xloc).values,1,2)
self['u'][2,:,:,:] = np.swapaxes(obj.data.w.sel(x=xloc).values,1,2)
self['t'] = obj.data.t.values

self['y'] = obj.data.y.values
self['z'] = obj.data.z.values
self['dt'] = obj.output_dt

self['ID'] = 7
ltime = time.strftime('%d-%b-%Y at %H:%M:%S', time.localtime())
self['info'] = 'Converted from AMRWind output file {0} {1:s}.'.format(filename,ltime)

iz = int(obj.nz/2)
self['zRef'] = float(obj.data.z[iz]) if zref is None else zref
if verbose:
print("Setting the TurbSim file reference height to z={0} m".format(self["zRef"]))

self['uRef'] = float(obj.data.u.sel(x=xloc).sel(y=0).sel(z=self["zRef"]).mean().values)
self['zRef'], self['uRef'], bHub = self.hubValues()

if fileout is not None:
filebase = os.path.splitext(filename)[1]
fileout = filebase+".bts"
if verbose:
print("===> {0}".format(fileout))
self.write(fileout)


def fromAMRWind_legacy(self, filename, dt, nt, y, z, sampling_identifier='p_sw2'):
"""
Convert current TurbSim file into one generated from AMR-Wind LES sampling data in .nc format
Assumes:
-- u, v, w (nt, nx * ny * nz)
-- u is aligned with x-axis (flow is not rotated) - this consideration needs to be added
INPUTS:
- filename: (string) full path to .nc sampling data file
- plane_label: (string) name of sampling plane group from .inp file (e.g. "p_sw2")
- sampling_identifier: (string) name of sampling plane group from .inp file (e.g. "p_sw2")
- dt: timestep size [s]
- nt: number of timesteps (sequential) you want to read in, starting at the first timestep available
INPUTS: TODO
- y: user-defined vector of coordinate positions in y
- z: user-defined vector of coordinate positions in z
- uref: (float) reference mean velocity (e.g. 8.0 hub height mean velocity from input file)
- zref: (float) hub height (e.t. 150.0)
"""
import xarray as xr

print('[TODO] fromAMRWind_legacy: function might be unfinished. Merge with fromAMRWind')
print('[TODO] fromAMRWind_legacy: figure out y, and z from data (see fromAMRWind)')

# read in sampling data plane
ds = xr.open_dataset(filename,
engine='netcdf4',
group=plane_label)
group=sampling_identifier)
ny, nz, _ = ds.attrs['ijk_dims']
noffsets = len(ds.attrs['offsets'])
t = np.arange(0, dt*(nt-0.5), dt)
Expand Down

0 comments on commit cd94e85

Please sign in to comment.