Introduction

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

What’s New?

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

About this guide

This document will teach you:

  1. setting initial conditions and boundary conditions
  2. setting mechanical properties for vessel walls (if deformable)
  3. setting parameters for flowsolver
  4. running the flow solver
  5. converting and analyzing the simualtion 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.

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 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 meshe. They are organized as shown in the example below.

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.

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.

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

Creating a Job for Simulation

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 alread 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 create input and data files yet.

Basic Parameters

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

Inlet Boundary Condition Specification

Let’s first create a simple file to provide inlet flow rates before seting 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 subfolder “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.

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.

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

Outlet Boundary Condition Specification

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.

Advanced Options for splitting BC

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"

Wall Property Specification

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

Go to "Wall Properties"
Type: Rigid

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

Solver Parameters

There are many parameters for flow solver, but only a few is required to set explicitly. Advanced paramters are optional. To check out the meannin 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.

Time Step Parameters Block

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.


Output Control Block

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


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

Creating Data Files for Simulation

Before running simualtion, 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 numberr 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 , ...


Creating Data Files

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)

Running Simualtion Jobs

Before running simualtion, you may need to click the button “Import Extra Files into Job” to add extra file into the simualtion job. In this example, we don’t need importint 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

At the same time, the solver copies all these files to a newly created folder called 8-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.

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

It shows the name of th running job, the percentage completed, and more info: [current time step] [time used(s)] [nonlinear residual] [logarithm of the redisual 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 “8-procs_case”. 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 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.

Running Simualtion Jobs at Clusters

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 requied 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 simualtion, 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.

Post-processing

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

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

    Goto "Export Results"
    Result Dir: select the running dir .../Simulations/steady/8-procs_case
    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 "Export..."
    Choose a directory for exporting

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

Creating VTU result files from svSolver restart files.

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

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.

Example 2

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/8-procs_case
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 (=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 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 displacment to a vtp file “displacement.vtp” to review the solution, and finaly append it to restar.0.1 we just copied from restart.500.0 of Example 2.

Initila Displacement from displacement.vtp

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

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.

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 displacment.vtp are created, which show the thickness and Young’s modulus, and initial displacement, respectively:

Variable Wall Thickness in varwallprop.vtp
Initila Displacement in displacement.vtp

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

Appendix

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

Configurating Solvers (Optional)

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.

Locations of Installed Solvers

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

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)

Format of BCT file

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.

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. Use 10 for pulsatible flow, or deformable wall casese
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, 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. $$