svMultiPhysics - Beta Version
The contents of the documentation are not finalized and could change in the future.
Some sections of the documentation are incomplete and are actively being updated.
All User Guide sections are being updated.
Introduction
svMultiPhysics is an open-source, parallel, finite element multi-physics solver providing capabilities to simulate the partial differential equations (PDEs) governing solid and fluid mechanics, diffusion, and electrophysiology. Equations can be solved coupled to simulate the interaction between multiple regions representing different physical systems. For example, in a coupled fluid-solid simulation the motion of the fluid can deform a solid region while the changing geometry of the solid region change the way the fluid flows. Equation coupling provides a framework for the computational modeling of whole heart dynamics.
svMultiPhysics implements boundary conditions useful for the simulation of blood flow in vascular models. It also has an interface to the svZeroDSolver that can be used for custom lumped parameter network (LPN) boundary conditions.
The svMultiPhysics solver is implemented in C++, a popular general-purpose object-oriented programming language that provides powerful features for producing modular software with a clear structure. C++ has many built-in data structures that can be used to efficiently organize and store data. The svMultiPhysics code is essentially a direct line-by-line translation of the Fortran svFSI solver code into C++.
The following sections describe the solver modeling capabilities, input data and format of the solver input parameters file used to run a simulation.
Documentation describing building, installing and testing svMultiPhysics can be found using the following links
Installing and building the solver
Testing Guide
Overview
The following sections provide an overview of the capabilities of the svMultiPhysics solver. In addition, they will provide a background for concepts and terminology that will be useful in later more detailed discussions of the solver.
Modeling capabilities
svMultiPhysics can numerically solve PDEs governing solid and fluid mechanics, diffusion, and electrophysiology. Some equations can be coupled and solved simultaneously to simulate the interaction between multiple regions representing different physical system.
svMultiPhysics supports the solution of the following equations
- Linear elastodynamics - Solves the linear elasticity problem.
- Linear elastodynamics modified for mesh motion used in fluid-solid interaction - Solves for the elastic deformation of a solid region as part of the Arbitrary Lagrangian-Eulerian (ALE) formulation in a fluid-solid interaction simulation.
- Nonlinear elastodynamics - Solves the nonlinear structural mechanics problem using a pure displacement-based formulation.
- Nonlinear elastodynamics using a mixed formulation - Solves the nonlinear structural dynamics problem where the structure's velocity and pressure are the unknown degrees of freedom.
- Nonlinear thin shell mechanics - Solves the nonlinear shell mechanics problem.
- Navier-Stokes equation - Solves the incompressible viscous fluid flow problem.
- Unsteady Stokes flow - Solves for unsteady flow where inertia is negligible (creeping flow).
- Unsteady diffusion - Solves the diffusion problem in a solid.
- Unsteady advection-diffusion - Solves the unsteady advection-diffusion problem in a fluid (e.g. dye transport).
- Coupled momentum method (CMM) - Solves the fluid-solid interaction problem using the coupled momentum method.
- Arbitrary Lagrangian-Eulerian (ALE) formulation - Solves the fluid-solid interaction problem using the Arbitrary Lagrangian-Eularian (ALE) method.
- Mono-domain model of cardiac electrophysiology - Solves for the electrical propagation in myocardial tissue.
Boundary conditions
Boundary conditions are an integral part of a problem definition and act as constraints necessary for the solution of a PDE.
Different types of boundary conditions can be imposed on the boundary of the domain
- Dirichlet boundary condition - Fixed solution values
- Neumann boundary condition - Derivative of solution values
- Robin boundary condition - Weighted combination of Dirichlet and Neumann boundary conditions
A temporal and spatial distribution of values can be applied to the boundary of a domain.
svMultiPhysics supports boundary conditions useful for the simulation of blood flow in vascular models
- Windkessel (RCR) boundary condition - Electric circuit analogue that has a proximal resistance R in series with a parallel arrangement of a capacitance C and a distal resistance Rd
- Custom lumped parameter network (LPN) boundary condition - Implemented using an interface to the svZeroDSolver
Finite element method
The finite element method (FEM) is used to numerically solve transient PDEs governing the behavior of a physical system in two or three space dimensions. FEM approximates the geometry of a physical domain by subdividing it into a collection of discrete finite elements called a finite element mesh. The finite elements use numerical interpolation to approximate field variables (e.g. velocity) within a geometric region of a physical system. Elements can be implemented using a combination of linear, quadratic, and cubic interpolating polynomials.
svMultiPhysics supports the following element types
- line - linear and quadratic
- triangle - linear and quadratic
- quadrilateral - bilinear, serendipity, biquadratic
- tetrahedron - linear and quadratic
- hexahedron - trilinear, quadratic/serendipity, triquadratic
- wedge - linear
The finite element mesh together with appropriate physical properties and boundary conditions generate a system of algebraic equations. Depending on the PDE being solved this system of equations can be linear or nonlinear. The solution of a linear system can be directly solved to produce an approximation to the actual solution to the PDE. A nonlinear system must be solved using an iterative procedure to generate and solve a series of linear systems of equations converging to an approximation to the actual solution to the nonlinear PDE.
Most automatic mesh generation software packages produce finite element meshes are composed primarily of tetrahedron.
The accuracy of a simulation depends on how well the mesh approximates
- Model geometry - changes in geometry (e.g. curvature) in a region
- Field variables - spatial variation in field values
Numerical linear algebra
Numerical linear algebra uses computer algorithms to solve the linear system generated by FEM. Linear algebra libraries provide access to specialized or general purpose routines implementing a significant number of computer algorithms useful when solving linear systems.
svMultiPhysics supports interfaces to the following numerical linear algebra packages
- Custom numerical linear algebra routines included in the svMultiPhysics implementation
- Portable, Extensible Toolkit for Scientific Computation (PETSc)
- Trilinos
Parallel computation
The system of equations generated by FEM can often be quite large, involving millions of unknowns for the accurate representation of the vasculature of an anatomical region. The numerical solution of such a large system typically requires it to be distributed onto on a high-performance computational (HPC) cluster and solved in parallel. The svMultiPhysics solver uses Message Passing Interface (MPI) library routines to automatically partition and distribute an FEM simulation to run in parallel on an HPC cluster.
Mathematical preliminaries
The following use index (Einstein) notation.
Fluid mechanics
Strong form
The incompressible Navier-Stokes equations governing fluid flow are
$$ \rho\left(\frac{d\boldsymbol{u}}{dt} + \boldsymbol{u} \cdot \boldsymbol{\nabla} \boldsymbol{u} - \boldsymbol{b}\right) = \boldsymbol{\nabla} \cdot \boldsymbol{\sigma}, $$
$$ \boldsymbol{\nabla} \cdot \boldsymbol{u} = 0, $$
where $\boldsymbol{u} = \boldsymbol{u}\left(\boldsymbol{x}, t\right)$ is the velocity, $p = p\left(\boldsymbol{x}, t\right)$ is the pressure, $\boldsymbol{b} = \boldsymbol{b}\left(\boldsymbol{x}, t\right)$ is the body force, and $\rho$ is the fluid density. The first equation corresponds to the momentum conservation in the flow and the second equation corresponds to mass conservation. The momentum equation can augmented with a Darcy permeability term, $-\frac{\mu}{K}\boldsymbol{u}$, on the right-hand side to yield the Navier-Stokes-Brinkman equation [6] [7], such that
$$ \rho\left(\frac{d\boldsymbol{u}}{dt} + \boldsymbol{u} \cdot \boldsymbol{\nabla} \boldsymbol{u} - \boldsymbol{b}\right) = \boldsymbol{\nabla} \cdot \boldsymbol{\sigma} - \frac{\mu}{K}\boldsymbol{u}, $$
This equation models incompressible fluid flow in porous media. Here $K$ is the permeability of the porous media. With this equation, we can of course recover the Navier-Stokes equations by simply removing the Darcy component (i.e., $K \rightarrow \infty$).
The Cauchy stress tensor is $\boldsymbol{\sigma} = \boldsymbol{\sigma}\left(\boldsymbol{x}, t\right) = -p\boldsymbol{I} + 2\mu\left(\boldsymbol{u}\right)\epsilon$, where $\epsilon = \epsilon\left(\boldsymbol{u}\right) = \nabla^{s} \boldsymbol{u} = \frac{1}{2}\left(\nabla \boldsymbol{u} + \left(\nabla\boldsymbol{u}\right)^{\text{T}} \right)$ is the strain rate tensor. The effective dynamic viscosity, $\mu\left(\boldsymbol{u}\right)$, is written generally as a function of velocity here to account for non-Newtonian fluids. For Newtonian fluids, $\mu$ is a simply constant. The divergence of the Cauchy stress tensor, written in both vector and index notation, is
$$ \boldsymbol{\nabla} \cdot \boldsymbol{\sigma} = -\boldsymbol{\nabla}p + 2\boldsymbol{\epsilon}\boldsymbol{\nabla}\mu + \mu\nabla^{2}\boldsymbol{u}, $$
$$ \sigma_{ij,j} = -p_{,i} + 2\epsilon_{ij}\frac{\partial \mu}{\partial x_{j}} + \mu u_{i,kk}, $$
The boundary conditions are
$$ \boldsymbol{u} = \boldsymbol{g}, $$
$$ \boldsymbol{\sigma} \cdot \boldsymbol{n} = \boldsymbol{h}, $$
where $\boldsymbol{g}$ is the prescribed velocity and $\boldsymbol{h}$ is the prescribed traction.
We will solve the Navier-Stokes and Navier-Stokes-Brinkman equations numerically, using the finite element method for spatial discretization [8].
Standard (Galerkin) weak form
For the finite element method, we will first derive the Galerkin weak form for the Navier-Stokes-Brinkman equations. We define our trial and weighting function spaces,
$$ u_{i} \in \tau_{i} : \{u_{i} \in H^{1}\left(\Omega\right) \mid u_{i} = g_{i} \text{ on } \Gamma_{g_{i}}\}, $$
$$ w_{i} \in \nu_{i} : \{w_{i} \in H^{1}\left(\Omega\right) \mid w_{i} = 0 \text{ on } \Gamma_{g_{i}}\}, $$
$$ p, q \in Q : \{p \in L^{2}\left(\Omega\right)\}, $$
where $\boldsymbol{w}$ is the weighting function for velocity and $q$ is the weighting function for pressure. We represent these weighting functions discretely on a per-element-basis as
$$ w_{i} = \sum_{a=1}^{n_{en}} N_{a}^{w}w_{ai}, $$
$$ q = \sum_{a=1}^{n_{en}} N_{a}^{q}q_{a}, $$
where $N_{a}^{w}$ and $N_{a}^{q}$ are the nodal shape (basis) functions for the velocity and pressure spaces, respectively, and $w_{ai}$ and $q_{a}$ are the associated arbitrary nodal coefficients. Similarly, the trial functions are represented by
$$ u_{i} = \sum_{a=1}^{n_{en}} N_{a}^{w}u_{ai}, $$
$$ p = \sum_{a=1}^{n_{en}} N_{a}^{q}p_{a}. $$
We then multiply the Navier-Stokes-Brinkman equations by $\boldsymbol{w}$ and $q$, respectively, and integrate by parts to obtain the standard Galerkin momentum and continuity weak forms [7] [5],
$$ \int_{\Omega} \rho w_{i}\frac{du_{i}}{dt} d\Omega + \int_{\Omega} \rho w_{i}u_{k}u_{i, k} d\Omega + \int_{\Omega} w_{i, j}\sigma_{ij} d\Omega + \int_{\Omega} \frac{\mu}{K}w_{i}u_{i} d\Omega - \int_{\Omega} w_{i}\rho b_{i} d\Omega - \int_{\Gamma_{h}} w_{i}h_{i} d\Gamma = 0, $$
$$ \int_{\Omega} qu_{i,i} d\Omega = 0. $$
These two equations can be added together to obtain
$$ \int_{\Omega} qu_{i,i} d\Omega + \int_{\Omega} \rho w_{i}\frac{du_{i}}{dt} d\Omega + \int_{\Omega} \rho w_{i}u_{k}u_{i, k} d\Omega + \int_{\Omega} w_{i, j}\sigma_{ij} d\Omega + \int_{\Omega} \frac{\mu}{K}w_{i}u_{i} d\Omega - \int_{\Omega} w_{i}\rho b_{i} d\Omega - \int_{\Gamma_{h}} w_{i}h_{i} d\Gamma = 0. $$
Stabilized weak form
The standard weak form is generally not stable. Additional terms must be added to stabilize it. We will apply the residual-based variational multiscale (RBVMS / VMS) method for stabilization [1] [5].
In VMS, the velocity and pressure terms are separated into coarse-scale and fine-scale components, such that
$$ \boldsymbol{u} = \boldsymbol{u}^{h} + \boldsymbol{u}', \ $$
$$ p = p^{h} + p', $$
where the $h$-superscript designates the coarse-scale components and the $'$-superscript denotes the fine-scale components. The fine-scale terms are defined as
$$ \boldsymbol{u}' = -\frac{\tau_{SUPS}}{\rho}\boldsymbol{r}_{M}\left(\boldsymbol{u}^{h}, p^{h}\right), $$
$$ p' = -\rho\nu_{LSIC}r_{C}\left(\boldsymbol{u}^{h}\right), $$
where the PDE residuals are
$$ \boldsymbol{r}_{M}\left(\boldsymbol{u}^{h}, p^{h}\right) = \rho\left(\frac{d\boldsymbol{u}^{h}}{dt} + \boldsymbol{u}^{h} \cdot \boldsymbol{\nabla} \boldsymbol{u}^{h} - \boldsymbol{b}\right) - \boldsymbol{\nabla} \cdot \boldsymbol{\sigma}^{h} + \frac{\mu}{K}\boldsymbol{u}^{h}, $$
$$ r_{C}\left(\boldsymbol{u}^{h}\right) = \boldsymbol{\nabla} \cdot \boldsymbol{u}^{h}. $$
The stabilization parameters are defined as
$$ \tau_{SUPS} = \tau_{M} = \left(\frac{4}{\Delta t^{2}} + \boldsymbol{u}^{h} \cdot \boldsymbol{G}\boldsymbol{u}^{h} + C_{1}\nu^{2}\boldsymbol{G}:\boldsymbol{G} + \left(\frac{\nu}{K}\right)^{2}\right)^{-1/2}, $$
$$ \nu_{LSIC} = \tau_{C} = \left(\tau_{SUPS} \text{tr}\boldsymbol{G} \right)^{-1}, $$
where $\boldsymbol{G}$ is the element metric tensor and $\text{tr}\boldsymbol{G}$ is the trace of the metric tensor [1].
Using standard Galerkin momentum and continuity weak forms, and removing the $h$-superscript from the coarse-scale components for notational simplicity (i.e., $\boldsymbol{u}^{h} \rightarrow \boldsymbol{u}$ and $p^{h} \rightarrow p$), we obtain
$$ \int_{\Omega} qu_{i,i} d\Omega + \int_{\Omega} \rho w_{i}\frac{du_{i}}{dt} d\Omega + \int_{\Omega} \rho w_{i}u_{k}u_{i, k} d\Omega + \int_{\Omega} w_{i, j}\sigma_{ij} d\Omega + \int_{\Omega} \frac{\mu}{K}w_{i}u_{i} d\Omega - \int_{\Omega} w_{i}\rho b_{i} d\Omega - \int_{\Gamma_{h}} w_{i}h_{i} d\Gamma + \int_{\Omega} \tau_{SUPS}\left(\frac{q_{,i}}{\rho} + w_{i,k}u_{k}\right)r_{Mi} d\Omega + \int_{\Omega} \rho \nu_{LSIC}r_{C}w_{i,i} d\Omega - \int_{\Omega} w_{i}\tau_{SUPS}r_{Mk}u_{i,k} d\Omega - \int_{\Omega} w_{i,k}\frac{\tau_{SUPS}^{2}}{\rho}r_{Mi}r_{Mk} d\Omega - \int_{\Omega} \frac{\nu}{K}w_{i}\tau_{SUPS}r_{Mi} d\Omega = 0. $$
This is the VMS-stabilized weak form for the Navier-Stokes-Brinkman equations [7] [5] [6]. The first seven terms on the left-hand side correspond to the standard Galerkin weak form. The last five terms are the stabilization terms obtained via VMS. In deriving this equation, we used the continuity equation to obtain $w_{i}u_{k}u_{i,k} = w_{i}\left(u_{k}u_{i}\right)_{,k}$. We also applied the following assumptions [5],
We then add an additional stabilization term, $\int_{\Omega} \frac{\bar{\tau}\tau_{SUPS}^{2}}{\rho} w_{i,k}r_{Mk}r_{Mj}u_{i,j} d\Omega$ [3] [4] [2], to obtain
$$ \int_{\Omega} qu_{i,i} d\Omega + \int_{\Omega} \rho w_{i}\frac{du_{i}}{dt} d\Omega + \int_{\Omega} \rho w_{i}u_{k}u_{i, k} d\Omega + \int_{\Omega} w_{i, j}\sigma_{ij} d\Omega + \int_{\Omega} \frac{\mu}{K}w_{i}u_{i} d\Omega - \int_{\Omega} w_{i}\rho b_{i} d\Omega - \int_{\Gamma_{h}} w_{i}h_{i} d\Gamma + \int_{\Omega} \tau_{SUPS}\left(\frac{q_{,i}}{\rho} + w_{i,k}u_{k}\right)r_{Mi} d\Omega + \int_{\Omega} \rho \nu_{LSIC}r_{C}w_{i,i} d\Omega - \int_{\Omega} w_{i}\tau_{SUPS}r_{Mk}u_{i,k} d\Omega - \int_{\Omega} w_{i,k}\frac{\tau_{SUPS}^{2}}{\rho}r_{Mi}r_{Mk} d\Omega - \int_{\Omega} \frac{\nu}{K}w_{i}\tau_{SUPS}r_{Mi} d\Omega + \int_{\Omega} \frac{\bar{\tau}\tau_{SUPS}^{2}}{\rho} w_{i,k}r_{Mk}r_{Mj}u_{i,j} d\Omega = 0. $$
This equation is the full stabilized weak form used in fluid.cpp in svMultiPhysics.
Residuals
We will temporally discretize the stabilized weak form using the generalized - $\alpha$ method. The resulting nonlinear equation will be linearized and solved iteratively using the Newton-Raphson (Newton) method.
To compute the residuals for each element in the mesh, we separate the stabilized weak form into momentum and continuity components,
$$ \int_{\Omega} \rho w_{i}\frac{du_{i}}{dt} d\Omega + \int_{\Omega} \rho w_{i}u_{k}u_{i, k} d\Omega + \int_{\Omega} w_{i, j}\sigma_{ij} d\Omega + \int_{\Omega} \frac{\mu}{K}w_{i}u_{i} d\Omega - \int_{\Omega} w_{i}\rho b_{i} d\Omega - \int_{\Gamma_{h}} w_{i}h_{i} d\Gamma + \int_{\Omega} \tau_{SUPS}w_{i,k}u_{k}r_{Mi} d\Omega + \int_{\Omega} \rho \nu_{LSIC}r_{C}w_{i,i} d\Omega - \int_{\Omega} w_{i}\tau_{SUPS}r_{Mk}u_{i,k} d\Omega - \int_{\Omega} w_{i,k}\frac{\tau_{SUPS}^{2}}{\rho}r_{Mi}r_{Mk} d\Omega - \int_{\Omega} \frac{\nu}{K}w_{i}\tau_{SUPS}r_{Mi} d\Omega + \int_{\Omega} \frac{\bar{\tau}\tau_{SUPS}^{2}}{\rho} w_{i,k}r_{Mk}r_{Mj}u_{i,j} d\Omega = 0, $$
$$ \int_{\Omega} qu_{i,i} d\Omega + \int_{\Omega} \tau_{SUPS}\frac{q_{,i}}{\rho}r_{Mi} d\Omega = 0. $$
Then, by plugging weighting functions into these equations and holding the results true for any arbitrary $w_{ai}$ and $q_{a}$, we obtain the momentum and continuity residuals,
$$ R_{ai}^{m} = \int_{\Omega} \rho N_{a}^{w}\frac{du_{i}}{dt} d\Omega + \int_{\Omega} \rho N_{a}^{w}u_{k}u_{i, k} d\Omega - \int_{\Omega} pN_{a, i}^{w} d\Omega + \int_{\Omega} N_{a, j}^{w}2\mu\epsilon_{ij} d\Omega + \int_{\Omega} \frac{\mu}{K}N_{a}^{w}u_{i} d\Omega - \int_{\Omega} N_{a}^{w} \rho b_{i} d\Omega + \int_{\Omega} \tau_{SUPS}N_{a, k}^{w}u_{k}r_{Mi} d\Omega + \int_{\Omega} \rho \nu_{LSIC}r_{C}N_{a, i}^{w} d\Omega - \int_{\Omega} N_{a}^{w}\tau_{SUPS}r_{Mk}u_{i,k} d\Omega - \int_{\Omega} N_{a, k}^{w}\frac{\tau_{SUPS}^{2}}{\rho}r_{Mi}r_{Mk} d\Omega - \int_{\Omega} \frac{\nu}{K}\tau_{SUPS}N_{a}^{w}r_{Mi} d\Omega + \int_{\Omega} \frac{\bar{\tau}\tau_{SUPS}^{2}}{\rho} N_{a, k}^{w}r_{Mk}r_{Mj}u_{i,j} d\Omega, $$
$$ R_{a}^{c} = \int_{\Omega} N_{a}^{q}u_{i,i} d\Omega + \int_{\Omega} \tau_{SUPS}\frac{N_{a, i}^{q}}{\rho}r_{Mi} d\Omega, $$
where, for the $a^{\text{th}}$ node in a given element, $R_{ai}^{m}$ is the momentum residual in the $i^{\text{th}}$ direction and $R_{a}^{c}$ is continuity residual. The full residual vector, as used in the generalized - $\alpha$ method, is $\boldsymbol{R} = \left[R_{ai}^{m}, R_{a}^{c}\right]^{T}$. $R_{ai}^{m}$ and $R_{a}^{c}$ are coded in the fluid_2d_m/fluid_3d_m and fluid_2d_c/fluid_3d_c functions, respectively, in fluid.cpp in svMultiPhysics.
Tangent matrices
To compute the elemental tangent matrices, as used in the Newton iterations in the generalized - $\alpha$ method, we plug trial functions into the residuals. We then differentiate the resulting equations with respect to $\frac{du_{n+1}}{dt}$ and $\frac{dp_{n+1}}{dt}$. This yields the tangent matrix,
$$ \boldsymbol{J} = \begin{bmatrix} \boldsymbol{K} & \boldsymbol{G} \ \cr \boldsymbol{D} & \boldsymbol{L} \end{bmatrix} = \begin{bmatrix} \left[K_{ab}^{ij}\right] & \left[G_{ac}^{i}\right] \ \cr \left[D_{ab}^{j}\right] & \left[L_{ac}\right] \end{bmatrix} , $$
where
$$ K_{ab}^{ij} = \frac{\partial R_{ai}^{m}}{\partial \dot{u}_{n+1,bj}}, $$
$$ G_{ac}^{i} = \frac{\partial R_{ai}^{m}}{\partial \dot{p}_{n+1,c}}, $$
$$ D_{ab}^{j} = \frac{\partial R_{a}^{c}}{\partial \dot{u}_{n+1,bj}}, $$
$$ L_{ac} = \frac{\partial R_{a}^{c}}{\partial \dot{p}_{n+1,c}}. $$
In fluid.cpp of svMultiPhysics, the following inconsistent tangent matrices are used,
$$ K_{ab}^{ij} = \alpha_{m} A_{ab}^{ij} + \alpha_{f}\gamma\Delta t B_{ab}^{ij} $$
$$ G_{ac}^{i} = \alpha_{f}\gamma\Delta t \left(-\int_{\Omega} N_{c}^{q}N_{a, i}^{w} d\Omega + \int_{\Omega} \tau_{SUPS} N_{a, g}^{w} u_{g} N_{c, i}^{q} d\Omega - \int_{\Omega} N_{a, k}^{w} \frac{\tau_{SUPS}^{2}}{\rho} N_{c, i}^{q} r_{Mk} d\Omega \right), $$
$$ D_{ab}^{j} = \alpha_{f}\gamma\Delta t \left(\int_{\Omega} N_{a}^{q}N_{b, j}^{w} d\Omega - \int_{\Omega} \tau_{SUPS}\frac{N_{a, i}^{q}}{\rho}\left(-\frac{\alpha_{m}}{\alpha_{f}\gamma\Delta t}\rho N_{b}^{w}\delta_{ij} - \frac{\partial r_{Mi}}{\partial u_{n+\alpha_f,bj}}\right) d\Omega\right), $$
$$ L_{ac} = \alpha_{f}\gamma\Delta t \int_{\Omega} \tau_{SUPS}\frac{N_{a, i}^{q}}{\rho}N_{c, i}^{q} d\Omega, $$
where
$$ A_{ab}^{ij} = \int_{\Omega} \left( \rho N_{a}^{w}N_{b}^{w} \delta_{ij} + \tau_{SUPS} N_{a,g}^{w} u_{g} \rho N_{b}^{w} \delta_{ij} - N_{a,k}^{w} \tau_{SUPS}^{2} N_{b}^{w} \delta_{ij} r_{Mk} \right) d\Omega , $$
$$ B_{ab}^{ij} = \int_{\Omega} \left( \rho N_{a}^{w} u_{k} N_{b, k}^{w} \delta_{ij} + N_{a, l}^{w} \mu N_{b, l}^{w} \delta_{ij} + N_{a, j}^{w} \mu N_{b, i}^{w} + \frac{\mu}{K} N_{a}^{w} N_{b}^{w} \delta_{ij} + \tau_{SUPS} N_{a,g}^{w} u_{g} \frac{\partial r_{Mi}}{\partial u_{n+\alpha_f,bj}} + \rho \nu_{LSIC} N_{b,j}^{w} N_{a,i}^{w} - N_{a}^{w} \tau_{SUPS} N_{b,k}^{w} \delta_{ij} r_{Mk} - N_{a,k}^{w} \frac{\tau_{SUPS}^{2}}{\rho} \frac{\partial r_{Mi}}{\partial u_{n+\alpha_f,bj}} r_{Mk} + \frac{4}{\gamma} \frac{\partial \mu}{\partial \gamma} \epsilon_{jk} N_{b,k}^{w} \epsilon_{il} N_{a,l}^{w} + \frac{\bar{\tau}\tau_{SUPS}^{2}}{\rho} N_{a,k}^{w} N_{b,z}^{w} r_{Mk} r_{Mz} \delta_{ij} \right) d\Omega , $$
and
$$ \frac{\partial r_{Mi}}{\partial u_{n+\alpha_f,bj}} = \left(\rho u_{k} N_{b,k}^{w} - \mu N_{b,kk}^{w} + \frac{\mu}{K} N_{b}^{w} - \frac{\partial \mu}{\partial x_{k}} N_{b,k}^{w} \right)\delta_{ij} - \frac{2}{\gamma} \frac{\partial \mu}{\partial \gamma} \epsilon_{il} N_{b,l}^{w} u_{j, kk} - \frac{\partial \mu}{\partial x_{j}} N_{b,i}^{w}. $$
These inconsistent tangent matrices were derived by using these assumptions:
$K_{ab}^{ij}$ and $G_{ac}^{i}$ are coded in the fluid_2d_m/fluid_3d_m functions, while $D_{ab}^{j}$ and $L_{ac}$ are coded in the fluid_2d_c/fluid_3d_c functions.
References
[1] Y. Bazilevs, V.M. Calo, J.A. Cottrell, T.J.R. Hughes, A. Reali, and G. Scovazzi (2007). Variational multiscale residual-based turbulence modeling for large eddy simulation of incompressible flows. Computer Methods in Applied Mechanics and Engineering, 197, 173-201, https://doi.org/10.1016/j.cma.2007.07.016
[2] C.A. Taylor, T.J.R. Hughes, and C.K. Zarins (1998). Finite element modeling of blood flow in arteries. Computer Methods in Applied Mechanics and Engineering, 158, 155-196, https://doi.org/10.1016/S0045-7825(98)80008-X
[3] M. Esmaily-Moghadam (2014). Development of multiscale modeling methods for clinical decision making in single ventricle heart patients. University of California, San Diego
[4] C.H. Whiting and K.E. Jansen (2001). A stabilized finite element method for the incompressible Navier–Stokes equations using a hierarchical basis. International Journal for Numerical Methods in Fluids, 35, 93-116, https://doi.org/10.1002/1097-0363(20010115)35:1<93::AID-FLD85>3.0.CO;2-G
[5] Y. Bazilevs, K. Takizawa, T.E. Tezduyar (2013). Computational Fluid–Structure Interaction: Methods and Applications. John Wiley & Sons, Ltd.
[6] F.M. Gerosa and A.L. Marsden (2024). A mechanically consistent unified formulation for fluid-porous-structure-contact interaction. Computer Methods in Applied Mechanics and Engineering, 425, 116942, https://doi.org/10.1016/j.cma.2024.116942
[7] J. Fuchsberger, P. Aigner, S. Niederer, G. Plank, H. Schima, G. Haase, and E. Karabelas (2022). On the incorporation of obstacles in a fluid flow problem using a Navier–Stokes–Brinkman penalization approach. Journal of Computational Science, 57, 101506, https://doi.org/10.1016/j.jocs.2021.101506
[8] T.J.R. Hughes (2003). CThe Finite Element Method: Linear Static and Dynamic Finite Element Analysis. Dover Publications
Solver Parameter Input File
The svMultiPhysics solver reads simulation parameters from a text file in Extensible Markup Language (XML) format. The XML file structurally organizes the svMultiPhysics solver parameters and input data (e.g. mesh files) needed to set up and execute a simulation.
What is XML
The XML format is made up of tags, elements and attributes. A tag begins with < and ends with >. There are two types of tag
- start-tag, such as <section>
- end-tag , such as </section>
An XML element is everything from (including) the element's start tag to (including) the element's end tag.
<Tolerance> 0.001 </Tolerance>
There are three types of elements
- text
- attributes - name-value pair within a start-tag: name="value" "
- other elements - elements nested within other elements
The following XML is an example of these different element types
<Add_mesh name="fluid"><Mesh_file_path> fluid_mesh.vtu </Mesh_file_path>
<Add_face name="inlet">
<Face_file_path> inlet.vtp </Face_file_path>
</Add_face>
<Domain> 1 </Domain>
<Mesh_scale_factor> 0.1 </Mesh_scale_facto>
</Add_mesh>
The above contains
- six elements: Add_mesh, Mesh_file_path, Add_face, Face_file_path, Domain and Mesh_scale_factor
- <Add_mesh name="fluid"> has an attribute name with value "fluid"
- <Add_mesh name="fluid"> element provides a context for the other elements under it: <Add_face name=inlet> associates the face named inlet to the mesh named fluid
- Mesh_file_path is a text element with text fluid_mesh.vtu
- Domain is a text element with text 1
- Mesh_scale_factor is a text element with text 0.1
How svMultiPhysics uses XML
For the svMultiPhysics solver input file the XML elements are used as the names of the solver parameters.
svMultiPhysics defines for each parameter
- Name - Case sensitive with the first letter capitalized
- Data type
- boolean - on, off
- integer - numeric with no decimal
- string - a contiguous list of characters
- real - numeric with decimal or scientific notation
- Context - parameters can only be found within the context of another parameter
If a value for a parameter is not valid svMultiPhysics will display an error message indicating where in the file the error occurred.
Some parameters are optional. There are two types of optional parameters
- Used if not set - a default value is used
- Not used if not set - no default value
The organization of the parameter input file
The parameter input file is organized into six sections used to provide context for the parameters defined under them
- <GeneralSimulationParameters> - General Simulation Parameters
- <Add_mesh> - Mesh Parameters
- <Add_equation> - Equation Section
- <LS> - Linear Solver Parameters
- <Output> - Output Parameters
<?xml version="1.0" encoding="UTF-8" ?>
<svMultiPhysicsFile version="0.1">
and end with the single line
</svMultiPhysicsFile>
where
- <?xml - identifies the file as an XML format file
- <svMultiPhysicsFile - identifies the file as an svMultiPhysics XML format file
The following outlines the basic structure of the parameter input XML file.
<GeneralSimulationParameters> [general parameters] </GeneralSimulationParameters><Add_mesh name=mesh_name_1> <Add_face name=mesh1_face_name_1> <Add_face name=mesh1_face_name_2> ... <Add_face name=mesh1_face_name_n> </Add_mesh> <Add_mesh name=mesh_name_2> <Add_face name=mesh2_face_name_1> <Add_face name=mesh2_face_name_2> ... <Add_face name=mesh2_face_name_n> </Add_mesh> ... <Add_mesh name=mesh_name_N> <Add_face name=meshN_face_name_1> <Add_face name=meshN_face_name_2> ... <Add_face name=meshN_face_name_M> </Add_mesh><Add_equation type=eq_1> [equation parameters] <Domain id=0> <Domain id=1> ... <Domain id=N> <LS type=ls_type> [linear solver parameters] </LS> <Output type=output_type> [solver output parameters] </Output> </Add_equation>
Parameter Documentation Conventions
Parameters are documented using the parameter Name and DataType. If the parameter is optional its default value is shown between square bracketsExample The optional Spectral_radius_of_infinite_time_step parameter with default value 0.5
<Spectral_radius_of_infinite_time_step> real [0.5] </Spectral_radius_of_infinite_time_step>
A string parameter will have a default value of none. If not set the parameter will have a value of an empty string ("" in C++) and will be treated as being unset and will not be used.
Example The optional Save_results_in_folder string parameter with default value none
<Save_results_in_folder> string [none] </Save_results_in_folder>
A required parameter will have no default value.
Example The required Number_of_time_steps parameter with no default value
<Number_of_time_steps> integer </Number_of_time_steps>
Solver Parameter Input File Examples
The svMultiPhysics GitHub repository Test Cases directory contains solver input XML files used to test svMultiPhysics. These XML files can be used as templates that can be customized for a particular application.General Simulation Parameters
The General Simulation Parameters section contains the basic parameters used to setup the simulation- spatial dimension
- time step control
- results output
Parameters
If true then restart the simulation with a restart file generated by svMultiPhysics.
<Number_of_spatial_dimensions> integer
The number of spatial dimensions the simulation is using. Acceptable values: 2 or 3
<Number_of_time_steps> integer
The number of time steps to run the simulation.
<Time_step_size> real
The value use to increment the simulation time.
<Starting_time_step> int [0]
The starting time step used for the simulation.
<Spectral_radius_of_infinite_time_step> real [0.5]
The spectral radius is used to compute parameters for the generalized alpha method. A value of 0.0 leads to an over-damped system while 1.0 leads to an undamped system.
If true then save simulation results to VTK-format files.
The name prefixed to a the file name when saving simulation results to VTK-format files.
Indicates how often to save simulation results to a VTK-format file.
Save simulation results in the specified directory.
If true then all available restart files will be converted to VTK format files.
<Start_saving_after_time_step> integer [1]
Save simulation results after the given time step.
<Save_averaged_results> boolean [false]
If true then compute time-averaged results from the entire simulation. The averaging is performed after the last time step of the simulation using VTU files written during the simulation.
<Start_averaging_from_zero> boolean [false]
If true then start everging results from time zero.
Indicates how often to save restart files.
Name of the file used to store simulation state data. This file is used for restarting simulations.
If true then the simulation results restart file is overwritten.
Initialize all state variables from a VTK VTU file.
Mesh Section
A Mesh Section of the solver parameters input file is primarily used to identify the finite element mesh files used in a simulation. The Add_mesh keyword defines the volumetric and surfaces meshes used in the simulation.A Mesh Section is organized as follows
<Add_face name=mesh_face_1_name>
[Face 1 parameters]
<Add_face>
<Add_face name=mesh_face_2_name>
[Face 2 parameters]
<Add_face>
...
<Add_face name=mesh_face_N_name>
[Face N parameters]
<Add_face>
</Add_mesh >
The <Add_mesh name=mesh_name> parameter adds a volumetric mesh named mesh_name to the simulation. Multiple <Add_mesh name=mesh_name> parameters can be given within a solver parameters input file.
The <Add_face name=mesh_face_name> parameters subsection under the <Add_mesh> parameter associates the surface mesh named mesh_face_name with the volumetric mesh. A surface mesh, referred to as a face, is used to apply boundary conditions.
Mesh Parameters
The path to a VTK VTU file defining a finite element volume mesh in the simulation.
<Mesh_domain> integer
The integer ID used to identify a domain in the simulation.
Creates IDs for multiple domains from a VTK VTU file. Domain IDs are assumed to be stored in a VTK Cell data array named DOMAIN_ID.
<Mesh_scale_factor> real [1.0]
The value used to scale mesh nodal coordinates.
The name of the VTK VTU file used to initialize simulation pressure values on a fluid domain. See Initial Condition Data.
The name of the VTK VTU file used to initialize simulation velocity values on a fluid domain. See Initial Condition Data.
The name of the VTK VTU file used to initialize simulation displacement values on a solid domain. See Initial Condition Data.
The name of the VTK VTU file used to initialize simulation stress values on a solid domain. See Initial Condition Data.
The name of the VTK VTU file used to define the fiber directions for a solid domain. See Fiber and Sheet Direction Data.
<Set_mesh_as_fibers> boolean [false]
If true then the mesh given in the Mesh_file_path keywork is treated as a mesh of one-dimensional elements in 3D. The mesh typically defines a network of Purkinje fibers used in cardiac electrophysiology simulations.
<Set_mesh_as_shell> boolean [false]
If true then the mesh given in the Mesh_file_path keywork is treated as a mesh of shell elements. The mesh can then be used for solving the shell equation type. A shell mesh is also needed for initializing a CMM equation using a prestress or displacements.
<quadrature_modifier_TET4> real [0.5854]
The quadrature used for tetrahedral elements.
Add_Face Parameters
<Face_file_path> string
The path to a VTK VTP file defining a finite element surface mesh for a face in the simulation.
<Quadrature_modifier_TRI3> real [0.6666]
The quadrature used for triangle elements.
</Add_face>
Contact Section
The Contact Section of the solver parameters input file is used to define the parameters controlling the contact computations that prevents the geometric interpenetration of element domains.A Contact Section is organized as follows
The Contact keyword is used to define the parameters used in a contact computation. The model attribute contact_model is the name of a method used for the contact computation. The "penalty" model is currently the only model supported.
Contact Parameters
The maximum depth of penetration between elements in contact.
<Desired_separation> real [0.05]
The minimum depth of penetration between elements in contact.
<Min_norm_of_face_normals> real [0.05]
The minimum norm of face normals in contact.
<Penalty_constant> real [1e5]
The penalty constant used to prevent element interpenetration.
Precomputed Solution Section
The Precomputed Solution Section of the solver parameters input file is used to define the parameters used to read in the data from a precomputed solution for the simulation state.A Precomputed Solution Section is organized as follows
The Precomputed_solution keyword is used to define the parameters to for reading in a precomputed solution data and using it in a simulation.
Precomputed Solution Parameters
If true then use a precomputed solution for a simulation.
<Field_name> string [none]
The name of the precomputed solution data field in the VTK VTU file.
<File_path> string [none]
The path to the VTK VTU file containing the precomputed solution data.
<Time_step> real [0.05]
The time step used to compute the precomputed solution data.
Projection Section
The Projection Section of the solver parameters input file is used to couple the interface between the fluid and the solid domains for a fluid-structure interaction (FSI) simulation. The shared boundary surfaces between the fluid and solid must have the same node IDs or coordinates.A Projection Section is organized as follows
The Projection keyword is used to create a projection. The name attribute project_from_face_name is the name of a face defined in the Mesh Section.
Projection Parameters
The name of a face defined in the Mesh Section that is to be coupled to the project_from_face_name face.
<Projection_tolerance> real [0.0]
Face node distances within this tolerance will be treated as equal.
Equation Section
The Equation Section of the solver parameters input file defines the properties of an equation- physics - fluid, structure, electrophysiology, etc.
- domains
The Equation Section is organized as follows
[ Equation Parameters ]
<Domain id=domain_id>
[ Domain Subsection ]
</Domain>
<ECGLeads>
[ ECG Leads Subsection ]
</ECGLeads>
<svZeroDSolver_interface>
[ svZeroDSolver Interface Subsection ]
</svZeroDSolver_interface>
<Variable_wall_properties mesh_name=mesh_name>
[ Variable wall properties Subsection ]
<Viscosity>
<Viscosity model=viscosity_model>
[ Viscosity Subsection ]
<Viscosity>
<LS type=linear_solver_type>
[ Linear Solver Subsection ]
<LS>
<Output type=output_type>
[ Output Subsection ]
<Output>
<svZeroDSolver_interface name=svzerod_interface_name>
[ svZeroDSolver Interface Subsection ]
<svZeroDSolver_interface>
<Add_BC name=bc_name>
[ Boundary Condition Subsection ]
<Add_BC>
<Add_equation>
The Add_equation keyword adds an equation of type equation_type to the simulation. Multiple Add_equation keywords can be given within a solver parameters input file.
The value of equation_type is a string identifying the equation type
- "advection_diffusion" - unsteady advection-diffusion
- "cardiac_electro_physiology" - solves the mono-domain model of cardiac electrophysiology
- "coupled_momentum" - coupled momentum method for fluid-structure interaction
- "fluid" - Navier-Stokes equations for unsteady viscous incompressible fluid flow
- "fluid-solid-interaction" - fluid-solid interaction using the arbitrary Lagrangian-Eulerian (ALE) formulation
- "linear_elasticity" - linear elastodynamics equation
- "mesh" - solves a modified linear elasticity equation for mesh motion in an fluid-solid interaction simulation
- "shell" - nonlinear thin shell mechanics (Kirchhoff-Love theory)
- "solid_heat" - unsteady diffusion equation
- "stokes" - unsteady Stokes equations
- "structural" - nonlinear elastodynamics equation
- "structural_velocity_pressure" - nonlinear elastodynamics equation using mixed VMS-stabilized formulation
Equation Parameters
If true then the convergence of a multi-equation system of equations is coupled: nonlinear iterations are performed within each time step on all the coupled system of equations until convergence is achieved. If false then the convergence of each equation is uncoupled and is achieved separately within each time step.
<Initialize> string [none]
Initialize the fluid wall for a CMM equation. Permissible values are ·inflate - Initialize the fluid wall using displacements read in from a file given in the Initial_displacements_file_path parameter. ·prestress - Initialize the fluid wall using traction values read in from a file given in the Prestress_file_path parameter.
<Initialize_RCR_from_flow> boolean [false]
If true then initialize RCR from flow data.
<Max_iterations> integer [1]
The maximum number of iterations used to solve a nonlinear system of equations.
<Min_iterations> integer [1]
The minimum number of iterations used to solve a nonlinear system of equations.
<Prestress> boolean [false]
If true then enable prestressing a solid domain to mimic physiological conditions. Valid for linear elastic and structural equations only
<Tolerance> real [0.5]
The solution of a nonlinear system of equations is considered to be converged (solved) if the nonlinear residual is less than this value.
<Use_taylor_hood_type_basis> boolean [false]
If true then use a Taylor-Hood element pair for increased stability.
ECG Leads Subsection
The ECG Leads Subsection of the Equation Section or Domain Subsection contains the parameters used to define ECG leads.The ECG Leads Subsection is organized as follows
The ECGLeads keyword defines a subsection for ECG Leads parameters.
ECG Leads Parameters
The path to the text file containing the x-coordinates for the ECG leads.
<Y_coords_file_path> string [none]
The path to the text file containing the y-coordinates for the ECG leads.
<Z_coords_file_path> string [none]
The path to the text file containing the z-coordinates for the ECG leads.
Variable wall properties Subsection
The Variable wall properties Subsection of the Equation Section or Domain Subsection defines the parameters used to define the node-based material properties used in a CMM equation.The Variable wall properties Subsection is organized as follows
[ Variable wall properties Parameters ]
</Variable_wall_properties>
The Variable_wall_properties keyword defines a subsection for variable wall properties parameters. The mesh_name attribute is the mesh name to assign the variable wall properties.
Variable wall properties Parameters
The path to the VTK VTP-format file used to set the node-wise elasticity modulus and shell thickness properties.
Viscosity Subsection
The Viscosity Subsection of the Equation Section or Domain Subsection defines the parameters for the fluid viscosity model used by fluid, CMM, and stokes equations, or the solid viscosity model used by struct and ustruct.The Viscosity Subsection is organized as follows
[ Viscosity Parameters ]
</Viscosity>
The Viscosity keyword defines a subsection for viscosity parameters.
For fluid, CMM, and stokes equations, the value of viscosity_model can be
- "newtonian" - Newtonian viscosity model
- "carreau-yasuda" - Carreau-Yasuda viscosity model
- "cassons" - Cassons viscosity model
Newtonian
<Value> real
The value of the viscosity constant.
</Viscosity>
Carreau-Yasuda
<Limiting_high_shear_rate_viscosity> real
The value of the limiting high shear rate viscosity parameter.
<Limiting_low_shear_rate_viscosity> real
The value of Limiting low shear rate viscosity parameter.
<Power_law_index> real
The value of the power law index parameter.
<Shear_rate_tensor_multiplier> real
The value of the shear rate tensor multiplier parameter.
<Shear_rate_tensor_exponent> real
The value of the shear rate tensor exponent parameter.
</Viscosity>
Cassons
<Asymptotic_viscosity_parameter> real
The value of the asymptotic viscosity parameter.
<Low_shear_rate_threshold> real
The value of the low shear rate threshold parameter.
<Yield_stress_parameter> real
The value of the yield stress parameter.
</Viscosity>
For struct and ustruct equations, the value of viscosity_model can be
- "Newtonian" - Newtonian viscosity model
- "Potential" - Pseudo-potential viscosity model
Newtonian
<Value> real
The value of the viscosity constant.
</Viscosity>
Potential
<Value> real
The value of the viscosity constant.
</Viscosity>
Domain Subsection
The Domain Subsection of the Equation Section defines the physical properties for a region of a finite element mesh used in solving the enclosing equation.The Domain Subsection is organized as follows
[ Domain Parameters ]
<Viscosity model=viscosity_model>
[ Viscosity Subsection ]
<Viscosity>
<Stimulus type=stimulus_type>
[ Stimulus Subsection ]
<Stimulus>
<Domain>
The Domain keyword adds a domain to the enclosing equation. domain_id is an integer ID defined in Mesh Section for a finite element mesh.
Domain Parameters by Equation
The following sections list the parameters associated with a domain for a specific equation.Advection Diffusion Domain Parameters
Conductivity Source_term
Cardiac Electrophysiology Domain Parameters
Anisotropic_conductivity Electrophysiology_model Myocardial_zone ODE_solver Stimulus Time_step_for_integration
Coupled Momentum Domain Parameters
Elasticity_modulus Fluid_density Solid_density Backflow_stabilization_coefficient Poisson_ratio Shell_thickness
Fluid Domain Parameters
Backflow_stabilization_coefficient Density Force_x Force_y Force_z Inverse_darcy_permeability
Fluid-Solid Interaction Domain Parameters
Fluid-solid interaction parameters will be a combination of fluid and linear elasticity, structural or structural velocity pressure parameters.Linear Elasticity Domain Parameters
Density Elasticity_modulus Force_x Force_y Force_z Poisson_ratio
Mesh Domain Parameters
Poisson_ratio
Shell Domain Parameters
Density Elasticity_modulus Force_x Force_y Force_z Poisson_ratio
Solid Heat Domain Parameters
Conductivity Density Source_term
Stokes Domain Parameters
Force_x Force_y Force_z Momentum_stabilization_coefficient
Structural Domain Parameters
Density Elasticity_modulus Force_x Force_y Force_z Poisson_ratio
Structural Velocity Pressure Domain Parameters
Density Elasticity_modulus Force_x Force_y Force_z Momentum_stabilization_coefficient Poisson_ratio
Domain Parameters
The following section lists all of the parameters associated with a domain. However, certain parameters only make sense when used in a domain defined for a specific type of equation.Convergence is deemed to occur when the size of the residual becomes sufficiently small compared with this value.
The anisotropic conductivity constant used to model the directional dependence of conduction velocity in cardiac tissue.
A parameter used to control flow reversal at outlet boundaries.
The thermal or diffusive conductivity within a domain.
<Continuity_stabilization_coefficient> real [0.0]
The stabilizing coefficient associated to the continuity equation in the Variational Multiscale Method.
The density property for a solid or fluid.
<Dilational_penalty_model> string [none]
The dilational or volume-preserving component of the hyperelastic strain energy function. The Shell equation uses ST91 dilational penalty model only Permissible values are
- Quadratic
- Simo-Taylor91
- Miehe94
The elastic modulus property for a solid or fluid.
The activation model used to simulation electrical activity in heart muscle. Permissible values are
- aliev-panfilov - Aliev-Panfilov model
- bueno-orovio - Bueno-Orovio model
- fitzhugh-nagumo - Fitzhugh-Nagumo model
- tentusscher-panfilov - tenTusscher-Panfilov model
<Feedback_parameter_for_stretch_activated_currents> real [0.5]
The constant for stretch-activated currents (SACs) used in cardiac electrophysiology simulations.
The density property for a fluid.
The x-component of a body force applied to a domain.
The y-component of a body force applied to a domain.
The z-component of a body force applied to a domain.
<G_Na> real [14.838]
The maximal sodium conductance used in the TenTusscher-Panfilov Ventricular Myocyte Model.
<G_CaL> real [3.98E-5]
The maximal calcium conductance used in the TenTusscher-Panfilov Ventricular Myocyte Model.
<G_Kr> real [0.153]
The maximal rapid delayed rectifier potassium current used in the TenTusscher-Panfilov Ventricular Myocyte Model.
<G_Ks> real [0.392]
The maximal rapid delayed rectifier sodium current used in the TenTusscher-Panfilov Ventricular Myocyte Model.
<G_to> real [0.294]
The maximal transient outward current conductance used in the TenTusscher-Panfilov Ventricular Myocyte Model.
Sets the inverse darcy permeability constant.
<Isotropic_conductivity> real [0.0]
Sets the isotropic conductivity constant.
<Mass_damping> real [0.0]
The mass damping coefficient used to dissipate energy.
The stabilizing coefficient associated to the momentum equation in the Variational Multiscale Method.
The region of the heart the domain is representing. This is needed for some activation models. Permissible values are
- endocardium
- epicardium
- myocardium
The time integration method to solve the ODEs used to integrate activation models.
- euler - first-order Euler method
- implicit - second-order Crank–Nicolson method
- runge - fourth-order Runge-Kutte method
<Penalty_parameter> real [0.0]
The penalty constant used to override the physical bulk modulus. This is needed if the Poisson ratio for a solid is close to 0.5.
The Poisson ratio property for a solid. Permissible values are between 0.0 and 0.5.
<Relative_tolerance> real [1e-4]
Convergence is deemed to occur when the size of the residual becomes sufficiently small compared with the first error estimate calculated.
The thickness of the shell elements in a domain.
The density property for a solid.
The volumetric source term (e.g. heat or dye) within a domain.
<tau_fi> real [0.110]
The fast inward current time scale used in the Bueno-Orovio cell activation model.
<tau_si> real [1.88750]
The slow inward current time scale used in the Bueno-Orovio cell activation model.
The time step size used to integrate activation models if different from simulation time step.
Stimulus Subsection
The Stimulus Subsection of the Domain Subsection defines the stimulus settings for pacemaker cells for a cardiac electrophysiology simulation.The Stimulus Subsection is organized as follows
The Stimulus keyword adds stimulus settings to the enclosing domain. The stimulus_type attribute can be
- "Istim" - the applied stimulus is a source current and not a voltage source
Parameters
The stimulus amplitude.
<Cycle_length> real [0.0]
The stimulus cycle length.
<Duration> real [0.0]
The stimulus duration.
<Start_time> real [0.0]
The stimulus start time.
Linear Solver Subsection
The Linear Solver Subsection of the Equation Section defines the linear solver to use to solve the linear system generated by FEM and several parameters used to control the solution process.The Linear Solver Subsection is organized as follows
<Linear_algebra type=linear_algebra_package>
[ Linear Algebra Parameters ]
<Linear_algebra>
[Linear Solver Parameters]
<LS>
The LS keyword adds a linear solver to the enclosing equation.
The value of linear_solver_type can be
- "BICG" - Biconjugate gradient method. This solver is used for the same equations used for GMRES.
- "CG" - Conjugate gradient. This solver is used for the following equations
- linear_elasticity
- mesh
- shell
- solid_heat
- structural
- "GMRES" - Generalized minimal residual method. This solver is used for the following equations
- cardiac_electro_physiology
- coupled_momentum
- fluid-solid-interaction
- advection_diffusion
- stokes
- structural_velocity_pressure
- "NS" - Navier-Stokes solver based on bipartition method for fluid equations.
Linear Solver Parameters
The maximum number of iterations to use for a linear solve.
<Krylov_space_dimension> integer [200]
The dimension of the Kyrlov matrix used to approximate a linear system.
<Tolerance> real [0.4]
The solution of a linear system of equations is considered to be converged (solved) if the linear residual is less than this value.
<Absolute_tolerance> real [1e-10]
The solution of a near system of equations is considered to be converged (solved) if the linear residual is less than this value.
<NS_GM_max_iterations> integer [5]
The number of maximum iterations used to the control the convergence of the GMRES portion of the bi-partitioned iterative algorithm.
<NS_GM_tolerance> real [1e-10]
The tolerance used to control the convergence of the GMRES portion of the bi-partitioned iterative algorithm.
<NS_CG_max_iterations> integer [5]
The number of maximum iterations used to the control the convergence of the conjugate gradient portion of the bi-partitioned iterative algorithm.
<NS_CG_tolerance> real [1e-10]
The tolerance used to the control the convergence of the conjugate gradient portion of the bi-partitioned iterative algorithm.
Linear Algebra Parameters
The Linear Algebra Parameters section of the solver parameters input file defines the linear algebra package to use to solve linear systems using the <Linear_algebra type=linear_algebra_type> parameter.The value of linear_algebra_type can be
- "fsils" - Custom numerical linear algebra routines included in the svMultiPhysics implementation
- "petsc" - Portable, Extensible Toolkit for Scientific Computation (PETSc)
- "trilinos" - The Trilinos numerical linear algebra package
The preconditioner applied to a linear system.
The value of the preconditioner can be
- fsils
- fsils
- rcs
- petsc
- jacobi
- rcs
- trilinos
- trilinos-diagonal
- trilinos-blockjacobi
- trilinos-ilu
- trilinos-ilut
- trilinos-ic
- trilinos-ict
- trilinos-ml
svZeroDSolver Interface Subsection
The svZeroDSolver Interface Subsection of the Equation Section defines the parameters needed by the svMultiPhysics solver to interface with the svZeroDSolver.The SimVascular svZeroDSolver simulates bulk cardiovascular flow rates and pressures using an arbitrary zero-dimensional (0D) lumped parameter model (LPM) of a discrete network of components analogous to electrical circuits. It provides an Application Programming Interface (API) that allows it to communicate and interact with external software applications directly using function calls to programmatically define custom inflow and outflow boundary conditions for a CFD simulation. The svMultiPhysics solver can directly access the svZeroDSolver API by loading the svZeroDSolver as a shared (dynamic) library available after installing the svZeroDSolver.
The svZeroDSolver Interface Subsection is organized as follows
The svZeroDSolver_interface keyword activates coupling for boundary conditions with the <Time_dependence> Coupled </Time_dependence> keyword for the enclosing equation coupled boundary condition.
The value of svzerodsolver_interface_name can be any string name and is not currently used.
svZeroDSolver Interface Parameters
The type of coupling used between the svMultiPhysics and svZeroDSolve solvers. The coupling type determines how contributions from the coupled svZeroDSolver are added to the svMultiPhysics solver tangent matrix.
The value of the coupling type can be
- explicit - No contribution is added. This may potentially cause poor nonlinear convergence.
- implicit - Contributions are added as outlet resistances for each svZeroDSolver and svMultiPhysics solver Newton iteration. This leads to theoretically quadratic nonlinear convergence but at the expense of additional computation.
- semi-implicit - Contributions obtained from the svZeroDSolver for outlet resistances are added only at the beginning of the simulation .
The name of the svZeroDSolver JSON configuration file containing the description of the components of a lumped parameter network and various simulation parameters.
<Shared_library> string [none]
The path and name of the svZeroDSolver shared library used to access the svZeroDSolver API. The name of the library is specific to the computer platform you are using
- Linux - libsvzero_interface.so
- MacOS - libsvzero_interface.dylib
- Windows - libsvzero_interface.dll
<Initial_flows> real [0.0]
Initialize the flows for all coupling blocks with the given value. If this parameter is not given then flows will not be initialized.
<Initial_pressures> real [0.0]
Initialize the pressures for all coupling blocks with the given value. If this parameter is not given then pressures will not be initialized.
Output Subsection
The Output Subsection section of the Equation Section defines the quantities written to simulation results files.The Output Subsection is organized as a list of quantity names with a boolean
[Quantity names]
<Output>
The Output keyword enables output of quantities for the enclosing equation. The value of output_type can be
- "Boundary_integral" - Quantities computed as a flux through boundary surfaces (faces) (e.g. velocity flux, energy flux)
- "Spatial" - Quantities written to VTK VTU files
- "Volume_integral" - Volume-averaged quantities
Quantity names are boolean parameters used to enable quantities for output by setting their value to true. The quantities available for output depends on the equation being solved. Below is a list of all quantity names available for output for the corresponding equation
- advection_diffusion - Temperature, Heat_flux
- cardiac_electro_physiology - Action_potential
- coupled_momentum - Velocity, Pressure, Displacement,WSS, Vorticity, Vortex, Energy_flux, Strain_invariants, Acceleration, Traction, Viscosity, Divergence
- fluid - Velocity, Pressure, WSS, Traction, Vorticity, Vortex, Energy_flux, Strain_invariants, Acceleration, Viscosity, Divergence
- fluid-solid-interaction - Velocity, Pressure, Displacement, VonMises_stress, WSS, Traction, Vorticity, Vortex, Energy_flux, Strain_invariants, Viscosity, Absolute_velocity, Stress, Strain, Cauchy_stress, Jacobian, Deg_grad, Area/Volume, Fiber_direction, Fiber_alignment, Divergence, Acceleration
- linear_elasticity - Displacement, VonMises_stress, Stress, Strain, Jacobian, Area/Volume, Velocity, Acceleration
- mesh - Displacement, Velocity, Acceleration
- shell - Displacement, Stress, Strain (Green-Lagrange in-plane), Def_grad, Area, Jacobian, Velocity, 1st invariant of Cauchy-Green strain (in-plane), Cauchy-Green strain tensor
- solid_heat - Temperature, Heat_flux
- stokes - Velocity, Pressure, WSS, Vorticity, Traction, Strain_invariants, Viscosity, Divergence
- structural - Displacement, VonMises_stress, Stress, Cauchy_stress, Strain, Jacobian, Def_grad, Area/Volume, Fiber_direction, Fiber_alignment, Velocity, Acceleration
- structural_velocity_pressure - Displacement, VonMises_stress, Stress, Cauchy_stress, Strain, Jacobian, Def_grad, Area/Volume, Fiber_direction, Fiber_alignment, Velocity, Pressure, Acceleration, Divergence
Example
Boundary Condition Parameters
The Boundary Condition Parameters section of the solver parameters input file adds boundary conditions to a simulation. A boundary condition is set on a mesh face (surface) using a face name added to the simulation using the <Add_face parameter.The Boundary Condition Parameters section is organized as follows
[ General Boundary Condition Parameters ]
<RCR_values>
[ RCR Parameters ]
<RCR_values>
<Add_BC>
The <Add_BC name=face_name> parameter adds a boundary condition to the enclosing equation.
Boundary condition types
Coupled momentum boundary condition
The coupled momentum boundary condition identifies the face to be treated using coupled momentum method.Dirichlet boundary condition
The Dirichlet boundary condition sets the values of a state variable on a face. A state variable refers to the unknown of the equation that is being solved for. For example, velocity is the state variable for a fluid equation.The following lists the state variable set for each type of equation
- velocity - stokes, fluid, structural_velocity_pressure, coupled_momentum, fluid-solid-interaction
- displacement - linear_elasticity, structural, shell, mesh
- temperature - advection_diffusion, solid_heat
- action_potential - cardiac_electro_physiology
Neumann boundary condition
The Neumann boundary condition imposes a normal force on the face.Robin boundary condition
The Robin boundary condition applies a force to a face. The force is represented as a spring-mass-damper system.Traction boundary condition
The Traction boundary condition applies a force on a face. The force can be along any direction.Temporal and spatial values for boundary conditions
A temporal distribution specifies the time-dependent values of a state variable. The values are stored as the value of a state variable integrated over a face stored in a text file.
A temporal and spatial distribution specifies time-dependent values for each node in a face. The values are stored in a text file.
General Parameters
Load a file used to set the spatial and temporal values for a boundary condition. The file can be a VTK VTP-format file or a text file.
<CST_shell_bc_type> string [none]
Constrain the degrees of freedom of a Constant Strain Triangle (CST) of a triangular shell element. condition. Permissible values are:
- fixed - Fix the edge
- hinged - Allow rotation about an edge
- free - Allow free motion about an edge
<Impose_flux> boolean [false]
If true then normalize the spatial profile with the area of the face so that the imposed flux value is converted into the state variable.
<Impose_on_state_variable_integral> boolean [false]
If true then used for applying a Dirichlet boundary condition on the displacement degrees of freedom when velocity is the state variable (e.g. fluid, CMM, FSI).
Use the VTK VTU-format file to initialize CMM using an inflation method resulting from a diastolic or time-averaged fluid traction.
<Penalty_parameter> real [0.0]
If the Poisson ratio for a given case is close to 0.5, then calculated bulk modulus used for dilational penalty model can be extremely high leading to poor linear solver convergence. The users may then override the physical bulk modulus with a penalty constant sufficiently large enough for the linear solver to converge.
Use the VTK VTP-format file to initialize CMM using a prestressed wall under equilibrium with fluid traction.
Set the spatial distribution of a state variable on the face. Acceptable values: Flat, Parabolic or User_defined
<Ramp_function> boolean [false]
If true then the first two entries in the file setting an unsteady boundary is used to linearly increment from the first value to the second value, and maintains a steady value thereafter.
Use the given text file to set the spatial distribution of a state variable for a User_defined profile.
The path to the VTK VTP-format file used to set the spatially varying body force to a face.
The name of the svZeroDSolver block this boundary face is coupled to.
The path to the text file containing temporal and spatial values.
The path to the text file containing temporal values.
<Time_dependence> string [steady]
The time dependence of a boundary condition. Permissible values are: ·general - Spatial and temporal variations are provided from a file ·spatial - Spatially varying values ·steady - A constant value is imposed ·unsteady - Time-dependent values are provide from a file
<Traction_values_file_path> string [none]
The path to the VTK VTP-format file containing nodally varying traction values.
<Traction_multiplier> real [1.0]
The value used to scale the traction values read from from a file.
<Type> string
The boundary condition type. Permissible values are: ·Coupled Momentum - Identifies the face to be treated using the coupled momentum method ·Dirichlet - Identifies the face to be treated as a Dirichlet boundary condition ·Neumann - Identifies the face to be treated as a Neumann boundary condition ·Robin - Identifies the face to be treated as a Robin boundary condition ·Traction - Identifies the face to be treated as a Traction boundary condition
<Value> real
The value of the state variable.
RCR Boundary Condition Parameters
<Capacitance> real
Capacitance.
<Distal_resistance> real
Distal resistance.
<Distal_pressure> real
The distal pressure used to initialize an RCR boundary condition.
<Initial_pressure> real
The initial pressure used to initialize an RCR boundary condition.
<Proximal_resistance> real
Proximal resistance.
</RCR_values>
Dirichlet Boundary Condition Parameters
Enforce a Dirichlet boundary condition along a given Cartesian coordinate 3-vector. The vector has the form (x,y,z).
<Penalty_parameter_normal> real [0.0]
If the Dirichlet boundary condition is weakly applied then it is applied in a normal direction.
<Penalty_parameter_tangential> real [0.0]
If the Dirichlet boundary condition is weakly applied then it is applied in a tangential direction.
<Weakly_applied> boolean [false]
If true then the Dirichlet boundary condition is applied weakly using augmented Lagrange-multiplier formulation. This setting is applied to fluid/FSI equations only.
<Zero_out_perimeter> boolean [true]
If true then the solver will zero out the nodes shared by two adjacent faces.
Neumann Boundary Condition Parameters
The distal pressure used to initialize an RCR boundary condition.
<Follower_pressure_load> boolean [false]
If true then the applied load follows the mesh deformation. This implies that the magnitude of the load is proportional to the surface area during the deformation.
A text file containing the Fourier coefficients interpolating a time-dependent state variable.
<Undeforming_neu_face> boolean [false]
If true then mimic clamped condition on a specimen routinely done in experiments. Clamping will not allow the surface, on which the load is applied, to deform. Used only for nonlinear elastodynamics.
Robin Boundary Condition Parameters
If true then when the Robin BC is applied along the normal direction the applied surface traction takes the form sigma.n = -(k u.n + c v.n - p)n
<Damping>> real [1.0]
The damping constant for the spring-mass-damper force used to implement a Robin boundary condition.
<Stiffness>> real [1.0]
The stiffness constant for the spring-mass-damper force used to implement a Robin boundary condition.
Boundary Condition Data
The data used for a boundary condition defines the spatial and/or the temporal distribution of values over a boundary surface. Data can be of five types- Fourier coefficients
- Spatial distribution
- Temporal distribution
- Temporal and spatial distribution
- User-defined flow profile
Data is stored using two file formats
Fourier coefficients Data
The data file contains the Fourier coefficients interpolating temporal distribution data for a boundary face. This data will be directly used to represent periodic data at intermediate time values.This file is set using the solver input file Fourier_coefficients_file_path keyword.
The file format is
InitialTime Period InitialValue_1 InitialValue_2 ... InitialValue_NDIM TimeDerivativeLinearPart1 TimeDerivativeLinearPart2 ... TimeDerivativeLinearPart_NDIM NumberOfFourierModes RealPart_1 ImaginaryPart_1 ...
where NDIM is the number of space dimensions used in the simulation.
Spatial Distribution Data
The data files for the spatial distribution of data are used for- Neumann boundary condition
- Traction boundary condition
- Coupled Momentum boundary condition
Neumann and Traction boundary conditions
For Neumann and Traction boundary conditions the data file is a VTK VTP format file containing a PointData array named- Pressure - Neumann boundary condition
- Traction - Traction boundary condition
A spatial data file is set using the solver input file Spatial_values_file keyword.
Coupled momentum boundary condition
For a coupled momentum boundary condition the data file is a VTK VTU format file used to initialized stresses (prestress) and displacements using PointData arrays named- Stress - Initialize stresses
- Displacement - Initialize displacements
These spatial data files are set using the solver input file Initial_displacements_file_path and Prestress_file_path and keywords.
Temporal Distribution Data
The data file contains a single time-varying value per time step for a boundary surface. The values (e.g. flow rate) will be distributed to the boundary surface nodes and interpolated in time using a given number of Fourier modes.Temporal data will be interpolated internally using Fourier series interpolation to obtain periodic data at intermediate time values.
This file is set using the solver input file Temporal_values_file_path keyword.
The file format is
NumberOfTimePoints NumberOfFourierModes Time_1 Value_1 Time_2 Value_2 ... Time_NumberOfTimePoints Value_NumberOfTimePoints
Example: Inlet fluid flow rate
Temporal and Spatial Distribution Data
The data file contains time-varying values per time step for each node on a boundary surface. The values are interpolated in time using a given number of Fourier modes.This file is set using the solver input file Temporal_and_spatial_values_file_path keyword.
The file format is
NumberDegreesOfFreedom NumberOfTimePoints NumberOfNodes Time_1 Time_2 ... Time_NumberOfTimePoints NodeID_1 Time_1 Value_1,1 Value_1,2 ... Value_1,NumberDegreesOfFreedom Time_2 Value_2,1 Value_2,2 ... Value_2,NumberDegreesOfFreedom .... NodeID_2 Time_1 Valu_1 Value_2 ... Value_NumberDegreesOfFreedom ... NodeID_NumberOfNodes ...
Example: Volumetric body force
User Profile Data
The data file contains values for each node on a boundary surface defining the shape of the profile for a boundary surface.This file is set using the solver input file Spatial_profile_file_path keyword.
The file format is
NodeID_1 Value_1 NodeID_2 Value_2 ... NodeID_NumberOfFaceNodes Value_NumberOfFaceNodes
where NumberOfFaceNodes is the number of nodes in the boundary surface.
Fiber and Sheet Direction Data
A fiber and sheet direction data file is used to define a local orthonormal coordinate system for each element in a finite element volume mesh representing cardiac muscle. The data file is a VTK VTU format file containing a CellData array of three-component vecors named FIB_DIR.The data file is set using the solver input file Fiber_direction_file_path keyword.
Initial Condition Data
Initial condition data for a finite element volume mesh can be initialized from values contained in a VTK VTU format file containing PointData arrays.Displacement, Pressure, Stress, and Velocity Data
Displacement, pressure, stress, and velocity data can be initialized using VTK VTU format file containing PointData arrays named Displacement, Pressure, Stress, and Velocity.The data files are set using the solver input file Initial_displacements_file_path, Initial_pressures_file_path and Prestress_file_path and Initial_velocities_file_path keywords.
State Variables
All state variables for a simulation can be initialized using VTK VTU format file containing the appropriate PointData arrays.The data file is set using the solver input file Simulation_initialization_file_path keyword.
Finite Element Mesh Data
Finite element mesh data is stored in the Visualization Toolkit (VTK) file format. Finite element mesh files are typically created using the SimVascular SV Meshing Tool .Finite Element Volume Mesh
A finite element volume mesh is stored in a VTK VTU-format file using the UnstructuredGrid data type. The file contains the following mesh data- number of nodes - NumberOfPoints
- nodal coordinates - DataArray type="Float32" Name="Points"
- number of elements - NumberOfCells
- element connectivity - DataArray type="Int64" Name="connectivity"
- user-defined arrays - CellData (element) or PointData (node) data specific to a SimVascular application
- GlobalElementID - CellData
- GlobalNodeID - PointData
- ModelRegionID - CellData
A finite element volume mesh file is set using the solver input file Mesh_file_path keyword.
Example: A finite element volume mesh comprising 8561 nodes and 45723 elements
<?xml version="1.0"?> <VTKFile type="UnstructuredGrid" version="0.1" byte_order="LittleEndian" header_type="UInt32" compressor="vtkZLibDataCompressor"> <UnstructuredGrid> <Piece NumberOfPoints="8561" NumberOfCells="45723" > <PointData> <DataArray type="Int32" Name="GlobalNodeID" format="appended" RangeMin="1" RangeMax="8561" offset="0" /> </PointData> <CellData Scalars="ModelRegionID"> <DataArray type="Int32" Name="ModelRegionID" format="appended" RangeMin="1" RangeMax="1" offset="15964" /> <DataArray type="Int32" Name="GlobalElementID" format="appended" RangeMin="1" RangeMax="45723" offset="16444" /> </CellData> <Points> <DataArray type="Float32" Name="Points" NumberOfComponents="3" format="appended" RangeMin="0.12837002719" RangeMax="30.066622254" offset="100980" /> </Points> <Cells> <DataArray type="Int64" Name="connectivity" format="appended" RangeMin="" RangeMax="" offset="227180" /> <DataArray type="Int64" Name="offsets" format="appended" RangeMin="" RangeMax="" offset="788628" /> <DataArray type="UInt8" Name="types" format="appended" RangeMin="" RangeMax="" offset="868152" /> </Cells> </Piece> </UnstructuredGrid> [COMPRESSED DATA]
The nodes and elements in the finite element mesh are identified by a unique ID in the GlobalNodeID and GlobalElementID arrays, respectively.
Finite Element Boundary Surface Mesh
A finite element boundary surface mesh is stored in a VTK VTP-format file using the PolyData data type. The file contains the following mesh data- number of nodes - NumberOfPoints
- nodal coordinates - DataArray type="Float32" Name="Points"
- number of face elements - NumberOfCells
- element face connectivity - DataArray type="Int64" Name="connectivity"
- user-defined arrays - CellData (element) or PointData (node) data specific to a SimVascular application
- GlobalElementID - CellData
- GlobalNodeID - PointData
- ModelFaceID - CellData
A finite element boundary surface mesh file is set using the solver input file Face_file_path keyword.
Example: A finite element boundary surface mesh comprising 2647 nodes and 5290 elements
<?xml version="1.0"?> <VTKFile type="PolyData" version="0.1" byte_order="LittleEndian" header_type="UInt32" compressor="vtkZLibDataCompressor"> <PolyData> <Piece NumberOfPoints="2647" NumberOfVerts="0" NumberOfLines="0" NumberOfStrips="0" NumberOfPolys="5290" > <PointData Scalars="GlobalNodeID"> <DataArray type="Int32" Name="GlobalNodeID" format="appended" RangeMin="1" RangeMax="2647" offset="0" /> </PointData> <CellData Scalars="ModelFaceID"> <DataArray type="Int32" Name="GlobalElementID" format="appended" RangeMin="2" RangeMax="45702" offset="4964" /> <DataArray type="Int32" Name="ModelFaceID" format="appended" RangeMin="1" RangeMax="3" offset="23712" /> </CellData> <Points> <DataArray type="Float32" Name="Points" NumberOfComponents="3" format="appended" RangeMin="0.12837002719" RangeMax="30.066622254" offset="23932" > </DataArray> </Points> <Verts> <DataArray type="Int64" Name="connectivity" format="appended" RangeMin="" RangeMax="" offset="62272" /> <DataArray type="Int64" Name="offsets" format="appended" RangeMin="" RangeMax="" offset="62288" /> </Verts> <Lines> <DataArray type="Int64" Name="connectivity" format="appended" RangeMin="" RangeMax="" offset="62304" /> <DataArray type="Int64" Name="offsets" format="appended" RangeMin="" RangeMax="" offset="62320" /> </Lines> <Strips> <DataArray type="Int64" Name="connectivity" format="appended" RangeMin="" RangeMax="" offset="62336" /> <DataArray type="Int64" Name="offsets" format="appended" RangeMin="" RangeMax="" offset="62352" /> </Strips> <Polys> <DataArray type="Int64" Name="connectivity" format="appended" RangeMin="" RangeMax="" offset="62368" /> <DataArray type="Int64" Name="offsets" format="appended" RangeMin="" RangeMax="" offset="99592" /> </Polys> </Piece> </PolyData> <AppendedData encoding="base64"> [COMPRESSED DATA]
The node IDs contained in the GlobalNodeID array identify the nodes stored in the boundary surface VTP with the nodes in the finite element volume mesh.
The element IDs contained in the GlobalElementID array identify the element in the finite element volume mesh that each polygon of the face belongs to.
Results Data
Simulation results are written to two types of filesSee the Simulation Results Output section for a description of the solver input file keywords used to control how results files are written.
Binary Results Restart File Format
The binary results restart file stores state variables and mesh data for a single time step in a custom binary format. A binary file is a series of sequential bytes, each of which is eight bits in length. The content of this file is must be accessed by a program (e.g. svMultiPhysics) that knows how that content is formatted and how to read the data programmatically using C++ or Python.Each restart file contains a header of seven integer (4-byte) values containing
- Number of processors used to run the simulation
- Number of equations
- Number of meshes
- Number of nodes
- Number of unknowns in the coupled 0D-3D problem
- Number of degrees of freedom
- An error flag
These values are followed by an integer (4-byte) and two double (8-byte) values containing
- Simulation time step
- Simulation time
- Wall clock time
Header values are checked when a simulation is restarted to ensure that the restart is valid for the the current simulation setup (e.g. same number of processors, same number of equations, etc.).
State variable data (e.g. velocity) is then written after the header. The results from all processors in a parallel simulation are written to the restart file..
VTK Results File Format
The data file is a VTK VTU format file. State variables and any output quantities given in the solver input file Output Subsection are stored using PointData arrays.A separate VTK file is used to store results for a single time step. Results can be visualized using the ParaView visualization program.
Data from all processors in a parallel simulation are written to the VTK files.
SimVascular MultiPhysics Tool
The SimVascular SV MultiPhysics Tool is used to interactively create the mesh and solver XML input files needed to run an svMultiPhysics simulation.
The MultiPhysics Tool stores data in the Project/svMultiPhysics/Data Node directory. The MultiPhysics Tool stores values input via the GUI controls into the ProjectMultiPhysics/Data Node.multiphysicsjob file.
The following sections describe how to use the SimVascular MultiPhysics Tool GUI to set the mesh, boundary conditions, material properties and solver parameters needed for a simulation.
Creating Mesh Files
The MultiPhysics Tool requires files containing the finite element volume mesh (VTK VTU) and boundary surface meshes (VTK VTP). These files are created by by right-clicking on the Meshes node under the SV Data Manager. Selecting the Export Mesh Complete menu item displays a file browser used to select to location to store the location for the finite element files.Creating an MultiPhysics Tool Instance
A MultiPhysics Tool instance is created by right-clicking on the MultiPhysics node under the SV Data Manager. Selecting the Create MultiPhysics job menu item displays a popup window.
Type fluid_simulation into the job Name DialogBox. This names is used to identify the MultiPhysics Tool instance and to name the files and directories stored under the SimVascular project's MultiPhysics directory. Selecting OK creates an MultiPhysics Tool instance node under the SV Data Manager. Selecting this instance displays the MultiPhysics Tool panel.

