Merge pull request #189 from keiyamamo/move_innerP

move InnerP to simulation_common and clean terminology
KVSlab · Nov 4, 2024 · 0b54e24 · 0b54e24
2 parents da7223d + d34352b
commit 0b54e24
Show file tree

Hide file tree

Showing 4 changed files with 109 additions and 359 deletions.
diff --git a/src/vasp/simulations/aneurysm.py b/src/vasp/simulations/aneurysm.py
@@ -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
@@ -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
@@ -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]
@@ -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)
@@ -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):
@@ -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)
@@ -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,
@@ -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)
diff --git a/src/vasp/simulations/offset_stenosis.py b/src/vasp/simulations/offset_stenosis.py
@@ -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
@@ -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
@@ -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
     ))
 
@@ -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
@@ -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)
@@ -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)
@@ -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)
@@ -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):