Skip to content

use Robin by default and remove FSI region #231

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

Merged
keiyamamo merged 3 commits into from
Aug 24, 2025
Merged

use Robin by default and remove FSI region #231

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
76d426b
use Robin by default and remove FSI region
keiyamamo Aug 13, 2025
39c6f76
avoid hardcoded numbers in the pytsest
keiyamamo Aug 20, 2025
b95337f
minor fix of aneurysm.py
keiyamamo Aug 20, 2025
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
33 changes: 8 additions & 25 deletions src/vasp/simulations/aneurysm.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -10,7 +10,7 @@
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, sqrt, FacetNormal, ds, \
from dolfin import HDF5File, Mesh, MeshFunction, assemble, FacetNormal, ds, \
DirichletBC, Measure, inner, parameters, VectorFunctionSpace, Function, XDMFFile

from vasp.simulations.simulation_common import load_probe_points, calculate_and_print_flow_properties, \
Expand Down Expand Up @@ -53,7 +53,6 @@ def set_problem_parameters(default_variables, **namespace):
inlet_id=2, # inlet id for the fluid
inlet_outlet_s_id=11, # inlet and outlet id for solid
fsi_id=22, # id for fsi surface
rigid_id=11, # "rigid wall" id for the fluid
outer_id=33, # id for the outer surface of the solid
# Fluid parameters
Q_mean=1.25E-06,
Expand All @@ -71,9 +70,10 @@ def set_problem_parameters(default_variables, **namespace):
nu_s=nu_s_val, # Solid Poisson ratio [-]
lambda_s=lambda_s_val, # Solid Young's modulus [Pa]
dx_s_id=2, # ID of marker in the solid domain
# FSI parameters
fsi_region=[0.123, 0.134, 0.063, 0.004], # x, y, and z coordinate of FSI region center,
# and radius of FSI region sphere
k_s=[1E5], # elastic constant of Robin [N/m^3]
c_s=[10], # viscous constant of Robin [N.s/m^3]
ds_s_id=[33], # ID of marker for the Robin boundary condition
robin_bc=True, # Use Robin boundary condition for the solid
# Simulation parameters
folder="aneurysm_results", # Folder name generated for the simulation
mesh_path="mesh/file_aneurysm.h5",
Expand All @@ -87,7 +87,7 @@ def set_problem_parameters(default_variables, **namespace):
return default_variables


def get_mesh_domain_and_boundaries(mesh_path, fsi_region, fsi_id, rigid_id, outer_id, **namespace):
def get_mesh_domain_and_boundaries(mesh_path, **namespace):

