Introduction

This tutorial demonstrates how to use the General Boundary Conditions (GenBC) framework to define custom inlet and outlet boundary conditions for a computational fluid dynamics (CFD) simulation in SimVascular. It will cover all of the settings and parameters that need to be specified. We will follow a simple example of a cylinder with a prescribed Dirichlet inflow and RCR outflow, which should cover the basics needed to extend the framework to most situations. Additional notes will be included to explain any additional steps needed to extend the framework for more complex models. If you have any questions or if anything is unclear, please do not hesitate to contact us through the public forums!

What is GenBC?

GenBC provides a framework to programmatically define custom inflow and outflow boundary conditions for a CFD simulation. The framework allows users to create an arbitrary lumped parameter network (LPN) layout suitable for their application. Some common examples include a lumped parameter heart model that models contraction of the heart chambers to use as an inlet boundary condition, sophisticated models of the downstream circulation for various areas of the body such as the legs and upper body, or a closed-loop formulation where all outflow of the SimVascular model returns back to the inflow after passing through the veins, heart, and pulmonary arteries.

The GenBC framework is implemented as a Fortran program called by the SimVascular flow solver svSolver. Users must derive the governing ordinary differential equations for the lumped parameter layout and implement them in Fortran inside the GenBC program. The Fortan program is then compiled to produce a GenBC executable program. This executable is called by svSolver during execution to provide values for custom boundary conditions.

SimVascular does not currently provide functionality to define GenBC boundary conditions using the GUI. The GUI is used to generate simulation files for built-in boundary conditions. These files are then manually edited to incorporate the commands needed to define data for the boundary conditions and to tell svSolver to use the GenBC framework.

Cylinder RCR example

This tutorial demonstrates how to use GenBC to define inflow and RCR boundary conditions for flow in a cylinder. Although these boundary conditions can be defined using the built-in features in SimVascular they are used here as a simple example of all the steps needed to implement GenBC boundary conditions. A graphical representation of the boundary conditions for the cylinder example is shown in the figure below.

Cylinder with sinusoidal inflow and RCR outflow.

The SimVascular project for the cylinder example can be downloaded using the Cylinder Example Project link in the menu on the left-hand side of this page. The project contains all of the data (image, model, and mesh) needed for the tutorial. You will need to manually load the mesh by selecting Meshes->cylinder_mesh->Load/Unload Mesh from the SV Data Manager menu. The project has two additional folders

1) GenBC-program - Contains the Fortran code implementing the GenBC framework 
2) flow-file - Contains the **inlet.flow** file used to to define the flow rate for Direchlet boundary conditions.

Preparing SimVascular solver files for GenBC

We will assume at this point that you will have already built the geometry of your model and meshed it using the model construction pipeline in SimVascular. This tutorial will start right at the Simulations module. Create a new Simulation job by right-clicking the Simulations item in the Data Manager, then click Create New Simulation Job. Select the appropriate Model from the drop down menu and give it the name cyl_sim. Navigate through the various tabs to set simulation parameters discussed in Quick Guide / Simulations. In particular, the values for Number of Time Steps, Time Step Size and Number of Timesteps between Restarts parameters under the Solver Parameters tab must be set.

During execution the solver will normalize the velocity profile of the Dirichlet inlets to have a flow of 1.0 and then multiply the velocity profile by the flow that is computed in GenBC. Caps with Dirichlet boundary conditions are typically connected to an inductor or have a specified flow waveform (e.g. measured or adapted from literature). For the cylinder example we want to specify the inlet face to have a Dirichlet boundary condition. To set a Dirichlet boundary condition for the inlet cap double-click the box in the row of the cap Name and the column of the BC Type and select Prescribed Velocities from the drop-down menu. Double-click the cap name under Name to open up a window with additional parameters. Set BC Type to Prescribed Velocities and set Analytic Shape to parabolic. For the flow rate select the inlet.flow file from the project’s simulation-files directory.

Prescribed inflow settings window.

A Neumann boundary condition is typically used for caps that interface directly with either a capacitor or a resistor. For the outlet cap we specify a Neumann boundary condition by selecting Resistance for the BC Type and assigning it a value of 0. These settings are used to create the appropriate entries in the solver.inp file we will later modify. The Inlet and Outlet BCs tab should now look like this

Inlet and outlet conditions.

To create the data files for the simulation go to the Create Files and Run Simulation tab, select the appropriate mesh from the drop down menu and click Create Data Files for Simulation. The simulation files are written to the project directory inside the Simulations folder. In this folder you will see a folder for the simulation job you created named cyl_sim. Inside this folder are the following files and folders:

