**Simulation Module in SimVascular** solves the three-dimensional incompressible Navier-Stokes equations in an arbitrary domain, generally a vascular model reconstructed from image data. This is a complex subject with extensive underlying theory, and therefore this document will focus mainly on the practical aspects of simulation and analysis. This module includes three parts: **Presolver(svPre), Flowsolver(svSolver), Postporcessing(svPost).**

The **svSolver** evolved from the academic finite element code PHASTA (Parallel, Hierarchical, Adaptive, Stabilized, Transient Analysis), developed mainly at RPI (Rensselaer Polytechnic Institute, Troy, NY) by Professor Kenneth Jansen. This code was in turned inspired by the Stabilized Finite Element theory developed by Professor Thomas J.R. Hughes during his Stanford years.

Building on the original PHASTA code, there have been a number of important additions and modifications. Professor Charles Taylor’s group at Stanford University developed key additions in the areas of Boundary Conditions and Fluid-Solid Interaction (FSI) coupling. These additions are crucial to represent with a high level of realism the way blood flows in arteries, since this flow is highly dependent on the characteristics of the vascular trees that are downstream of our three-dimensional model, and the compliance of the three-dimensional vascular tree.

Building on the above features, the Marsden lab at UCSD has added additional key functionality enabling efficient and stable solutions with models of the circulatory physiology:

**Backflow stabilization.**This problem usually arises in large vessels that are exposed to backflow in 3D and 2D flow simulations. This phenomenon may be a cause of divergence of the numerical scheme due to bulk reversal of the flow through an outlet, localized areas of flow reversal or use of a boundary 0D circulation model.Custom and efficient

**linear solver.**Accurate simulation of blood flow in vessels require the repeated solution of linear systems of equations with millions of unknowns. Moreover, use of closed-loop boundary models significantly increases the degree of coupling between boundary degrees of freedoms. The**svLS**linear solver is designed to efficiently handle large cardiovascular simulations with arbitrary boundary conditions and reduce solution times.Multiscale Coupling for

**closed loop boundary conditions.**Coupling a three-dimensional finite element solver with a 0D lumped circulation model drastically improves the possibility of realistically simulate patient-specific hemodynamics and phisiology.

This document will teach you the fundamentals of:

**svPre:**Preparing the necessary svSolver input files**svSolver:**Running a flow analysis**svPost:**Looking and providing interpretation to the results generated by the code

In addition, this tutorial will show you a number of good practices that are important to observe during the simulation process. We will do this considering very simple geometry (a straight cylinder) to illustrate different points in a simple way.

The theory and implementation details are not covered in this document. For more information about those details, please refer to our publications page.

The following figure contains a schematic representation of the processes involved in running a simulation using SimVascular.

We start off with the files coming from the meshing of the analysis: these files contain nodal and connectivity information for the finite element mesh, located in the the *mesh-complete/mesh-surfaces/* folder.

We then run **Presolver(svPre)** using the *.svpre_ file. The *.svpre file contains the set of instructions that define the boundary conditions, initial conditions, and geometrical configuration of our problem. The output of **svPre** is a set of files (**bct.dat, geombc.dat.1, restart.0.1, numstart.dat**) that are ready to be processed by **svSolver** to run a blood flow analysis. Running svSolver also need **solver.inp**, which provide further info for flowsolver.

Once the analysis is finished, the solver outputs files that characterize the finite element solution over a defined time period, typically several cardiac cycles. These files are taken by **svPost** to generate visualization files (typically *.vtu and *.vtp files) that are used to post-process and analyze the desired hemodynamic results.

In the following sections the components of this flow chart will be discussed in detail.

**svSolver**, just like many other Finite Element Programs, does not enforce a consistent set of physical units in the analysis, but it is up to the analyst to make sure that input data are dimensionally consistent.

To have a consistent set of units, users are advised to either work in CGS, MGS, or SI units; see the following table.

Quantity | CGS Unit | MGS Unit | SI Unit |
---|---|---|---|

Length | cm | mm | m |

Mass | gr | gr | Kg |

Time | s | s | s |

The following table gathers several important physical constants of blood given in different unit systems.

Property | CGS Unit | MGS Unit | SI Unit |
---|---|---|---|

Dynamic viscosity $\mu$ [M· L -1 · T -1 ] | 0.04 poise [gr· cm -1 · s -1 ] | 0.004 [gr· mm -1 · s -1 ] | 0.004 [Pa· s -1 ] |

Blood density $\rho$ [M· L -3 ] | 1.06 [gr· cm -3 ] | 0.00106 [gr· mm -3 ] | 1060 [Kg· m -3 ] |

Boundary conditions are crucial to obtaining high quality cardiovascular simulation results. It is essential that boundary conditions accurately capture the physiology of vascular networks outside of the 3D domain of the model. **SimVascular** provides several different options for boundary condition assignment at inlets and outlets that are described in this section. Typically, we begin with some physiologic information about our problem, for instance:

- Flow rate info coming from
**MRI**or**ultrasound**measurements. - Pressure values in the model obtained with a
**catheter**transducer or a pressure cuff. - Vessel wall elastic properties (if we are modeling the vessel walls as deformable).

From a conceptual standpoint, no matter how geometrically complex a vascular model is (we’ll refer to this as $\Omega$), its boundaries can be classified into three groups (see figure above):

- An
**inflow**boundary $\Gamma_g$. This is the set of face(s) of the model where we will prescribe a flow wave as obtained from a clinical measurement. - A vessel
**wall**boundary $\Gamma_s$. This boundary represents the interface between the fluid domain and the vessel wall. In the physical world, this boundary is lined by a layer of endothelial cells, the treatment of which can be complex. Many blood flow simulations have traditionally used a**rigid wall assumption**. Under these circumstances, a zero velocity condition is applied on these surfaces.**SimVascular**also offers options for fluid structure interaction as discussed below. - An
**outflow**boundary $\Gamma_h$. On this boundary, we will typically prescribe a pressure value that is uniform over the face (i.e. spatially not temporally constant) in a**weak manner**. A**weakly applied**pressure means that we are not prescribing nodal values of pressure at the nodes of the outlet face as Dirichlet boundary conditions. Instead, we apply this pressure by enforcing that the integral of the pressure field on that face must be a certain value.

These boundaries have an absolutely critical impact on the numerical simulation results. The SimVascular package contains several options for boundary condition assignment. All of these use a weakly prescribed pressure formulation, with the purpose of accounting for effects of the downstream vasculature on the vascular model (see figure below). These boundary conditions include:

**Resistance Boundary condition**. Here, the condition prescribed on the face is a relationship between flow and pressure of the form $p = p_0 + R\,Q$, where $R$ is the resistance parameter that characterizes the downstream vasculature, $p$ is the weakly prescribed pressure, $Q$ is the flow rate passing through the face and $p_0$ is a “flag” that sets the boundary as a “weakly-prescribed pressure boundary”. This flag has a “zero” numerical value, so the total value of the pressure on that face is simply given by $R\,Q$.**Impedance Boundary condition**. Here, the condition prescribed on the face is a relationship of the form:

$$ p(t)=p_0 + \frac{1}{T}\,\int_{t-T}^{t} Z(t−\tau) \, Q(\tau) \, d\tau $$

where $Z$ is the Inverse Fourier Transform of an impedance function that characterizes the downstream vasculature, $p$ is the weakly prescribed pressure, $Q$ is the flow rate passing through the face, and $T$ is the cardiac cycle.

The mathematical definition of an impedance function is:

$$ Z(\omega)=P(\omega)\,Q(\omega),\,\omega=0,1,2,\dots $$

That is, a relationship between pressure and flow modes for different frequencies. The figure above shows the typical shape of these impedances as a function of the frequency in the human iliac artery under rest and exercise conditions. Getting a good characterization of these impedance functions is once again critical to accurately represent the way blood flows in our computational domain.

**RCR lumped parameters condition**. In this boundary condition type, we use a*reduced-order*model of the downstream vasculature, considering an electric circuit analog. In this theory, the behavior of the vasculature is represented by three parameters: a proximal resistance $R_p$, a capacitance $C$, and a distal resistance $R_d$.

**Coronary boundary conditions**. These conditions simulate what happens at the coronary outlets. The implementation in the**svSolver**follows the derivations in this paper.**Closed-loop boundary circulation model**. The capability of coupling a 3D finite element model with a lumped parameter model is built into the**svSolver**. Documentation on this feature will be available with later releases of the code.

Before we move on, let us recap the type of *physical problem* (**Example 1**) we are solving: the geometry used in this problem consists of an idealized blood vessel, represented by a cylindrical segment with a radius $r=2$ cm and length $L=30$ cm. We prescribe an idealized constant inlet volumetric flow rate $Q$ of $100$ cc/sec to a parabolic profile at the inlet face of the model ($\Gamma_g$). The dynamic viscosity $\mu$ and density $\rho$ of the blood are 0.04 poise and 1.06 gr/cm$^3$, respectively. The lateral surface of the vessel ($\Gamma_{s}$) is considered to be rigid (therefore, we will apply a zero-velocity condition there). For the outlet boundary ($\Gamma_h$), a spatially-constant pressure boundary condition is weakly enforced via a resistance $R$.
In this problem, we will consider a resistance of $R = 1333.0$ dynes·s/cm$^5$.

This resistance will give a (weakly-applied) pressure at the outlet face of

$$ p=p_0 + R\,Q = 0.0 + 1333.0 \cdot 100.0=133300.00 \approx 100\,\text{mmHg} $$

(recall that $1.0$ mmHg = $1333.2$ dyn/cm$^2$). For steady flow in a long tube with a circular cross section, the flow will develop a flow profile known as the *Poiseuille* flow profile assuming the flow remains laminar. The flow will remain laminar in a circular tube assuming that the non-dimensional parameter given by the *Reynolds* number $Re$ is $Re < 2100$.

The definition of the Reynolds number is given by:

$$ Re = \frac{\rho\,D\,V}{\mu} $$

where $V$ is a representative velocity of the flow, $D$ is a characteristic dimension of the vessel where the flow is happening (in this case, the diameter), and $\mu$ and $\rho$ are the dynamic viscosity and density, respectively.

For this problem, the Reynolds number is about $884$, so it is well within the laminar flow regime. For a fully developed flow, the axisymmetric profile is parabolic and is given by:

$$ v_z(a) = v_z^{max}\left[1-\left(\frac{a}{r}\right)^2\right] $$

where $v_z^{max}$ is the maximum velocity at the center of the vessel, a is the radial coordinate from center of the vessel $0\le a \le r$ . The expression for maximum velocity is given by:

$$ v_z^{max} = 2\frac{Q}{\pi\,r^2} $$

where $Q$ is the volumetric flow rate. The wall shear stress $\tau$, is given by

$$ \tau = \frac{2\,\mu\,v_z^{max}}{r} $$

For this example, the maximum velocity is $v_z^{max} = 16.68$ cm/s and the wall shear stress is $\tau$ = $0.67$ dynes/cm$^2$.

**svPre** is the preprocessor for **svSolver**. The input files to svPre contain a complete description of the discrete model: nodal coordinates, element connectivity, element adjacency information and connectivity of boundary nodes and elements. Running **svPre** with an input *.svpre file with the appropriate commands, generates the input files for **svSolver**.

The **svPre** program is called either from the command line (in terminal) or using the **SimVascular** GUI. The input files for **svPresolver** are those generated by **Simvascular Meshing Module**. We will review this process briefly with a simple example of steady flow through a cylinder (Example 1). Before we start, first set the project folder as the example folder (…/cylinder)

These prerequisite files for svPre are generate by the output from Meshing Module (Click *Write Files* button in Mesh tab after meshing).

These files are:

in the **mesh-complete/** folder:

**mesh-complete.mesh.vtu**, vtu file containing the solid mesh generated with TetGen.**mesh-complete.exterior.vtp**, vtp file containig all the exterior elements of the mesh generated with TetGen.**walls_combined.vtp**, vtp file containing all surface elements assigned to the wall, possibily combined from various surfaces.

in the **mesh-complete/mesh-surfaces/** folder:

**inflow.vtp**, vtp file containing the meshed inlet surface.**outlet.vtp**, vtp file containing the meshed outlet surface.**wall.vtp**, vtp file containing the meshed wall surface.

The files for Example 1 can be found here. Create an empty folder on your hard drive to unzip the content of the archive. The following files are contained:

**HINT** - It is advisable that you set the project folder as *cylinder*. **SimVascular** will use this folder name as the default when creating new files. Using a meaningful folder name will make sure that your model files are named consistently. Also store the files containing the inlet flow rates in a folder called *flow-files*. Your problem may have more that one inflow wave form file. In this case, we only have a single flow file (called inflow.flow).

The format of the **steady.flow** file is as follows:

```
# Time (sec) Flow (cc/sec)
0 -100.0
1 -100.0
```

The first line is a comment line that you don’t need to include, but it helps to remember what units you used in the analysis. Then, two columns of numbers follow. The first column contains time values, and the second column flow values.

**WARNING**: please note that flow coming **into** the model (forward flow) will have a negative sign, (like in the example considered here), whereas flow coming **out of** the model (backflow) will be positive. A good way to remember that is that in the case of forward flow, the vector that gives you the direction of the flow and the normal to the face of the model point in opposite directions, and therefore their dot product will be negative.

On the other hand, in a situation of back flow, the numerical value in the *.flow file with be positive.

In this problem, since we are running a steady case, our physical time goes from 0.0 to 1.0 seconds, and the flow is constant with a value of 100.0 cc/sec.

**HINT**: it is very important that you are absolutely sure about the physical dimensions of your model: every unit (length, time, flow, density, etc.) in your analysis must be dimensionally consistent. You can easily check the size of your model in Paraview before importing it into **SimVascular**.

In this case, our cylinder has a radius $r=2.0$ cm and length $L=30$ cm.

First we need to create a file called **bct.dat** (its vtp format file is also created, called **bct.vtp**) that specifies the boundary conditions at the inlet defined by inflow.vtp

In the **SimVascular** GUI window, go to the *Inflow BC* subtab under *Simulations*. You will have to enter the following values in the various boxes/buttons of the GUI (see figure below):

Under

**Analytic Shape of Profile**, select the**parabolic**radio button. This options allows a parabolic velocity profile to be applied at the inlet. Other options are**plug**, which applies a constant velocity profile throughout the inlet face and Womerseley, that uses a closed form solution for the velocity profile of pulsatile flow in arteries.In the

**Mesh Face File (vtp)**box, click on the**Browse**button and select the**inflow.vtp**file in the*mesh-complete/mesh-surfaces/*folder.In the

**Flow Rate File**box, click on**Browse**and select the desired *.flow file. In this case, the file is**steady.flow**.

Under the **Parameters** menu, enter the following values:

**Period**: $1.0$ sec. For a steady flow problem like ours, the concept of*period*is somewhat vague. In this case, $1.0$ means the amount of physical time that we are going to run our simulation for.**viscosity**: $0.04$ Poise (gr/cm/s).

**WARNING**: Be very careful with the units! This is one the points where it is easy to make a mistake with the dimensions. For consistency, besides entering the right numerical value, try to also modify the units listed next to it (see figure below).

**density**: $1.06$ gr/cm$^3$. Same as before, be very careful with these units!

Under the **Output Parameters** menu, enter the following values:

**num of periods**: always enter 1 here. If you want to run your simulation for multiple cardiac cycles, it is possible to ask**svSolver**to loop over the information given by the**bct.dat**file for just one cycle (using the**solver.inp**file). By doing this, you will minimize the size of this file that can potentially be very large.**num pts in period**: 2. This is the number of*temporal*data points that you want to have in your bct.dat. This is not necessarily the number of time points in the *.flow file. In this case, they match (2 in both cases), but this is because this is a very simple example using steady flow and two time points is all we need to characterize a constant flow. In general, your *.flow file will have in the order of $20$ data points over the cardiac cycle (that’s about how many points you will be able to reconstruct using**PC-MRI**, for example), and your bct.dat will have on the order of $100$-$200$ points. Whatever is enough to have a smooth representation of the inflow wave mapping to velocity vectors at the inlet face.**num fourier modes**: 1. Fourier smoothing allows to smooth your inlet flow curve and to make sure that you have a periodic function in the specified interval.

**WARNING**: Be careful with this! **SimVascular** is doing a Fourier Series approximation of the data that you provide in the *.flow file. Since in this case, our data is constant flow, we only need one Fourier mode to capture this appropriately. For pulsatile flow problems, we will need more Fourier Modes to accurately represent the *.flow data (usually, $10$ Fourier modes is enough for a pulsatile problem).

After you are done entering all these parameters, click on the *CREATE 3-D FLOW SOLVER BC FILE* button to generate the **bct.dat** file. The format of this file is as follows:

```
np nl
x1 y1 z1 nl nn
vx1 vy1 vz1 t1
. . . .
. . . .
. . . .
vxnl vynl vznl tnl
. . . .
. . . .
. . . .
xnp ynp znp nl
vx1 vy1 vz1 t1
. . . .
. . . .
. . . .
```

The file defines the inflow boundary condition both spatially and in time. The spatial definition is obtained through $n_p$ point-wise velocity input blocks. In this case, $n_p = 102$, this is the total number of nodes on the **inflow.vtp** face. The temporal definition is given by $n_l$ input lines of the values at a certain position for $n_l$ time points, $t_1$ to $t_{n_l}$. In this case, $n_l = 2$ points (this is the value we entered in the *num pts in period* box.

Each block of data has, for each of the $n_p = 102$ spatial points, the following info:

The coordinates of the point: $x_1$ $x_2$ $x_3$ and the number of time points $n_l$.

The list of velocity vectors $v_x$ $v_y$ $v_z$ at time t.

A vtp file **bct.vtp** can be written using this option **Create Vtp** to graphically visualize the velocity distribution at the inlet surface with Paraview, as shown in the picture below.

To define the initial and boundary conditions of this problem, svPre need a script file (*.svpre) file. We go to the *Simulations->Create 3-D Solver Files* tab. In the “Create PreSolver script file” menu (see figure below), make sure that the right *.svpre file is loaded in the box (in this case, it should be cylinder.svpre . Click on the “**Load PreSolver scriptfile**” button. The following screen will appear:

The contents of the cylinder.svpre script file are:

```
# Read Mesh info
mesh_and_adjncy_vtu mesh-complete/mesh-complete.mesh.vtu
# Assign IDs to the surfaces
set_surface_id_vtp mesh-complete/mesh-complete.exterior.vtp 1
set_surface_id_vtp mesh-complete/mesh-surfaces/inflow.vtp 2
set_surface_id_vtp mesh-complete/mesh-surfaces/outlet.vtp 3
# Set Inlet BC
prescribed_velocities_vtp mesh-complete/mesh-surfaces/inflow.vtp
# Set BCT for Inlet
fluid_density 1.06
fluid_viscosity 0.04
bct_analytical_shape parabolic
bct_period 1.0
bct_point_number 2
bct_fourier_mode_number 1
bct_create mesh-complete/mesh-surfaces/inflow.vtp flow-files/steady.flow
bct_write_dat
bct_write_vtp
# Set Outlet BC
zero_pressure_vtp mesh-complete/mesh-surfaces/outlet.vtp
# Set Wall BC
noslip_vtp mesh-complete/walls_combined.vtp
# Write geometry and property data to geombc.dat.1
write_geombc
# Write initial values of velocity, pressure, etc to restart.0.1
write_restart
```

As we said before, each line of this *.svpre file represents a command that will be executed by **svPre**. This file needs to be edited to incorporate the right parameters/conditions for this problem. A complete list of svpre commands available is this section.Here is a description of each line.

The first line is used to define the topology of the finite element mesh. The first command **mesh_and_adjncy_vtu** used a vtk unstructured mesh file to define node coordinates, element connectivities and an adjacency relationship between elements.

```
mesh_and_adjncy_vtu mesh-complete/mesh-complete.mesh.vtu
```

The following command is used to assign an ID to all model surfaces.

```
set_surface_id_vtp cylinder.exterior.vtp 1
```

**HINT**: This line tags all the exterior faces in the model with an identifier (a Suface ID) , in this case, the number one. The assignment can tell flowsolver later which faces needed for wall stress calculation. We also need to introduce a new command if we want to activate the resistance boundary condition at the outlet face. We had previously determined that a resistance equal to $R = 1333.0$ dynes$\,$s/cm$^5$ needs to be applied at that outlet.

In order to do this, we need to assign a Surface ID that will help us later to identify the face and assign the correct resistance value. This is a trivial case, because we only have one single outflow face, and therefore one single resistance. But imagine one case where many more are needed. In this case, it is very important to meticulously label all the outlet faces with a meaningful name, and to make a sketch that helps you remember the list of Surface IDs that you considered in the model. Each surface ID will have a corresponding Resistance value (or impedance function, or set of RCR parameters etc.).

Going back to our *.svpre file, we need to add a line specifying the surface ID for the outlet face. We also number all other model surfaces in case we need to apply different boundary conditions through the **solver.inp** file.

```
set_surface_id_vtp mesh-complete/mesh-surfaces/inflow.vtp 2
set_surface_id_vtp mesh-complete/mesh-surfaces/outlet.vtp 3
```

**HINT**: Since a face can only have one ID, now ID 1 just represents cylnder wall because the inlet and outlet are assigned with ID 2 and 3.

The following line uses the existing **inflow.vtp** file to define a boundary subregion with applied velocities.

```
prescribed_velocities_vtp mesh-complete/mesh-surfaces/inflow.vtp
```

Note that we are just pointing to the right file (inflow.vtp) in the mesh-surfaces folder where we want some velocity vectors to be prescribed. These velocity vectors are given by the **bct.dat** file, already created from GUI as shown above. We must use the command prescribed*velocities*vtp in each surface where we prescribe some flow wave via a Dirichlet condition on the velocity vectors of that face. Instead, we can also create bct.dat (and bct.vtp) here using script commands as below. Similar to the GUI example above, these commands need to provide info about fluid density, fluid viscosiyt, velocity profile shape, period length, number of points in one period, number of fourier modes, inlet face, and flow file.

```
fluid_density 1.06
fluid_viscosity 0.04
bct_analytical_shape parabolic
bct_period 1.0
bct_point_number 2
bct_fourier_mode_number 1
bct_create mesh-complete/mesh-surfaces/inflow.vtp flow-files/steady.flow
bct_write_dat
bct_write_vtp
```

Like before, we are only pointing to the right path of the surface file where we want to prescribe the nonslip (rigid wall) condition. This is also a Dirichlet condition that makes all the velocity vectors of the nodes of the surface **wall** to be zero.

```
zero_pressure_vtp mesh-complete/mesh-surfaces/outlet.vtp
```

By using this condition, we are making the face **outlet** into a **weakly-pressure** face. This is mathematically expressed by the expressions we saw before:

$$ p = p_0 + R\,Q $$

$$ p(t)=p_0 + \frac{1}{T}\,\int_{t-T}^{t} Z(t−\tau) \, Q(\tau) \, d\tau $$

This expression sets $p_0 = 0$ for the face **outlet**. The total pressure set on that face will be the result of the flow-pressure operator considered (i.e., resistance, impedance, RCR, Coronary etc.).

The following line set no-slip boundary conditions for all walls. Since for this simple cylinder, there is only one wall and We can also use mesh-complete/mesh-surfaces/wall.vtp instead of mesh-complete/walls_combined.vtp.

```
noslip_vtp mesh-complete/walls_combined.vtp
```

The last two lines are the culmination of all the work you have been doing in **SimVascular** to this point!

```
write_geombc
write_restart
```

They generate two binary files (geombc.dat.1) and (restart.0.1) that are used as inputs for **svSolver** and are used to run the flow analysis.

**geombc.dat.1**contains the combination of geometry, material properties and boundary conditions specified in the problem.**restart.0.1**contains the set of initial conditions for our problem. We haven’t said anything about initial conditions till now. If you do not do something different,**SimVascular**will specify an almost-zero velocity initial condition for all the nodes of the mesh and a zero pressure**initial condition**. Here, the number**0**refers to**time step zero**, as it corresponds to the first file of the simulation.

Now, click on **Run PreSolver**. This command will actually launch the **svPre** application. A window will pop up, and you will see the list of commands of your *.svpre file being executed. After a few seconds (or minutes, depending on the size of the problem), the files **geombc.dat.1** and **restart.0.1** will be generated.

You can do the same if, instead of using the **SimVascular** GUI, you edit the *.svpre file like shown above, and then, from the command line, type:

```
%svpre cylinder.svpre
```

**HINT**: In both files, the number “.1” refers to the **partition number** of the file. Remember **svSolver** has the ability of running a problem *in parallel* (i.e., using multiple processors or computing cores), using MPI (message-passing interface). When we run a job using multiple processors, the first thing that happens is the “splitting” of these two files into as many processors we are going to use in our analysis, so the calculations can be performed faster. For example, if we use $4$ processors later in svSolver, these files will be split as follows:

```
geombc.dat.1 => geombc.dat.1 , geombc.dat.2 , geombc.dat.3 , geombc.dat.4
restart.0.1 => restart.0.1 , restart.0.2 , restart.0.3 , restart.0.4
```

Roughly speaking, each of the four files is $1⁄4$ of the size of the original un-split file. For a generic time step **n**, the solution will be given by the following files:

```
restart.n.1 , restart.n.2 , restart.n.3 , restart.n.4 , ...
```

We are almost done. There is only one thing left in the **svPresolver** part: to generate the numstart.dat dat. To do this, click on the **Create File** button under **Analysis Files**.

This file is really simple: it contains the scalar **0**. This number is used by the solver to identify the restart file that should be used as initial condition. In this case, since this file is **restart.0.1**, the file **numstart.dat** should contain a **0**. If for whatever reason, the initial file of our simulation was **restart.300.1**, the numstart.dat file should have a **300** entry. The value stored in this file gets updated as the simulation advances in time (we will see this later one).

At this point, we are almost ready to run the flow solver. Using **svPre**, we have generated the following files that are direct inputs to the solver:

**geombc.dat.1**: this file contains the combination of geometry and boundary conditions specified in the problem.**restart.0.1**: this file contains the set of initial conditions for our problem.**numstart.dat**: this file contains the scalar**0**. This number is used by the solver to identify the restart file that should be used as initial condition.**bct.dat**: this file contains the history of velocity vectors at the inflow face of the model according to a prescribed flow wave coming from a *.flow file.

**svSolver** is the flowsolver for SimVascular simulation. The **svSolver** program is called either using the **SimVascular** GUI or from the command line (in terminal). The input files to svSolver contain a complete description of model geometry,material properties, initial condition, boundary condition, and various parameter to control simulation. We will review this process briefly by keeping using the above cylinder example. First make sure the project folder is the example folder (…/cylinder)

Besides bct.dat, geombc.dat.1,restart.0.1 and numstart.dat, we are only missing one file in order to be able to run our analysis. This file is another input file for the solver that controls the actual flow of the numerical simulation, specifying parameters such as time step size, number of time steps, number of nonlinear iterations, boundary condition control, etc. This file needs to have the name **solver.inp** (input file for the solver), and we will characterize it in detail in the following section. A detailed description is also presented in this section.

These five files are the bare minimum we need to run an analysis. However, if we want to perform more complicated simulations, involving more complex boundary conditions, we will need more input files.

In addition to the five standard files (geombc.dat.1, restart.0.1, numstart.dat, bct.dat, solver.inp), we will need to provide impedance functions in the time domain for each impedance outlet, as well as a history of flow rates for each outlet. We will have therefore two additional ascii files: **impt.dat** (containing the impedance functions for each of the outlets), and **Qhistor.dat** (containing the flow rate history). A detailed description is here.

In addition to the five standard files (geombc.dat.1, restart.0.1, numstart.dat, bct.dat, solver.inp), we will need to provide the RCR parameters in a ascii file that will set the relationship between flow and pressure on each outflow face. This is done by defining a file named **rcrt.dat** containing such parameters. A detailed description is here.

In addition to the five standard files (geombc.dat.1, restart.0.1, numstart.dat, bct.dat, solver.inp), we will need to provide the coronary model parameters in a ascii file that will set the relationship between flow and pressure on each outflow face. This is done by defining a file named **cort.dat** containing such parameters. A detailed description is here.

This will required an executable that implements a lumped parameter network model for the patient circulation. This will be covered in a later version of this tutorial. Stay tuned!

We have completed the section on preprocessing the model. Let’s move on to **svSolver**, define the solver.inp file and run the analysis.

The main goal of this section is to define the file we are missing to run the analysis. This is the **solver.inp** file (i.e., input parameters for the solver). Most parameters are already assigned default values for cardiovascular simulation. Only a very small number of parameters must be set up in solver.inp. For this problem, the file we need will look like this:

```
# ================
# SOLUTION CONTROL
# ================
Number of Timesteps: 200
Time Step Size: 0.03
# ==============
# OUTPUT CONTROL
# ==============
Number of Timesteps between Restarts: 10
Number of Force Surfaces: 1
Surface ID's for Force Calculation: 1
# ===================
# MATERIAL PROPERTIES
# ===================
Viscosity: 0.04
Density: 1.06
# ==================================
# CARDIOVASCULAR MODELING PARAMETERS
# ==================================
Number of Coupled Surfaces: 1
Number of Resistance Surfaces: 1
List of Resistance Surfaces: 3
Resistance Values : 1333
# =============
# STEP SEQUENCE
# =============
Step Construction: 0 1 0 1
```

The file consists of a number of blocks, each block containing a number of lines that are instructions for the solver.

**WARNING**: it is very important that the wording of each line is exactly as presented here. Even a slight change (for instance, a change from uppercase to lowercase) will make the solver not understand the command!

The lines preceded by a **#** sign are comments and are ignored by the solver. Likewise, anything placed after a **#** on a given line is also ignored.

In this block, the different commands are:

**Number of Timesteps: 200** and **Time Step Size: 0.03** - These two lines control the amount of physical time that you run your problem for. In this case,

$$ \text{Total physical time} = \text{N. time steps} \times \text{Time Step Size} = T = N \times \Delta t = 200 \times 0.03 = 6.0\,\text{sec} $$

Note that this doesn’t match the **period** options we specified to generate the **bct.dat**. In this case, like we mentioned before, it does not really make sense to talk about a *cardiac cycle* (this is a steady flow), but if we wanted to run this analysis for *six* cardiac cycles, we would have to run the problem for $6.0$ seconds of physical time. If we kept our choice of time step size the same ( $\Delta t = 0.03$ sec), we will need a total number of time steps of $N = 200$.

**WARNING**: Note that this $N$ is the total number of time steps you need in your numerical simulation to model a certain physical time, given a prescribed $\Delta t$. This is not to be confused with the previous number of time steps you used to generate the bct.dat!

**WARNING**: Now the question is: how do you come up with a reasonable value for $\Delta t$? There is not a straightforward answer for this. $\Delta t$ is the parameter that controls your **temporal discretization**, which is something that works in a similar fashion to the **spatial discretization** given by your mesh: the finer, the more accurate the results, but also the bigger the size of the problem and the time to solve it! We don’t want to get into a lot of theoretical details now, so we will just provide you with a reasonable **recipe** to evaluate this parameter. The **recipe** to estimate a reasonable $\Delta t$ is based on a dimensionless parameter called the **CFL** number. The CFL number relates the velocities happening in the fluid domain ($v$), a temporal discretization parameter ($\Delta t$), and a mesh discretization parameter (i.e. mesh size) ($h$) as follows:

$$ \text{CFL} = \frac{v\,\Delta t}{h} $$

We want this **CFL** number to be around $1.0$. This will mean that, for the velocities present in our fluid domain, the temporal and spatial discretizations are *balanced*. In our problem, it can be shown that the average expected velocity is about $v = 16.7$ cm/s; the spatial discretization parameter or finite element size is $h = 0.5$. Therefore, if we shoot for a CFL number close to one, we have:

$$ \Delta t = \frac{h}{v} = \frac{0.5}{16.7} = 0.03 $$

Of course, you can imagine that in a real-world problem things are way more complicated to evaluate: it will be much harder to estimate where your model will have the largest velocities, what the mesh element size will be there, etc. The time step size $\Delta t$ is a parameter that will have a very important impact on the performance of the linear solver of equations. The smaller you make it, the *easier* you will be for the solver to find a solution, but the longer it will take you to reach a certain point in time.

This block contains the values for density and dynamic viscosity of blood: nothing really new here. Be careful though and make sure that you use the same units you have been using through the simulation process!

In this block, the meaning of the command is:

**Number of Timesteps between Restarts: 10** - This line tells the solver how often it should save solution files. In this problem, you are really calculating $200$ solutions to the problem at $200$ different time points, but in general you do not want to save a solution file for every single time step. Keep in mind that two consecutive solutions are only $\Delta t = 0.03$ seconds apart! In this line, we are asking the solver to save every other $20$ files. Therefore, the output files of the solver will look like this: restart.0.*, restart.10.*, restart.20.*, restart.30.*, …., restart.190.*, restart.200.*

**Number of Force Surfaces: 1** - This is the number of surfaces of the model where we are calculating the wall stress.

**Surface ID’s for Force Calculation: 1** - This line is list of surface ID’s considered for walls stress calculation. In our case, we only defined one surface ID (the number 1, assigned to the cylinder in svPre).

This is the block that controls the Boundary Conditions and the other features such as deformable wall parameters. The meaning of each command is:

**Number of Coupled Surfaces: 1**- This is the number of surfaces of the model where we are specifying a relation that couples Flow and Pressure. In order words, this number is the total number of**Resistance, Impedance, RCR and coronary surfaces**we have in our problem. In this case, since we only have one outlet with a resistance boundary condition, we enter a 1 in this line.**WARNING**: this line refers to coupled surfaces. A constant pressure outlet with no coupling between flow and pressure**does not**qualify as a coupled surface!**Number of Resistance Surfaces: 1**- This line sets the number of resistance surfaces in the model. In our case, we have one resistance surface.**List of Resistance Surfaces: 3**- This line the list of surface ID’s considered in the model for Boundary Condition specification. In our case, we only defined one surface ID (the number 3), at the outlet face of the model. It is very important that this number matches what you used in your *.svpre file. Otherwise, things will not work!**Resistance Values : 1333.0**- This line the list of resistancese considered in the outlets of the model. In our case, this resistance is 1333.0.

**WARNING**: Be very careful with the units! It is also very important that ordering of the resistance values in this line and the surface ID’s you provided in the previous line is consistent. This is a very common place to make a mistake. It is also very important that whatever you enter in these last two lines is consistent with want you entered in the *.svpre file.

Let us illustrate this with a more complex problem with 4 outlets (see figure below)

The *.svpre file should read something like this:

```
.
.
.
zero_pressure_vtp mesh-surfaces/outlet1.vtp
zero_pressure_vtp mesh-surfaces/outlet2.vtp
zero_pressure_vtp mesh-surfaces/outlet3.vtp
zero_pressure_vtp mesh-surfaces/outlet4.vtp
#
set_surface_id_vtp exterior_faces.vtp 1
set_surface_id_vtp mesh-surfaces/outlet1.vtp 2
set_surface_id_vtp mesh-surfaces/outlet2.vtp 3
set_surface_id_vtp mesh-surfaces/outlet3.vtp 4
set_surface_id_vtp mesh-surfaces/outlet4.vtp 5
.
.
.
```

And the solver.inp file:

```
.
.
.
Number of Resistance Surfaces: 4
List of Resistance Surfaces: 2 3 4 5
Resistance Values : 20000 10000 15000 21000
.
.
.
```

This line controls the non-linear iteration loop within the time step. The syntax of the line represents a two nonlinear iteration process for each time step. The **0** tells the solver to make a solve, the **1** to make an update of the solution. Since this sequence 0 1 is repeated, the two iterations are defined.

**WARNING**: Deciding on the adequate number of non-linear iterations for a problem is also a non-trivial problem. In principle, we need to iterate until the residual (i.e., the *error*) of our numerical solution is small enough. But doing many non-linear iterations on each time step is very costly. In general, for steady flow problems, 1 or 2 non-linear iterations are enough. For pulsatile problems, at least three non-linear iterations are needed. For deformable wall problems, 4 or more non-linear iterations are required. This parameter, together with the time step size $\Delta t$ and the quality of the spatial discretization given by the finite element mesh, will completely determine the performance of the linear solver of equations. The better chosen these parameters are, the faster and more accurately our simulation will run. We will talk more about this later.

The set of instructions explained here constitute a very small sample of all the possible instructions the **svSolver** can take via a solver.inp file. A more detailed discussion can be found in this section.

At this point we have generated all the files we need for this problem. Start simulation from the GUI.

```
Tab: Simulation -> Run Solver -> localhost
Select the project dir for Run Dir
Select Log Dir
Starting Step Number: 0
Click the button "whoami" to set username
Choose the number of processors - localhost num procs:4
Click "Run Simulation"
Wait a few seconds, Click the last "Start Trail" button to track the simulation progress
```

You cana also run simulation by a command lines.

```
%mpiexec -np 4 svsolver
```

This will launch a four-processor job in your computer. Therefore, the input file are split as follows:

```
geombc.dat.1 => geombc.dat.1, geombc.dat.2, geombc.dat.3, geombc.dat.4
restart.0.1 => restart.0.1, restart.0.2, restart.0.3, restart.0.4
```

At the same time, the solver copies all these files to a newly created folder called **4-procs_case/** under the project folder, and this is where all the output files of the analysis will be written to. In general, if you launch a **n** processor job, all the files will be copied to a **n-procs_case/** folder.

You can check the simulation progress in tab Console. It contains containing details that allows to evaluate how well your numerical simulation is doing. Here’s a brief description of what each of those columns means.

a | b | c | d | e | f | g | h |
---|---|---|---|---|---|---|---|

1 | 1.30E+001 | 1.16E-002 | (0) | 2.10E+002 | 2.62E+028 | < -10474 1 | 15 > | [199-190] |

1 | 2.50E+001 | 7.35E-003 | (-1) | 2.93E-001 | 5.15E+000 | < -3237 1 | 13> | [117-200] |

2 | 2.80E+001 | 5.13E-001 | (-16) | 1.75E-001 | 1.69E-001 | < -1357 1 | 5> | [63-1] |

2 | 3.00E+001 | 2.05E-002 | (-2) | 8.07E-002 | 2.67E-002 | < -3286 1 | 11> | [21-13] |

3 | 3.20E+001 | 1.20E-001 | (-10) | 8.75E-002 | 2.44E-002 | < -2342 1 | 7> | [36-1] |

3 | 3.40E+001 | 5.18E-003 | (-3) | 2.13E-002 | 3.59E-003 | < -3277 1 | 10> | [6-6] |

4 | 3.60E+001 | 2.14E-002 | (-2) | 5.57E-002 | 6.13E-003 | < -3146 1 | 9> | [24-2] |

4 | 3.80E+001 | 2.18E-003 | (-7) | 7.33E-003 | 3.15E-004 | < -3233 1 | 11> | [9-5] |

5 | 4.00E+001 | 1.52E-002 | (-1) | 4.22E-002 | 3.45E-004 | < -3141 1 | 10> | [27-3] |

5 | 4.20E+001 | 1.11E-003 | (-10) | 3.57E-003 | 1.24E-004 | < -3237 1 | 12> | [7-5] |

6 | 4.40E+001 | 1.24E-002 | (0) | 3.31E-002 | 3.16E-004 | < -3024 1 | 10> | [16-2] |

6 | 4.60E+001 | 1.01E-003 | (-10) | 3.98E-003 | 3.30E-005 | < -3220 1 | 11> | [6-5] |

7 | 4.80E+001 | 1.05E-002 | (0) | 2.70E-002 | 1.87E-004 | < -3019 1 | 10> | [8-2] |

7 | 5.00E+001 | 8.74E-004 | (-11) | 3.50E-003 | 4.35E-005 | < -2993 1 | 11> | [5-5] |

8 | 5.20E+001 | 9.09E-003 | (-1) | 2.25E-002 | 1.16E-004 | < -2993 1 | 10> | [8-2] |

8 | 5.30E+001 | 7.42E-004 | (-11) | 2.92E-003 | 5.63E-005 | < -2871 1 | 12> | [5-5] |

9 | 5.50E+001 | 7.95E-003 | (-1) | 1.89E-002 | 8.19E-005 | < -2876 1 | 10> | [7-2] |

9 | 5.70E+001 | 6.58E-004 | (-12) | 2.47E-003 | 7.07E-005 | < -2850 1 | 12> | [7-5] |

10 | 5.90E+001 | 6.99E-003 | (-2) | 1.65E-002 | 8.52E-005 | < -2871 1 | 10> | [11-3] |

10 | 6.10E+001 | 4.26E-004 | (-14) | 1.39E-003 | 3.55E-005 | < -6284 2 | 11> | [7-5] |

11 | 6.30E+001 | 6.17E-003 | (-2) | 1.42E-002 | 5.24E-005 | < -2845 1 | 10> | [6-3] |

11 | 6.50E+001 | 3.86E-004 | (-14) | 1.28E-003 | 3.41E-005 | < -6284 2 | 11> | [8-5] |

12 | 6.70E+001 | 5.45E-003 | (-3) | 1.23E-002 | 4.06E-005 | < -2728 1 | 10> | [7-3] |

12 | 6.90E+001 | 3.48E-004 | (-15) | 1.18E-003 | 3.95E-005 | < -6284 2 | 11> | [9-6] |

**a: Time step number**- We see that each time step number appears twice: this is because we are considering two non-linear iterations. This column will go from 1 to the total number of time steps in our problem (in this case, 100).**b: CPU time in seconds**- This is counted since you launched the analysis.**c: Measure of the nonlinear residual**- This gives you and idea of how accurate your solution is. Note that for each time step, the second entry is smaller than the first entry. This is a good sign that shows that the nonlinear iteration loop of the solver is doing a good job in improving the solution. You should always aim to a nonlinear residual at the end of the time step of at most $1\times10^{-3}$.**d: Logarithm of the residual change**- This provides you a very good way of quickly evaluating how well the solution is doing. If this number is very small and negative, then it is a good sign. An entry with the value “-10” means that you have reduced the residual by an order of magnitude from the beginning of the analysis, a -20, by**2 orders of magnitude**, and so on.**e: Entropy norm of the residual for the velocity field (max $\Delta u/ u$)****f: entropy norm of the residual for the pressure field (max $\Delta p/ p$)****g: > a – b | c>**a: block where the maximum residual happens (each block has 255 elements by default).

b: node in the block with the maximum error.

c: logarithmic measure of the ratio between the maximum residusl and the average residual: want this number to be as small as possible: it will be an indicator of the spatial uniformity of the residual.

**h: [a-b]**a: number of Krylov vectors used in the pre-conditioner solver.

b: number of Krylov vectors used in the GMRES solver.

**NOTE** - these numbers are zero when using the default **svLS**.

Once the analysis is done, you will see the collection of restart files corresponding to the different points in the time you decided to save.

Sometimes you want to restart the finished simulation.

```
Select .../cylinder/4-procs_case for Run Dir
Starting Step Number: 200 (the last time steop from previous simulation)
Choose the same number of processors as the previous
Click "Run Simulation"
```

Not the solver starts another 200-step simulation from time step 200.

As the simulation is completed, we are now ready to look at the restart files containing the solution. You will read these files in **svPost** to generate the visualization files. We explain that process in the following section.

In order to generate the visualization files (*.vip) and (*.vtu) files:

```
Tab: Simulations->Create VTU Files
Select a restart file in the running dir .../cylinder/4-procs_case for Input Files/Dir
start: the initial restart file number (0)
stop: the final restart file number (200)
increment: the increment between restart files (10).
Toggle on "Volume Mesh" and "Surface Mesh"
Click the button "Convert Files Only"
```

*.vtp or *.vtu files are generated in the project folder, containing the results for the whole finite element model using partitioned restart files as inputs.

Other options are also provided:

**Single File**, this options combines all the results at different time steps into a single *.vtp or *.vtu file.**Last Step to restart.x.0**, this options combines all the resart files of the last step (200) to a single restart file restart.200.0. which can be used to start a new simulation after renamed as restart.0.1.**Sim Units**, enter**cm**as we used**GCS**units throughout this tutorial.**Use Wall Options**, toggle on to apply more options for postprocessing.

**Wall File**, this limit the data output only on the specified wall.**Recalculate Wall Stresses**, this enables postprocessing to calculate the wall stress again.**Apply Wall Deformation**, this enables the update wall coordinates for deformable wall.**Viscosity**, this sets a new viscosity value for wall stress calculation.

If Single File option is on when postprocessing, all_results.vtp or all_results.vtu will be produced. In this case, svPost can calculate pressure and flow rate for outlets. Just click “Calculate Flows Only”.

To visualize the time dependent results we use **ParaView**.

Launch

**Paraview**.Select

*File->Open…*. The*Open File*dialog should appear. Go to the project folder , select the entry**all_results..vtu**and click**OK**.Under

*Properties*click**Apply**. The solid model will show up on the screen.At this point, you can interact with the model by rotating it using the rotation or translation buttons. Use the

**Surface with Edges**option to visualize the finite element mesh.

- First you should increase your current result time from
**0**to**20**(the last available time step).

- You should now see the available result quantities for your model, i.e., cellsNormals, GlobalElementID, GlobalNodeID, pressure, timeDeriv, traction, velocity, WSS.

Go to

*Properties-Coloring*click on**Show**and**Rescale**. The default color map of the model is showing the pressure results in dynes/cm$^{2}$. Let’s transform the pressure scale to mmHg (1 mmHg = 1333.2 dyn/cm$^{2}$). To do this, use the**calculator**filter.Add a calculator filter, choose pressure in the

*scalars*drop-down list and divide*pressure*by the conversion factor 1333.2.Enter a meaningful name in the “Results Array Name” box (for example, Pressure in mmHg)

You should now see the following contour plot.

The first example shown above is a cylinder with parabolic steady flow at the inlet and resistance at the outlet. We continue to use the cylinder in the following examples to explore plug flow, RCR, deformable wall, variable wall properties.

**Example 2** - the cylinder with plug steady flow at the inlet and RCR at the outlet.

**Example 3** - the cylinder with plug steady flow at the inlet, RCR at the outlet and deformable wall with uniform wall properties.

**Example 4** - the cylinder with plug steady flow at the inlet, RCR at the outlet and deformable wall with variable wall properties.

All the examples are included in here. Example 1 *cylinder* includes all the files from presolver to postprocessing; checkout readme.txt in the *cylinder* for more details about the files. Example 2-3 only include the necessary files to start simulation from presolver.

This example shows how to simulate a cylinder with plug steady flow at the inlet and RCR at the outlet. The script file cylinder.svpre used is as below:

```
# Read Mesh info
mesh_and_adjncy_vtu mesh-complete/mesh-complete.mesh.vtu
# Assign IDs to the surfaces
set_surface_id_vtp mesh-complete/mesh-complete.exterior.vtp 1
set_surface_id_vtp mesh-complete/mesh-surfaces/inflow.vtp 2
set_surface_id_vtp mesh-complete/mesh-surfaces/outlet.vtp 3
# Set Inlet BC
prescribed_velocities_vtp mesh-complete/mesh-surfaces/inflow.vtp
# Set BCT for Inlet
fluid_density 1.06
fluid_viscosity 0.04
bct_analytical_shape plug
bct_period 1.0
bct_point_number 2
bct_fourier_mode_number 1
bct_create mesh-complete/mesh-surfaces/inflow.vtp flow-files/steady.flow
bct_write_dat
bct_write_vtp
# Set Outlet BC
zero_pressure_vtp mesh-complete/mesh-surfaces/outlet.vtp
# Set Wall BC
noslip_vtp mesh-complete/walls_combined.vtp
# Write geometry and property data to geombc.dat.1
write_geombc
# Write initial values of velocity, pressure, etc to restart.0.1
write_restart
```

The only difference from Example 1 in the presolver script is using plug flow:

bct_analytical_shape plug

When the inlet of a vessel is close to the heart, the velocity profile at the inlet has a plug shape, instead of a parabolic shape. The solver.inp used for this example is as below:

```
# ================
# SOLUTION CONTROL
# ================
Number of Timesteps: 500
Time Step Size: 0.001
# ==============
# OUTPUT CONTROL
# ==============
Number of Timesteps between Restarts: 20
Number of Force Surfaces: 1
Surface ID's for Force Calculation: 1
# ===================
# MATERIAL PROPERTIES
# ===================
Viscosity: 0.04
Density: 1.06
# ==================================
# CARDIOVASCULAR MODELING PARAMETERS
# ==================================
Number of Coupled Surfaces: 1
Number of RCR Surfaces: 1
List of RCR Surfaces: 3
RCR Values From File: True
# =============
# STEP SEQUENCE
# =============
Step Construction: 0 1 0 1
```

There are a few differences compared to Example 1 in solver.inp:

Number of RCR Surfaces: 1

List of RCR Surfaces: 3

RCR Values From File: True

The outlet boundary condition is RCR type, which needs a file rcrt.dat to define for flowsolver. The file format is described in here.

*rcrt.dat*

```
2
2
121
0.000015
1212
0.0 0.0
1.0 0.0
```

In the rcrt.dat:

Line 1: the maximum number (2) of time points for all the outlets

Line 2: the number (2) of time points for the first outlet (this example only has one outlet)

Line 3: the proximal resistance (121) for the first outlet

Line 4: the capacitance (0.000015) for the first outlet

Line 5: the distal resistance (1212) for the first outlet

Line 6: the distal pressure (0.0) at time (0.0) for the first outlet

Line 7: the distal pressure (0.0) at time (1.0) for the first outlet

Number of Timesteps: 500

Time Step Size: 0.001

In this example, a smaller time step size is used, because the effect of the time constant in RCR type BC needs to be considered. For more accurate numerical computation, the time step size should be much smaller than the time constant (=R*C=(121+1212)*0.000015=0.02). With the small step size, the number of time steps is also increased.

This example shows how to simulate a cylinder with plug steady flow at the inlet, RCR at the outlet and deformable wall with uniform wall properties. The most robust way to simulate a deformable wall case is to first simulate a rigid wall case with the same boundary conditions. For this example, Example 2 is the rigid case we need to run. After the simulation of Example 2 is complete, use svPost to reduce the restart files for the last step to a single file, which would be restart.500.0. Copy this file to the project folder of Example 3 and rename it as restart.0.1.

The script file cylinder.svpre used is as below:

```
# Read Mesh info
mesh_and_adjncy_vtu mesh-complete/mesh-complete.mesh.vtu
# Assign IDs to the surfaces
set_surface_id_vtp mesh-complete/mesh-complete.exterior.vtp 1
set_surface_id_vtp mesh-complete/mesh-surfaces/inflow.vtp 2
set_surface_id_vtp mesh-complete/mesh-surfaces/outlet.vtp 3
# Set Inlet BC
prescribed_velocities_vtp mesh-complete/mesh-surfaces/inflow.vtp
# Set BCT for Inlet
fluid_density 1.06
fluid_viscosity 0.04
bct_analytical_shape plug
bct_period 1.0
bct_point_number 2
bct_fourier_mode_number 1
bct_create mesh-complete/mesh-surfaces/inflow.vtp flow-files/steady.flow
bct_write_dat
bct_write_vtp
# Set Outlet BC
zero_pressure_vtp mesh-complete/mesh-surfaces/outlet.vtp
# Set Wall BC
deformable_wall_vtp mesh-complete/walls_combined.vtp
deformable_thickness 0.2
deformable_E 4000000.0
deformable_nu 0.5
deformable_kcons 0.833333
deformable_pressure 133300
deformable_solve_displacements
wall_displacements_write_vtp
# Write Geometry and Property data to geombc.dat.1
write_geombc
# Append displacement to restart.0.1
append_displacements
```

The difference from Example 2 in the presolver script is to indicate which surfaces are deformable, set up uniform wall properties, provide an initial pressure (which is estimated form the simulation result of the rigid case Example 2), solve initial displacement, write the initial displacment to a vtp for review the solution, and finaly append it to restar.0.1 we just copied and renamed.

deformable_wall_vtp mesh-complete/walls_combined.vtp

deformable_thickness 0.2

deformable_E 4000000.0

deformable_nu 0.5

deformable_kcons 0.833333

deformable_pressure 133300

deformable_solve_displacements

wall_displacements_write_vtp

append_displacements

After running the script, displacment.vtp is create in the project folder, which shows the initial displacement:

The solver.inp used for this example is as below:

```
# ================
# SOLUTION CONTROL
# ================
Number of Timesteps: 500
Time Step Size: 0.001
# ==============
# OUTPUT CONTROL
# ==============
Number of Timesteps between Restarts: 20
Number of Force Surfaces: 1
Surface ID's for Force Calculation: 1
# ===================
# MATERIAL PROPERTIES
# ===================
Viscosity: 0.04
Density: 1.06
# ==========================
# DEFORMABLE WALL PROPERTIES
# ==========================
Deformable Wall: True
Thickness of Vessel Wall: 0.2
Young Mod of Vessel Wall: 4000000.0
Density of Vessel Wall: 1.0
Poisson Ratio of Vessel Wall: 0.5
Shear Constant of Vessel Wall: 0.833333
# ==================================
# CARDIOVASCULAR MODELING PARAMETERS
# ==================================
Number of Coupled Surfaces: 1
Number of RCR Surfaces: 1
List of RCR Surfaces: 3
RCR Values From File: True
# =============
# STEP SEQUENCE
# =============
Step Construction: 0 1 0 1 0 1 0 1
# =============
# LINEAR SOLVER
# =============
svLS Type: GMRES
```

There are a few differences compared to Example 2 in solver.inp:

Deformable Wall: True

Thickness of Vessel Wall: 0.2

Young Mod of Vessel Wall: 4000000.0

Density of Vessel Wall: 1.0

Poisson Ratio of Vessel Wall: 0.5

Shear Constant of Vessel Wall: 0.833333

The solver.inp indicates the wall is deformable and assign uniform wall properties for flowsolver simulation.

Step Construction: 0 1 0 1 0 1 0 1

For deformable cases, the flowsolver needs more iterations of step sequence to derive accurate solution.

svLS Type: GMRES

For deformable cases, linear solver has better performance with GMRES.

This example shows how to simulate a cylinder with plug steady flow at the inlet, RCR at the outlet and deformable wall with variable wall properties. Similar to Example 3, we use the simulation result of Example 2 as the initial point.

The script file cylinder.svpre used is as below:

```
# Read Mesh info
mesh_and_adjncy_vtu mesh-complete/mesh-complete.mesh.vtu
# Assign IDs to the surfaces
set_surface_id_vtp mesh-complete/mesh-complete.exterior.vtp 1
set_surface_id_vtp mesh-complete/mesh-surfaces/inflow.vtp 2
set_surface_id_vtp mesh-complete/mesh-surfaces/outlet.vtp 3
# Set Inlet BC
prescribed_velocities_vtp mesh-complete/mesh-surfaces/inflow.vtp
fluid_density 1.06
fluid_viscosity 0.04
bct_analytical_shape plug
bct_period 1.0
bct_point_number 2
bct_fourier_mode_number 1
bct_create mesh-complete/mesh-surfaces/inflow.vtp flow-files/steady.flow
bct_write_dat
bct_write_vtp
# Set Outlet BC
zero_pressure_vtp mesh-complete/mesh-surfaces/outlet.vtp
# Set Wall BC
deformable_wall_vtp mesh-complete/walls_combined.vtp
set_surface_thickness_vtp mesh-complete/mesh-surfaces/inflow.vtp 0.2
set_surface_thickness_vtp mesh-complete/mesh-surfaces/outlet.vtp 0.1
solve_varwall_thickness
set_surface_E_vtp mesh-complete/walls_combined.vtp 4000000
solve_varwall_E
varwallprop_write_vtp
deformable_nu 0.5
deformable_kcons 0.833333
deformable_pressure 133300
deformable_solve_varwall_displacements
wall_displacements_write_vtp
# Write Geometry and Property data to geombc.dat.1
write_geombc
# Append displacement to restart.0.1
append_displacements
```

The difference from Example 3 in the presolver script is to solve variable thickness and Young’s modulus, and assign them to the wall, instead of giving uniform thickness and Young’s modulus.

set_surface_thickness_vtp mesh-complete/mesh-surfaces/inflow.vtp 0.2

set_surface_thickness_vtp mesh-complete/mesh-surfaces/outlet.vtp 0.1

solve_varwall_thickness

set_surface_E_vtp mesh-complete/walls_combined.vtp 4000000

solve_varwall_E

varwallprop_write_vtp

After running the script, varwallprop.vtp and displacment.vtp are create in the project folder, which shows the thickness and Young’s modulus, and initial displacement, respectively:

The solver.inp used for this example is as below:

```
# ================
# SOLUTION CONTROL
# ================
Number of Timesteps: 1000
Time Step Size: 0.001
# ==============
# OUTPUT CONTROL
# ==============
Number of Timesteps between Restarts: 20
Number of Force Surfaces: 1
Surface ID's for Force Calculation: 1
# ===================
# MATERIAL PROPERTIES
# ===================
Viscosity: 0.04
Density: 1.06
# ==========================
# DEFORMABLE WALL PROPERTIES
# ==========================
Deformable Wall: True
Variable Wall Thickness and Young Mod: True
Density of Vessel Wall: 1.0
Poisson Ratio of Vessel Wall: 0.5
Shear Constant of Vessel Wall: 0.833333
# ==================================
# CARDIOVASCULAR MODELING PARAMETERS
# ==================================
Number of Coupled Surfaces: 1
Number of RCR Surfaces: 1
List of RCR Surfaces: 3
RCR Values From File: True
# =============
# STEP SEQUENCE
# =============
Step Construction : 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1
# =============
# LINEAR SOLVER
# =============
svLS Type: GMRES
```

There are a few differences compared to Example 3 in solver.inp:

Variable Wall Thickness and Young Mod: True

Step Construction : 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1

The solver.inp indicates the wall properties are variable, and uses the derived solution form svPre, instead of assigning uniform thickness and Young’s modulus. For variable wall cases, the flowsolver need much more iterations of step sequence to get accurate solution.

In this section, more details are provided for presolver commands, solver.inp and other files for running various simulations.

This section lists the available **svPre** commands, the associated parameters and what they do.

svPre Command | Argument Format | Description |
---|---|---|

mesh_and_adjncy_vtu | (file name) | Read node coordinates,element connectivities and adjacencies from the give vtu file |

svPre Command | Argument Format | Description |
---|---|---|

set_surface_id_vtp | (file name) (integer) | Set a ID for the element faces provided by the given vtp file. Mostly used to tag exterior surfaces for Boundary Condition attributes, and also to compute tractions at the boundary |

svPre Command | Argument Format | Description |
---|---|---|

noslip_vtp | (file name) | Noslip condition will be applied on the nodes provided by the given vtp file |

prescribed_velocities_vtp | (file name) | Same as “noslip_vtp” |

zero_pressure_vtp | (file name) | Zero pressure boundary condition will be appied on the element faces provided by the given vtp file |

pressure_vtp | (file name) (double) | Prescribed pressure boundary condition will be appied on the element faces provided by the given vtp file |

svPre Command | Argument Format | Description |
---|---|---|

fluid_density | (double) | Set the fluid density to the given value |

fluid_viscosity | (double) | Set the fluid viscosity to the given value |

bct_analytical_shape | (parabolic, plug, womersley) | Indicate the velocity profile for the inlet |

bct_period | (double) | Set the period to the given value |

bct_point_number | (integer) | Set the number of points in one period |

bct_fourier_mode_number | (integer) | Set the mode number for Fourier smoothing of inlet flow |

bct_flip | Flip the flow direction at the inlet | |

bct_create | (vtp file name) (flow file name) | Calculate the velocity for the specified inlet with the given flow file |

bct_merge_on | Merge all bct files to one single file if there are multiple inlets. Use before write_bct_dat and write_bct_vtp | |

bct_write_dat | Write calculated BCT data to bct.dat | |

bct_write_vtp | Write calculated BCT data to bct.vtp |

svPre Command | Argument Format | Description |
---|---|---|

initial_pressure | (double) | Set the initial pressure as the given value in the model if not read from other files. The default is $p_0$ = $0$ |

initial_velocity | (double) (double) (double) | Set the initial velocity as the given values in the model if not read from other files. The default is $v_0 = 0.001,\,0.001,\,0.001$ |

svPre Command | Argument Format | Description |
---|---|---|

write_restart | (file name) | Write the specified file (restart.0.1 if a file name not provided) for svSolver. It contains initial conditions:pressure, velocity (and displacement, time derivative of solution if applicable). |

write_geombc | (file name) | Write the specified file (geombc.dat.1 if a file name not provided) for svSolver. It contains the info of geometry, boundary conditions and material properties. |

write_numstart | (integer) | Write numstart.dat with the specified time step number (0 if a number not provided) for svSolver. |

svPre Command | Argument Format | Description |
---|---|---|

read_restart_solution | (restart file name) | Read a previously computed restart solution (velocity and pressure fields). This could be used in deformable wall simulations |

read_restart_displacements | (restart file name) | Read a previously computed restart solution (displacement field). This is used in deformable wall simulations |

read_restart_accelerations | (restart file name) | Read a previously computed restart solution (acceleration field) |

read_restart_varwallprop | (restart file name) | Read a previously computed restart solution (variable wall properties) |

read_geombc_varwallprop | (geombc file name) | Read a previously computed geombc file (variable wall properties) |

svPre Command | Argument Format | Description |
---|---|---|

read_all_variables_vtu | (file name) | Read the values for all variables: pressure, velocity (and displacement, time derivative of solution, wall property for variable wall, if available) from the given vtu file. |

read_pressure_velocity_vtu | (file name) | Read the values for pressure and velocity from the given vtu file. |

read_pressure_vtu | (file name) (variable name in vtu) | Read the values for pressure from the given vtu file using an optional variable name. |

read_velocity_vtu | (file name) (variable name in vtu) | Read the values for velocity from the given vtu file using an optional variable name. |

read_displacements_vtu | (file name) (variable name in vtu) | Read the values for displacements from the given vtu file using an optional variable name. |

read_accelerations_vtu | (file name) (variable name in vtu) | Read the values for time derivative of solution from the given vtu file using an optional variable name. |

read_varwallprop_vtu | (file name) (variable name in vtu) | Read the values for variable wall properties from the given vtu file using an optional variable name. |

svPre Command | Argument Format | Description |
---|---|---|

deformable_wall_vtp_simple | (file name) | Set the face specified by the given vtp file as deformable wall |

fix_free_edge_nodes_vtp | (file name) | Fix the free edges of the face specified by the given vtp file with a zero-displacement condition |

deformable_create_mesh_vtp | (file name) | Generate a new finite element mesh with the nodes and elements on the face specified by the given vtp file (filename), with the purpose of solving the solid mechanics finite element problem of the pressurization of that surface with the fluid’s internal pressure |

deformable_wall_vtp | (file name) | For convenience, the command is the combination of the above three commands: deformable_wall_vtp_simple, fix_free_edge_nodes_vtp and deformable_create_mesh_vtp. |

deformable_thickness | (double) | Set the thickness of the vessel wall to the given value, assumed uniform. |

deformable_E | (double) | Set the elastic modulus of the vessel wall to the given value, assumed uniform |

deformable_nu | (double) | Set the Poisson’s ratio of the vessel wall to the given value, assumed uniform. |

deformable_kcons | (double) | Set the Shear correction factor of the vessel wall to the given value, assumed uniform. |

deformable_pressure | (double) | Pressure used to load the vessel wall structure. This pressure should be representative of pressure of the internal fluid. The result of this pressurization is a “loaded” vessel wall in equilibrium with the internal blood |

deformable_direct_solve_displacements | (none) | Use the direct, sparse solver for the initial displacements of the CMM-FSI method (deformable wall model) |

deformable_solve_displacements | (none) | Use the iterative solver for the initial displacements of the CMM-FSI method (deformable wall model) |

wall_displacements_write_vtp | (file name) | Write the displacment to the specified vtp file (displacement.vtp if a file name not provided) for review. |

append_displacements | (file name) | Append the displacement field to the specified existing file (restart.0.1 if a file name not provided). This command does not need a posterior “write_restart” command. |

svPre Command | Argument Format | Description |
---|---|---|

set_surface_thickness_vtp | (file name) (double) | Set the thickness of the face specified by the given vtp file (file name) to a prescribed value. This can be used to assign the thickness to various outlet surfaces and then interpolate the wall thickness from these values using the Laplace_Thickness command. |

solve_varwall_thickness | (none) | Solve the Laplace problem and determines a variable thinkness distribution in the wall. |

set_surface_E_vtp | (file name) (double) | Set the elastic modulus of the face specified by the given vtp file(file name) to a prescribed value. This can be used to assign the elastic modulus to various outlet surfaces and then interpolate the wall elastic modulus from these values using the Laplace_Evw command. |

solve_varwall_E | (none) | Solve the Laplace problem and determines a variable elastic modulus distribution in the wall. |

set_surface_initial_E_vtp | (file name) (double) | Set the initial elastic modulus of the face specified by the given vtp file (file name) to a prescribed value for posterior transient Laplace solving (command Transient_Laplace_Evw). |

solve_transient_varwall_E | (none) | Solve the Transient Laplace problem and determines a variable elastic modulus distribution in the wall. |

varwallprop_write_vtp | (file name) | Write the variable thickness and elastic modulus to the specified vtp file |

deformable_solve_varwall_displacements | (none) | Use the iterative solver for the initial displacements of the CMM-FSI method (deformable wall model) with variable wall properties |

append_varwallprop | (file name) | Append the variable wall thickness and elastic modulus to the specified file (geombc.dat.1 if a file name not provided) |

This section discusses the options available in the **solver.inp** file.

Command | Default | Possible Values | Description |
---|---|---|---|

Default Input File | File name with relative or absolute path | Most parameters are already assigned default values by hard-coding for cardiovascular simulation as shown in the following tables. Only a very small number of parameters must be set up in solver.inp. If the user needs different default values for a few parameters, the new values can also be assigned for them in solver.inp. But if the user needs different default values for many parameters, a default input file can be created and the new default values are put in this file. Users can use the file default.inp under SimVascular home(installation) folder as a reference to create solver.inp or a custom default input file. |

Command | Default | Possible Values | Description |
---|---|---|---|

Time Varying Boundary Conditions From File | (True) | True,False | If the bct.dat file was created containg prescribed velocity at the inlet (this will be the case for most simulations), this option should be set to True. |

BCT File Type | (DAT) | DAT,VTP | This entry tells the solver to read inflwo boundary conditions from bct.dat or bct.vtp. |

Number of BCT Files | (1) | (integer) | This entry tells the solver how many bct files need to be read. If there is only one, name it as bct.dat (or bct.vtp); if there are multiple, name them as bct1.dat, bct2.dat…(or bct1.vtp, bct2.vtp…) |

BCT Matching Type | (Global Node ID) | Global Node ID,Coordinates | This entry tells the solver to match bct nodes by global node id or coordinates. |

BCT Time Scale Factor | (1.0) | (double) | Defines an amplification factor for the velocity data contained in the bct.dat file. It allows scaling the time history given in the bct.dat file by the factor specified in this line. For example, if your original bct.dat has a period of $0.8$ seconds, and if you wanted to simulate a problem with the same inflow wave shape, but with a period of $0.4$ seconds, you would have to enter a BCT Time Scale Factor of 0.5 in this line. For most cases, it should be 1.0. |

Command | Default | Possible Values | Description |
---|---|---|---|

Equation of State | (Incompressible) | Incompressible, Compressible | This entry tells the solver to solve the INCOMPRESSIBLE or COMPRESSIBLE Navier-Stokes equations. |

Viscous Control | (Viscous) | Viscous | This entry activates the viscous terms in the Navier-Stokes equations. A value different from the default will neglect the contribution of the viscosity and solver the inviscid Navier-Stokes equations. |

Number of Timesteps | NO DEFAULT | (integer) | Total number of timesteps in the simulation |

Time Step Size | NO DEFAULT | (double) | Time size of each step |

Command | Default | Possible Values | Description |
---|---|---|---|

Viscosity | NO DEFAULT | (double) | Fluid viscosity |

Density | NO DEFAULT | (double) | Fluid density |

Command | Default | Possible Values | Description |
---|---|---|---|

Body Force Option | (None) | (Vector,User e3source.f) | The vector option applies a constant body force. A user-defined body force can also be specified through the e3source.f Fortran subroutine |

Body Force | (0.0 0.0 0.0) | (double double double) | The vector valued constant body force |

Command | Default | Possible Values | Description |
---|---|---|---|

Number of timesteps between restarts | NO DEFAULT | (integer) | Number of time steps between a new restart. |

Print Average Solution | (True) | False,True | Print time-averaged pressure and velocities in restart files. These values can be used by the svAdapt application to refine the computational mesh increasing the accuracy of the resulting pressure/velocity distribution. |

Print Error Indicators | (False) | False,True | Print time-accumulated errors of pressure and velocities to their averaged values in restart files. |

Number of Surfaces which Output Pressure and Flow | (0) | (integer) | Total number of surfaces where average pressure and total flow rate need to be integrated. The results will be written in the simulation folder to the files PressHist.dat and FlowHist.dat |

List of Output Surfaces | NO DEFAULT | (space separated integer list) | List of ID for surfaces where average pressure and total flow rate need to be integrated |

Number of Force Surfaces | (0) | (integer) | This line specifies the number of surfaces where we are going to collect tractions (i.e., wall shear stress) in our model. Tractions are collected as a ‘post-processing’ step after the velocity and pressure fields are obtained. |

Surface ID’s for Force Calculation | NO DEFAULT | (space separated integer list) | List of ID for surfaces tagged for Force Calculation (i.e., shear stress calculation). |

Force Calculation Method | (Velocity Based) | Velocity Based,Residual Based | Specify a method for surface force calculation (i.e., traction, shear stress). |

Apply Wall Deformation | (False) | True,False | Decide whether to udpate wall coordinates for deformable wall during wall stress calculation |

Command | Default | Possible Values | Description |
---|---|---|---|

Solver Type | (svLS) | svLS | sv linear solver selected by default |

svLS Type | (NS) | CG,GMRES,NS | Type selected for svLS. For deformable wall cases, GMRES is recommended to improve computing performance. |

Maximum Number of Iterations for svLS NS Solver | (1) | (integer) | Maximum number of iterations for the full Navier-Stokes solver |

Maximum Number of Iterations for svLS Momentum Loop | (2) | (integer) | Maximum number of iterations for the Momentum equation solver |

Maximum Number of Iterations for svLS Continuity Loop | (400) | (integer) | Maximum number of iterations for the continuity equation solver |

Number of Solves per Left-hand-side Formation | (1) | (integer) | It tells the solver how often the left-hand-side of the system of equations should be updated in the non-linear iteration loop. 1 is a good value, we recommend you to use it. |

**WARNING:** For simulations of deformable vessels these defaults may need to be changed to 10,20,400, respectvely.

Command | Default | Possible Values | Description |
---|---|---|---|

Time Integration Rule | (Second Order) | First Order,Second Order | This option allows the user to select the parameters in the generalized $\alpha$ time integration method for systems of first order ordinary differential equations in time. Setting a first order scheme is equivalent to adopt a backward (or implicit) Euler scheme. If the user selects a second order scheme then the value of the generalized alpha parameters are determined as $\alpha_m = \frac{1}{2}\frac{3-\rho_{\infty}}{1+\rho_{\infty}}$, $\alpha_f = \frac{1}{1+\rho_{\infty}}$, and $\gamma = \frac{1}{2} + \alpha_m - \alpha_f$ |

Time Integration Rho Infinity | (0.5) | (double in [0,1]) | Value of $\rho_{\infty}$ for the generalized $\alpha$ method. It sets the parameter that controls the amount of user-defined numerical dissipation. This parameter takes values in the range (0,1). The value 0 corresponds to maximal numerical dissipation, where the under-resolved frequencies are annihilated within one time step. The value 1 implies no numerical dissipation: all the frequencies present in the simulation are maintained through the time steps. These can be dangerous if the temporal and spatial discretizations are not adequate. |

Flow Advection Form | (Convective) | Convective,Conservative | It sets the solver to use the advective or conservative form of the Navier-Stokes equations. |

Quadrature Rule on Interior | (2) | 1,2,3,4 | This option sets the integration order for elements that do not contain boundary faces. Orders 1,2,3,4 correspond to 1,4,16,29 integration points, respectively. |

Quadrature Rule on Boundary | (3) | 1,2,3,4 | This option sets the integration order for boundary elements. Orders 1,2,3,4 correspond to 1,3,6,12 integration points, respectively. |

Number of Elements Per Block | (255) | (integer) | To improve the efficiency, the elements belonging to each processor are stored in separate groups. The size of these groups can be set by the user through this option |

Command | Default | Possible Values | Description |
---|---|---|---|

Number of Coupled Surfaces | (0) | (integer) | Total number of boundary surfaces with flow-pressure coupling |

Pressure Coupling | (Implicit) | None,Explicit,Implicit,P-Implicit | This line sets the implicit implementation of the Coupled-Multidomain Method in the code. |

Number of Resistance Surfaces | (0) | (integer) | Total number of surfaces with assigned resistance boundary condition |

List of Resistance Surfaces | NO DEFAULT | (space separated integer list) | List of IDs (same ID defined in the model.svpre file) of surfaces with resistance boundary condition applied |

Resistance Values | NO DEFAULT | (space separated double list) | Resistance values for the previously specified list of faces IDs |

Number of Impedance Surfaces | (0) | (integer) | Total number of faces with impedance boundary condition |

List of Impedance Surfaces | NO DEFAULT | (space-separated integer list) | List of IDs (same ID defined in the model.svpre file) of surfaces with impedance boundary condition applied |

Impedance From File | (False) | False,True | Must be set to True to read the impedance information from file |

Number of RCR Surfaces | (0) | (integer) | Total number of faces with boundary RCR blocks |

List of RCR Surfaces | NO DEFAULT | (space-separated integer list) | List of IDs (same ID defined in the model.svpre file) of surfaces with boundary RCR block applied |

RCR Values From File | (False) | False,True | Must be set to True to read the impedance information from file |

Number of COR Surfaces | (0) | (integer) | Total number of faces with coronary boundary conditions |

List of COR Surfaces | NO DEFAULT | (space-separated integer list) | List of IDs (same ID defined in the model.svpre file) of surfaces with applied coronary boundary condition |

COR Values From File | (False) | False,True | Must be set to True to read the impedance information from file |

Command | Default | Possible Values | Description |
---|---|---|---|

Backflow Stabilization Coefficient | (0.2) | (double in [0,1]) | Backflow stabilization coefficient. For the definition of these coefficient, see the [following publications](docsRefs.html#refSec2) |

Number of Surfaces which zero out in-plane tractions | (0) | (integer) | Number of surfaces with prescribed zero in-plane tractions. |

List of Surfaces which zero out in-plane tractions | NO DEFAULT | (space-separated integer list) | List of surfaces which zero out in-plane tractions. This list should adopt the same numbering strategy defined in the svPre input file |

Lagrange Multipliers | (False) | False,True | Activate/deactivate Lagrange multipliers |

Number of Constrained Surfaces | (0) | (integer) | Total number of face with applied Lagrange multiplier constraints |

List of Constrained Surfaces | NO DEFAULT | (space-separated integer list) | List of surfaces with applied Lagrange multipliers. This list should adopt the same numbering strategy defined in the svPre input file |

Constrained Surface Information From File | (False) | False,True | This flag needs to be set to true if a file is available with the information on how to set Lagrange multipliers at selected outlets |

Command | Default | Possible Values | Description |
---|---|---|---|

Find the GenBC Inside the Running Directory | (True) | False,True | Look for a GenBC executable inside the current simulation folder. This executable implements a 0D lumped circulation model providing the boundary conditions to the 3D Finite Element solver |

Number of Timesteps for GenBC Initialization | (0) | (integer) | This defines the number of steps to be performed before activation the the closed-loop boundary conditions. For timesteps smaller than this value, the GenBC application will provide fixed boundary pressures at each outlet equal to the initial values provided by the user. |

Number of Dirichlet Surfaces | (0) | (integer) | This is the total number of surfaces where the 3D model exchanges Pressure information with the 0D lumped parameter network model. For details the reader should read [this publication](docsRefs.html#refSec2) |

List of Dirichlet Surfaces | NO DEFAULT | (space-separated integer list) | This is the list of surface IDs where the 3D model exchanges pressure information with the 0D lumped parameter network model |

Number of Neumann Surfaces | (0) | (integer) | This is the total number of surfaces where the 3D model exchanges flow rate information with the 0D lumped parameter network model. For details the reader should read [this publication](docsRefs.html#refSec2) |

List of Neumann Surfaces | NO DEFAULT | (space-separated integer list) | This is the list of surface IDs where the 3D model exchanges flow rate information with the 0D lumped parameter network model |

Command | Default | Possible Values | Description |
---|---|---|---|

Residual Control | (True) | False,True | Activates the possibility of adjusting the number of iterations based on the value of the residual norm at the current iteration |

Residual Criteria | (0.01) | (double) | Lower bound on the residual norm that triggers convergence. In other words, if the residual norm is lower than this value the solver will consider the current time step converged and continue to the next step |

Minimum Required Iterations | (3) | (integer) | This are the number of iterations that are performed independently on the value of the residual norm in the current time step |

Step Construction | (0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1) | (Sequence of 0 and 1) | Non linear iteration sequence. It controls the non-linear iteration loop within the time step . The Navier-Stokes equations constitute a no-linear system of partial differential equations. Like any non-linear system, in order to find a solution for a given time step, we must undergo an iteration process to obtain a solution that reduces the residual (i.e., the error) as we iterate more and more. The 0 tells the solver to make a solve, the 1 to make an update of the solution. For example, a ten iterations solution is specified as default as a sequence of solve/update operations |

Command | Default | Possible Values | Description |
---|---|---|---|

Deformable Wall | (False) | False,True | Activates/deactivates the deformability of vessel walls using the coupled-momentum method **LINK** |

Number of Wall Properties per Node | (10) | 10,21 | Defines the material model for the vessel wall. Option “10” (default option) means that an isotropic material model needs to be used. Option “21” means that an orthotropic material model needs to be used |

Density of Vessel Wall | NO DEFAULT | (double) | Mass density of the vessel wall material |

Thickness of Vessel Wall | NO DEFAULT | (double) | Uniform thinkness of the vessel wall material |

Young Mod of Vessel Wall | NO DEFAULT | (double) | Uniform elastic modulus for the vessel walls |

Poisson Ratio of Vessel Wall | (0.5) | (double in [0.0,0.5]) | Uniform Poisson ratio for the vessel walls |

Shear Constant of Vessel Wall | NO DEFAULT | (double) | Uniform Shear constant for the vessel walls |

Wall Mass Matrix for LHS | (True) | False,True | Assembles the contribution of the wall mass matrix in the global LHS matrix |

Wall Stiffness Matrix for LHS | (True) | False,True | Assembles the contribution of the wall stiffness matrix in the global LHS matrix |

Variable Wall Thickness and Young Mod | (False) | False,True | Activates/deactivates the possibility of specifying a variable thickness/elastic modulus for the vessel wall material |

Command | Default | Possible Values | Description |
---|---|---|---|

Solver Task | (Full Simulation) | Full Simulation,Only Partition | Decide whether to do full simulation or just partition restart.01 and geombc.dat.1 |

Impedance boundary conditions are defined through the **Qhistor.dat** file where the flow rate time history at each face is specifiedd and the **impt.dat** files containing impedance data. These two files must be present in the folder with all others solver files when executing **svSolver**.

The **Qhistor.dat** file has the following format:

```
totTS
Q0,0 Q0,1 ... Q0,n
Q1,0 Q1,1 ... Q1,n
...
...
QtotTS,0 QtotTS,1 ... QtotTS,n
```

where **n** denotes the total number of faces with impedance boundary condition applied and **totTS** the total number of time steps.

The **impt.dat** file has the following format:

```
nptsImpmax
numDataImp,1
t0 impVal0
...
t_numDataImp impVal_numDataImp
...
...
numDataImp,n
t0 impVal0
...
t_numDataImp impVal_numDataImp
```

where **nptsImpmax** is the maximum number of time steps defined on all the boundary faces with impedance boundary condition.
There are **n** blocks in the file, each defining impedance data for each face. Every block is defined by two columns, the first denoting time and the second impedance values.

RCR boundary conditions are defined through the **rcrt.dat** file using the following format:

```
nptsRCRmax
numDataRCR_1
Rp_1
C_1
Rd_1
time_1_1 Pdist_1_1
...
...
time_1_numDataRCR Pdist_1_numDataRCR
...
...
numDataRCR_i
Rp_i
C_i
Rd_i
time_i_1 Pdist_i_1
...
...
time_i_numDataRCR Pdist_i_numDataRCR
```

The first quantity **nptsRCRmax** defines the maximum number of time points for the curves defined at each outlet.
This quantity is followed by one block for each outlet, containing **numDataRCR_i**, i.e., the number of time point for RCR outlet i.
The three values **Rp_i**, **C_i**, **Rd_i** are defined for the proximal resistance, compliance and distal vessel resistance, respectively.
A time series follows, defining the evolution in time of the reference pressure at the distal end of the RCR block.

Coronary boundary conditions are defined through the **cort.dat** file using the following format

```
nptsCORmax
numDataCOR_1
q0_1
q1_1
q2_1
p0_1
p1_1
p2_1
b0_1(=0)
b1_1
b2_1(=0)
time_1_1 Plv_1_1 (/Prv_1_1)
...
...
time_1_numDataCOR Plv_1_numDataCOR
...
...
...
numDataCOR_i
q0_i
q1_i
q2_i
p0_i
p1_i
p2_i
b0_i(=0)
b1_i
b2_i(=0)
time_i_1 Plv_i_1 (/Prv_i_1)
...
...
time_i_numDataCOR Plv_i_numDataCOR
```

The first quantity **nptsCORmax** defines the maximum number of time points for the curves defined at each outlet defining the ventricular pressures.
This quantity is followed by one block for each outlet, containing **numDataCOR_i**, i.e., the number of time point for Coronary outlet i.
Nine constants need also to be defined for each coronary outlet, i.e., $q_0$, $q_1$, $q_2$, $p_0$, $p_1$, $p_2$, $b_0$, $b_1$, $b_2$.
The physical meaning of these constants is related to the resistances and capacitances used to simulated each coronary outlet, as shown in the picture below:

The following expressions are used in this paper.

$$ p_0 = 1,\quad p_1 = R_{a-micro}C_a + (R_v + R_{v-micro})(C_a + C_{im}),\quad p_2 = C_{a}\,C_{im}\,R_{a-micro}(R_v + R_{v-micro}). $$

$$ q_0 = R_{a} + R_{a-micro} + R_{v} + R_{v-micro},\quad q_1 = R_{a}C_{a}(R_{a-micro} + R_{v} + R_{v-micro}) + C_{im}(R_{a} + R_{a-micro})(R_{v} + R_{v-micro}). $$

$$ q_2 = C_{a}C_{im}R_{a}R_{a-micro}(R_v + R_{v-micro}),\quad b_0 = 0,\quad b_1 = C_{im}(R_v + R_{v-micro}),\quad b_2 = 0. $$