Merge pull request #208 from keiyamamo/update_stenosis_tutorial

Update aneurysm turotial

docs/aneurysm.md

@@ -1,3 +1,157 @@
 (tutorial:aneurysm)=
 
-# Cerebral aneurysm simulation
+# Cerebral aneurysm simulation 
+
+A cerebral aneurysm is a pathological dilation of an artery in the brain. `VaSP` originated from research on cerebral aneurysms via high-fidelity fluid-structure interaction (FSI) simulations. For more context, please refer to the previous publications by Souche et al. {cite}`Souche2022` and Bruneau et al. {cite}`Bruneau2023`. In this tutorial, we focus on meshing a cerebral aneurysm with interactive and local refinement and how to set up a high-fidelity FSI simulation under physiological conditions.
+
+## Meshing with interactive and local refinement ##
+
+```{note}
+In this meshing section, we will use the data from [AneuriskWeb](http://ecm2.mathcs.emory.edu/aneuriskweb/index). Specifiacally, we will use the geometry `C0002` which is a 3D model of a cerebral aneurysm and can be downloaded from [here](http://ecm2.mathcs.emory.edu/aneuriskweb/repository#C0002).
+```
+
+In the previous tutorial, we used an objective method for determining the local mesh density, namely diameter-based meshing. In this tutorial, we will introduce a more interactive method for refining the mesh locally. This method is particularly useful when you want to refine the mesh in specific regions of the geometry, such as the aneurysm sac. To interactively refine the mesh, run:
+
+```console
+vasp-generate-mesh -i C0002/surface/model.vtp -sch5 0.001 -st constant -m distancetospheres  -mp 0 0.2 0.4 0.7 0 0.1 0.6 0.7 -fli 5 -flo 2 
+```
+
+The flag `-sch5` or (`--scale-factor-h5`) specifies the scaling factor from the original surface mesh to the volumetric mesh with `.h5` extension. The flag `-st constant` (or `--solid-thickness constant`) specifies that the wall thickness of the mesh should be constant across the geometry. The flag `-m distancetospheres` (or `--meshing-method distancetospheres`) specifies the method for meshing the geometry. The flag `-mp` (or `--meshing-parameters`) specifies the parameters for the distance-to-sphere scaling function. When `-m distancetospheres` is selected, you must define multiple sets of four parameters for the distance-to-sphere scaling function: offset, scale, min, and max. `--meshing-parameters 0 0.2 0.4 0.7 0 0.1 0.6 0.7` will run the function twice and the render window will pop up twice, where you can adjust the local mesh density interactively. The following gif shows an example of how to interactively refine the mesh:
+
+```{figure} figures/aneurysm_meshing.gif
+---
+width: 600px
+align: center
+name: aneurysm_meshing
+---
+An interactive method for specifying local mesh refinement, using VMTK. Press `d` on the keyboard to display the distances to the sphere, which will define the local mesh density.
+```
+
+```{figure} figures/aneurysm_mesh.png
+---
+align: center
+name: aneurysm_mesh
+---
+The mesh of a cerebral aneurysm generated with local refinement as shown in the previous gif.
+On the left, you can see the mesh with edges. On the right, you can see the mesh color coded by the local mesh size.
+``` 
+
+## Running a FSI simulation ##
+
+```{attention}
+To properly run high-fidelity FSI simulations of cerebral aneurysms, it is necessary to use a high-performance computing (HPC) cluster.
+```
+
+The simulation setup of a cerebral aneurysm is essentially the same as the offset stenosis case from the previous tutorial. Here, we focus on physiological boundary conditions and how they are defined in `aneurysm.py`.
+
+## Fluid boundary conditions ##
+
+First and foremost, we use Womersley velocity profile for the fluid inlet boundary condition where the mathematical expression is given by:
+
+```{math}
+u(r,t) = \frac{2C_0}{\pi R^2}\left[1-\left(\frac{r}{R}\right)^2 \right] + \sum_{n=1}^{N}\frac{C_n}{\pi R^2} \left[\frac{J_0(\alpha_n i^{\frac{3}{2}})-J_0(\alpha_n \frac{r}{R}i^{\frac{3}{2}})}{J_0(\alpha_n i^{\frac{3}{2}})-\frac{2}{\alpha_n i^{\frac{3}{2}}}J_1(\alpha_n i^{\frac{3}{2}})} \right]e^{in\omega t}
+```
+
+Here, $r$ and $t$ denote the cylindrical coordinates and time, respectively, and $R$ is the radius of the inlet cross-section. $J_{0}$ and $J_{1}$ are Bessel functions of the first kind of orders 0 and 1, and $\alpha_n$ is the Womersley number defined as $\alpha_n = R\sqrt{n \omega/\nu}$ where $\omega$ is the angular frequency of one cardiac cycle and $\nu$ is the kinematic viscosity of the fluid. Finally, $C_n$ are complex Fourier coefficients derived from the time-dependent flow rate as:
+
+```{math}
+    Q(t) = \sum_{n=0}^{N}C_n e^{in\omega t}
+```
+
+In `VaSP`, the coefficients $C_n$ are obtained from the internal carotid arteries of older adults {cite}`Hoi2010`, located at `src/vasp/simulations/FC_MCA10`
+
+In the `aneurysm.py` file, we define the Womersley velocity profile through `VaMPy` as follows:
+
+```python
+from vampy.simulation.Womersley import make_womersley_bcs, compute_boundary_geometry_acrn
+... # omitted code
+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(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)
+
+    # Create Womersley boundary condition at inlet
+    tmp_element = DVP.sub(1).sub(0).ufl_element()
+    inlet = make_womersley_bcs(T_Cycle, None, mu_f, tmp_center, tmp_radius, tmp_normal, tmp_element, Cn=Cn)
+    # Initialize inlet expressions with initial time
+    for uc in inlet:
+        uc.set_t(t)
+
+    # Create Boundary conditions for the velocity
+    u_inlet = [DirichletBC(DVP.sub(1).sub(i), inlet[i], boundaries, inlet_id) for i in range(3)]
+```
+`FC_file` is the file containing the Fourier coefficients, `Q_mean` is the mean flow rate, and `T_Cycle` is the cardiac cycle period (mostly 0.951 s). The function `compute_boundary_geometry_acrn` computes the center, radius, and normal of the fluid inlet. The function `make_womersley_bcs` returns `FEniCS` `Expression` objects for the Womersley velocity profile, which are then used to define the fluid inlet boundary condition using `DirichletBC`.
+
+Additionally, to avoid the sudden increase of flow at the beginning of the simulation, we use a ramp function to gradually increase the flow rate to the desired value. The following code shows how it is implemented in `aneurysm.py`:
+
+```python
+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)
+
+        # Multiply by cosine function to ramp up smoothly over time interval 0-250 ms
+        if t < 0.25:
+            uc.scale_value = -0.5 * np.cos(np.pi * t / 0.25) + 0.5
+        else:
+            uc.scale_value = 1.0
+```
+
+## Fluid-solid interface boundary conditions ##
+
+Secondly, we also apply arterial pressure at the fluid-solid interface. The waveform is the same as the flow rate waveform, but the pressure is scaled to vary between 70 and 110 mmHg. This boundary condition is implemented weakly by modifying the `FEniCS` variational form as follows:
+
+```python
+from vasp.simulations.simulation_common import InterfacePressure
+... # omitted code
+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 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)
+    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(J_(d_["n"]("+")) * inv(F_(d_["n"]("+"))).T * n("+")) * dSS(fsi_id)
+```
+
+Here, `P_FC_File` is the file containing the Fourier coefficients for the pressure waveform, and `P_mean` is the mean pressure value. The `InterfacePressure` class is defined in `VaSP` and returns an floating value (Pa) for the pulsatile pressure waveform. The pressure is applied weakly to the fluid-solid interface by modifying the variational form `F_solid_linear`. Note that the normal vector `n` needs to be updated as the fluid-solid interface moves during the simulation. To do that, we use Nanson's formula:
+
+```{math}
+n^{'} = J F^{-T} n
+```
+
+where $F$ is the deformation gradient, $J$ is the determinant of $F$, and $n^{'}$ is the updated normal vector.
+
+## Incorporating perivascular damping ##
+
+Finally, to incorporate viscoelastic damping from the perivascular environment, we use Robin boundary conditions at the solid outer wall, as introduced by Moireau et al. {cite}`Moireau2012`. This has been already implemented in `turtFSI` and can be activated by setting the respective parameter as follows:
+
+```python
+def set_problem_parameters(default_variables, **namespace):
+    default_variables.update(dict(
+    ... # omitted other parameters
+    robin_bc=True,
+    k_s=[1E5], # elastic coefficient [N*s/m^3]
+    c_s=[10], # viscous coefficient [N*s/m^3]
+    ds_s_id=[33] # surface ID for robin boundary condition
+    ))
+    return default_variables
+```
+
+Here, `k_s` and `c_s` are the elastic and viscous coefficients, respectively, and `ds_s_id` is the surface ID for the Robin boundary condition.
+
+```{bibliography}
+:filter: docname in docnames
+```

docs/avf.md

@@ -28,7 +28,7 @@ name: avf_meshing
 An interactive method for specifying variable wall thickness, using VMTK. Press `d` on the keyboard to display the distances to the sphere, which will define the local wall thickness. In this case, the red regions will have a wall thickness of 0.15, while the blue regions will have a wall thickness of 0.3.
 ```
 
-The third option `-eb` (or `--extract-branch`) option enables the extraction of specific branches from the mesh, such as the artery or vein in an arteriovenous fistula (AVF). When this option is enabled, `VaSP` assigns a unique solid mesh ID to the extracted branch, with an optional ID offset, controlled by `--branch-ids-offset`. This offset allows for the modification of the original solid ID number. For instance, if the original ID for a solid is 2 and the offset is set to 1000, the extracted branch will be assigned a new solid ID of 1002. This feature is particularly useful when you want to separately treat specific branches, like the artery and vein, with different wall properties. Again, with this option enabled, a render window will pop up, allowing you to specify the branch interactively.
+The third option `-eb` (or `--extract-branch`) option enables the extraction of specific branches from the mesh, such as the artery or vein in an AVF. When this option is enabled, `VaSP` assigns a unique solid mesh ID to the extracted branch, with an optional ID offset, controlled by `--branch-ids-offset`. This offset allows for the modification of the original solid ID number. For instance, if the original ID for a solid is 2 and the offset is set to 1000, the extracted branch will be assigned a new solid ID of 1002. This feature is particularly useful when you want to separately treat specific branches, like the artery and vein, with different wall properties. Again, with this option enabled, a render window will pop up, allowing you to specify the branch interactively.
 
 ```{figure} figures/avf_branch.gif
 ---
@@ -73,21 +73,23 @@ def set_problem_parameters(default_variables, **namespace):
 Moreover, we use patient-specific boundary conditions for the AVF model. The following code snippet shows how to define the boundary conditions for the AVF model:
 
 ```python
-# read patient-specific data
-patient_data = np.loadtxt(patient_data_path, skiprows=1, delimiter=",", usecols=(0, 1, 2))
-v_PA = patient_data[:, 0]
-v_DA = patient_data[:, 1]
-PV = patient_data[:, 2]
-
-len_v = len(v_PA)
-t_v = np.arange(len(v_PA))
-num_t = int(T / dt)  # 30.000 timesteps = 3s (T) / 0.0001s (dt)
-tnew = np.linspace(0, len_v, num=num_t)
-
-interp_DA = np.array(np.interp(tnew, t_v, v_DA))
-interp_PA = np.array(np.interp(tnew, t_v, v_PA))
-# pressure interpolation (velocity and pressure waveforms must be syncronized)
-interp_P = np.array(np.interp(tnew, t_v, PV))
+def create_bcs(DVP, mesh, boundaries, T, dt, fsi_id, inlet_id1, inlet_id2, rigid_id, psi, F_solid_linear,
+               vel_t_ramp, p_t_ramp_start, p_t_ramp_end, p_deg, v_deg, patient_data_path, **namespace):
+    # read patient-specific data
+    patient_data = np.loadtxt(patient_data_path, skiprows=1, delimiter=",", usecols=(0, 1, 2))
+    v_PA = patient_data[:, 0]
+    v_DA = patient_data[:, 1]
+    PV = patient_data[:, 2]
+
+    len_v = len(v_PA)
+    t_v = np.arange(len(v_PA))
+    num_t = int(T / dt)  # 30.000 timesteps = 3s (T) / 0.0001s (dt)
+    tnew = np.linspace(0, len_v, num=num_t)
+
+    interp_DA = np.array(np.interp(tnew, t_v, v_DA))
+    interp_PA = np.array(np.interp(tnew, t_v, v_PA))
+    # pressure interpolation (velocity and pressure waveforms must be syncronized)
+    interp_P = np.array(np.interp(tnew, t_v, PV))
 ```
 
 In this code snippet, patient-specific boundary conditions are defined by interpolating measured data to match the simulation's temporal resolution. The patient data is read from a CSV file, which extracts the velocity values for the proximal artery (`v_PA`), distal artery (`v_DA`), and the pressure values (`PV`). The time indices for the original dataset are then generated as t_v.

docs/installation.md

@@ -25,3 +25,6 @@ with `conda` can be found [here](install:conda).
 To install `VaSP` and all its dependencies to an *isolated environment*, we recommend using the dedicated Docker
 container. Instructions for installing `VaSP` with Docker can be found [here](install:docker).
 
+## Installing in high-performance computing (HPC) clusters
+
+To install and use `VaSP`, it is recommended to first install `FEniCS` on the cluster and then install `turtleFSI`, `VaMPy`, and `VaSP` separately. To install `FEniCS`, it is recommended to build it from source and instructions can be found [here](https://fenics.readthedocs.io/en/latest/installation.html). After installing `FEniCS`, you can install `turtleFSI`, `VaMPy`, and `VaSP` via `pip`.

docs/references.bib

@@ -55,4 +55,54 @@ @article{Soliveri2024
 publisher = {Springer Berlin Heidelberg},
 title = {{Toward a physiological model of vascular wall vibrations in the arteriovenous fistula}},
 year = {2024}
+}
+
+@article{Hoi2010,
+author = {Hoi, Yiemeng and Wasserman, Bruce A and Xie, Yuanyuan J and Najjar, Samer S and Ferruci, Luigi and Lakatta, Edward G and Gerstenblith, Gary and Steinman, David A},
+doi = {10.1088/0967-3334/31/3/002},
+issn = {0967-3334},
+journal = {Physiological Measurement},
+month = {mar},
+number = {3},
+pages = {291--302},
+title = {{Characterization of volumetric flow rate waveforms at the carotid bifurcations of older adults}},
+volume = {31},
+year = {2010}
+}
+
+@article{Moireau2012,
+author = {Moireau, P. and Xiao, N. and Astorino, M. and Figueroa, C. A. and Chapelle, D. and Taylor, C. A. and Gerbeau, J. F.},
+doi = {10.1007/s10237-011-0289-z},
+issn = {16177959},
+journal = {Biomechanics and Modeling in Mechanobiology},
+keywords = {Aorta modeling,Blood flows,Boundary conditions,Fluid-structure interaction,Patient-specific geometry,Tissue support},
+number = {1-2},
+pages = {1--18},
+pmid = {21308393},
+title = {{External tissue support and fluid-structure simulation in blood flows}},
+volume = {11},
+year = {2012}
+}
+
+@article{Souche2022,
+author = {Souche, Alban and Valen-Sendstad, Kristian},
+doi = {10.1016/j.jbiomech.2022.111369},
+issn = {00219290},
+journal = {Journal of Biomechanics},
+month = {dec},
+pages = {111369},
+title = {{High-fidelity fluid structure interaction simulations of turbulent-like aneurysm flows reveals high-frequency narrowband wall vibrations: A stimulus of mechanobiological relevance?}},
+volume = {145},
+year = {2022}
+}
+
+@article{Bruneau2023,
+author = {Bruneau, David A. and Steinman, David A. and Valen-Sendstad, Kristian},
+doi = {10.1038/s43856-023-00396-5},
+journal = {Communications Medicine},
+number = {1},
+pages = {1--11},
+title = {{Understanding intracranial aneurysm sounds via high-fidelity fluid-structure-interaction modelling}},
+volume = {3},
+year = {2023}
 }