Introduction

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.

What’s New?

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.

About this guide

This document will teach you the fundamentals of:

  1. svPre: Preparing the necessary svSolver input files
  2. svSolver: Running a flow analysis
  3. 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.

Theory and Implementation

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

Overview

Process Flow of SimVascular Simulation

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

Workflow for generating hemodynamic results of a cylindrical model starting from a stereolithography of its exterior surface

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.

Units in Simulation

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

Useful constants in Simulation

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 Condition Specification: the Physical Side of the Problem

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).
Inflow, outflow and wall characterization in a simple cylindrical vessel

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.

Cardiovascular model with various boundary conditions
Frequency content of impedance from the literature for the left external iliac and a canine femoral artery

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$.
Circuit representation of RCR block
  • 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.

Boundary conditions considered in Example 1

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

Presolver(svPre)

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)

Prerequisite Files for svPre

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

Folder structure and file created after clicking on Write Files

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.

Cylinder with negative inflow

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

Cylinder with positive inflow

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.

Inlet Boundary Condition Specification (From GUI)

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

Creating a bct.dat file through the GUI
  • 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.

Visualizing the inlet velocity profile in Paraview

Running Script

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:

Running svPre through the GUI

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

Final recap of the files generated by svPre

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.

Flowsolver (svSolver)

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)

Prerequisite Files for svSolver

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.

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!

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

solver.inp

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.


Solution Control Block

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.


Material Properties Block

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!


Output Control Block

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


Cardiovascular Modeling Parameters Block

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

  1. 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!

  2. Number of Resistance Surfaces: 1 - This line sets the number of resistance surfaces in the model. In our case, we have one resistance surface.

  3. 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!

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

Schematic representation of a model with four outlets

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


Step sequence Block

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.

Running Simulation

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
Running simulation through the GUI

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"
Restart simulation through the GUI

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.

Post-processing (svPost)

Generating the *.vip and *.vtu files in SimVascular

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"
Creating VTU result files from svSolver restart files.

*.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 Options
  • 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”.

Calculate Pressure and Flow rate

Visualizing the results in ParaView

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.

By showing edges, the finite element mesh will become apparent
Visualizing the cylinder mesh in Paraview

Visualizing Pressure Results

  • First you should increase your current result time from 0 to 20 (the last available time step).
Choosing the last time step results in Paraview
  • You should now see the available result quantities for your model, i.e., cellsNormals, GlobalElementID, GlobalNodeID, pressure, timeDeriv, traction, velocity, WSS.
Available model results
  • 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.

Final contour of pressures in mmHg

Examples

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.

Example 2

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 (=RC=(121+1212)0.000015=0.02). With the small step size, the number of time steps is also increased.

Example 3

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:

Initila Displacement Calculated form svPre

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.

Example 4

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:

Variable Wall ThicknessInitila Calculated form svPre
Initila Displacement Calculated form svPre

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.

Appendix

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

svPre Commands

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

Read Mesh info

svPre Command Argument Format Description
mesh_and_adjncy_vtu (file name) Read node coordinates,element connectivities and adjacencies from the give vtu file

Assign IDs to Surfaces

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

Set Boundary Conditions for Walls, Inlets, Outlets

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

Creat BCT file for Inlet

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

Initial Conditions

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$

Write simulation files

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.

Read variable values from non-vtk files

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)

Read variable values from vtk files

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.

Deformable walls with uniform thickness and elastic modulus

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.

Deformable walls with variable thickness and elastic modulus

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)

Solver.inp File Guide

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

Input Control

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.

BCT Control

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.

Solution Control

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

Material Properties

Command Default Possible Values Description
Viscosity NO DEFAULT (double) Fluid viscosity
Density NO DEFAULT (double) Fluid density

Body Forces

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

Output Control

Command Default Possible Values Description
Number of timesteps between restarts NO DEFAULT (integer) Number of time steps between a new restart.. is saved. This values decides how often the complete status of the model is saved to disk.
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

Linear Solver Options

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.

Discretization Control

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

Cardiovascular Modelling Parameters

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

Backflow Control

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

Closed Loop Boundary Conditions

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

Non-linear Iteration Control

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

Deformable Wall Simulations

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

Task Control

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



Format of Impedance boundary condition file

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.

Format of RCR boundary condition file

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.

Format of COR boundary condition file

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:

Circuit visualization for coronary boundary condition

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