Commit 1f66ed0b authored by Mads M. Pedersen's avatar Mads M. Pedersen
Browse files
parents e8e82a29 23f27eaf
......@@ -326,130 +326,4 @@ class HTCSensor(HTCLine):
("", "\t" + self.str_values())[bool(self.values)],
("", "\t" + self.comments)[bool(self.comments.strip())])
class HTCDefaults(object):
empty_htc = """begin simulation;
time_stop 600;
solvertype 1; (newmark)
on_no_convergence continue;
convergence_limits 1E3 1.0 1E-7; ; . to run again, changed 07/11
begin newmark;
deltat 0.02;
end newmark;
end simulation;
;
;----------------------------------------------------------------------------------------------------------------------------------------------------------------
;
begin new_htc_structure;
begin orientation;
end orientation;
begin constraint;
end constraint;
end new_htc_structure;
;
;----------------------------------------------------------------------------------------------------------------------------------------------------------------
;
begin wind ;
density 1.225 ;
wsp 10 ;
tint 1;
horizontal_input 1 ; 0=false, 1=true
windfield_rotations 0 0.0 0.0 ; yaw, tilt, rotation
center_pos0 0 0 -30 ; hub heigth
shear_format 1 0;0=none,1=constant,2=log,3=power,4=linear
turb_format 0 ; 0=none, 1=mann,2=flex
tower_shadow_method 0 ; 0=none, 1=potential flow, 2=jet
end wind;
;
;----------------------------------------------------------------------------------------------------------------------------------------------------------------
;
begin dll;
end dll;
;
;----------------------------------------------------------------------------------------------------------------------------------------------------------------
;
begin output;
general time;
end output;
exit;"""
def add_mann_turbulence(self, L=29.4, ae23=1, Gamma=3.9, seed=1001, high_frq_compensation=True,
filenames=None,
no_grid_points=(16384, 32, 32), box_dimension=(6000, 100, 100),
dont_scale=False,
std_scaling=None):
wind = self.add_section('wind')
wind.turb_format = 1
mann = wind.add_section('mann')
if 'create_turb_parameters' in mann:
mann.create_turb_parameters.values = [L, ae23, Gamma, seed, int(high_frq_compensation)]
else:
mann.add_line('create_turb_parameters', [L, ae23, Gamma, seed, int(high_frq_compensation)], "L, alfaeps, gamma, seed, highfrq compensation")
if filenames is None:
fmt = "mann_l%.1f_ae%.2f_g%.1f_h%d_%dx%dx%d_%.3fx%.2fx%.2f_s%04d%c.turb"
import numpy as np
dxyz = tuple(np.array(box_dimension) / no_grid_points)
filenames = ["./turb/" + fmt % ((L, ae23, Gamma, high_frq_compensation) + no_grid_points + dxyz + (seed, uvw)) for uvw in ['u', 'v', 'w']]
if isinstance(filenames, str):
filenames = ["./turb/%s_s%04d%s.bin" % (filenames, seed, c) for c in ['u', 'v', 'w']]
for filename, c in zip(filenames, ['u', 'v', 'w']):
setattr(mann, 'filename_%s' % c, filename)
for c, n, dim in zip(['u', 'v', 'w'], no_grid_points, box_dimension):
setattr(mann, 'box_dim_%s' % c, "%d %.4f" % (n, dim / (n - 1)))
if dont_scale:
mann.dont_scale = 1
else:
try:
del mann.dont_scale
except KeyError:
pass
if std_scaling is not None:
mann.std_scaling = "%f %f %f" % std_scaling
else:
try:
del mann.std_scaling
except KeyError:
pass
def add_turb_export(self, filename="export_%s.turb", samplefrq = None):
exp = self.wind.add_section('turb_export', allow_duplicate=True)
for uvw in 'uvw':
exp.add_line('filename_%s'%uvw, [filename%uvw])
sf = samplefrq or max(1,int( self.wind.mann.box_dim_u[1]/(self.wind.wsp[0] * self.deltat())))
exp.samplefrq = sf
if "time" in self.output:
exp.time_start = self.output.time[0]
else:
exp.time_start = 0
exp.nsteps = (self.simulation.time_stop[0]-exp.time_start[0]) / self.deltat()
for vw in 'vw':
exp.add_line('box_dim_%s'%vw, self.wind.mann['box_dim_%s'%vw].values)
def import_dtu_we_controller_input(self, filename):
dtu_we_controller = [dll for dll in self.dll if dll.name[0] == 'dtu_we_controller'][0]
with open (filename) as fid:
lines = fid.readlines()
K_r1 = float(lines[1].replace("K = ", '').replace("[Nm/(rad/s)^2]", ''))
Kp_r2 = float(lines[4].replace("Kp = ", '').replace("[Nm/(rad/s)]", ''))
Ki_r2 = float(lines[5].replace("Ki = ", '').replace("[Nm/rad]", ''))
Kp_r3 = float(lines[7].replace("Kp = ", '').replace("[rad/(rad/s)]", ''))
Ki_r3 = float(lines[8].replace("Ki = ", '').replace("[rad/rad]", ''))
KK = lines[9].split("]")
KK1 = float(KK[0].replace("K1 = ", '').replace("[deg", ''))
KK2 = float(KK[1].replace(", K2 = ", '').replace("[deg^2", ''))
cs = dtu_we_controller.init
cs.constant__11.values[1] = "%.6E" % K_r1
cs.constant__12.values[1] = "%.6E" % Kp_r2
cs.constant__13.values[1] = "%.6E" % Ki_r2
cs.constant__16.values[1] = "%.6E" % Kp_r3
cs.constant__17.values[1] = "%.6E" % Ki_r3
cs.constant__21.values[1] = "%.6E" % KK1
cs.constant__22.values[1] = "%.6E" % KK2
'''
Created on 20/01/2014
@author: MMPE
See documentation of HTCFile below
'''
from __future__ import division
from __future__ import unicode_literals
from __future__ import print_function
from __future__ import absolute_import
from builtins import zip
from builtins import int
from builtins import str
from future import standard_library
import os
from wetb.wind.shear import log_shear, power_shear
standard_library.install_aliases()
class HTCDefaults(object):
empty_htc = """begin simulation;
time_stop 600;
solvertype 1; (newmark)
on_no_convergence continue;
convergence_limits 1E3 1.0 1E-7; ; . to run again, changed 07/11
begin newmark;
deltat 0.02;
end newmark;
end simulation;
;
;----------------------------------------------------------------------------------------------------------------------------------------------------------------
;
begin new_htc_structure;
begin orientation;
end orientation;
begin constraint;
end constraint;
end new_htc_structure;
;
;----------------------------------------------------------------------------------------------------------------------------------------------------------------
;
begin wind ;
density 1.225 ;
wsp 10 ;
tint 1;
horizontal_input 1 ; 0=false, 1=true
windfield_rotations 0 0.0 0.0 ; yaw, tilt, rotation
center_pos0 0 0 -30 ; hub heigth
shear_format 1 0;0=none,1=constant,2=log,3=power,4=linear
turb_format 0 ; 0=none, 1=mann,2=flex
tower_shadow_method 0 ; 0=none, 1=potential flow, 2=jet
end wind;
;
;----------------------------------------------------------------------------------------------------------------------------------------------------------------
;
begin dll;
end dll;
;
;----------------------------------------------------------------------------------------------------------------------------------------------------------------
;
begin output;
general time;
end output;
exit;"""
def add_mann_turbulence(self, L=29.4, ae23=1, Gamma=3.9, seed=1001, high_frq_compensation=True,
filenames=None,
no_grid_points=(16384, 32, 32), box_dimension=(6000, 100, 100),
dont_scale=False,
std_scaling=None):
wind = self.add_section('wind')
wind.turb_format = 1
mann = wind.add_section('mann')
if 'create_turb_parameters' in mann:
mann.create_turb_parameters.values = [L, ae23, Gamma, seed, int(high_frq_compensation)]
else:
mann.add_line('create_turb_parameters', [L, ae23, Gamma, seed, int(high_frq_compensation)], "L, alfaeps, gamma, seed, highfrq compensation")
if filenames is None:
fmt = "mann_l%.1f_ae%.2f_g%.1f_h%d_%dx%dx%d_%.3fx%.2fx%.2f_s%04d%c.turb"
import numpy as np
dxyz = tuple(np.array(box_dimension) / no_grid_points)
filenames = ["./turb/" + fmt % ((L, ae23, Gamma, high_frq_compensation) + no_grid_points + dxyz + (seed, uvw)) for uvw in ['u', 'v', 'w']]
if isinstance(filenames, str):
filenames = ["./turb/%s_s%04d%s.bin" % (filenames, seed, c) for c in ['u', 'v', 'w']]
for filename, c in zip(filenames, ['u', 'v', 'w']):
setattr(mann, 'filename_%s' % c, filename)
for c, n, dim in zip(['u', 'v', 'w'], no_grid_points, box_dimension):
setattr(mann, 'box_dim_%s' % c, "%d %.4f" % (n, dim / (n - 1)))
if dont_scale:
mann.dont_scale = 1
else:
try:
del mann.dont_scale
except KeyError:
pass
if std_scaling is not None:
mann.std_scaling = "%f %f %f" % std_scaling
else:
try:
del mann.std_scaling
except KeyError:
pass
def add_turb_export(self, filename="export_%s.turb", samplefrq = None):
exp = self.wind.add_section('turb_export', allow_duplicate=True)
for uvw in 'uvw':
exp.add_line('filename_%s'%uvw, [filename%uvw])
sf = samplefrq or max(1,int( self.wind.mann.box_dim_u[1]/(self.wind.wsp[0] * self.deltat())))
exp.samplefrq = sf
if "time" in self.output:
exp.time_start = self.output.time[0]
else:
exp.time_start = 0
exp.nsteps = (self.simulation.time_stop[0]-exp.time_start[0]) / self.deltat()
for vw in 'vw':
exp.add_line('box_dim_%s'%vw, self.wind.mann['box_dim_%s'%vw].values)
def import_dtu_we_controller_input(self, filename):
dtu_we_controller = [dll for dll in self.dll if dll.name[0] == 'dtu_we_controller'][0]
with open (filename) as fid:
lines = fid.readlines()
K_r1 = float(lines[1].replace("K = ", '').replace("[Nm/(rad/s)^2]", ''))
Kp_r2 = float(lines[4].replace("Kp = ", '').replace("[Nm/(rad/s)]", ''))
Ki_r2 = float(lines[5].replace("Ki = ", '').replace("[Nm/rad]", ''))
Kp_r3 = float(lines[7].replace("Kp = ", '').replace("[rad/(rad/s)]", ''))
Ki_r3 = float(lines[8].replace("Ki = ", '').replace("[rad/rad]", ''))
KK = lines[9].split("]")
KK1 = float(KK[0].replace("K1 = ", '').replace("[deg", ''))
KK2 = float(KK[1].replace(", K2 = ", '').replace("[deg^2", ''))
cs = dtu_we_controller.init
cs.constant__11.values[1] = "%.6E" % K_r1
cs.constant__12.values[1] = "%.6E" % Kp_r2
cs.constant__13.values[1] = "%.6E" % Ki_r2
cs.constant__16.values[1] = "%.6E" % Kp_r3
cs.constant__17.values[1] = "%.6E" % Ki_r3
cs.constant__21.values[1] = "%.6E" % KK1
cs.constant__22.values[1] = "%.6E" % KK2
class HTCExtensions(object):
def get_shear(self):
shear_type, parameter = self.wind.shear_format.values
z0 = -self.wind.center_pos0[2]
wsp = self.wind.wsp[0]
if shear_type==1: #constant
return lambda z : parameter
elif shear_type==3:
return power_shear(parameter, z0, wsp)
else:
raise NotImplementedError
\ No newline at end of file
......@@ -18,8 +18,8 @@ from wetb.utils.cluster_tools.cluster_resource import unix_path_old
standard_library.install_aliases()
from collections import OrderedDict
from wetb.hawc2.htc_contents import HTCContents, HTCSection, HTCLine, \
HTCDefaults
from wetb.hawc2.htc_contents import HTCContents, HTCSection, HTCLine
from wetb.hawc2.htc_extensions import HTCDefaults, HTCExtensions
import os
from copy import copy
......@@ -27,7 +27,7 @@ from copy import copy
def fmt_path(path):
return path.lower().replace("\\","/")
class HTCFile(HTCContents, HTCDefaults):
class HTCFile(HTCContents, HTCDefaults, HTCExtensions):
"""Wrapper for HTC files
Examples:
......
......@@ -10,62 +10,168 @@ from __future__ import absolute_import
from builtins import range
from io import open
from future import standard_library
from wetb.hawc2.htc_file import HTCFile
standard_library.install_aliases()
import numpy as np
import os
def save(filename, x_coordinates, z_coordinates, u=None, v=None, w=None):
"""
Parameters
----------
filename : str
filename
x_coordinates : array_like
lateral coordinates
z_coordinates : array_like
vertical coordinates
u : array_like, optional
shear_u component, normalized with U_mean\n
shape must be (#z_coordinates, #x_coordinates) or (#z_coordinates,)
v : array_like, optional
shear_v component, normalized with U_mean\n
shape must be (#z_coordinates, #x_coordinates) or (#z_coordinates,)
w : array_like, optional
shear_w component, normalized with U_mean\n
shape must be (#z_coordinates, #x_coordinates) or (#z_coordinates,)
"""
shape = (len(z_coordinates), len(x_coordinates))
vuw = [v, u, w]
for i in range(3):
if vuw[i] is None:
if i == 1:
vuw[i] = np.ones((shape))
class ShearFile(object):
"""HAWC2 user defined shear file
Examples:
---------
>>> sf = ShearFile([-55, 55], [30, 100, 160] , u=np.array([[0.7, 1, 1.3], [0.7, 1, 1.3]]).T)
>>> print (sf.uvw([-55,0],[65,135])) #uvw factors
[array([ 0.85 , 1.175]), array([ 0., 0.]), array([ 0., 0.])]
>>> from wetb.wind.shear import power_shear
>>> print (sf.uvw([-55,0],[65,135], shear=power_shear(.2,100,10))) # uvw wind speeds
[array([ 7.79832978, 12.47684042]), array([ 0., 0.]), array([ 0., 0.])]
>>> sf.save('test.dat')
"""
def __init__(self,v_positions, w_positions, u=None, v=None, w=None, shear=None):
"""
Parameters
----------
v_positions : array_like
lateral coordinates
w_positions : array_like
vertical coordinates
u : array_like, optional
shear_u component, normalized with U_mean\n
shape must be (#w_positions, #v_positions) or (#w_positions,)
v : array_like, optional
shear_v component, normalized with U_mean\n
shape must be (#w_positions, #v_positions) or (#w_positions,)
w : array_like, optional
shear_w component, normalized with U_mean\n
shape must be (#w_positions, #v_positions) or (#w_positions,)
"""
self.v_positions = v_positions
self.w_positions = w_positions
shape = (len(w_positions), len(v_positions))
uvw = [u, v, w]
for i in range(3):
if uvw[i] is None:
if i == 0:
uvw[i] = np.ones((shape))
else:
uvw[i] = np.zeros((shape))
else:
vuw[i] = np.zeros((shape))
else:
vuw[i] = np.array(vuw[i])
if len(vuw[i].shape) == 1 and vuw[i].shape[0] == shape[0]:
vuw[i] = np.repeat(np.atleast_2d(vuw[i]).T, shape[1], 1)
uvw[i] = np.array(uvw[i])
if len(uvw[i].shape) == 1 and uvw[i].shape[0] == shape[0]:
uvw[i] = np.repeat(np.atleast_2d(uvw[i]).T, shape[1], 1)
assert uvw[i].shape == shape, (i, uvw[i].shape, shape)
self.u, self.v, self.w = uvw
self.shear = shear
def uvw(self, v,w,shear=None):
"""Calculate u,v,w wind speeds at position(s) (v,w)
Parameters
----------
v : int, float or array_like
v-coordinate(s)
w : int, float or array_like
w-coordinates(s)
shear : function or None
if function: f(height)->wsp
if None: self.shear is used, if not None.
Otherwise wind speed factors instead of absolute wind speeds are returned
Returns
-------
u,v,w
wind speed(s) or wind speed factor(s) if shear not defined
"""
shear = shear or self.shear or (lambda z: 1)
from scipy.interpolate import RegularGridInterpolator
wv = np.array([w,v]).T
return [RegularGridInterpolator((self.w_positions, self.v_positions), uvw)(wv)*shear(w)
for uvw in [self.u, self.v, self.w] if uvw is not None]
def save(self, filename):
"""Save user defined shear file
Parameters
----------
filename : str:
Filename
"""
# exist_ok does not exist in Python27
filename = os.path.abspath(filename)
if not os.path.exists(os.path.dirname(filename)):
os.makedirs(os.path.dirname(filename))#, exist_ok=True)
with open(filename, 'w', encoding='utf-8') as fid:
fid.write(" # autogenerated shear file\n")
fid.write(" %d %d\n" % (len(self.v_positions), len(self.w_positions)))
for i, (l,vuw) in enumerate(zip(['v', 'u', 'w'],[self.v,self.u,self.w])):
fid.write(" # shear %s component\n " % l)
fid.write("\n ".join([" ".join(["%.10f" % v for v in r ]) for r in vuw]))
fid.write("\n")
for yz, coor in (['v', self.v_positions], ['w', self.w_positions]):
fid.write(" # %s coordinates\n " % yz)
fid.write("\n ".join("%.10f" % v for v in coor))
fid.write("\n")
@staticmethod
def load(filename):
"""Load shear file
Parameters
----------
filename : str
Filename
Returns
-------
shear file : ShearFile-object
"""
with open(filename) as fid:
lines = fid.readlines()
no_V,no_W=map(int,lines[1].split())
v,u,w = [np.array([row.split() for row in lines[3+(no_W+1)*i:3+(no_W+1)*(i+1)-1] ],dtype=float) for i in range(3)]
v_positions = np.array(lines[3+(no_W+1)*3:3+(no_W+1)*3+no_V],dtype=float)
w_positions = np.array(lines[3+(no_W+1)*3+(no_V+1):3+(no_W+1)*3+(no_V+1)+no_W],dtype=float)
return ShearFile(v_positions,w_positions,u,v,w)
@staticmethod
def load_from_htc(htc_file):
"""Load shear file from HTC file including shear function
Parameters
----------
htc_file : str or HTCFile
Filename or HTCFile
Returns
-------
shear file : ShearFile-object
"""
if isinstance(htc_file,str):
htc_file = HTCFile(htc_file)
user_defined_shear_filename = os.path.join(htc_file.modelpath, htc_file.wind.user_defined_shear[0])
shear_file = ShearFile.load(user_defined_shear_filename)
shear_file.shear = htc_file.get_shear()
return shear_file
def save(filename, v_coordinates, w_coordinates, u=None, v=None, w=None):
"""Save shear file (deprecated)"""
ShearFile(v_coordinates, w_coordinates,u,v,w).save(filename)
assert vuw[i].shape == shape, (i, vuw[i].shape, shape)
# exist_ok does not exist in Python27
if not os.path.exists(os.path.dirname(filename)):
os.makedirs(os.path.dirname(filename))#, exist_ok=True)
with open(filename, 'w', encoding='utf-8') as fid:
fid.write(" # autogenerated shear file\n")
fid.write(" %d %d\n" % (len(x_coordinates), len(z_coordinates)))
for i, l in enumerate(['v', 'u', 'w']):
fid.write(" # shear %s component\n " % l)
fid.write("\n ".join([" ".join(["%.10f" % v for v in r ]) for r in vuw[i]]))
fid.write("\n")
for yz, coor in (['v', x_coordinates], ['w', z_coordinates]):
fid.write(" # %s coordinates\n " % yz)
fid.write("\n ".join("%.10f" % v for v in coor))
fid.write("\n")
if __name__ == "__main__":
save("test.dat", [-55, 55], [30, 100, 160] , u=np.array([[0.7, 1, 1.3], [0.7, 1, 1.3]]).T)
from wetb.wind.shear import power_shear
sf = ShearFile([-55, 55], [30, 100, 160] , u=np.array([[0.7, 1, 1.3], [0.7, 1, 1.3]]).T)
print (sf.uvw([-55,0],[65,135])) #uvw factors
#[array([ 0.85 , 1.175]), array([ 0., 0.]), array([ 0., 0.])]
print (sf.uvw([-55,0],[65,135], shear=power_shear(.2,100,10))) # uvw wind speeds
#[array([ 7.79832978, 12.47684042]), array([ 0., 0.]), array([ 0., 0.])]
sf.save('test.dat')
# autogenerated shear file
2 3
# shear v component
0.0000000000 0.0000000000
0.0000000000 0.0000000000
0.0000000000 0.0000000000
# shear u component
0.8000000000 0.6000000000
1.0000000000 1.0000000000
1.2000000000 1.4000000000
# shear w component
0.0000000000 0.0000000000
0.0000000000 0.0000000000
0.0000000000 0.0000000000
# v coordinates
-55.0000000000
55.0000000000
# w coordinates
30.0000000000
100.0000000000
160.0000000000
......@@ -260,6 +260,7 @@ begin wind;
tower_shadow_method 3; 0=none, 1=potential flow, 2=jet
scale_time_start 0;
wind_ramp_factor 0.0 100 0.8 1.0;
user_defined_shear ./data/user_shear.dat;
; iec_gust;
;
; wind_ramp_abs 400.0 401.0 0.0 1.0; wsp. after the step: 5.0
......
'''
Created on 17/07/2014
@author: MMPE
'''
from __future__ import print_function
from __future__ import unicode_literals
from __future__ import division
from __future__ import absolute_import
from io import open
from builtins import str
from builtins import zip
from future import standard_library
standard_library.install_aliases()
import os
import unittest
from datetime import datetime
from wetb.hawc2.htc_file import HTCFile, HTCLine
import numpy as np
tfp = os.path.join(os.path.dirname(__file__), 'test_files/htcfiles/') # test file path
class TestHtcFile(unittest.TestCase):
def test_get_shear(self):
htc = HTCFile(tfp+'test.htc')
self.assertEqual(htc.get_shear()(100), 10*(100/119)**.2)
if __name__ == "__main__":
#import sys;sys.argv = ['', 'Test.testName']
unittest.main()