1. mesh-complete – Contains mesh information and discretization in VTK format files
2. bct.dat and bct.vtp – Contains velocity profile information on your Dirichlet caps
3. cyl_sim.svpre file – Pre-solver script which instructs SimVascular how to prepare the simulation files
           from your mesh. Contains the ID numbers for mesh surfaces (used later)
4. geombc.dat.1 – Simulation file that contains geometry information
5. inlet.flow - .flow file which has the template flow information at your inlet face(s)
6. numstart.dat – Dummy file which tells the solver on which timestep to start on
7. restart.0.1 – Contains initial conditions for pressure and velocity in your simulations
8. solver.inp – Solver settings file

Surfaces are identified using integer IDs automatically generated by SimVascular and assigned to surfaces in the .svpre file using the set_surface_id_vtp command

Surface IDs in .svpre file.

We need to manually modify the solver.inp file to use GenBC. Edit the solver.inp file with a text editor and scroll down to the section right below Step Construction. You will see following three commands for Resistance boundary conditions

Number of Resistance Surfaces: 1
List of Resistance Surfaces: 3
Resistance Values: 0

These commands were created for the outlet cap Neumann boundary condition but are not needed and so we delete them.

We next need to add commands to instruct the solver to look for GenBC and all the inlet/outlet surfaces. Below the Step Construction command add

Find the GenBC Inside the Running Directory: True
Number of Timesteps for GenBC Initialization: 0

We next add commands identifying the inlet Dirichlet surfaces

Number of Dirichlet Surfaces: 1
List of Dirichlet Surfaces: 2

The Number of Dirichlet Surfaces command specifies the total number of caps in your model where Dirichlet boundary condition are defined for GenBC. The List of Dirichlet Surfaces is a space-delimited list surface IDs for the Dirichlet boundary conditions. For the cylinder model the inlet is assigned ID 2.

We now add a similar commands for surfaces with Neumann boundary conditions

Number of Neumann Surfaces: 1    
List of Neumann Surfaces: 3 

We now change the value of Number of Coupled Surfaces to the total number of Dirichlet and Neumann surfaces in your model. This is to indicate that these surfaces are coupled in the sense that the flow and pressure have the potential to be coupled in the GenBC formulation. For this example the number of coupled surfaces is 2 because we have 1 Dirichlet surface and 1 Neumann surface.

After making these changes the solver.inp file should look like

Configured solver.inp for GenBC simulation of an RCR cylinder.

Setting up GenBC

The Fortan files implementing GenBC are found in the project directory under the GenBC-program folder. This folder contains three Fortran source files (GenBC.f, Modules.f and USER.f) and a Makefile.

All of the implementation will be done in USER.f. Open this file with a text editor. We will proceed down this code and explain it in blocks. The first block that the user can change is the “unit conversion” block:

Unit conversion block of GenBC.

In this block, you can specify conversion factors between pressure and flow in the GenBC code, and pressure and flow in the SimVascular flowsolver. In general, SimVascular models and flowsolvers use the CGS unit system (centimeter-gram-second). The units for pressure and flow in this system are dyne/cm^2 and cm^3/s, respectively. In the GenBC framework, it is often convenient to work with “clinical” units of pressure (i.e. mmHg). The conversion factor between mmHg and dyne/cm^2 is 1334 (dyne/cm^2) / mmHg, hence the value of 1.334D3 for the pConv variable. The clinical unit for flow is mL/s, which is equivalent to cm^3/s, so we leave the qConv variable at 1.0. If you want to develop your LPN model using CGS units for the parameter values, then you can change both Conv variables to 1.0. Note that the system you will choose to work with will affect the values of the parameters in your LPN network, and the values for the initial conditions in your LPN network as will be described below.

The pCorr and qCorr variables are seldom used and can be left as .FALSE. The next block is the inputs block:

Inputs block of GenBC.

In this block, you must specify the number of Dirichlet and Neumann surfaces that match you specified earlier in your solver.inp. Next, you need to specify the number of unknowns. This is the total number of state variables that you will solve for in your boundary condition model. In this case, we have just a simple model with 2 unknowns for our 2 caps, so we set this number to 2. But more sophisticated models will often involve solving for more unknowns than there are caps. This is common for heart model LPNs, where the volumes, pressures, and flows need to be computed as state variables for each chamber of the heart. But these “internal variables” do not necessarily get applied as boundary conditions to the 3D model.

The next parameter is the number of sub time steps used in time-integrating the system of ordinary differential equations to the next time step. Because the boundary condition model is often significantly cheaper than computational fluid dynamics, we can afford to iterate the boundary condition model at a much smaller time step size. This parameter specifies how many timesteps will be used to integrate the boundary condition model IN BETWEEN time steps of the 3D model. 100 is a good number here and you will likely never have to change this.

Nxprint specifies the number of extra information you would like to output. This will be further explained later, but for now we can leave this as 3.

