-
Notifications
You must be signed in to change notification settings - Fork 158
/
call_capytaine.py
347 lines (294 loc) · 12.4 KB
/
call_capytaine.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
'''
Function to run Capytaine for hydrodynamics and Meshmagick for hydrostatics.
'''
import os
import numpy as np
from multiprocessing import Process
capy_v2 = False
import capytaine as cpt
# Get the version of capytaine
capytaine_version = cpt.__version__
# If the version is less than 2.0, import meshmagick:
if int(capytaine_version.split('.')[0]) >= 2:
capy_v2 = True
if not capy_v2:
import meshmagick.mesh as mmm
try:
# Latest version on github removed the previous
# meshmagick.hydrostatics.Hydrostatics() method. Use old module w/ new version
import meshmagick.hydrostatics_old as mmhs
except ModuleNotFoundError:
# Older versions of meshmagick should have meshmagick.hydrostatics.Hydrostatics() method
import meshmagick.hydrostatics as mmhs
import xarray as xr
import logging as LOG
from glob import glob
import shutil
import platform
import sys
# Set the affinity back on all the cores solving the single core issue.
# 0xf is essentially a hexadecimal bitmask, corresponding to 4 cores
# 0xffffffffff is used here for a maximum of 40 cores.
# This is only used for Linux machines
if platform == "linux" or platform == "linux2":
os.system("taskset -p 0xffffffffff %d" % os.getpid())
def __init__(self):
LOG.info("Capytaine imported.")
def hydrostatics(myBodies, savepath=''):
'''
use Capytaine/Meshmagick functions to calculate and output the hydrostatic
stiffness, inertia, center of gravity, center of buoyancy and displaced
volume of a capytaine bodies. Output is saved to Hydrostatics.dat and KH.dat
files in the same manner as Nemoh
Example of output format:
Hydrostatics.dat:
XF = 0.000 - XG = 0.000
YF = 0.000 - YG = 0.000
ZF = -2.500 - ZG = -2.500
Displacement = 0.4999997E+03
Waterplane area = 0.1000002E+03
KH.dat:
0.0000000E+00 0.0000000E+00 0.0000000E+00 0.0000000E+00 0.0000000E+00 0.0000000E+00
0.0000000E+00 0.0000000E+00 0.0000000E+00 0.0000000E+00 0.0000000E+00 0.0000000E+00
0.0000000E+00 0.0000000E+00 0.9810053E+06 -0.1464844E-01 -0.5859375E-02 0.0000000E+00
0.0000000E+00 0.0000000E+00 -0.1464844E-01 0.8160803E+07 0.0000000E+00 0.0000000E+00
0.0000000E+00 0.0000000E+00 -0.5859375E-02 0.0000000E+00 0.8160810E+07 0.0000000E+00
0.0000000E+00 0.0000000E+00 0.0000000E+00 0.0000000E+00 0.0000000E+00 0.0000000E+00
Parameters
----------
myBodies : List
A list of capytaine floating bodies.
Returns
-------
None
'''
nbod = len(myBodies)
for i, body in enumerate(myBodies):
cg = body.center_of_mass
# Set file index
fileind = '' if nbod == 1 else '_' + str(i)
# use meshmagick to compute hydrostatic stiffness matrix
# NOTE: meshmagick currently has issue if a body is copmletely submerged (OSWEC base)
# use try-except statement to catch this error use alternate function for cb/vo
try:
if capy_v2:
# Capytaine version is >= 2.0
body_hs = body.compute_hydrostatics(rho=1023.0)
vo = body_hs['disp_volume']
cb = body_hs['center_of_buoyancy']
khs = body_hs['hydrostatic_stiffness']
else:
# Capytaine version is < 2.0; use meshmagick
body_mesh = mmm.Mesh(body.mesh.vertices, body.mesh.faces, name= body.mesh.name)
body_hs = mmhs.Hydrostatics(working_mesh=body_mesh, cog=body.center_of_mass, rho_water=1023.0, grav=9.81)
vo = body_hs.displacement_volume
cb = body_hs.buoyancy_center
khs = body_hs.hydrostatic_stiffness_matrix
except:
# Exception if body is fully submerged
vo = body.volume if capy_v2 else body_mesh.volume
cb = body.center_of_buoyancy if capy_v2 else cg
khs = np.zeros((3,3))
# Write hydrostatic stiffness to KH.dat file
khs_full = np.zeros((6,6))
if capy_v2:
khs_full[2:5, 2:5] += khs[2:5, 2:5]
else:
khs_full[2:5, 2:5] += khs
tmp = savepath + 'KH' + fileind +'.dat'
np.savetxt(tmp, khs_full)
# Write the other hydrostatics data to Hydrostatics.dat file
tmp = savepath + 'Hydrostatics' + fileind + '.dat'
# f = open(rf'{tmp}','w')
f = open(tmp,'w')
for j in [0,1,2]:
line = f'XF = {cb[j]:7.3f} - XG = {cg[j]:7.3f} \n'
f.write(line)
line = f'Displacement = {vo:E}'
f.write(line)
f.close()
def call_capy(meshFName, wCapy, CoG=([0,0,0],), headings=[0.0],ncFName=None,
wDes=None, body_name=('',), depth=np.infty, density=1025.0,
additional_dofs_dir=None,num_threads=1):
'''
Setup the problem and call the capytine solver.
Setup parallel computing for different frequencies and combine the data
after the parallel simulation is completed.
May be called with multiple bodies (automatically implements B2B).
In this case, the meshFName, CoG, body_name should be a tuple of the
values for each body.
Also has the ability to add generalized body modes by inputting the path to
a 'gbm_dofs.py' file (see RM3 example).
Parameters
----------
meshFName : tuple of strings
Tuple containing a string for the path to each body's hydrodynamic mesh.
mesh must be cropped at waterline (OXY plane) and have no lid
wCapy: array
array of frequency points to be computed by Capytaine
CoG: tuple of lists
tuple contains a 3x1 list of each body's CoG
headings: list
list of wave headings to compute [rad]
saveNc: Bool
save results to .nc file
ncFName: str
name of .nc file
wDes: array
array of desired frequency points
(for interpolation of wCapy-based Capytaine data)
body_name: tuple of strings
Tuple containing strings. Strings are the names of each body.
Prevent the body name from being a long file path
depth: float
Water depth. Should be positive downwards. Use decimal value to prevent
Capytaine outputting int32 types. Default is -np.infty
density: float
Water density. Use decimal value to prevent Capytaine outputting int32
types. Default 1025.0
additional_dofs: string
path to a gbm_dofs.py file that returns GBM dofs to this function
Returns
-------
capyData: xarray Dataset
Hydrodynamic coefficients as computed
problems: list
capytaine Problems that were solved
'''
# check that old output is not being overwritten (runs take awhile)
if os.path.isfile(ncFName):
print(f'Output ({ncFName}) file already exists and will be overwritten. '
'Do you wish to proceed? (y/n)')
try:
ans = input()
except EOFError:
# Catch error that occurs when this script is run in a
# non-interactive way ('python CASE.py' in run_cases.py, etc) and
# default to overwriting the output file
ans = 'y'
pass
if ans.lower() != 'y':
print('\nEnding simulation. file not overwritten')
sys.exit(0)
bodies = []
for i in np.arange(0, len(meshFName)):
bodies.append(cpt.FloatingBody.from_file(meshFName[i]))
bodies[i].center_of_mass = CoG[i]
bodies[i].keep_immersed_part()
if body_name[i] != '':
bodies[i].name = body_name[i]
bodies[i].add_all_rigid_body_dofs()
# calculate hydrostatics and output to KH.dat and Hydrostatics.dat files
path,tmp = os.path.split(ncFName)
path += os.path.sep
hydrostatics(bodies,path)
# add gbm dofs
# 1. pass gbm_dofs.py path to call_capy
# 2. Here: pass mesh to a local gbm_dofs.py script
# 3. in the gbm_dofs.py script, create dofs based on body mesh that is passed in
# 4. pass dof dict back to call_capy and continue adding to body
if additional_dofs_dir is not None:
old_dir = os.getcwd()
os.chdir(additional_dofs_dir)
import gbm_dofs
additional_dofs = gbm_dofs.new_dofs(bodies)
for i in np.arange(0, len(meshFName)):
if body_name[i] in additional_dofs:
for k,v in additional_dofs[body_name[i]].items():
bodies[i].dofs[k] = v
os.chdir(old_dir)
# combine all bodies to account for B2B interactions
combo = bodies[0]
for i in np.arange(1,len(bodies),1):
combo += bodies[i]
# call Capytaine solver
print(f'\n-------------------------------\n'
f'Calling Capytaine BEM solver...\n'
f'-------------------------------\n'
f'mesh = {meshFName}\n'
f'w range = {wCapy[0]:.3f} - {wCapy[-1]:.3f} rad/s\n'
f'dw = {(wCapy[1]-wCapy[0]):.3f} rad/s\n'
f'no of headings = {len(headings)}\n'
f'no of radiation & diffraction problems = {len(wCapy)*(len(headings) + len(combo.dofs))}\n'
f'-------------------------------\n')
wCapy_threads = np.array_split(np.array(wCapy),num_threads)
if num_threads != 1:
try:
shutil.rmtree('capyParallelFolder')
except OSError as e:
pass
os.mkdir('capyParallelFolder')
# An array for the processes.
processing_jobs = []
for i in range(num_threads):
if num_threads == 1:
ncFName_each_thread = ncFName
else:
os.chdir("./capyParallelFolder")
ncFName_each_thread = os.getcwd() + os.path.sep + "capyParallel_{}.nc".format(i+1)
os.chdir("../")
p = Process(target=capy_solver, args= (wCapy_threads[i], CoG, headings,ncFName_each_thread,wDes, body_name, depth, density,
combo,additional_dofs_dir))
processing_jobs.append(p)
p.start()
# Wait for all processes to finish.
for proc in processing_jobs:
proc.join()
if num_threads == 1:
capyData = read_netcdfs(ncFName, dim='omega')
else:
os.chdir("./capyParallelFolder")
ncFName_thread = os.getcwd() + os.path.sep + 'capyParallel_*.nc'
os.chdir("../")
capyData = read_netcdfs(ncFName_thread, dim='omega')
print('\nCombine Capytaine data and saved to \n' + ncFName +'\n\n')
capyData.to_netcdf(ncFName)
# Remove saved Capytaine data from each thread.
try:
shutil.rmtree('capyParallelFolder')
except OSError as e:
pass
# Create a dataset of parameters.
# 'fill_dataset()' automatically creates problems and solves them.
problems = xr.Dataset(coords={
'omega': wCapy,
'wave_direction': headings,
'radiating_dof': list(combo.dofs),
'water_depth': [depth],
'rho': [density],
})
print('\nCapytaine call complete. \n\n')
return capyData, problems
def read_netcdfs(files, dim):
# glob expands paths with * to a list of files, like the unix shell
paths = sorted(glob(files))
datasets = [xr.open_dataset(p) for p in paths]
combined = xr.concat(datasets, dim)
return combined
def capy_solver(wCapy, CoG, headings,ncFName,wDes, body_name, depth, density,
combo,additional_dofs_dir):
'''
call Capytaine for a given mesh, frequency range and wave headings
This function is modified from David Ogden's work
(see https://github.com/mattEhall/FrequencyDomain/blob/b89dd4f4a732fbe4afde56efe2b52c3e32e22d53/FrequencyDomain.py#L842 for the original function).
save the results to a file in Network Common Data Form.
Returns
-------
None
'''
problems = xr.Dataset(coords={
'omega': wCapy,
'wave_direction': headings,
'radiating_dof': list(combo.dofs),
'water_depth': [depth],
'rho': [density],
})
solver = cpt.BEMSolver()
capyData = solver.fill_dataset(problems, [combo], hydrostatics=False)
# # add kochin diffraction results
# kochin = cpt.io.xarray.kochin_data_array(results, headings)
# capyData.update(kochin)
# save to .nc file
cpt.io.xarray.separate_complex_values(capyData).to_netcdf(ncFName, encoding={'radiating_dof': {'dtype': 'U'}, 'influenced_dof': {'dtype': 'U'}})
print('\nCapytaine call complete. Data saved to \n' + ncFName +'\n\n')
return