The panel contains four buttons at the top used to primarily to create new svMultiPhysics simulations (jobs)
New job Button - Selecting to create a new MultiPhysics instance
Save job Button - Select to write the values of GUI controls into the SVPROJECT/MultiPhysicsJob Name.multiphysicsjob file.
Load job Button - Select to display a file browser used to select a .multiphysicsjob file. The values of GUI controls are then set from the values in this file.
Load mesh - Reads in the finite element mesh stored in the Project/Meshes directory. At startup SimVascular does not read in the finite element mesh so this must be initiated manually.
The panel contains four sub-panels used to input a specific category of data needed to create the files needed to run a svMultiPhysics solver simulation
- Domains
- Physics
- Simulation Parameters
- Run Simulation
A selecting a sub-panel name brings up the sub-panel's controls. The following sections describe how each of the sub-panels are used.
Physics Panel
The Physics panel is used to set the equation(s), its physical properties, boundary conditions and linear solver solution strategy needed to numerically solve it in a simulation.Panel Layout

Usage
Add or remove equations - This GUI control is used to add equations to the simulation. Selecting the equation name in the left panel and selecting the > Button moves it to the right panel and adds that equation to the simulation. Use the < Button to remove that equation from the simulation.
When an equation is added the Buttons at the bottom of this control (Properties, Output, etc.) are used to set the parameters for the selected equation. The GUI controls for setting parameters for each of these are displayed directly beneath them.
Properties Button - Select to display the GUI controls used to set the physicals properties (e.g. density) needed by the equation.
Output Button - Select to display the GUI controls to add/remove names of the quantities for output.
BCs Button - Select to display the GUI controls to add/remove names of the quantities for output.
Advanced Button - Select to display a GUI panel to set advanced solver options.
Linear Solver Button - Select to display a GUI panel to set linear solver options.
Remesher Button - Select to display a GUI panel to set remeshing options.
GUI controls for setting parameters
Properties Button - Setting physical properties
The physical properties displayed in the panel depend of the selected equation.Fluid equation properties