The next block is the initial conditions block:

Initial conditions block of GenBC.

Because the boundary condition model is comprised of a system of ordinary differential equations, we need to specify the initial values for all the state variables. In this example, we have defined the first state variable to be the inflow to the model, and give it an initial condition of -10 mL/s (the negative is to denote inflow, rather than outflow of the model). The second state variable is the pressure on the RCR capacitor, which we initialize to 90 mmHg.

The next block is the surface mapping block.

Surface mapping block of GenBC.

In this block, we specify the mapping between state variables in the boundary condition model and surface IDs in the SimVascular simulation. Recall that in our solver.inp we had to specify a list of the surface IDs for the Dirichlet and Neumann surfaces. These two vectors (srfToXdPtr and srfToXPtr) are the corresponding lists in the boundary condition model. The elements of these vectors are the state variable indices that will be applied as boundary conditions on the model.

Take the first vector as an example. This vector is called “srfToXdPtr” and specifies the mapping between state variables and the Dirichlet surfaces of the SimVascular model. Recall from the section on the initial conditions that we defined the first state variable (i.e. index 1) to correspond to the flow at the cylinder inlet. Hence, the vector “srfToXdPtr” has a 1 to indicate this.

The second vector “srfToXPtr” specifies the mapping between state variables and the Neumann surfaces in our SimVascular model. Since we are using the second state variable for the RCR capacitance pressure, we put a 2 in this list.

Additional entries can be added to these vectors for more complex models separated by commas:

srfToXdPts = (/1,4,5,6/)

A good safety check at this point is the make sure the vectors in this block have the exact same length as the corresponding lists of surface IDs in solver.inp. SrfToXdPtr should have the same length as “List of Dirichlet Surfaces” in solver.inp, and srfToXPtr should have the same length as “List of Neumann Surfaces”. Each element in the USER.f lists should have a matching pair in the solver.inp lists. To recap for this example, state variable 1 gets applied to surface ID 2 as a Dirichlet boundary condition, and state variable 2 gets applied to surface ID 3 as a Neumann boundary condition.

The final block in USER.f is the equations block:

Equations block of GenBC.

In this block, we specify parameter values, write out the governing differential equations, and specify any extra outputs we would like to be printed out. We first define the parameters of our RCR outlet circuit. We need to specify three parameters: proximal resistance Rp, capacitance C, and distal resistance Rd. As is typical of Windkessel models, the distal resistance is significantly higher. Here, we use a simple 10:1 ratio, which is fairly standard. Typically, the values of these parameters are specified or tuned such that the simulation results match either population average data or patient-specific data. Here, we just use nominal values for this example. Note the Fortran method for specifying variables. First the type and name of the variables must be declared:

REAL(8) Rp, C, Rd

Then, the values can be set on subsequent lines. Also note how all non-comment lines of code start in column 7, as is required by the Fortran programming language.

Next, we specify the governing differential equations for our LPN model. These are saved in the output vector “f”. Note that “f” has length equal to the number of unknowns. Since we have 2 unknowns in this cylinder model, we need to specify a differential equation for both. For our inflow (i.e. state variable 1) recall that we set an initial condition of 10 mL/s. Just for the sake of this example, let’s apply a simple sinusoidal oscillation such that we get a little unsteadiness in our result. We apply a time derivative of cos(15t) to the inflow. Note that if a steady inflow is desired, simply set the time derivative to 0. This will cause the boundary condition model to maintain the initial condition.

For the RCR circuit, we apply the differential equation governing the pressure evolution on the capacitor. Assuming an inflow of Q into the capacitor, P is the current pressure on the capacitor, and applying Kirchoff’s laws of current and voltage conservation, the time derivative of the pressure is:

RCR equation for pressure.

This is reflected in the code in the expression for f(2). Finally, we need to add an offset for the proximal resistance. The offset serves a purpose of adding to or subtracting from a state variable before it gets applied as a boundary condition for the fluid dynamics simulation. In this case, the pressure that is applied on the simulation needs to be offset by the action of the proximal resistance. This offset pressure is simply the flow into the capacitor multiplied by the proximal resistance.

A quick explanation about the P and Q variables. Note how these are inputs to the FINDF subroutine we have been modifying. These are pressures and flows that are computed by the fluid dynamics simulation and passed to the boundary condition model as inputs. For Dirichlet boundary conditions, SimVascular computes an area-averaged pressure across that face and passes it to P vector. For Neumann boundary conditions, SimVascular computes the flow out of that outlet and passes it to Q vector. For most coupled boundary conditions, these P and Q values are needed as inputs. Take our RCR outlet condition for example. It needed the flow out of the outlet face (i.e. Q(1)) in order to specify the differential equation for the pressure. This is what is meant by “coupled” boundary conditions, where the pressure and flow at the caps of our model are intimately linked through the boundary condition model. We did not make use of the P vector in this case since we are using a special case of a completely prescribed inflow. But this will be used for more general Dirichlet boundary conditions.

