Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

move InnerP to simulation_common and clean terminology #189

Merged
merged 5 commits into from
Nov 4, 2024
Merged

move InnerP to simulation_common and clean terminology #189

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
5 commits
Select commit Hold shift + click to select a range
a7b8b7a
move interface pressure to simulation_common
keiyamamo Sep 18, 2024
570df6d
remove innerP from offset stenosis
keiyamamo Sep 18, 2024
a72d2d5
remove innerP from aneurysm
keiyamamo Sep 18, 2024
143bf0e
switch from os to pathlib
keiyamamo Nov 1, 2024
d34352b
remove offset_stenosis_box and merge with the original offset_stenosis
keiyamamo Nov 1, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
76 changes: 20 additions & 56 deletions src/vasp/simulations/aneurysm.py
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,16 +1,17 @@
"""
Problem file for aneurysm FSI simulation
"""
import os
from pathlib import Path
import numpy as np

from vampy.simulation.Womersley import make_womersley_bcs, compute_boundary_geometry_acrn
from vampy.simulation.simulation_common import print_mesh_information
from turtleFSI.problems import *
from dolfin import HDF5File, Mesh, MeshFunction, facets, assemble, UserExpression, sqrt, FacetNormal, ds, \
from dolfin import HDF5File, Mesh, MeshFunction, facets, assemble, sqrt, FacetNormal, ds, \
DirichletBC, Measure, inner, parameters, VectorFunctionSpace, Function, XDMFFile

from vasp.simulations.simulation_common import load_probe_points, calculate_and_print_flow_properties, \
print_probe_points
print_probe_points, InterfacePressure

# set compiler arguments
parameters["form_compiler"]["quadrature_degree"] = 6
Expand All @@ -33,7 +34,7 @@ def set_problem_parameters(default_variables, **namespace):