Viscosity TextBox - Fluid viscosity
Linear elasticity equation properties

Elastic modulus TextBox - Elastic modulus
Poisson ratio TextBox - Poisson ratio. Value <= 0.5
Structure equation properties

Elastic modulus TextBox - Elastic modulus
Poisson ratio TextBox - Poisson ratio. Value <= 0.5
Constitutive model ComboBox - Constitutive model selected from
- stVK - St. Venant-Kirchhoff material model
- nHK - NeoHookean material model
FSI equation properties

Viscosity TextBox - Fluid viscosity
Solid Density TextBox - Solid density
Poisson ratio TextBox - Poisson ratio. Value <= 0.5
Constitutive model ComboBox - Constitutive model selected from
- stVK - St. Venant-Kirchhoff material model
- nHK - NeoHookean material model
Output Button - Adding output quantities
The quantities output by default are displayed in the right sub-panel.Selecting the output quantity name in the left panel and selecting the > Button moves it to the right panel and adds that quantity for output.

BCs Button - Adding Boundary Conditions


Boundary Condition Setup Parameters
Face list: ListTable - Lists the boundary surface face names defined for a domain.
Selecting a face name adds a boundary condition for that face.
BC Type: ComboBox - Select boundary condition type from the list [ Neumann Dirichlet]
Directional BC CheckBox - If selected
Steady CheckBox - If selected
Value TextBox -
Unsteady CheckBox - If selected
Temporal values file FileBrowser - Select to identify the time values
Value TextBox -
Coupled - CheckBox - If selected
Temporal and spatial values file FileBrowser - Select to identify the time values
Projection - CheckBox - If selected
Zero out perimeter - CheckBox - If selected
Impose flux - CheckBox - If selected
Flat - CheckBox - If selected
Parabolic - CheckBox - If selected
User defined - CheckBox - If selected
Spatial profile file path FileBrowser - Select to identify the time values
Impose on state variable integral - CheckBox - If selected
Effective direction TextBox -
Reset Button -
Cancel Button -
OK Button -
Advanced Button - Setting advanced solver parameters