Another quick explanation on the “x” variable. This vector contains the current values for the state variables. In this case, x(1) contains the current value of flow at the inflow, and x(2) contains the current value of pressure on the RCR capacitor. The current values of the state variables are often used in numerical solutions of systems of ordinary differential equations, so this x variable is how you access them. In this example, we only utilized the current pressure on the capacitor in specifying f(2), but more complex models will make more use of the “x” variable.

We recommend writing out the governing equations using a pencil and paper before implementing them in the Fortran code. This will keep you organized and ensure that the equations are being derived and implemented in the code correctly.

Finally, we can specify additional outputs that we would like to see in the Xprint vector. The variables specified here will appear in the outputs. This part is completely optional, but it can be useful to print out additional variables when debugging. Recall in the input block that we specified 3 Xprints, which will allocate 3 slots in the Xprint vector. You can raise this number if you would like to print out more. For this case, we have printed out three variables: the current time (t), the flow at the Neumann outlet (Q(1)), and the value of the proximal resistance offset.

Compiling GenBC

Once all of the settings and equations have been set in USER.f you can compile the source code. Navigate to the GenBC-program folder on your Linux or Mac machine using the terminal and type make. This creates an executable named GenBC in the current folder. Copy the executable to the project Simulations cyl_sim folder.

Running a flow simulation and checking results

Now navigate to the Simulations cyl_sim folder. You can now run the svSolver program. Once the simulation starts it will create a file named AllData. This is an output file produced by GenBC which contains the values of the state variables as a function of time and the additional Xprint outputs we specified earlier.

Prescribing a pressure waveform as an inlet boundary condition

This section of the tutorial describes how to use the GenBC framework to define a Neumann boundary condition as a pressure waveform at an inlet surface. The pressure boundary condition is implemented in Fortran using the GenBC framework. The pressure at the cylinder inlet is obtained by interpolating a pressure waveform defined in the GenBC Fortan code.

Pressure waveform used for the inlet Neumann boundary condition.

The pressure waveform inlet boundary condition simulation is similar to the RCR simulation we created above except for changes to the boundary conditions and some changes in the simulation parameters needed for solver stability. Create a new Simulation job by right-clicking the Simulations item in the Data Manager, then click Create New Simulation Job. Select the appropriate Model from the drop down menu and give it the name cyl_pres. Navigate to the Solver Parameters tab and set

Number of Time Steps: 2000 

Time Step Size: 0.000588 

Number of Timesteps between Restarts: 20 

We will define Neumann boundary conditions for both the inlet and outlet caps selecting Resistance for the BC Type and assigning it a value of 0. These settings are used to create the appropriate entries in the solver.inp file we will later modify. The Inlet and Outlet BCs tab should now look like this

To create the data files for the simulation go to the Create Files and Run Simulation tab, select the appropriate mesh from the drop down menu and click Create Data Files for Simulation. The simulation files are written to the project directory inside the Simulations folder. In this folder you will see a folder for the simulation job you created named cyl_pres. Now navigate to the cyl_pres folder. Inside this folder are the following files and folders:

    1. mesh-complete – Contains mesh information and discretization in VTK format files
    2. cyl_pres.svpre file – Pre-solver script which instructs SimVascular how to prepare the simulation files
           from your mesh. Contains the ID numbers for mesh surfaces (used later)
    3. geombc.dat.1 – Simulation file that contains geometry information
    4. numstart.dat – Dummy file which tells the solver on which timestep to start on
    5. restart.0.1 – Contains initial conditions for pressure and velocity in your simulations
    6. solver.inp – Solver settings file

We need an additional rcrt.dat file containing the RCR parameters that set the relationship between flow and pressure on the outflow face. Copt this file from the simulation-files folder to here.

the presover solver script (svpre) is needed to be modified by treating the inlet as a Neunmann face (see steady_rcr.svpre). Thus, to prepare a simulation with a prescribed pressure waveform, we run the presolver with the modified presolver script to generate geombc.dat and restart.0.1 files with surfaces properly tagged. Then we compile the GenBC files (user defined Fortran files) before calling the flow solver.

Note that apart from potentially increased computational cost and numerical instability (which might be insignificant), prescribing a pressure waveform doesn’t necessarily make flow simulations better agree with clinical target values. So extra adjustments are often needed in order to match target pressures and flow. Don’t forget to check resulting inflow when the Neumann BC (pre-defined pressure or a lumped parameter model defined pressure) is applied to the inlet.