From b233dd6db069de33509acc33397318479beb4dd1 Mon Sep 17 00:00:00 2001
From: "Mads M. Pedersen" <mmpe@dtu.dk>
Date: Fri, 28 Jul 2017 08:53:47 +0200
Subject: [PATCH] added constraint_file
---
wetb/wind/turbulence/constraint_file.py | 149 ++++++++++++++++++
.../turbulence/tests/test_constraint_file.py | 114 ++++++++++++++
.../tests/test_files/constraints_test.con | 3 +
3 files changed, 266 insertions(+)
create mode 100644 wetb/wind/turbulence/constraint_file.py
create mode 100644 wetb/wind/turbulence/tests/test_constraint_file.py
create mode 100644 wetb/wind/turbulence/tests/test_files/constraints_test.con
diff --git a/wetb/wind/turbulence/constraint_file.py b/wetb/wind/turbulence/constraint_file.py
new file mode 100644
index 0000000..0d20024
--- /dev/null
+++ b/wetb/wind/turbulence/constraint_file.py
@@ -0,0 +1,149 @@
+import numpy as np
+import os
+from wetb.hawc2.htc_file import HTCFile
+class ConstraintFile(object):
+
+ def __init__(self, center_gl_xyz, box_transport_speed, no_grid_points=(4096, 32, 32), box_size=(6000, 100, 100)):
+ """Generate list of constraints for turbulence constraint simulator
+
+ Parameters
+ ----------
+ center_gl_xyz : (float, float, float)
+ Box front plane center position in global coordinates
+ box_transport_speed : float or int
+ Reference wind speed [m/s] (box transportation speed)
+ no_grid_points : [int, int, int]
+ Number of grid points in x,y,z direction of the turbulence box
+ box_size : [float, float, float]
+ Dimension(length) of box in x,y,z direction [m]
+ """
+ self.center_gl_xyz = center_gl_xyz
+ self.box_transport_speed = box_transport_speed
+ self.no_grid_points = no_grid_points
+ self.box_size = box_size
+ self.dxyz = [d / (n - 1) for d, n in zip(self.box_size, self.no_grid_points)]
+
+ self.constraints = {'u':[], 'v':[], 'w':[]}
+
+ def load(self, filename):
+ """Load constraint from existing constraint file (usable for plotting constraints via time_series()"""
+ with open(filename) as fid:
+ lines = fid.readlines()
+ for l in lines:
+ values = l.split(";")
+ mx,my,mz= map(int, values[:3])
+ comp = values[3]#{"100":'u','010':'v','001':'w'}["".join(values[3:6])]
+ self.constraints[comp].append((mx,my,mz,float(values[-1])))
+
+
+ def add_constraints(self, glpos, tuvw, subtract_mean=True, fail_outside_box=True, nearest=True):
+ """Add constraints to constraint file
+
+ parameters
+ ----------
+ glpos : (float or array_like, float or array_like, float or array_like)
+ global position(s), (x,y,z) of measurement point point(s)\n
+ x: horizontal left seen in direction of mean wind, y: direction of mean wind, z: vertical down\n
+ tuvw: array_like (shape: no_obs x (1+no_components))
+ time and u[, v[, w]] components of wind at measurement point\n
+ v and w components are optional\n
+ subtract_mean : boolean, optional
+ if True, default, the mean values are subtracted from u,v,w
+ fail_outside_box : boolean, optional
+ if True, default, an error is raised if any positions are outside the box\n
+ if False, mann coordinates modulo N is used
+ nearest : boolean, optional
+ if True, the wsp at the position closest to the mann grid point are applied\n
+ if False, the average of all the wsp, that map to the mann grid point are applied
+ """
+ nx, ny, nz = self.no_grid_points
+ dx, dy, dz = self.dxyz
+ center_x, center_y, center_z = self.center_gl_xyz
+
+ time, u, v, w = (list(tuvw.T) + [None, None])[:4]
+ x, y, z = glpos
+ mxs = np.round((time * self.box_transport_speed + (center_y - y)) / dx).astype(np.int)
+ if not np.array(y).shape==mxs.shape:
+ y = np.zeros_like(mxs)+y
+ if not np.array(x).shape==mxs.shape:
+ x = np.zeros_like(mxs)+x
+ if not np.array(z).shape==mxs.shape:
+ z = np.zeros_like(mxs)+z
+ x_err = y-(-(mxs*dx-time*self.box_transport_speed-center_y))
+ mys = np.round(((ny - 1) / 2. + (-x+center_x) / dy)).astype(np.int)
+ mzs = np.round(((nz - 1) / 2. + (center_z - z) / dz)).astype(np.int)
+
+ y_err = x-(-((mys-(ny-1)/2)*dy-center_x))
+ z_err = z-(-(mzs-(nz-1)/2)*dz+center_z)
+ pos_err = np.sqrt(x_err**2+y_err**2+z_err**2)
+ if fail_outside_box:
+ for ms, n in zip([mxs,mys,mzs],self.no_grid_points):
+ if ms.min() + 1 < 1:
+ i = np.argmin(ms)
+ raise ValueError("At time, t=%s, global position (%s,%s,%s) maps to grid point (%d,%d,%d) whichs is outside turbulence box, range (1..%d,1..%d,1..%d)" % (time[i], x[i], y[i], z[i], mxs[i] + 1, mys[i] + 1, mzs[i] + 1, nx , ny , nz))
+ if ms.max() + 1 > n:
+ i = np.argmax(ms)
+ raise ValueError("At time, t=%s, global position (%s,%s,%s) maps to grid point (%d,%d,%d) whichs is outside turbulence box, range (1..%d,1..%d,1..%d)" % (time[i], x[i], y[i], z[i], mxs[i] + 1, mys[i] + 1, mzs[i] + 1, nx , ny , nz))
+ else:
+ mxs%=nx
+ mys%=ny
+ mzs%=nz
+ mann_pos = ["%06d%02d%02d"%(mx,my,mz) for mx,my,mz in zip(mxs,mys,mzs)]
+ i = 0
+ N = len(mann_pos)
+ uvw_lst = [uvw for uvw in [u,v,w] if uvw is not None]
+ if subtract_mean:
+ uvw_lst = [uvw - uvw.mean() for uvw in uvw_lst]
+ while i<N:
+ for j in range(i+1,N+1):
+ if j==N or mann_pos[i]!=mann_pos[j]:
+ break
+ for comp, wsp in zip(['u', 'v', 'w'], uvw_lst):
+ if nearest:
+ value = wsp[i:j][np.argmin(pos_err[i:j])]
+ else: #mean
+ value = wsp[i:j].mean()
+ self.constraints[comp].append((mxs[i] + 1, mys[i] + 1, mzs[i] + 1, value))
+ i=j
+
+
+ def __str__(self):
+ return "\n".join(["%d;%d;%d;%s;%.10f" % (mx, my, mz, comp, value) for comp in ['u', 'v', 'w']
+ for mx, my, mz, value in self.constraints[comp]])
+
+ def save(self, path, name, folder="./constraints/"):
+ path = os.path.join(path, folder)
+ if not os.path.isdir(path):
+ os.makedirs(path, exist_ok=True)
+ with open(os.path.join(path, "%s.con" % name), 'w') as fid:
+ fid.write(str(self))
+
+ def time_series(self, y_position=0):
+ time = (np.arange(self.no_grid_points[0])*float(self.dxyz[0]) + y_position)/ self.box_transport_speed
+ u,v,w = [(np.zeros_like(time)+ np.nan)]*3
+ for uvw, constr in zip([u,v,w], [np.array(self.constraints[uvw]) for uvw in 'uvw']):
+ if constr.shape[0]>0:
+ uvw[constr[:,0].astype(np.int)-1] =constr[:,3]
+ return time, np.array([u,v,w]).T
+
+
+ def simulation_cmd(self, ae23, L, Gamma, seed, name, folder="./constraints/"):
+ assert isinstance(seed, int) and seed > 0, "seed must be a positive integer"
+ cmd = "./csimu2.exe %d %d %d %.3f %.3f %.3f " % (self.no_grid_points + self.box_size)
+ cmd += "%.6f %.2f %.2f %d " % (ae23, L, Gamma, seed)
+ cmd += "./turb/%s_s%s %s/%s.con" % (name, seed, folder, name)
+ return cmd
+
+ def hawc2_mann_section(self, name, seed):
+ htc = HTCFile()
+ mann = htc.wind.add_section("mann")
+ for uvw in 'uvw':
+ setattr(mann, 'filename_%s'%uvw, "./turb/%s_s%04d%s.bin" % ( name, int(seed),uvw))
+ for uvw, l,n in zip('uvw', self.box_size, self.no_grid_points):
+ setattr(mann, 'box_dim_%s'%uvw, [n,l/n])
+ mann.dont_scale = 1
+ return mann
+
+
+
+
diff --git a/wetb/wind/turbulence/tests/test_constraint_file.py b/wetb/wind/turbulence/tests/test_constraint_file.py
new file mode 100644
index 0000000..40233f6
--- /dev/null
+++ b/wetb/wind/turbulence/tests/test_constraint_file.py
@@ -0,0 +1,114 @@
+'''
+Created on 02/10/2015
+
+@author: MMPE
+'''
+import unittest
+
+import numpy as np
+from wetb.wind.turbulence.constraint_file import ConstraintFile
+from wetb.utils.test_files import get_test_file
+class TestConstraintFile(unittest.TestCase):
+
+
+ def testConstraintFile(self):
+ time = [0, 1, 2]
+ u = [5, 7, 9]
+ v = [1.5, 1, .5]
+ w = [3, 1, 5]
+ tu = np.array([time, u]).T
+ tuvw = np.array([time, u, v, w]).T
+ glpos = (0, 0, -85)
+ constraint_file = ConstraintFile(center_gl_xyz=(-5, 0, -90), box_transport_speed=10, no_grid_points=(16, 8, 8), box_dimension=(30, 70, 70))
+ constraint_file.add_constraints(glpos, tu)
+ np.testing.assert_array_equal(constraint_file.constraints['u'][1], (6, 4, 4, 0))
+ constraint_file = ConstraintFile(center_gl_xyz=(-5, 0, -90), box_transport_speed=5, no_grid_points=(16, 8, 8), box_dimension=(30, 70, 70))
+ constraint_file.add_constraints(glpos, tuvw)
+ np.testing.assert_array_equal(constraint_file.constraints['u'][2], (6, 4, 4, 2))
+ np.testing.assert_array_equal(constraint_file.constraints['v'][2], (6, 4, 4, -.5))
+ np.testing.assert_array_equal(constraint_file.constraints['w'][2], (6, 4, 4, 2))
+ constraint_file = ConstraintFile(center_gl_xyz=(-5, 0, -90), box_transport_speed=10, no_grid_points=(16, 8, 8), box_dimension=(30, 70, 70))
+ constraint_file.add_constraints((-10, -4, -95), tu)
+ np.testing.assert_array_equal(constraint_file.constraints['u'][1], (8, 5, 5, 0))
+ constraint_file = ConstraintFile(center_gl_xyz=(-5, 4, -90), box_transport_speed=10, no_grid_points=(16, 8, 8), box_dimension=(30, 70, 70))
+ constraint_file.add_constraints((-10, 0, -95), tu)
+ np.testing.assert_array_equal(constraint_file.constraints['u'][1], (8, 5, 5, 0))
+ #constraint_file.save("./test_files", "constraints_test",'')
+
+ def testConstraintFile_Mxyz(self):
+ time = [0]
+ u = [5]
+ tu = np.array([time, u]).T
+
+ glpos = (0, 0, -85)
+ constraint_file = ConstraintFile(center_gl_xyz=(0, 0, 0), box_transport_speed=10, no_grid_points=(16, 8, 8), box_dimension=(30, 70, 70))
+ constraint_file.add_constraints([0,0,0], tu, subtract_mean=False)
+ np.testing.assert_array_equal(constraint_file.constraints['u'][-1], (1, 5, 5, 5))
+ constraint_file.add_constraints([5,0,0], tu, subtract_mean=False)
+ np.testing.assert_array_equal(constraint_file.constraints['u'][-1], (1, 4, 5, 5))
+ constraint_file.add_constraints([0,0,5], tu, subtract_mean=False)
+ np.testing.assert_array_equal(constraint_file.constraints['u'][-1], (1, 5, 4, 5))
+ constraint_file.add_constraints([0,-10,0], tu, subtract_mean=False)
+ np.testing.assert_array_equal(constraint_file.constraints['u'][-1], (6, 5, 5, 5))
+
+ def testConstraintFile_center_gl_xyz(self):
+ time = [0]
+ u = [5]
+ tu = np.array([time, u]).T
+
+ constraint_file = ConstraintFile(center_gl_xyz=(0, 0, 0), box_transport_speed=10, no_grid_points=(16, 8, 8), box_dimension=(30, 70, 70))
+ constraint_file.add_constraints([0,0,0], tu, subtract_mean=False)
+ np.testing.assert_array_equal(constraint_file.constraints['u'][-1], (1, 5, 5, 5))
+
+ constraint_file = ConstraintFile(center_gl_xyz=(0, 0, -5), box_transport_speed=10, no_grid_points=(16, 8, 8), box_dimension=(30, 70, 70))
+ constraint_file.add_constraints([0,0,0], tu, subtract_mean=False)
+ np.testing.assert_array_equal(constraint_file.constraints['u'][-1], (1, 5, 4, 5))
+
+ constraint_file = ConstraintFile(center_gl_xyz=(-5, 0, 0), box_transport_speed=10, no_grid_points=(16, 8, 8), box_dimension=(30, 70, 70))
+ constraint_file.add_constraints([0,0,0], tu, subtract_mean=False)
+ np.testing.assert_array_equal(constraint_file.constraints['u'][-1], (1, 4, 5, 5))
+
+ constraint_file = ConstraintFile(center_gl_xyz=(0, 10, 0), box_transport_speed=10, no_grid_points=(16, 8, 8), box_dimension=(30, 70, 70))
+ constraint_file.add_constraints([0,0,0], tu, subtract_mean=False)
+ np.testing.assert_array_equal(constraint_file.constraints['u'][-1], (6, 5, 5, 5))
+
+
+ def test_outputbox_errors(self):
+ time = [0, 1, 2]
+ u = [5, 7, 9]
+ v = [1.5, 1, .5]
+ w = [3, 1, 5]
+ tu = np.array([time, u]).T
+ tuvw = np.array([time, u, v, w]).T
+ glpos = (0, 0, -85)
+ constraint_file = ConstraintFile(center_gl_xyz=(0, 0, -85), box_transport_speed=10, no_grid_points=(16, 8, 8), box_dimension=(30, 70, 70))
+ self.assertRaisesRegex(ValueError, "At time, t=0, global position \(0,0,-44\)", constraint_file.add_constraints, (0, 0, -85 + 35 + 6), tu)
+ self.assertRaisesRegex(ValueError, "At time, t=0, global position \(0,0,-126\)", constraint_file.add_constraints, (0, 0, -85 - 35 - 6), tu)
+ self.assertRaisesRegex(ValueError, "At time, t=0, global position \(-41,0,-85\)", constraint_file.add_constraints, (-35 - 6, 0, -85), tu)
+ self.assertRaisesRegex(ValueError, "At time, t=0, global position \(41,0,-85\)", constraint_file.add_constraints, (35 + 6, 0, -85), tu)
+ self.assertRaisesRegex(ValueError, "At time, t=0, global position \(0,2,-85\)", constraint_file.add_constraints, (0, 2, -85), tu)
+ self.assertRaisesRegex(ValueError, "At time, t=2, global position \(0,-13,-85\)", constraint_file.add_constraints, (0, -13, -85), tu)
+
+ def test_hawc2_cmd(self):
+ constraint_file = ConstraintFile(center_gl_xyz=(0, 0, -85), box_transport_speed=10, no_grid_points=(16, 8, 8), box_dimension=(30, 70, 70))
+ mann = constraint_file.hawc2_mann_section("test",1)
+ self.assertEqual(mann.box_dim_u.values, [16,1.875])
+
+ def test_load(self):
+ constraint_file = ConstraintFile(center_gl_xyz=(-5, 4, -90), box_transport_speed=10, no_grid_points=(16, 8, 8), box_dimension=(30, 70, 70))
+ constraint_file.load(get_test_file('constraints_test.con'))
+ with open(get_test_file('constraints_test.con')) as fid:
+ self.assertEqual(str(constraint_file), fid.read())
+
+ def test_time_series(self):
+ constraint_file = ConstraintFile(center_gl_xyz=(-5, 4, -90), box_transport_speed=10, no_grid_points=(16, 8, 8), box_dimension=(30, 70, 70))
+ constraint_file.load(get_test_file('constraints_test.con'))
+ time,uvw = constraint_file.time_series(-4)
+ u,v,w = uvw.T
+ m = ~np.isnan(u)
+ np.testing.assert_array_equal(time[m],[0,1,2])
+ np.testing.assert_array_equal(u[m],[-2,0,2])
+
+if __name__ == "__main__":
+ #import sys;sys.argv = ['', 'Test.testConstraintFile']
+ unittest.main()
diff --git a/wetb/wind/turbulence/tests/test_files/constraints_test.con b/wetb/wind/turbulence/tests/test_files/constraints_test.con
new file mode 100644
index 0000000..600a34d
--- /dev/null
+++ b/wetb/wind/turbulence/tests/test_files/constraints_test.con
@@ -0,0 +1,3 @@
+3;5;5;u;-2.0000000000
+8;5;5;u;0.0000000000
+13;5;5;u;2.0000000000
\ No newline at end of file
--
GitLab