Coupled - CheckBox - If selected
Tolerance TextBox -
Residual dB reduction SpinBox -
Min iterations SpinBox -
Max iterations SpinBox -
Backflow stabilization SpinBox -
Reset Button -
Linear solverButton - Setting linear solver parameters

Coupled - CheckBox - If selected
Type ComboBox - Select the linear solver to use.
Max iterations TextBox -
Tolerance TextBox -
Krylov dimension TextBox -
Absolute Tolerance TextBox -
Preconditioner ComboBox - Select the linear solver preconditioner to use.
Simulation Parameters Panel
The Simulation Parameters panel is used to set general solver parameters.Panel Layout

Usage
Time control
Start from previous simulation CheckBox - If selected then start the simulation from the saved state of a previous simulation. Number of time steps TextBox - The number of time steps to run the simulation. Time step size TextBox - The simulation time step.
Restart files
Name TextBox - The name prefix of the files used to store simulation state data. Save increment TextBox - Indicates how often to save restart files. Save results after time step TextBox - Start saving simulation results after this time.
Save as VTK CheckBox - If selected then save simulation output to VTK files. Name TextBox - The name prefix of the VTK files used to store output data. Save increment TextBox - Indicates how often to save the VTK files.
Time-averaged results CheckBox - If selected then compute time-averaged results for
the entire simulation using the using VTU files written during the simulation.
Advanced options CheckBox - If selected then enable setting advanced parameters. Spectral radius of infinite time step TextBox - The spectral radius is used to compute parameters for the generalized alpha method. A value of 0.0 leads to an over-damped system while 1.0 leads to an undamped system. Remeshing CheckBox - If selected then the finite element mesh for FSI simulations will be automatically remeshed during a simulation. Warning CheckBox - If selected then print out warning messages indicating unexpected inputs or simulation behavior. Verbose CheckBox - If selected print detailed messages when processing input. Debug CheckBox - If selected the print additional information that may be used for debugging purposes.
Run Simulation Panel
The Run Simulation panel is used to create the solver input files needed to run a simulation.Panel Layout

Usage
Create Input File Button - Select to create solver input files.
Number of processes Slider - Select the number of processors (cores) for MPI to use for a parallel simulation.
Run Button - Select to start a simulation.
Stop Button - Select to stop a simulation.
Fluid-Structure Interaction Modeling using svMultiPhysics
Cardiovascular applications that involve significant deformation are difficult to simulate and model. Some examples include vascular simulations where the wall dilates to more than 10% of the vessel diameter, cardiac simulations where the heart wall displaces as it contracts and expands, or coronary simulations where the coronary vessels displace with the motion of the heart.
One way to computationally model situations like these is to use a fluid-structure interaction (FSI) solver. In FSI, separate domains are defined for the fluid part and solid parts of the computational geometry. The respective equations governing fluid flow (typically Navier-Stokes) are solved in the fluid domain, while the equations governing solid mechanics are solved in the solid domain. The two domains then interact through their interface where the solution variables (displacements, velocities, pressures, stresses) are required to match. This interface acts as a coupled boundary condition for both domains, as the solution of one domain will affect the solution in the other, and vice-versa.
svMultiPhysics utilizes Arbitrary Lagrangian-Eulerian (ALE) method to perform FSI. As the name implies, this method is a combination of the Eulerian and Lagrangian descriptions of motion that is particularly well suited for FSI problems. In the Eulerian description, the computational mesh stays fixed and the motion is characterized as velocities flowing past the grid nodes. The Eulerian description is typical for most fluid dynamics problems without FSI. While in the Lagrangian description, the mesh nodes move exactly with the motion of the fluid or solid and thus are often characterized as moving domain problems. Lagrangian problems are also often described in terms of the reference configuration (the initial computational geometry before any motion) and the current configuration (the current state of the geometry and mesh nodes).ALE combines these descriptions into a formulation that is convenient for describing FSI problems, where the fluid mesh motion does not follow the motion of the fluid, but instead the mesh motion is determined by some moving scheme, which usually includes the solution of equations of elasticity. More information can be found in the literature.
References
[1] Ju Liu, Alison L. Marsden. A unified continuum and variational multiscale formulation for fluids, solids, and fluid–structure interaction. Computer Methods in Applied Mechanics and Engineering; 1 August 2018.
[2] Jie Cheng and Lucy T. Zhang. A General Approach to Derive Stress and Elasticity Tensors for Hyperelastic Isotropic and Anisotropic Biomaterials. International Journal of Computational Methods; Vol. 15, No. 04, 1850028 (2018).
[3] T. Christian Gasser, Ray W Ogden and Gerhard A Holzapfel.Hyperelastic modelling of arterial layers with distributed collagen fibre orientations. Royal Society. 28 September 2005.
[4] J. M. Guccione, A. D. McCulloch, L. K. Waldman. Passive Material Properties of Intact Ventricular Myocardium Determined From a Cylindrical Model.J Biomech Eng. Feb 1991, 113(1): 42-55.
[5] Gerhard A. Holzapfel and Ray W. Ogden. Constitutive modelling of passive myocardium: a structurally based framework for material characterization. Royal Society. 13 September 2009.
[6] Lei Shi, Ian Y. Chen, Hiroo Takayama, Vijay Vedula. An optimization framework to personalize passive cardiac mechanics.Computer Methods in Applied Mechanics and Engineering. 1 December 2024.