Skip to content
Snippets Groups Projects
Commit 182c4d61 authored by aschubig's avatar aschubig
Browse files

Merge remote-tracking branch 'origin/development' into development

parents b31be69b a8e10008
Branches
No related tags found
1 merge request!1Development
Hide whitespace changes
Inline Side-by-side
main.py
+ 98
− 93
View file @ 182c4d61
...@@ -8,98 +8,103 @@ import west.pre as pre ...@@ -8,98 +8,103 @@ import west.pre as pre
import west.utils as utl import west.utils as utl
import west.fluent as fl import west.fluent as fl
# read input config
with open('config.yaml') as f:
config = yaml.load(f, Loader=yaml.FullLoader)
# print(config)
# store input variables
u_ref = config['boundary_condition']['u_ref']
z_ref = config['boundary_condition']['z_ref']
z_0 = config['boundary_condition']['z_0']
rho = config['setup']['fluid']['density']
nu = config['setup']['fluid']['viscosity']
zmax = config['mesh']['domain_height']
nx = config['mesh']['nx']
nz = config['mesh']['nz']
# --- create some directories ---------------------------------------------------------------------------------------- #
dirs = ['post', 'pw', 'stl']
for d in dirs:
(Path.cwd() / 'simulation' / d).mkdir(parents=True, exist_ok=True)
wd = Path.cwd()
for case in config['setup']['wind_directions']:
p = wd / 'simulation' / 'rans' / f'winddir-{case:03.0f}'
p.mkdir(parents=True, exist_ok=True)
# read terrain data
terrain = utl.load_terrain_data(wd / 'data' / config['site']['terrain_data'])
terrain = utl.filter_domain(terrain, config['site']['center'].values(), config['site']['domain_radius'])
dZ = np.nanmax(terrain['z']) - np.nanmin(terrain['z'])
Re = utl.get_Re(config, dZ)
print(f'Re = {int(Re)}')
# get first element height
h_0 = utl.get_first_element_height_by_yplus(u_ref, rho, nu, dZ, 75)
# create stl
if config['options']['do_stl']:
print('Creating STL')
pre.create_stl(terrain)
# create mesh
if config['options']['do_mesh']:
print('Creating Mesh')
wpw.create_structured_mesh(Path.cwd() / 'simulation' / 'stl' / 'terrain.stl', nx=nx, nz=nz, nspokes=24, num_seg=12,
d_h=1, h_min=1500, solve_doms=True) # ,
# path=r'D:\home\aschubig\terrain.cas')
# read terrain data
if config['options']['do_z0']:
print('Creating Roughness Map')
roughness = utl.load_terrain_data(wd / 'data' / config['site']['roughness_data'])
roughness = utl.filter_domain(roughness, config['site']['center'].values(), config['site']['domain_radius'])
if roughness.iloc[:, 2].mean() > 100:
# convert clc to z0
roughness.iloc[:, 2] = roughness.iloc[:, 2].map(lambda x: utl.clc_index_to_z0(x))
# convert z0 to ks
roughness.iloc[:, 2] = roughness.iloc[:, 2].map(lambda x: utl.convert_z0_to_ks(x, cs=1))
# write roughness profile
utl.write_roughness_profile(roughness)
poi = []
poix = []
poiy = []
poiz = []
if config['site']['poi']['poi_1']['x']:
for p in config['site']['poi']:
poi.append(p)
poix.append(str(config['site']['poi'][p]['x']))
poiy.append(str(config['site']['poi'][p]['y']))
poiz.append(str(config['site']['poi'][p]['z'] + config['site']['poi'][p]['h']))
poi_str = utl.convert_to_str(poi)
poix_str = utl.convert_to_str(poix)
poiy_str = utl.convert_to_str(poiy)
poiz_str = utl.convert_to_str(poiz)
# create inlet profiles
if config['options']['do_bc']:
print('Creating Inlet Profiles')
for case in config['setup']['wind_directions']:
utl.create_inlet_profile(config, case, u_ref=u_ref, z_ref=z_ref, z0=z_0, zmin=0, zmax=500)
for case in config['setup']['wind_directions']: if __name__ == '__main__':
print('Creating BCs')
inlet = utl.get_inletzones(case) # read input config
outlet = utl.get_outletzones(case) with open('config.yaml') as f:
inlet_zones = [f'zone-{x}' for x in inlet] config = yaml.load(f, Loader=yaml.FullLoader)
outlet_zones = [f'zone-{x}' for x in outlet] # print(config)
fl.write_jrnl(config, case, poi_str, poix_str, poiy_str, poiz_str, inlet_zones, outlet_zones)
# store input variables
u_ref = config['boundary_condition']['u_ref']
z_ref = config['boundary_condition']['z_ref']
z_0 = config['boundary_condition']['z_0']
rho = config['setup']['fluid']['density']
nu = config['setup']['fluid']['viscosity']
zmax = config['mesh']['domain_height']
nx = config['mesh']['nx']
nz = config['mesh']['nz']
# --- create some directories ---------------------------------------------------------------------------------------- #
dirs = ['post', 'pw', 'stl']
for d in dirs:
(Path.cwd() / 'simulation' / d).mkdir(parents=True, exist_ok=True)
wd = Path.cwd()
for case in config['setup']['wind_directions']:
p = wd / 'simulation' / 'rans' / f'winddir-{case:03.0f}'
p.mkdir(parents=True, exist_ok=True)
# read terrain data
terrain = utl.load_terrain_data(wd / 'data' / config['site']['terrain_data'])
terrain = utl.filter_domain(terrain, config['site']['center'].values(), config['site']['domain_radius'])
dZ = np.nanmax(terrain['z']) - np.nanmin(terrain['z'])
Re = utl.get_Re(config, dZ)
print(f'Re = {int(Re)}')
# get first element height
h_0 = utl.get_first_element_height_by_yplus(u_ref, rho, nu, dZ, 75)
# create stl
if config['options']['do_stl']:
print('Creating STL')
pre.create_stl(terrain)
# create mesh
if config['options']['do_mesh']:
print('Creating Mesh')
wpw.create_structured_mesh(Path.cwd() / 'simulation' / 'stl' / 'terrain.stl', nx=nx, nz=nz, nspokes=24, num_seg=12,
d_h=1, h_min=1500, solve_doms=True, export_cae=True)
# path=r'D:\home\aschubig\terrain.cas')
# read terrain data
if config['options']['do_z0']:
print('Creating Roughness Map')
roughness = utl.load_roughness_data(wd / 'data' / config['site']['roughness_data'])
roughness = utl.filter_domain(roughness, config['site']['center'].values(), config['site']['domain_radius'])
if roughness.iloc[:, 2].mean() > 100:
# convert clc to z0
roughness.iloc[:, 2] = roughness.iloc[:, 2].map(lambda x: utl.clc_index_to_z0(x))
elif 100 > roughness.iloc[:, 2].mean() > 1.2:
roughness.iloc[:, 2] = roughness.iloc[:, 2].map(lambda x: utl.clc_to_z0(x))
# convert z0 to ks
roughness.iloc[:, 2] = roughness.iloc[:, 2].map(lambda x: utl.convert_z0_to_ks(x, cs=1))
# write roughness profile
utl.write_roughness_profile(roughness)
poi = []
poix = []
poiy = []
poiz = []
if config['site']['poi']['poi_1']['x']:
for p in config['site']['poi']:
poi.append(p)
poix.append(str(config['site']['poi'][p]['x']))
poiy.append(str(config['site']['poi'][p]['y']))
poiz.append(str(config['site']['poi'][p]['z'] + config['site']['poi'][p]['h']))
poi_str = utl.convert_to_str(poi)
poix_str = utl.convert_to_str(poix)
poiy_str = utl.convert_to_str(poiy)
poiz_str = utl.convert_to_str(poiz)
# create inlet profiles
if config['options']['do_bc']:
print('Creating Inlet Profiles')
for case in config['setup']['wind_directions']:
utl.create_inlet_profile(config, case, u_ref=u_ref, z_ref=z_ref, z0=z_0, zmin=0, zmax=500)
print('Preprocessing done') for case in config['setup']['wind_directions']:
print('Creating BCs')
inlet = utl.get_inletzones(case)
outlet = utl.get_outletzones(case)
inlet_zones = [f'zone-{x}' for x in inlet]
outlet_zones = [f'zone-{x}' for x in outlet]
fl.write_jrnl(config, case, poi_str, poix_str, poiy_str, poiz_str, inlet_zones, outlet_zones)
print('Preprocessing done')
west/pointwisemesh.py
+ 17
− 6
View file @ 182c4d61
from pathlib import Path from pathlib import Path
from pointwise.glyphapi import * from pointwise.glyphapi import *
import yaml
# Connect to the Pointwise server listening on localhost at the default port # Connect to the Pointwise server listening on localhost at the default port
# with no authentication token... # with no authentication token...
...@@ -68,6 +68,10 @@ def create_structured_mesh(stl_file, nx=256, nz=32, nspokes=24, num_seg=12, d_h= ...@@ -68,6 +68,10 @@ def create_structured_mesh(stl_file, nx=256, nz=32, nspokes=24, num_seg=12, d_h=
# solve_doms = True # solve_doms = True
# solve_blks = True # solve_blks = True
# -------------------------------------------------------------------------------------------------------------------- # # -------------------------------------------------------------------------------------------------------------------- #
# read input config
with open('config.yaml') as f:
config = yaml.load(f, Loader=yaml.FullLoader)
print('Starting Pointwise') print('Starting Pointwise')
# Run in batch mode # Run in batch mode
# Pointwise directory has to be in PATH e.g. C:\Program Files\Pointwise\PointwiseV18.4R4\win64\bin # Pointwise directory has to be in PATH e.g. C:\Program Files\Pointwise\PointwiseV18.4R4\win64\bin
...@@ -116,6 +120,8 @@ def create_structured_mesh(stl_file, nx=256, nz=32, nspokes=24, num_seg=12, d_h= ...@@ -116,6 +120,8 @@ def create_structured_mesh(stl_file, nx=256, nz=32, nspokes=24, num_seg=12, d_h=
z0 = extents[0][2] z0 = extents[0][2]
z1 = extents[1][2] z1 = extents[1][2]
# print(x0, x1, y0, y1, z0, z1)
center = [(x1 + x0) / 2, (y1 + y0) / 2, z0] center = [(x1 + x0) / 2, (y1 + y0) / 2, z0]
# set new center of domain if coordinates are given in c_domain # set new center of domain if coordinates are given in c_domain
...@@ -131,17 +137,22 @@ def create_structured_mesh(stl_file, nx=256, nz=32, nspokes=24, num_seg=12, d_h= ...@@ -131,17 +137,22 @@ def create_structured_mesh(stl_file, nx=256, nz=32, nspokes=24, num_seg=12, d_h=
d = x1 - x0 d = x1 - x0
dx = d / nx dx = d / nx
print(f'dx: {dx}')
# set domain radius to x*diameter # set domain radius to x*diameter
r_outter = 3 * d r_outter = 3 * d
# calculate essential coordinates # calculate essential coordinates
p_outter = get_outer_coordinates(num_seg, r_outter, center) p_outter = get_outer_coordinates(num_seg, r_outter, center)
corners_start = [(x0, y0, z0), (x0, y1, z0), (x1, y0, z0), (x0, y0, z0)] # corners_start = [(x0, y0, z0), (x0, y1, z0), (x1, y0, z0), (x0, y0, z0)]
corners_end = [(x0, y1, z0), (x1, y1, z0), (x1, y1, z0), (x1, y0, z0)] # corners_end = [(x0, y1, z0), (x1, y1, z0), (x1, y1, z0), (x1, y0, z0)]
offset = 2*dx # move boundary inside
corners_start = [(x0+offset, y0+offset, z0), (x0+offset, y1-offset, z0), (x1-offset, y0+offset, z0), (x0+offset, y0+offset, z0)] # temporary workaround to fix holes in the stl boundary
corners_end = [(x0+offset, y1-offset, z0), (x1-offset, y1-offset, z0), (x1-offset, y1-offset, z0), (x1-offset, y0+offset, z0)] # temporary workaround to fix holes in the stl boundary
z_max = z1 + h_min z_max = z1 + h_min
corners_top_start = [(x0, y0, z_max), (x0, y1, z_max), (x1, y0, z_max), (x0, y0, z_max)] # corners_top_start = [(x0, y0, z_max), (x0, y1, z_max), (x1, y0, z_max), (x0, y0, z_max)]
corners_top_end = [(x0, y1, z_max), (x1, y1, z_max), (x1, y1, z_max), (x1, y0, z_max)] # corners_top_end = [(x0, y1, z_max), (x1, y1, z_max), (x1, y1, z_max), (x1, y0, z_max)]
corners_top_start = [(x0+offset, y0+offset, z_max), (x0+offset, y1-offset, z_max), (x1-offset, y0+offset, z_max), (x0+offset, y0+offset, z_max)] # temporary workaround to fix holes in the stl boundary
corners_top_end = [(x0+offset, y1-offset, z_max), (x1-offset, y1-offset, z_max), (x1-offset, y1-offset, z_max), (x1-offset, y0+offset, z_max)] # temporary workaround to fix holes in the stl boundary
# -------------------------------------------------------------------------------------------------------------------- # # -------------------------------------------------------------------------------------------------------------------- #
# inner block ########################################################################################################## # inner block ##########################################################################################################
# create connectors on ground stl # create connectors on ground stl
......
west/pre.py
+ 3
− 2
View file @ 182c4d61
...@@ -5,12 +5,12 @@ import open3d as o3d ...@@ -5,12 +5,12 @@ import open3d as o3d
import pandas as pd import pandas as pd
from pyproj import Transformer from pyproj import Transformer
import west.utils as wu import west.utils as utl
def create_stl(df, method='Poisson', downsample=False, cleanup=False, export_stl=True, def create_stl(df, method='Poisson', downsample=False, cleanup=False, export_stl=True,
path=str(Path.cwd() / 'simulation' / 'stl' / 'terrain.stl')): path=str(Path.cwd() / 'simulation' / 'stl' / 'terrain.stl')):
extents = wu.get_extents(df) extents = utl.get_extents(df)
# create open3d object # create open3d object
point_cloud = np.array(df) point_cloud = np.array(df)
...@@ -30,6 +30,7 @@ def create_stl(df, method='Poisson', downsample=False, cleanup=False, export_stl ...@@ -30,6 +30,7 @@ def create_stl(df, method='Poisson', downsample=False, cleanup=False, export_stl
# o3d.visualization.draw_geometries([pcd]) # o3d.visualization.draw_geometries([pcd])
bbox = pcd.get_axis_aligned_bounding_box() bbox = pcd.get_axis_aligned_bounding_box()
print(bbox)
# bbox = o3d.geometry.AxisAlignedBoundingBox(min_bound=(extents[0]+1, extents[2]+1, extents[4]+1), # bbox = o3d.geometry.AxisAlignedBoundingBox(min_bound=(extents[0]+1, extents[2]+1, extents[4]+1),
# max_bound=(extents[1]-1, extents[3]-1, extents[5]-1)) # max_bound=(extents[1]-1, extents[3]-1, extents[5]-1))
......
west/utils.py
+ 19
− 4
View file @ 182c4d61
...@@ -75,6 +75,10 @@ def load_terrain_data(path): ...@@ -75,6 +75,10 @@ def load_terrain_data(path):
data = pd.read_csv(path, sep=',', header=None) data = pd.read_csv(path, sep=',', header=None)
data.columns = ['x', 'y', 'z'] data.columns = ['x', 'y', 'z']
print('Data imported') print('Data imported')
print(f'min\n{data.min()}')
print(f'mean\n{data.mean()}')
print(f'max\n{data.max()}')
print(f'delta\n{abs(data.max() - data.min())}')
return data return data
...@@ -82,6 +86,10 @@ def load_roughness_data(path): ...@@ -82,6 +86,10 @@ def load_roughness_data(path):
data = pd.read_csv(path, sep=',', header=None) data = pd.read_csv(path, sep=',', header=None)
data.columns = ['x', 'y', 'z0'] data.columns = ['x', 'y', 'z0']
print('Data imported') print('Data imported')
print(f'min\n{data.min()}')
print(f'mean\n{data.mean()}')
print(f'max\n{data.max()}')
print(f'delta\n{abs(data.max() - data.min())}')
return data return data
...@@ -91,6 +99,10 @@ def filter_domain(df, center, radius): ...@@ -91,6 +99,10 @@ def filter_domain(df, center, radius):
df.iloc[:, 1].between(list(center)[1] - radius, list(center)[1] + radius) & df.iloc[:, 1].between(list(center)[1] - radius, list(center)[1] + radius) &
df.iloc[:, 2].between(0, 10000)] df.iloc[:, 2].between(0, 10000)]
print('Data filtered') print('Data filtered')
print(f'min\n{data.min()}')
print(f'mean\n{data.mean()}')
print(f'max\n{data.max()}')
print(f'delta\n{abs(data.max() - data.min())}')
return data return data
...@@ -201,10 +213,10 @@ def convert_ks_to_z0(ks, cs=0.5, solver='fluent'): ...@@ -201,10 +213,10 @@ def convert_ks_to_z0(ks, cs=0.5, solver='fluent'):
return ks / 30 return ks / 30
def write_roughness_profile(site_data, data): # def write_roughness_profile(site_data, data):
with open(Path.cwd() / 'rans' / 'roughness.csv', 'w', newline='\n') as f: # with open(Path.cwd() / 'rans' / 'roughness.csv', 'w', newline='\n') as f:
f.write('[Name] \nroughness \n\n[Data] \nx,y,z0 \n') # f.write('[Name] \nroughness \n\n[Data] \nx,y,z0 \n')
data.to_csv(f, index=False, header=False) # data.to_csv(f, index=False, header=False)
def get_u_abl(u_ref, z_ref, z0, c=0.09, kappa=0.42): def get_u_abl(u_ref, z_ref, z0, c=0.09, kappa=0.42):
...@@ -462,5 +474,8 @@ def linear_regression(x, y, title=[]): ...@@ -462,5 +474,8 @@ def linear_regression(x, y, title=[]):
with open(Path.cwd() / 'data/plots' / f'linear_regression.txt', 'w', newline='\n') as f: with open(Path.cwd() / 'data/plots' / f'linear_regression.txt', 'w', newline='\n') as f:
f.write(str(results.summary())) f.write(str(results.summary()))
with open(Path.cwd() / 'data/plots' / f'linear_regression.txt', 'w', newline='\n') as f:
f.write(str(results.summary()))
print(results.summary()) print(results.summary())
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment