**SV Simulation** tool can solve 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.

**Notice**: To get full functions from the Simulation tool, it uses three solvers: **Presolver (svPre), Flowsolver (svSolver), Postsolver (svPost)**. Normally, SimVascular already includes the solvers and can find them automatically. Users don’t need to set up the solvers. However, in case SimVascular can’t find them while users are using the Simulation tool, refer to Solver Configuration.

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 Stanford 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 physiology.

This document will teach you:

- setting initial conditions and boundary conditions
- setting mechanical properties for vessel walls (if deformable)
- setting parameters for flowsolver
- running the flow solver
- converting and analyzing the simulation results

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 folder *mesh-complete/mesh-surfaces/*.

**svPre** is the preprocessor for **svSolver**. The input data 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. The **svPre** program can be called either from the command line (in terminal) or the Simulation tool (in GUI). The input data files for **svPresolver** are created from the mesh. They are organized as shown in the example below.

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 containing all the exterior elements of the mesh generated with TetGen.**walls_combined.vtp**, vtp file containing all surface elements assigned to the wall, possibly 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.

SimVascular runs **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.

**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 [Kg · m^{-1}· 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} = 15.92$ cm/s and the wall shear stress is $\tau$ = $0.64$ dynes/cm$^2$.

Let’s try Example 1. The files for Example 1 can be found from the link on the left “Cylinder Example”. Download and open the project. Now we assume you already have the model and mesh for it.

**To create a simulation job (empty):**

```
Right click the data node "Simulations" in the SV project in Data Manager
Click "Create Simulation Job" in the popup menu
Select Model: cylinder
Job Name: steady
Click "OK"
```

Now a new data node “steady” for the job is created under the data node “Simulations” in Data Manager. Double click the data node “steady" and the tool “SV Simulation” automatically shows up. The new job is empty and has not created input and data files yet.

Basic parameters include fluid material properties, and initial conditions, which are some basic information to be used for Presolver and Flowsolver. Here we use the default values (in CSG units). Initial pressure is “0”. Initial velocities (vx, vy, vz) is “0 0 0”.

Let’s first create a simple file to provide inlet flow rates before setting inlet boundary condition. Your problem may have more that one inflow wave form file. In this case, we only have a single flow file (called steady.flow). We put the file in a sub-folder “flow-files” under the project folder in the disk.

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.

Let’s start setting inlet boundary condition.

```
Go to "Inlet and Outlet BCs"
All the faces of type "cap" are already listed in the table
Double click the row for the face "inflow"
A dialog pops up and also the face is highlighted in the 3D-view window.
In the dialog:
BC Type: Prescribed Velocities
Analytic Shape: parabolic
Point Number: 2
Fourier Modes: 1
Click the tool button "..." to select the file "steady.flow" we just created as above
Period: (filled automatically based on the data from the file)
Click "OK"
```

**Help Hints:** Please make use you specify face types when you create the corresponding model;otherwise, the cap faces won’t appear.

**Point Number**: This is the number of*temporal*data points that you want to have in one cycle. 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 interpolated data from it 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.**Fourier Modes**: 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).

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

Let’s start setting outlet boundary condition.

```
Double click the row for the face "outlet"
A dialog pops up and also the face is highlighted in the 3D-view window.
In the dialog:
BC Type: Resistance
Resistance: 1333 (resistance value)
Distal Pressure: (0 by default)
Click "OK"
```

**Help Hints:** Please make use you specify face types when you create the corresponding model;otherwise, the cap faces won’t appear.

In many clinical cases, we would like to split a total resistance or capacitance among a group of outlets, instead of explicitly assigning resistance or capacitance for each individual outlet.

To use this feature, first set up the group of outlets with initial zero resistance or capacitance.

For resistance BC:

```
Select a group of outlets, right click
Click "Set BC"
BC Type: Resistance
Resistance: 0
Click "OK"
```

For RCR BC:

```
Select a group of outlets, right click
Click "Set BC"
BC Type: RCR
Rp, C, Rd: 0 0 0
Click "OK"
```

For Coronary BC:

```
Select a group of outlets, right click
Click "Set BC"
BC Type: Coronary BC
Rp, Ca, Ra-micro,Cim, RV: 0 0 0 0 0
Pim (from file): Click the button "..." to select file for Pim
Pim Period: filled automatically using the data from the file. You can change the period, SimVascular will automatically scale the time from the file
Pim Scaling: use this value to scale presure values from the file
Click "OK"
```

After you set up the group, let’s split the total resistance or capacitance.

To split total resistance:

```
Select the group of outlets, right click
Click "Split Resistance"
Total Resistance: total resistance for the group of outlets
Murray's Law Coefficient: (2~3), 2 in general, 2.6 for coronary arteries, 2.3 for pulmonary arteries
Ratio of resistance for each outlet:
Rp:Rd (for RCR BC):
Ra:Ra-micro:Rv (for Coronary BC):
Click "OK"
```

To split total capacitance:

```
Select the group of outlets, right click
Click "Split Capacitance"
Total Capacitance: total capacitance for the group of outlets
Ratio of capacitance for each outlet:
Ca:Cim (for Coronary BC):
Click "OK"
```

In Example 1, the wall is rigid and set no-slip boundary condition.

```
Go to "Wall Properties"
Type: Rigid
```

SimVascular also has options like “Deformable(Constance)” and “Deformable(Variable)”. We will explain them in Examples 3 and 4.

There are many parameters for flow solver, but you are only required to set explicitly a few of them. Advanced parameters are optional. To check out the meaning of all the parameters, refer to appendix.

```
Number of Timesteps: 200
Time Step Size: 0.03
Number of Timesteps between Restarts: 10
Step Construction : 2 # standard two iterations, enough for constant problems.
```

**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 before. 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$.

**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.

**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 one in every $10$ files. Therefore, the output files of the solver will look like this: restart.0.*, restart.10.*, restart.20.*, restart.30.*, …., restart.190.*, restart.200.*

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. Each iteration tells the solver to make a solve and an update of the solution.

**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 parameters the **svSolver** can take. A more detailed discussion can be found in this section.

Before running simulation, you need to create some required data files for presolver and flowsolver.

For Presolver:

**.svpre**: script input file

For Flowsolver:

**geombc.dat.1**: containing mesh info and boundary conditions specified in the problem, created by presolver**restart.0.1**: containing initial conditions for our problem, , created by presolver**numstart.dat**: contains a number**0**. This number is used by the solver to identify the restart file that should be used as initial condition.**bct.dat**: containing time-dependent velocity vectors at the inflow face of the model according to a prescribed flow wave coming from a *.flow file. See this section.**solver.inp**: providing further info for flowsolver, specifying parameters such as time step size, number of time steps, number of nonlinear iterations, boundary condition control, etc. A detailed description is 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.

**Impedance Boundary Condition simulations:**

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.

**RCR Boundary Condition simulations:**

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.

**Coronary Boundary Condition simulations:**

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.

**Closed-loop boundary conditions:**

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!

**HINT**: In both files geombc.data.1 and restart.0.1, 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 , ...
```

This step creates files inside the job folder in the SV project

```
Choose Mesh: cylinder
Click the button "Create Data Files for Simulation"
```

A new directory using the job name “steady” is created under the folder “Simulations” and contains the following folder/files:

```
For presolver:
mesh-complete (folder)
inflow.flow
steady.svpre
For flowsolver
geombc.dat.1
restart.0.1
bct.dat (and bct.vtp)
solver.inp
numstart.dat
rcrt.dat (if applicable)
cort.dat (if applicable)
```

Before running simulation, you may need to click the button “Import Extra Files into Job” to add extra file into the simulation job. In this example, we don’t need important extra files

**To run a simulation job:**

```
Number of Processes: 8 (select the maximum number in your case)
Starting Step Number: (be default, use 0 or the last time step from previous simulation)
Click the button "Run Job"
```

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

```
geombc.dat.1 => geombc.dat.1, geombc.dat.2, ..., geombc.dat.8
restart.0.1 => restart.0.1, restart.0.2, ..., restart.0.8
```

During the job running, the progress of simulation is shown at “Job Status” and the status bar at the the bottom of the main window:

It shows the name of the running job, the percentage completed, and more info: [current time step] [time used(s)] [nonlinear residual] [logarithm of the residual change] …

During the job running, you can terminate it before it finishes. To stop the simulation:

```
Right click the job data node "steady" at Data Manager
Click "Stop Simulation" in the popup menu
```

After the simulation is finished. A dialog pops up to inform that the job is done. You can click the button “Show Details…” to get more info about the simulation progress. You can also check the simulation progress by open the file “histor.dat” in the folder “steady”. 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 residual 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 continue the finished simulation job. Make use the same number of processes as the previous one and just click the button “Run Job”. Now 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. We’ll convert these files to generate the visualization files. We explain that process in the following section.

You can also export the required data files and upload to a computer cluster to run the simulation. Make sure SimVascular flow solver is available in the cluster.

To export the files:

```
Make sure you have created data files for this job.
Right click the job node "steady" in Data Manager
Click "Export Data Files"
Select a directory for exporting.
```

A new folder “steady-sim-files” is created, which includes the required files by flowsolver

To run the simulation at the cluster:

```
Make sure SimVascular flow solver is available on the cluster.
Upload the folder "steady-sim-files" to the cluster
Login the cluster
Go to the folder "steady-sim-files"
Run "mpiexec -n [number of processes] [solver path]/svsolver" or you need create/submit a job file as required by the cluster to run the simulation
```

During the simulation, result files are saved at the folder “steady-sim-files/[number of processes]]-procs_case/”. You can download the files back to your computer and convert them to vtp/vtu files.

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

```
Go to "Convert Results"
Result Dir: select the running dir .../Simulations/steady/
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..."
Choose a directory for exporting
```

A new folder “steady-converted-results” is created, which contains all the converted vtp and vtu files.

After converting is done, there are also four text files created in the exported folder if “Calculate Flows” is toggled on.

```
all_results-flows.txt: flowrate for each face with time steps
all_results-pressures.txt: pressure for each face with time steps
all_results-averages.txt: the average, maximum, minimum values of presure, flowrate for each face
all_results-averages-from_cm-to-mmHg-L_per_min.txt: same info as all_results-averages.txt, but pressure is in mmHg, flowrate is L/min.
```

Other options are also provided:

**As 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 restart 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.

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.

This example shows how to simulate a cylinder with plug steady flow at the inlet and RCR at the outlet.

Let create a new job by copying/pasting the job of Example 1.

```
Right click the data node "steady" in Data Manager
Click "Copy"
Right click the data ndoe "Simulations" in Data Manager
Click "Paste"
A new data node "steady_copy" is created.
Rename it as "steady_rcr"
```

Let modify some BCs and parameters.

```
Go to "Inlet and Outlet BCs"
Double click the row of "inflow", Analytic Type: plug
Double click the row of "outlet", BC Type: RCR, Rp,C,Rd: 121 0.000015 1212
```

```
Go to "Solver Parameters"
Number of Timesteps: 500
Time Step Size: 0.001
Number of Timesteps between Restarts: 20
```

Similarly to Example 1, to run the job

```
Go to "Create Files and Run Job"
Choose Mesh: cylinder
Click the button "Create Data Files for Simulation"
Number of Processes: 8
Staring Step Number: (leave it empty)
Click the button "Run Job"
```

To export results:

```
Goto "Convert Results"
Result Dir: select the running dir .../Simulations/steady\_rcr/
Start: the initial restart file number (0)
Stop: the final restart file number (500)
Increment: the increment between restart files (20).
Toggle on "Volume Mesh" and "Surface Mesh"
Toggle on "Last Step to restart.x.0"
Click the button "Convert ..."
Choose a directory for exporting
```

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 and exported, the restart files for the last step is reduced to a single file, which would be restart.500.0.

Let create a new job by copying/pasting the job of Example 2.

```
Right click the data node "steady_rcr" in Data Manager
Click "Copy"
Right click the data ndoe "Simulations" in Data Manager
Click "Paste"
A new data node "steady_rcr_copy" is created.
Rename it as "steady_rcr_deformable"
```

Let modify some IC, BCs and parameters.

```
Go to "Basic Parameters"
Double click "IC File" and a dialog pops up
Select the file "restart.500.0" exported from Example 2
Go to "Wall Properties"
Type: Deformable(Constant)
Thickness: 0.2
Elastic Modulus: 4000000
Poisson Ratio: 0.5
Shear Constant: 0.833333
Density: 1.0
Pressure: 133300 (initial pressure, estimated form the simulation result of the rigid case Example 2)
```

```
Go to "Solver Parameters"
Step Construction: 4
```

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

To create data files:

```
Go to "Create Files and Run Job"
Choose Mesh: cylinder
Click the button "Create Data Files for Simulation"
```

Different from Example 2, the step also solves initial displacement, write the initial displacement to a vtp file “displacement.vtp” to review the solution, and finally append it to restar.0.1 we just copied from restart.500.0 of Example 2.

Similarly to Example 2, run the job and convert the results.

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.

Let create a new job by copying/pasting the job of Example 3.

```
Right click the data node "steady_rcr_deformable" in Data Manager
Click "Copy"
Right click the data ndoe "Simulations" in Data Manager
Click "Paste"
A new data node "steady_rcr_deformable" is created.
Rename it as "steady_rcr_variable"
```

Let modify wall properties.

```
Go to "Wall Properties"
Type: Deformable(Variable)
Poisson Ratio: 0.5
Shear Constant: 0.833333
Density: 1.0
Pressure: 133300 (initial pressure, estimated form the simulation result of the rigid case Example 2)
wall: E. Modulus: 4000000
inflow: Thickness: 0.2
outlet: Thickness: 0.1
```

```
Go to "Solver Parameters"
Step Construction: 10
```

For variable wall cases, the flowsolver needs much more iterations of step sequence to get accurate solution.

To create data files:

```
Go to "Create Files and Run Job"
Choose Mesh: cylinder
Click the button "Create Data Files for Simulation"
```

Different from Example 3, the step also solves variable thickness or Young’s modulus, and assign them to the wall, instead of giving uniform thickness or Young’s modulus. varwallprop.vtp and displacement.vtp are created, which show the thickness and Young’s modulus, and initial displacement, respectively:

Similarly to Example 3, run the job and convert the results.

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

SimVascular solvers(presolver, flowsolver, postsolver) are used as individual executable programs. They are called by SV Simulation tool to create data files, run simulation, and convert results to files in VTK formats. Normally those solvers are installed separately. However, SimVascular can automatically find those solvers. In case your SimVascular can’t find them, or you want to use different solvers, you need to explicitly tell SimVascular where to find or how to use them.

```
Goto Menu: Preferences -> SimVascular Simulation
Use MPI: whether to use mpi to run flowsolver. It depends on if your flowsolver supports MPI. By default, we assume using MPI.
MPIExec: (mpiexec, by default) It's opional. Specify the full path and file name for the mpiexec to use.
Use Solver Input Custom Template: whether to use a custom tempalte file for the table of "Solver Parameters". Only for advanced users.
Custom Tempalte: the full path and file name for your custom template file.
Presolver: Specify the full path and file name for the presolver to use.
Flowsolver: Specify the full path and file name for the flowsolver to use.
Postsolver: Specify the full path and file name for the postsolver to use.
```

**Linux:**

```
Presolver: /usr/local/sv/svsolver/yyyy-mm-dd/bin/svpre
Flowsolve: /usr/local/sv/svsolver/yyyy-mm-dd/bin/svsolver
Postsolver: /usr/local/sv/svsolver/yyyy-mm-dd/bin/svpost
```

In case you can’t run mpiexec when using flowsolver, please make sure mpi is installed.

To install MPI:

```
For Ubuntu 14:
sudo apt-get install libmpich2-dev
For Ubuntu 16:
sudo apt-get install libmpich-dev
```

**Mac OS X**

```
MPIExec: /usr/local/sv/svsolver/yyyy-mm-dd/bin/mpiexec
Presolver: /usr/local/sv/svsolver/yyyy-mm-dd/svpre
Flowsolve: /usr/local/sv/svsolver/yyyy-mm-dd/svsolver
Postsolver: /usr/local/sv/svsolver/yyyy-mm-dd/bin/svpost
```

**Windows**

For the versions from 2017-04-09

```
The solvers are located at C:\Program Files\SimVascular\svSolver\yyyy-mm-dd\
Presolver: C:\Program Files\SimVascular\svSolver\yyyy-mm-dd\svpre-bin.exe
Flowsolver: C:\Program Files\SimVascular\svSolver\yyyy-mm-dd\svsolver-msmpi-bin.exe
Postsolver: C:\Program Files\SimVascular\svSolver\yyyy-mm-dd\svpost-bin.exe
Flowsolver (without mpi): C:\Program Files\SimVascular\svSolver\yyyy-mm-dd\svsolver-nompi-bin.exe
```

For the versions before 2017-04-09

```
The solvers are located at C:\Program Files (x86)\SimVascular\svSolver\Release\xxxxxxxxxx\
Presolver: C:\Program Files (x86)\SimVascular\svSolver\Release\xxxxxxxxxx\svpre-bin.exe
Flowsolver: C:\Program Files (x86)\SimVascular\svSolver\Release\xxxxxxxxxx\svsolver-msmpi-bin.exe
Postsolver: C:\Program Files (x86)\SimVascular\svSolver\Release\xxxxxxxxxx\svpost-bin.exe
Flowsolver (without mpi): C:\Program Files (x86)\SimVascular\svSolver\Release\xxxxxxxxxx\svsolver-nompi-bin.exe
```

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 applied on the element faces provided by the given vtp file |

pressure_vtp | (file name) (double) | Prescribed pressure boundary condition will be applied 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 displacement 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 thinness 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) |

Inlet prescribed velocity profile is defined through the **bct.dat** file using the following format:

```
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.

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 containing 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 inflow 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 update 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. Use 10 for pulsatile flow, or deformable wall cases |

Maximum Number of Iterations for svLS Momentum Loop | (2) | (integer) | Maximum number of iterations for the Momentum equation solver. Use 20 for deformable wall cases |

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, respectively.

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 thickness 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 specified 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. $$