# Read mesh
mesh = Mesh()
Expand All @@ -100,28 +100,12 @@ def get_mesh_domain_and_boundaries(mesh_path, fsi_region, fsi_id, rigid_id, oute

print_mesh_information(mesh)

# Only consider FSI in domain within this sphere
sph_x = fsi_region[0]
sph_y = fsi_region[1]
sph_z = fsi_region[2]
sph_rad = fsi_region[3]

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()
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

return mesh, domains, boundaries


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,
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
Expand All @@ -144,10 +128,9 @@ def create_bcs(t, DVP, mesh, boundaries, mu_f,
# Solid Displacement BCs
d_inlet = DirichletBC(DVP.sub(0), (0.0, 0.0, 0.0), boundaries, inlet_id)
d_inlet_s = DirichletBC(DVP.sub(0), (0.0, 0.0, 0.0), boundaries, inlet_outlet_s_id)
d_rigid = DirichletBC(DVP.sub(0), (0.0, 0.0, 0.0), boundaries, rigid_id)

# Assemble boundary conditions
bcs = u_inlet + [d_inlet, u_inlet_s, d_inlet_s, d_rigid]
bcs = u_inlet + [d_inlet, u_inlet_s, d_inlet_s]

# Load Fourier coefficients for the pressure
An_P, Bn_P = np.loadtxt(Path(__file__).parent / P_FC_File).T
Expand Down
112 changes: 13 additions & 99 deletions tests/test_simulations.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -65,18 +65,7 @@ def test_predeform_problem(input_mesh, tmpdir):
cmd = ("turtleFSI -p predeform -dt 0.01 -T 0.03 --verbose True" +
f" --folder {tmpdir} --sub-folder 1 --new-arguments mesh_path={input_mesh}")
result = subprocess.check_output(cmd, shell=True, cwd="src/vasp/simulations/")

output_re = r"v \(centerline, at inlet\) = (\d+\.\d+|\d+) m/s"
output_match = re.findall(output_re, str(result))

assert output_match is not None, "Regular expression did not match the output."

expected_velocity = 0.009549150281252628
velocity_at_inlet = float(output_match[-1])

print("Velocity: {}".format(velocity_at_inlet))

assert np.isclose(velocity_at_inlet, expected_velocity), "Velocity does not match expected value."
check_velocity_cfl_reynolds(result)


@pytest.mark.parametrize("input_mesh", [input_data_paths[1]])
Expand All @@ -87,62 +76,7 @@ def test_cylinder_problem(input_mesh, tmpdir):
cmd = ("turtleFSI -p cylinder -dt 0.001 -T 0.004 --verbose True " +
f"--folder {tmpdir} --sub-folder 1 --new-arguments mesh_path={input_mesh}")
result = subprocess.check_output(cmd, shell=True, cwd="src/vasp/simulations/")
# check flow rate at inlet
output_flow_rate = r"Flow Rate at Inlet: (\d+(?:\.\d+)?(?:e-\d+)?)"

output_match_flow_rate = re.findall(output_flow_rate, str(result))
assert output_match_flow_rate is not None, "Regular expression did not match the output."

expected_flow_rate = 1.6913532412047225e-09
flow_rate_at_inlet = float(output_match_flow_rate[-1])

print(f"Flow rate: {flow_rate_at_inlet}")
assert np.isclose(flow_rate_at_inlet, expected_flow_rate), "Flow rate does not match expected value."

# check velocity mean, min, max in the domain
ourput_velocity = (r"Velocity \(mean, min, max\): (\d+(?:\.\d+)?(?:e-\d+)?)\s*,\s*(\d+(?:\.\d+)?(?:e-\d+)?)\s*,"
r"\s*(\d+(?:\.\d+)?(?:e-\d+)?)")

output_match_velocity = re.findall(ourput_velocity, str(result))
assert output_match_velocity is not None, "Regular expression did not match the output."

expected_velocity_mean_min_max = [0.0015175903693111237, 2.83149082127162e-06, 0.004025814882456499]
velocity_mean_min_max = [float(output_match_velocity[-1][0]), float(output_match_velocity[-1][1]),
float(output_match_velocity[-1][2])]

print(f"Velocity mean, min, max: {velocity_mean_min_max}")
assert np.isclose(velocity_mean_min_max, expected_velocity_mean_min_max).all(), \
"Velocity mean, min, max does not match expected value."

# check CFL number
output_cfl_number = (r"CFL \(mean, min, max\): (\d+(?:\.\d+)?(?:e-\d+)?)\s*,\s*(\d+(?:\.\d+)?(?:e-\d+)?)\s*,"
r"\s*(\d+(?:\.\d+)?(?:e-\d+)?)")

output_match_cfl_number = re.findall(output_cfl_number, str(result))
assert output_match_cfl_number is not None, "Regular expression did not match the output."

expected_cfl_number = [0.016930040421859752, 3.1587742666035394e-05, 0.044911466275218616]
cfl_number_mean_min_max = [float(output_match_cfl_number[-1][0]), float(output_match_cfl_number[-1][1]),
float(output_match_cfl_number[-1][2])]

print(f"CFL number mean, min, max: {cfl_number_mean_min_max}")
assert np.isclose(cfl_number_mean_min_max, expected_cfl_number).all(), \
"CFL number mean, min, max does not match expected value."

# Check Reynolds number
output_re_number = (r"Reynolds Numbers \(mean, min, max\): (\d+(?:\.\d+)?(?:e-\d+)?)\s*,"
r"\s*(\d+(?:\.\d+)?(?:e-\d+)?)\s*,\s*(\d+(?:\.\d+)?(?:e-\d+)?)")

output_match_re_number = re.findall(output_re_number, str(result))
assert output_match_re_number is not None, "Regular expression did not match the output."

expected_re_number = [0.4304434011992387, 0.0008031129903162388, 1.1418663992904048]
reynolds_number_mean_min_max = [float(output_match_re_number[-1][0]), float(output_match_re_number[-1][1]),
float(output_match_re_number[-1][2])]

print(f"Reynolds number mean, min, max: {reynolds_number_mean_min_max}")
assert np.isclose(reynolds_number_mean_min_max, expected_re_number).all(), \
"Reynolds number mean, min, max does not match expected value."
check_velocity_cfl_reynolds(result)


@pytest.mark.parametrize("input_mesh", [input_data_paths[2]])
Expand All @@ -153,32 +87,26 @@ def test_aneurysm_problem(input_mesh, tmpdir):
cmd = ("turtleFSI -p aneurysm -dt 0.001 -T 0.004 --verbose True " +
f"--folder {tmpdir} --sub-folder 1 --new-arguments inlet_id=4 mesh_path={input_mesh}")
result = subprocess.check_output(cmd, shell=True, cwd="src/vasp/simulations/")
# check flow rate at inlet
output_flow_rate = r"Flow Rate at Inlet: (\d+(?:\.\d+)?(?:e-\d+)?)"
check_velocity_cfl_reynolds(result)

output_match_flow_rate = re.findall(output_flow_rate, str(result))
assert output_match_flow_rate is not None, "Regular expression did not match the output."

expected_flow_rate = 7.297633240079062e-10
flow_rate_at_inlet = float(output_match_flow_rate[-1])

print(f"Flow rate: {flow_rate_at_inlet}")
assert np.isclose(flow_rate_at_inlet, expected_flow_rate), "Flow rate does not match expected value."

def check_velocity_cfl_reynolds(result):
"""
Sanity check of the velocity, CFL number, and Reynolds number in the result.
"""
# check velocity mean, min, max in the domain
output_velocity = (r"Velocity \(mean, min, max\): (\d+(?:\.\d+)?(?:e-\d+)?)\s*,\s*(\d+(?:\.\d+)?(?:e-\d+)?)\s*,"
r"\s*(\d+(?:\.\d+)?(?:e-\d+)?)")

output_match_velocity = re.findall(output_velocity, str(result))
assert output_match_velocity is not None, "Regular expression did not match the output."

expected_velocity_mean_min_max = [0.0007154906300607233, 6.665204824466191e-18, 0.002775071833322646]
velocity_mean_min_max = [float(output_match_velocity[-1][0]), float(output_match_velocity[-1][1]),
float(output_match_velocity[-1][2])]

print(f"Velocity mean, min, max: {velocity_mean_min_max}")
assert np.isclose(velocity_mean_min_max, expected_velocity_mean_min_max).all(), \
"Velocity mean, min, max does not match expected value."
assert all(np.isfinite(v) for v in velocity_mean_min_max), "Velocity mean, min, max should be finite values."
assert all(v >= 0 for v in velocity_mean_min_max), "Velocity mean, min, max should be non-negative."

# check CFL number
output_cfl_number = (r"CFL \(mean, min, max\): (\d+(?:\.\d+)?(?:e-\d+)?)\s*,\s*(\d+(?:\.\d+)?(?:e-\d+)?)\s*,"
Expand All @@ -187,25 +115,11 @@ def test_aneurysm_problem(input_mesh, tmpdir):
output_match_cfl_number = re.findall(output_cfl_number, str(result))
assert output_match_cfl_number is not None, "Regular expression did not match the output."

expected_cfl_number = [0.004760513375812616, 4.434690740353818e-17, 0.01846392674667467]
cfl_number_mean_min_max = [float(output_match_cfl_number[-1][0]), float(output_match_cfl_number[-1][1]),
float(output_match_cfl_number[-1][2])]

print(f"CFL number mean, min, max: {cfl_number_mean_min_max}")
assert np.isclose(cfl_number_mean_min_max, expected_cfl_number).all(), \
"CFL number mean, min, max does not match expected value."

# Check Reynolds number
output_re_number = (r"Reynolds Numbers \(mean, min, max\): (\d+(?:\.\d+)?(?:e-\d+)?)\s*,"
r"\s*(\d+(?:\.\d+)?(?:e-\d+)?)\s*,\s*(\d+(?:\.\d+)?(?:e-\d+)?)")

output_match_re_number = re.findall(output_re_number, str(result))
assert output_match_re_number is not None, "Regular expression did not match the output."

expected_re_number = [0.4637715859370615, 4.32029782385228e-15, 1.7987649469568674]
reynolds_number_mean_min_max = [float(output_match_re_number[-1][0]), float(output_match_re_number[-1][1]),
float(output_match_re_number[-1][2])]

print(f"Reynolds number mean, min, max: {reynolds_number_mean_min_max}")
assert np.isclose(reynolds_number_mean_min_max, expected_re_number).all(), \
"Reynolds number mean, min, max does not match expected value."
assert all(np.isfinite(cfl) for cfl in cfl_number_mean_min_max), \
"CFL number mean, min, max should be finite values."
assert all(cfl >= 0 for cfl in cfl_number_mean_min_max), \
"CFL number mean, min, max should be non-negative."
Loading