-
Notifications
You must be signed in to change notification settings - Fork 0
/
Artery.py
457 lines (372 loc) · 18.6 KB
/
Artery.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
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
from os import path, makedirs
import json
import pickle
from pprint import pprint
import numpy as np
#from Probe import Probes
from Womersley import make_womersley_bcs, compute_boundary_geometry_acrn
from oasis.problems.NSfracStep import *
"""
Problem file for running CFD simulation in arterial models consisting of one inlet, and two or more outlets.
A Womersley velocity profile is applied at the inlet, and a flow split pressure condition is applied at the outlets,
following [1]. Flow rate for the inlet condition, and flow split values for the outlets are computed from the
pre-processing script automatedPreProcessing.py. The simulation is run for two cycles (adjustable), but only the
results/solutions from the second cycle are stored to avoid non-physiological effects from the first cycle.
One cardiac cycle is set to 0.951 s from [2], and scaled by a factor of 1000, hence all parameters are in [mm] or [ms].
[1] Gin, Ron, Anthony G. Straatman, and David A. Steinman. "A dual-pressure boundary condition for use in simulations
of bifurcating conduits." J. Biomech. Eng. 124.5 (2002): 617-619.
[2] Hoi, Yiemeng, et al. "Characterization of volumetric flow rate waveforms at the carotid bifurcations of older
adults." Physiological measurement 31.3 (2010): 291.
"""
# FEniCS specific command to control the desired level of logging, here set to critical errors
set_log_level(50)
def problem_parameters(commandline_kwargs, NS_parameters, NS_expressions, **NS_namespace):
if "restart_folder" in commandline_kwargs.keys():
restart_folder = commandline_kwargs["restart_folder"]
f = open(path.join(restart_folder, 'params.dat'), 'rb')
NS_parameters.update(pickle.load(f))
NS_parameters['restart_folder'] = restart_folder
else:
# Parameters are in mm and ms
t_values, Q_ = np.loadtxt(path.join(path.dirname(path.abspath(__file__)), "inflow.flow")).T
cardiac_cycle = t_values[-1]*1000
number_of_cycles = 2
NS_parameters.update(
# Fluid parameters
nu=3.3018e-3, # Kinematic viscosity: 0.0035 Pa-s / 1060 kg/m^3 = 3.3018E-6 m^2/s = 3.3018-3 mm^2/ms
# Geometry parameters
id_in=[], # Inlet boundary ID
id_out=[], # Outlet boundary IDs
area_ratio=[], # Area ratio for the flow outlets
area_inlet=[], # Area of inlet in [mm^2]
# Simulation parameters
cardiac_cycle=cardiac_cycle, # Duration of cardiac cycle [ms]
T=cardiac_cycle * number_of_cycles, # Simulation end time [ms]
dt=cardiac_cycle/1000, # Time step size [ms]
dump_probe_frequency=100, # Dump frequency for sampling velocity & pressure at probes along the centerline
save_solution_frequency=5, # Save frequency for velocity and pressure field
save_solution_after_cycle=1, # Store solution after 1 cardiac cycle
# Oasis specific parameters
checkpoint=500, # Checkpoint frequency
print_intermediate_info=100, # Frequency for printing solver statistics
folder="results_artery", # Preferred results folder name
mesh_path=commandline_kwargs["mesh_path"], # Path to the mesh
# Solver parameters
velocity_degree=1, # Polynomial order of finite element for velocity. Normally linear (1) or quadratic (2)
pressure_degree=1, # Polynomial order of finite element for pressure. Normally linear (1)
use_krylov_solvers=True,
krylov_solvers=dict(monitor_convergence=False)
)
mesh_file = NS_parameters["mesh_path"].split("/")[-1]
case_name = mesh_file.split(".")[0]
NS_parameters["folder"] = path.join(NS_parameters["folder"], case_name)
if MPI.rank(MPI.comm_world) == 0:
print("=== Starting simulation for case: {} ===".format(case_name))
print("Running with the following parameters:")
pprint(NS_parameters)
def mesh(mesh_path, **NS_namespace):
# Read mesh and print mesh information
mesh = Mesh(mesh_path)
print_mesh_information(mesh)
return mesh
def create_bcs(t, NS_expressions, V, Q, area_ratio, area_inlet, mesh, mesh_path, nu, id_in, id_out, pressure_degree,
**NS_namespace):
# Mesh function
boundary = MeshFunction("size_t", mesh, mesh.geometry().dim() - 1, mesh.domains())
# Read case parameters
info = mesh_path.split(".xml")[0] + "_info.json"
with open(info) as f:
info = json.load(f)
id_in[:] = info['inlet_id']
id_out[:] = info['outlet_ids']
id_wall = int(min(id_in + id_out)) - 1
Q_mean = info['mean_flow_rate']
area_ratio[:] = info['area_ratio']
area_inlet.append(info['inflow_area'])
# Load normalized time and flow rate values
t_values, Q_ = np.loadtxt(path.join(path.dirname(path.abspath(__file__)), "inflow.flow")).T
Q_values = Q_mean * Q_ # Specific flow rate = Normalized flow wave form * Prescribed flow rate
t_values *= 1000 # Scale time in normalised flow wave form to [ms]
tmp_area, tmp_center, tmp_radius, tmp_normal = compute_boundary_geometry_acrn(mesh, id_in[0], boundary)
# Create Womersley boundary condition at inlet
inlet = make_womersley_bcs(t_values, Q_values, mesh, nu, tmp_area, tmp_center, tmp_radius, tmp_normal,
V.ufl_element())
NS_expressions["inflow"] = inlet
# Initialize inlet expressions with initial time
for uc in inlet:
uc.set_t(t)
# Create pressure boundary condition
area_out = []
for i, ind in enumerate(id_out):
dsi = ds(ind, domain=mesh, subdomain_data=boundary)
area_out.append(assemble(Constant(1.0, name="one") * dsi))
bc_p = []
if MPI.rank(MPI.comm_world) == 0:
print("=== Initial pressure and area fraction ===")
for i, ID in enumerate(id_out):
p_initial = area_out[i] / sum(area_out)
outflow = Expression("p", p=p_initial, degree=pressure_degree)
bc = DirichletBC(Q, outflow, boundary, ID)
bc_p.append(bc)
NS_expressions[ID] = outflow
if MPI.rank(MPI.comm_world) == 0:
print(("Boundary ID={:d}, pressure: {:0.6f}, area fraction: {:0.4f}".format(ID, p_initial, area_ratio[i])))
# No slip condition at wall
wall = Constant(0.0)
# Create Boundary conditions for the velocity
bc_wall = DirichletBC(V, wall, boundary, id_wall)
bc_inlet = [DirichletBC(V, inlet[i], boundary, id_in[0]) for i in range(3)]
# Return boundary conditions in dictionary
return dict(u0=[bc_inlet[0], bc_wall],
u1=[bc_inlet[1], bc_wall],
u2=[bc_inlet[2], bc_wall],
p=bc_p)
# Oasis hook called before simulation start
def pre_solve_hook(mesh, V, Q, newfolder, mesh_path, restart_folder, velocity_degree, cardiac_cycle,
save_solution_after_cycle, dt, **NS_namespace):
# Mesh function
boundary = MeshFunction("size_t", mesh, mesh.geometry().dim() - 1, mesh.domains())
# Create point for evaluation
n = FacetNormal(mesh)
eval_dict = {}
""" rel_path = mesh_path.split(".xml")[0] + "_probe_point"
probe_points = np.load(rel_path, encoding='latin1', fix_imports=True, allow_pickle=True)
# Store points file in checkpoint
if MPI.rank(MPI.comm_world) == 0:
probe_points.dump(path.join(newfolder, "Checkpoint", "points"))
eval_dict["centerline_u_x_probes"] = Probes(probe_points.flatten(), V)
eval_dict["centerline_u_y_probes"] = Probes(probe_points.flatten(), V)
eval_dict["centerline_u_z_probes"] = Probes(probe_points.flatten(), V)
eval_dict["centerline_p_probes"] = Probes(probe_points.flatten(), Q)"""
if restart_folder is None:
# Get files to store results
files = get_file_paths(newfolder)
NS_parameters.update(dict(files=files))
else:
files = NS_namespace["files"]
# Save mesh as HDF5 file for post processing
with HDF5File(MPI.comm_world, files["mesh"], "w") as mesh_file:
mesh_file.write(mesh, "mesh")
# Create vector function for storing velocity
Vv = VectorFunctionSpace(mesh, "CG", velocity_degree)
U = Function(Vv)
u_mean = Function(Vv)
u_mean0 = Function(V)
u_mean1 = Function(V)
u_mean2 = Function(V)
# Tstep when solutions for post processing should start being saved
save_solution_at_tstep = int(cardiac_cycle * save_solution_after_cycle / dt)
return dict(eval_dict=eval_dict, boundary=boundary, n=n, U=U, u_mean=u_mean, u_mean0=u_mean0, u_mean1=u_mean1,
u_mean2=u_mean2, save_solution_at_tstep=save_solution_at_tstep)
#Backflow stabilization
def velocity_tentative_hook(mesh, mesh_path, id_in, id_out, ui, A, u_ab, u, v, q_1, b, x_1, boundary, **NS_namespace):
backflow_beta = 0.2
for i, ind in enumerate(id_out):
ds = Measure("ds", domain=mesh, subdomain_data=boundary)
n = FacetNormal(mesh)
K = assemble(u_dot_n(u_ab, n) * dot(u, v) * ds(ind))
A.axpy(-backflow_beta * 0.5, K, True)
b[ui].axpy(backflow_beta * 0.5, K * x_1[ui])
# Oasis hook called after each time step
def temporal_hook(u_, p_, mesh, tstep, dump_probe_frequency, eval_dict, newfolder, id_in, id_out, boundary, n,
save_solution_frequency, NS_parameters, NS_expressions, area_ratio, t, save_solution_at_tstep,
U, area_inlet, nu, u_mean0, u_mean1, u_mean2, mesh_path, ui, A, u_ab, u, v, q_1, b, x_1, **NS_namespace):
# Update boundary condition to current time
for uc in NS_expressions["inflow"]:
uc.set_t(t)
# Compute flux and update pressure condition
if tstep > 2:
Q_ideals, Q_in, Q_outs = update_pressure_condition(NS_expressions, area_ratio, boundary, id_in, id_out, mesh, n,
tstep, u_, mesh_path, ui, A, u_ab, u, v, q_1, b, x_1, **NS_namespace)
# Compute flow rates and updated pressure at outlets, and mean velocity and Reynolds number at inlet
if MPI.rank(MPI.comm_world) == 0 and tstep % 10 == 0:
U_mean = Q_in / area_inlet[0]
diam_inlet = np.sqrt(4 * area_inlet[0] / np.pi)
Re = U_mean * diam_inlet / nu
print("=" * 10, "Time step " + str(tstep), "=" * 10)
print("Sum of Q_out = {:0.4f}, Q_in = {:0.4f}, mean velocity (inlet): {:0.4f}, Reynolds number (inlet): {:0.4f}"
.format(sum(Q_outs), Q_in, U_mean, Re))
for i, out_id in enumerate(id_out):
print(("For outlet with boundary ID={:d}: target flow rate: {:0.4f} mL/s, " +
"computed flow rate: {:0.4f} mL/s, pressure updated to: {:0.4f}")
.format(out_id, Q_ideals[i], Q_outs[i], NS_expressions[out_id].p))
print()
"""# Sample velocity and pressure in points/probes
eval_dict["centerline_u_x_probes"](u_[0])
eval_dict["centerline_u_y_probes"](u_[1])
eval_dict["centerline_u_z_probes"](u_[2])
eval_dict["centerline_p_probes"](p_)
# Store sampled velocity and pressure
if tstep % dump_probe_frequency == 0:
# Save variables along the centerline for CFD simulation
# diagnostics and light-weight post processing
filepath = path.join(newfolder, "Probes")
if MPI.rank(MPI.comm_world) == 0:
if not path.exists(filepath):
makedirs(filepath)
arr_u_x = eval_dict["centerline_u_x_probes"].array()
arr_u_y = eval_dict["centerline_u_y_probes"].array()
arr_u_z = eval_dict["centerline_u_z_probes"].array()
arr_p = eval_dict["centerline_p_probes"].array()
# Dump stats
if MPI.rank(MPI.comm_world) == 0:
arr_u_x.dump(path.join(filepath, "u_x_%s.probes" % str(tstep)))
arr_u_y.dump(path.join(filepath, "u_y_%s.probes" % str(tstep)))
arr_u_z.dump(path.join(filepath, "u_z_%s.probes" % str(tstep)))
arr_p.dump(path.join(filepath, "p_%s.probes" % str(tstep)))
# Clear stats
MPI.barrier(MPI.comm_world)
eval_dict["centerline_u_x_probes"].clear()
eval_dict["centerline_u_y_probes"].clear()
eval_dict["centerline_u_z_probes"].clear()
eval_dict["centerline_p_probes"].clear()"""
# Save velocity and pressure for post processing
if tstep % save_solution_frequency == 0 and tstep >= save_solution_at_tstep:
# Assign velocity components to vector solution
assign(U.sub(0), u_[0])
assign(U.sub(1), u_[1])
assign(U.sub(2), u_[2])
# Get save paths
files = NS_parameters['files']
p_path = files['p']
u_path = files['u']
file_mode = "w" if not path.exists(p_path) else "a"
# Save pressure
viz_p = HDF5File(MPI.comm_world, p_path, file_mode=file_mode)
viz_p.write(p_, "/pressure", tstep)
viz_p.close()
# Save velocity
viz_u = HDF5File(MPI.comm_world, u_path, file_mode=file_mode)
viz_u.write(U, "/velocity", tstep)
viz_u.close()
# Accumulate velocity
u_mean0.vector().axpy(1, u_[0].vector())
u_mean1.vector().axpy(1, u_[1].vector())
u_mean2.vector().axpy(1, u_[2].vector())
# Oasis hook called after the simulation has finished
def theend_hook(u_mean, u_mean0, u_mean1, u_mean2, T, dt, save_solution_at_tstep, save_solution_frequency,
**NS_namespace):
# get the file path
files = NS_parameters['files']
u_mean_path = files["u_mean"]
# divide the accumlated veloicty by the number of steps
NumTStepForAverage = (T / dt - save_solution_at_tstep) / save_solution_frequency + 1
u_mean0.vector()[:] = u_mean0.vector()[:] / NumTStepForAverage
u_mean1.vector()[:] = u_mean1.vector()[:] / NumTStepForAverage
u_mean2.vector()[:] = u_mean2.vector()[:] / NumTStepForAverage
assign(u_mean.sub(0), u_mean0)
assign(u_mean.sub(1), u_mean1)
assign(u_mean.sub(2), u_mean2)
# Save u_mean
with HDF5File(MPI.comm_world, u_mean_path, "w") as u_mean_file:
u_mean_file.write(u_mean, "u_mean")
def beta(err, p):
"""
Adjusted choice of beta for the dual-pressure boundary condition.
Ramped up to desired value if flow rate error (err) increases
Args:
err (float): Flow split error
p (float): Pressure value
Returns:
beta (float): Variable factor in flow split method
"""
if p < 0:
if err >= 0.1:
return 0.5
else:
return 1 - 5 * err ** 2
else:
if err >= 0.1:
return 1.5
else:
return 1 + 5 * err ** 2
def update_pressure_condition(NS_expressions, area_ratio, boundary, id_in, id_out, mesh, n, tstep, u_, mesh_path, ui, A, u_ab, u, v, q_1, b, x_1, **NS_namespace):
"""
Use a dual-pressure boundary condition as pressure condition at outlet.
"""
# velocity_tentative_hook(mesh, mesh_path, id_in, id_out, ui, A, u_ab, u, v, q_1, b, x_1, boundary, **NS_namespace)
Q_in = abs(assemble(dot(u_, n) * ds(id_in[0], domain=mesh, subdomain_data=boundary)))
Q_outs = []
Q_ideals = []
for i, out_id in enumerate(id_out):
Q_out = abs(assemble(dot(u_, n) * ds(out_id, domain=mesh, subdomain_data=boundary)))
Q_outs.append(Q_out)
Q_ideal = area_ratio[i] * Q_in
Q_ideals.append(Q_ideal)
#NS_expressions[out_id].p=NS_expressions[out_id].p*(Q_ideal/Q_out)
p_old = NS_expressions[out_id].p
R_optimal = area_ratio[i]
R_actual = Q_out / Q_in
M_err = abs(R_optimal / R_actual)
R_err = abs(R_optimal - R_actual)
if p_old < 0:
E = 1 + R_err / R_optimal
else:
E = -1 * (1 + R_err / R_optimal)
# 1) Linear update to converge first 100 tsteps of first cycle
delta = (R_optimal - R_actual) / R_optimal
if tstep < 100:
h = 0.1
if p_old > 1 and delta < 0:
NS_expressions[out_id].p = p_old
else:
NS_expressions[out_id].p = p_old * (1 - delta * h)
# 2) Dual pressure BC
else:
if p_old > 2 and delta < 0:
NS_expressions[out_id].p = p_old
else:
NS_expressions[out_id].p = p_old * beta(R_err, p_old) * M_err ** E
return Q_ideals, Q_in, Q_outs
def print_mesh_information(mesh):
comm = MPI.comm_world
local_xmin = mesh.coordinates()[:, 0].min()
local_xmax = mesh.coordinates()[:, 0].max()
local_ymin = mesh.coordinates()[:, 1].min()
local_ymax = mesh.coordinates()[:, 1].max()
local_zmin = mesh.coordinates()[:, 2].min()
local_zmax = mesh.coordinates()[:, 2].max()
xmin = comm.gather(local_xmin, 0)
xmax = comm.gather(local_xmax, 0)
ymin = comm.gather(local_ymin, 0)
ymax = comm.gather(local_ymax, 0)
zmin = comm.gather(local_zmin, 0)
zmax = comm.gather(local_zmax, 0)
local_num_cells = mesh.num_cells()
local_num_edges = mesh.num_edges()
local_num_faces = mesh.num_faces()
local_num_facets = mesh.num_facets()
local_num_vertices = mesh.num_vertices()
num_cells = comm.gather(local_num_cells, 0)
num_edges = comm.gather(local_num_edges, 0)
num_faces = comm.gather(local_num_faces, 0)
num_facets = comm.gather(local_num_facets, 0)
num_vertices = comm.gather(local_num_vertices, 0)
volume = assemble(Constant(1) * dx(mesh))
if MPI.rank(MPI.comm_world) == 0:
print("=== Mesh information ===")
print("X range: {} to {} (delta: {:.4f})".format(min(xmin), max(xmax), max(xmax) - min(xmin)))
print("Y range: {} to {} (delta: {:.4f})".format(min(ymin), max(ymax), max(ymax) - min(ymin)))
print("Z range: {} to {} (delta: {:.4f})".format(min(zmin), max(zmax), max(zmax) - min(zmin)))
print("Number of cells: {}".format(sum(num_cells)))
print("Number of cells per processor: {}".format(int(np.mean(num_cells))))
print("Number of edges: {}".format(sum(num_edges)))
print("Number of faces: {}".format(sum(num_faces)))
print("Number of facets: {}".format(sum(num_facets)))
print("Number of vertices: {}".format(sum(num_vertices)))
print("Volume: {:.4f}".format(volume))
print("Number of cells per volume: {:.4f}".format(sum(num_cells) / volume))
def get_file_paths(folder):
# Create folder where data and solutions (velocity, mesh, pressure) is stored
common_path = path.join(folder, "Solutions")
if MPI.rank(MPI.comm_world) == 0:
if not path.exists(common_path):
makedirs(common_path)
file_p = path.join(common_path, "p.h5")
file_u = path.join(common_path, "u.h5")
file_u_mean = path.join(common_path, "u_mean.h5")
file_mesh = path.join(common_path, "mesh.h5")
files = {"u": file_u, "u_mean": file_u_mean, "p": file_p, "mesh": file_mesh}
return files
def u_dot_n(u, n):
return (dot(u, n) - abs(dot(u, n))) / 2