default_variables.update(dict(
# Temporal parameters
T=1.902, # Simulation end time
T=0.002, # Simulation end time
dt=0.001, # Timne step size
theta=0.501, # Theta scheme parameter
save_step=1, # Save frequency of files for visualisation
Expand Down Expand Up @@ -94,6 +95,8 @@ def get_mesh_domain_and_boundaries(mesh_path, fsi_region, fsi_id, rigid_id, oute
domains = MeshFunction("size_t", mesh, 3)
hdf.read(domains, "/domains")

print_mesh_information(mesh)

# Only consider FSI in domain within this sphere
sph_x = fsi_region[0]
sph_y = fsi_region[1]
Expand All @@ -113,53 +116,13 @@ def get_mesh_domain_and_boundaries(mesh_path, fsi_region, fsi_id, rigid_id, oute
return mesh, domains, boundaries


class InnerP(UserExpression):
def __init__(self, t, t_ramp, An, Bn, period, P_mean, **kwargs):
self.t = t
self.t_ramp = t_ramp
self.An = An
self.Bn = Bn
self.omega = (2.0 * np.pi / period)
self.P_mean = P_mean
self.p_0 = 0.0 # Initial pressure
self.P = self.p_0 # Apply initial pressure to inner pressure variable
super().__init__(**kwargs)

def update(self, t):
self.t = t
# apply a sigmoid ramp to the pressure
if self.t < self.t_ramp:
ramp_factor = -0.5 * np.cos(np.pi * self.t / self.t_ramp) + 0.5
else:
ramp_factor = 1.0
if MPI.rank(MPI.comm_world) == 0:
print("ramp_factor = {} m^3/s".format(ramp_factor))

# Caclulate Pn (normalized pressure)from Fourier Coefficients
Pn = 0 + 0j
for i in range(len(self.An)):
Pn = Pn + (self.An[i] - self.Bn[i] * 1j) * np.exp(1j * i * self.omega * self.t)
Pn = abs(Pn)

# Multiply by mean pressure and ramp factor
self.P = ramp_factor * Pn * self.P_mean
if MPI.rank(MPI.comm_world) == 0:
print("P = {} Pa".format(self.P))

def eval(self, value, x):
value[0] = self.P

def value_shape(self):
return ()


def create_bcs(t, DVP, mesh, boundaries, mu_f,
fsi_id, inlet_id, inlet_outlet_s_id,
rigid_id, psi, F_solid_linear, p_deg, FC_file,
Q_mean, P_FC_File, P_mean, T_Cycle, **namespace):

# Load fourier coefficients for the velocity and scale by flow rate
An, Bn = np.loadtxt(os.path.join(os.path.dirname(os.path.abspath(__file__)), FC_file)).T
An, Bn = np.loadtxt(Path(__file__).parent / FC_file).T
# Convert to complex fourier coefficients
Cn = (An - Bn * 1j) * Q_mean
_, tmp_center, tmp_radius, tmp_normal = compute_boundary_geometry_acrn(mesh, inlet_id, boundaries)
Expand All @@ -183,20 +146,21 @@ def create_bcs(t, DVP, mesh, boundaries, mu_f,
# Assemble boundary conditions
bcs = u_inlet + [d_inlet, u_inlet_s, d_inlet_s, d_rigid]

# Load Fourier coefficients for the pressure and scale by flow rate
An_P, Bn_P = np.loadtxt(os.path.join(os.path.dirname(os.path.abspath(__file__)), P_FC_File)).T
# Load Fourier coefficients for the pressure
An_P, Bn_P = np.loadtxt(Path(__file__).parent / P_FC_File).T

# Apply pulsatile pressure at the fsi interface by modifying the variational form
n = FacetNormal(mesh)
dSS = Measure("dS", domain=mesh, subdomain_data=boundaries)
p_out_bc_val = InnerP(t=0.0, t_ramp=0.2, An=An_P, Bn=Bn_P, period=T_Cycle, P_mean=P_mean, degree=p_deg)
F_solid_linear += p_out_bc_val * inner(n('+'), psi('+')) * dSS(fsi_id)
interface_pressure = InterfacePressure(t=0.0, t_ramp_start=0.0, t_ramp_end=0.2, An=An_P,
Bn=Bn_P, period=T_Cycle, P_mean=P_mean, degree=p_deg)
F_solid_linear += interface_pressure * inner(n('+'), psi('+')) * dSS(fsi_id)

# Create inlet subdomain for computing the flow rate inside post_solve
dsi = ds(inlet_id, domain=mesh, subdomain_data=boundaries)
inlet_area = assemble(1.0 * dsi)
return dict(bcs=bcs, inlet=inlet, p_out_bc_val=p_out_bc_val, F_solid_linear=F_solid_linear, n=n, dsi=dsi,
inlet_area=inlet_area)
return dict(bcs=bcs, inlet=inlet, interface_pressure=interface_pressure, F_solid_linear=F_solid_linear, n=n,
dsi=dsi, inlet_area=inlet_area)


def initiate(mesh_path, scale_probe, mesh, v_deg, p_deg, **namespace):
Expand All @@ -215,7 +179,7 @@ def initiate(mesh_path, scale_probe, mesh, v_deg, p_deg, **namespace):
return dict(probe_points=probe_points, d_mean=d_mean, u_mean=u_mean, p_mean=p_mean)


def pre_solve(t, inlet, p_out_bc_val, **namespace):
def pre_solve(t, inlet, interface_pressure, **namespace):
for uc in inlet:
# Update the time variable used for the inlet boundary condition
uc.set_t(t)
Expand All @@ -227,9 +191,9 @@ def pre_solve(t, inlet, p_out_bc_val, **namespace):
uc.scale_value = 1.0

# Update pressure condition
p_out_bc_val.update(t)
interface_pressure.update(t)

return dict(inlet=inlet, p_out_bc_val=p_out_bc_val)
return dict(inlet=inlet, interface_pressure=interface_pressure)


def post_solve(dvp_, n, dsi, dt, mesh, inlet_area, mu_f, rho_f, probe_points, t,
Expand Down Expand Up @@ -265,6 +229,6 @@ def finished(d_mean, u_mean, p_mean, visualization_folder, save_solution_after_t
]

for vector, data_name in data_names:
file_path = os.path.join(visualization_folder, data_name)
with XDMFFile(MPI.comm_world, file_path) as f:
file_path = Path(visualization_folder) / data_name
with XDMFFile(MPI.comm_world, str(file_path)) as f:
f.write_checkpoint(vector, data_name, 0, XDMFFile.Encoding.HDF5)
94 changes: 35 additions & 59 deletions src/vasp/simulations/offset_stenosis.py
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,17 +1,17 @@
"""
Problem file for offset stenosis FSI simulation
"""
import os
from pathlib import Path
import numpy as np

from vampy.simulation.Womersley import make_womersley_bcs, compute_boundary_geometry_acrn
from vampy.simulation.simulation_common import print_mesh_information
from turtleFSI.problems import *
from dolfin import HDF5File, Mesh, MeshFunction, facets, cells, UserExpression, FacetNormal, ds, \
DirichletBC, Measure, inner, parameters, assemble
from turtleFSI.problems import * # noqa: F401
from dolfin import HDF5File, Mesh, MeshFunction, facets, assemble, sqrt, cells, FacetNormal, ds, \
DirichletBC, Measure, inner, parameters

from vasp.simulations.simulation_common import load_probe_points, print_probe_points, \
calculate_and_print_flow_properties, load_solid_probe_points, print_solid_probe_points
calculate_and_print_flow_properties, load_solid_probe_points, print_solid_probe_points, InterfacePressure

# set compiler arguments
parameters["form_compiler"]["quadrature_degree"] = 6
Expand Down Expand Up @@ -43,8 +43,8 @@ def set_problem_parameters(default_variables, **namespace):
linear_solver="mumps",
atol=1e-6, # Absolute tolerance in the Newton solver
rtol=1e-6, # Relative tolerance in the Newton solver
recompute=20, # Recompute the Jacobian matix within time steps
recompute_tstep=20, # Recompute the Jacobian matix over time steps
recompute=20, # Recompute the Jacobian matrix within time steps
recompute_tstep=20, # Recompute the Jacobian matrix over time steps
# boundary condition parameters
inlet_id=3, # inlet id for the fluid
inlet_outlet_s_id=11, # inlet and outlet id for solid
Expand Down Expand Up @@ -74,7 +74,7 @@ def set_problem_parameters(default_variables, **namespace):
mesh_path="mesh/file_stenosis.h5",
FC_file="FC_MCA_10", # File name containing the Fourier coefficients for the flow waveform
P_FC_File="FC_Pressure", # File name containing the Fourier coefficients for the pressure waveform
compiler_parameters=_compiler_parameters, # Update the defaul values of the compiler arguments (FEniCS)
compiler_parameters=_compiler_parameters, # Update the default values of the compiler arguments (FEniCS)
save_deg=2, # Degree of the functions saved for visualisation
))

Expand Down Expand Up @@ -105,11 +105,26 @@ def get_mesh_domain_and_boundaries(mesh_path, fsi_region, dx_f_id, fsi_id, rigid
idx_facet = boundaries.array()[i]
if idx_facet == fsi_id or idx_facet == outer_id:
mid = submesh_facet.midpoint()
dist_sph_center = np.sqrt((mid.x() - sph_x) ** 2 + (mid.y() - sph_y) ** 2 + (mid.z() - sph_z) ** 2)
dist_sph_center = sqrt((mid.x() - sph_x) ** 2 + (mid.y() - sph_y) ** 2 + (mid.z() - sph_z) ** 2)
if dist_sph_center > sph_rad:
boundaries.array()[i] = rigid_id # changed "fsi" idx to "rigid wall" idx
i += 1

# NOTE: Instead of using a sphere, we can also use a box to define the FSI region as follows
# Only consider FSI in domain within fsi_x_range, fsi_y_range, fsi_z_range
# i = 0
# for submesh_facet in facets(mesh):
# idx_facet = boundaries.array()[i]
# if idx_facet == fsi_id or idx_facet == outer_id:
# mid = submesh_facet.midpoint()
# if mid.x() < fsi_region[0] or mid.x() > fsi_region[1]:
# boundaries.array()[i] = rigid_id # changed "fsi" id to "rigid wall" id
# elif mid.y() < fsi_region[2] or mid.y() > fsi_region[3]:
# boundaries.array()[i] = rigid_id
# elif mid.z() < fsi_region[4] or mid.z() > fsi_region[5]:
# boundaries.array()[i] = rigid_id
# i += 1

# In this region, make fluid more viscous
x_min = 0.024
i = 0
Expand All @@ -124,46 +139,6 @@ def get_mesh_domain_and_boundaries(mesh_path, fsi_region, dx_f_id, fsi_id, rigid
return mesh, domains, boundaries


class InnerP(UserExpression):
def __init__(self, t, t_ramp, An, Bn, period, P_mean, **kwargs):
self.t = t
self.t_ramp = t_ramp
self.An = An
self.Bn = Bn
self.omega = (2.0 * np.pi / period)
self.P_mean = P_mean
self.p_0 = 0.0 # Initial pressure
self.P = self.p_0 # Apply initial pressure to inner pressure variable
super().__init__(**kwargs)

def update(self, t):
self.t = t
# apply a sigmoid ramp to the pressure
if self.t < self.t_ramp:
ramp_factor = -0.5 * np.cos(np.pi * self.t / self.t_ramp) + 0.5
else:
ramp_factor = 1.0
if MPI.rank(MPI.comm_world) == 0:
print("ramp_factor = {} m^3/s".format(ramp_factor))

# Caclulate Pn (normalized pressure)from Fourier Coefficients
Pn = 0 + 0j
for i in range(len(self.An)):
Pn = Pn + (self.An[i] - self.Bn[i] * 1j) * np.exp(1j * i * self.omega * self.t)
Pn = abs(Pn)

# Multiply by mean pressure and ramp factor
self.P = ramp_factor * Pn * self.P_mean
if MPI.rank(MPI.comm_world) == 0:
print("P = {} Pa".format(self.P))

def eval(self, value, x):
value[0] = self.P

def value_shape(self):
return ()


def initiate(mesh_path, **namespace):

probe_points = load_probe_points(mesh_path)
Expand All @@ -178,7 +153,7 @@ def create_bcs(t, DVP, mesh, boundaries, mu_f,
Q_mean, P_FC_File, P_mean, T_Cycle, **namespace):

# Load Fourier coefficients for the velocity and scale by flow rate
An, Bn = np.loadtxt(os.path.join(os.path.dirname(os.path.abspath(__file__)), FC_file)).T
An, Bn = np.loadtxt(Path(__file__).parent / FC_file).T
# Convert to complex fourier coefficients
Cn = (An - Bn * 1j) * Q_mean
_, tmp_center, tmp_radius, tmp_normal = compute_boundary_geometry_acrn(mesh, inlet_id, boundaries)
Expand All @@ -202,23 +177,24 @@ def create_bcs(t, DVP, mesh, boundaries, mu_f,
# Assemble boundary conditions
bcs = u_inlet + [d_inlet, u_inlet_s, d_inlet_s, d_rigid]

# Load Fourier coefficients for the pressure and scale by flow rate
An_P, Bn_P = np.loadtxt(os.path.join(os.path.dirname(os.path.abspath(__file__)), P_FC_File)).T
# Load Fourier coefficients for the pressure
An_P, Bn_P = np.loadtxt(Path(__file__).parent / P_FC_File).T

# Apply pulsatile pressure at the fsi interface by modifying the variational form
n = FacetNormal(mesh)
dSS = Measure("dS", domain=mesh, subdomain_data=boundaries)
p_out_bc_val = InnerP(t=0.0, t_ramp=0.2, An=An_P, Bn=Bn_P, period=T_Cycle, P_mean=P_mean, degree=p_deg)
F_solid_linear += p_out_bc_val * inner(n('+'), psi('+')) * dSS(fsi_id)
interface_pressure = InterfacePressure(t=0.0, t_ramp_start=0.0, t_ramp_end=0.2, An=An_P,
Bn=Bn_P, period=T_Cycle, P_mean=P_mean, degree=p_deg)
F_solid_linear += interface_pressure * inner(n('+'), psi('+')) * dSS(fsi_id)

# Create inlet subdomain for computing the flow rate inside post_solve
dsi = ds(inlet_id, domain=mesh, subdomain_data=boundaries)
inlet_area = assemble(1.0 * dsi)
return dict(bcs=bcs, inlet=inlet, p_out_bc_val=p_out_bc_val, F_solid_linear=F_solid_linear, n=n, dsi=dsi,
inlet_area=inlet_area)
return dict(bcs=bcs, inlet=inlet, interface_pressure=interface_pressure, F_solid_linear=F_solid_linear, n=n,
dsi=dsi, inlet_area=inlet_area)


def pre_solve(t, inlet, p_out_bc_val, **namespace):
def pre_solve(t, inlet, interface_pressure, **namespace):
for uc in inlet:
# Update the time variable used for the inlet boundary condition
uc.set_t(t)
Expand All @@ -230,9 +206,9 @@ def pre_solve(t, inlet, p_out_bc_val, **namespace):
uc.scale_value = 1.0

# Update pressure condition
p_out_bc_val.update(t)
interface_pressure.update(t)

return dict(inlet=inlet, p_out_bc_val=p_out_bc_val)
return dict(inlet=inlet, interface_pressure=interface_pressure)


def post_solve(probe_points, solid_probe_points, dvp_, dt, mesh, inlet_area, dsi, mu_f, rho_f, n, **namespace):
Expand Down
Loading
Loading