svMultiPhysics - Beta Version

The svMultiPhysics documentation is in a Beta phase of release.

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

  • Solid mechanics 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.
  • Fluid mechanics equations
    • Navier-Stokes equation - Solves the incompressible viscous fluid flow problem.
    • Unsteady Stokes flow - Solves for unsteady flow where inertia is negligible (creeping flow).
  • Diffusion equations
    • 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 fluid-solid equations
    • 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.
  • Cardiac electrophysiology
    • 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

    1. Model geometry - changes in geometry (e.g. curvature) in a region
    2. 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
    Numerical linear algebra packages typically provide a number of preconditioners used to accelerate the solution of a linear system.

    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.

  • Time derivative: $\frac{d\boldsymbol{a}}{dt} = \dot{\boldsymbol{a}}$
  • Dot product: $\boldsymbol{a} \cdot \boldsymbol{b} = a_{i}b_{i}$
  • Gradient of vector: $\boldsymbol{\nabla} \boldsymbol{a} = \frac{da_{i}}{dx_{j}} = a_{i,j}$
  • Tensor-vector product: $\boldsymbol{A}\boldsymbol{a} = A_{ij}a_{j}$
  • Double dot product: $\boldsymbol{A}:\boldsymbol{B} = A_{ij}B_{ij}$
  • Tensor (dyadic) product: $\boldsymbol{C} = C_{ij} = \boldsymbol{a} \otimes \boldsymbol{b} = a_{i}b_{j}$
  • Divergence of a tensor: $\boldsymbol{\nabla} \cdot \boldsymbol{\sigma} = \sigma_{ij,j}$
  • Gradient of vector double dotted with a tensor: $\boldsymbol{\nabla} \boldsymbol{a} : \left(\boldsymbol{b} \otimes \boldsymbol{c}\right) = a_{i,j}b_{i}c_{j}$
  • 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],

  • $\frac{du'}{dt} = 0$,
  • $u' = 0$ on $\Gamma_{g}$ and $\Gamma_{h}$,
  • $\nabla^{s}\boldsymbol{w}:2\mu\nabla^{s}\boldsymbol{u}' = 0$.
  • 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:

  • convective velocities (the $\boldsymbol{u}$ in $\boldsymbol{u} \cdot \boldsymbol{\nabla} u_{i}$) are constant,
  • stabilization parameters, $\tau_{SUPS}$ and $\nu_{LSIC}$, are constant.
  • $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

    1. <GeneralSimulationParameters> - General Simulation Parameters
    2. <Add_mesh> - Mesh Parameters
    3. <Add_equation> - Equation Section
    4. <LS> - Linear Solver Parameters
    5. <Output> - Output Parameters
    The parameter input XML file must start with the following two lines
    <?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 brackets

    Example 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

    <Continue_previous_simulation> boolean [false] </Continue_previous_simulation>
    If true then restart the simulation with a restart file generated by svMultiPhysics.
    <Number_of_spatial_dimensions> integer </Number_of_spatial_dimensions>
    The number of spatial dimensions the simulation is using. Acceptable values: 2 or 3
    <Number_of_time_steps> integer </Number_of_time_steps>
    The number of time steps to run the simulation.
    <Time_step_size> real </Time_step_size>
    The value use to increment the simulation time.
    <Starting_time_step> int [0] </Starting_time_step>
    The starting time step used for the simulation.
    <Spectral_radius_of_infinite_time_step> real [0.5] </Spectral_radius_of_infinite_time_step>
    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.
    <Save_results_to_VTK_format> boolean [false] </Save_results_to_VTK_format>
    If true then save simulation results to VTK-format files.
    <Name_prefix_of_saved_VTK_files> string [result] </Name_prefix_of_saved_VTK_files>
    The name prefixed to a the file name when saving simulation results to VTK-format files.
    <Increment_in_saving_VTK_files> integer [1] </Increment_in_saving_VTK_files>
    Indicates how often to save simulation results to a VTK-format file.
    <Save_results_in_folder> string [none] </Save_results_in_folder>
    Save simulation results in the specified directory.
    <Convert_BIN_to_VTK_format> boolean [false] </Convert_BIN_to_VTK_format>
    If true then all available restart files will be converted to VTK format files.
    <Start_saving_after_time_step> integer [1] </Start_saving_after_time_step>
    Save simulation results after the given time step.
    <Save_averaged_results> boolean [false] </Save_averaged_results>
    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] </Start_averaging_from_zero>
    If true then start everging results from time zero.
    <Increment_in_saving_restart_files> integer [10] </Increment_in_saving_restart_files>
    Indicates how often to save restart files.
    <Restart_file_name> string [stFile.bin] </Restart_file_name>
    Name of the file used to store simulation state data. This file is used for restarting simulations.
    <Overwrite_restart_file> boolean [false] </Overwrite_restart_file>
    If true then the simulation results restart file is overwritten.
    <Simulation_initialization_file_path> string [none] </Simulation_initialization_file_path>
    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_mesh name=mesh_name>

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

    [ Mesh Parameters ]

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

    <Mesh_file_path> string </Mesh_file_path>
    The path to a VTK VTU file defining a finite element volume mesh in the simulation.
    <Mesh_domain> integer </Mesh_domain>
    The integer ID used to identify a domain in the simulation.
    <Domain_file_path> string </Domain_file_path>
    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] </Mesh_scale_factor>
    The value used to scale mesh nodal coordinates.
    <Initial_pressures_file_path> string [none] </Initial_pressures_file_path>
    The name of the VTK VTU file used to initialize simulation pressure values on a fluid domain. See Initial Condition Data.
    <Initial_velocities_file_path> string [none] </Initial_velocities_file_path>
    The name of the VTK VTU file used to initialize simulation velocity values on a fluid domain. See Initial Condition Data.
    <Initial_displacements_file_path> string [none] </Initial_displacements_file_path>
    The name of the VTK VTU file used to initialize simulation displacement values on a solid domain. See Initial Condition Data.
    <Prestress_file_path> string [none] </Prestress_file_path>
    The name of the VTK VTU file used to initialize simulation stress values on a solid domain. See Initial Condition Data.
    <Fiber_direction_file_path> string [none] </Fiber_direction_file_path>
    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] </Set_mesh_as_fibers>
    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] </Set_mesh_as_shell>
    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] </quadrature_modifier_TET4>
    The quadrature used for tetrahedral elements.

    Add_Face Parameters

    <Add_face>
    <Face_file_path> string </Face_file_path>
    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] </Quadrature_modifier_TRI3>
    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

    <Contact model=contact_model>
    [ Contact Parameters ]
    </Contact>

    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

    <Closest_gap_to_activate_penalty> real [1.0] </Closest_gap_to_activate_penalty>
    The maximum depth of penetration between elements in contact.
    <Desired_separation> real [0.05] </Desired_separation>
    The minimum depth of penetration between elements in contact.
    <Min_norm_of_face_normals> real [0.05] </Min_norm_of_face_normals>
    The minimum norm of face normals in contact.
    <Penalty_constant> real [1e5] </Penalty_constant>
    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

    <Precomputed_solution>
    [ Precomputed Solution Parameters ]
    </Precomputed_solution>

    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

    <Use_precomputed_solution> boolean [false] </Use_precomputed_solution>
    If true then use a precomputed solution for a simulation.
    <Field_name> string [none] </Field_name>
    The name of the precomputed solution data field in the VTK VTU file.
    <File_path> string [none] </File_path>
    The path to the VTK VTU file containing the precomputed solution data.
    <Time_step> real [0.05] </Time_step>
    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

    <Projection name=project_to_face_name>
    [ Projection Parameters ]
    </Projection>

    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

    <Project_from_face> string </Project_from_face>
    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] </Projection_tolerance>
    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

    <Add_equation type=equation_type>

    [ 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

    <Coupled> boolean [true] </Coupled>
    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>
    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] </Initialize_RCR_from_flow>
    If true then initialize RCR from flow data.
    <Max_iterations> integer [1] </Max_iterations>
    The maximum number of iterations used to solve a nonlinear system of equations.
    <Min_iterations> integer [1] </Min_iterations>
    The minimum number of iterations used to solve a nonlinear system of equations.
    <Prestress> boolean [false] </Prestress>
    If true then enable prestressing a solid domain to mimic physiological conditions. Valid for linear elastic and structural equations only
    <Tolerance> real [0.5] </Tolerance>
    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] </Use_taylor_hood_type_basis>
    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

    <ECGLeads>

    [ ECG Leads Parameters ]

    </ECGLeads>

    The ECGLeads keyword defines a subsection for ECG Leads parameters.

    ECG Leads Parameters
    <X_coords_file_path> string [none] </X_coords_file_path>
    The path to the text file containing the x-coordinates for the ECG leads.
    <Y_coords_file_path> string [none] </Y_coords_file_path>
    The path to the text file containing the y-coordinates for the ECG leads.
    <Z_coords_file_path> string [none] </Z_coords_file_path>
    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 mesh_name=mesh_name>

    [ 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
    <Wall_properties_file_path> string [none] </Wall_properties_file_path>
    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 model=viscosity_model>

    [ 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
    <Viscosity model="newtonian">
    <Value> real </Value>
    The value of the viscosity constant.
    </Viscosity>
    Carreau-Yasuda
    <Viscosity model="carreau-yasuda">
    <Limiting_high_shear_rate_viscosity> real </Limiting_high_shear_rate_viscosity>
    The value of the limiting high shear rate viscosity parameter.
    <Limiting_low_shear_rate_viscosity> real </Limiting_low_shear_rate_viscosity>
    The value of Limiting low shear rate viscosity parameter.
    <Power_law_index> real </Power_law_index>
    The value of the power law index parameter.
    <Shear_rate_tensor_multiplier> real </Shear_rate_tensor_multiplier>
    The value of the shear rate tensor multiplier parameter.
    <Shear_rate_tensor_exponent> real </Shear_rate_tensor_exponent>
    The value of the shear rate tensor exponent parameter.
    </Viscosity>
    Cassons
    <Viscosity model="cassons">
    <Asymptotic_viscosity_parameter> real </Asymptotic_viscosity_parameter>
    The value of the asymptotic viscosity parameter.
    <Low_shear_rate_threshold> real </Low_shear_rate_threshold>
    The value of the low shear rate threshold parameter.
    <Yield_stress_parameter> real </Yield_stress_parameter>
    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
    <Viscosity model="Newtonian">
    <Value> real </Value>
    The value of the viscosity constant.
    </Viscosity>
    Potential
    <Viscosity model="Potential">
    <Value> real </Value>
    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 id=domain_id>

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

    If the simulation comprises a single domain with uniform properties then the domain parameters don't need to be defined within a Domain Subsection, they can just be given in the enclosing equation.

    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.
    <Absolute_tolerance> real [1e-6] </Absolute_tolerance>
    Convergence is deemed to occur when the size of the residual becomes sufficiently small compared with this value.
    <Anisotropic_conductivity> vector </Anisotropic_conductivity>
    The anisotropic conductivity constant used to model the directional dependence of conduction velocity in cardiac tissue.
    <Backflow_stabilization_coefficient> real [0.2] </Backflow_stabilization_coefficient>
    A parameter used to control flow reversal at outlet boundaries.
    <Conductivity> real [0.0] </Conductivity>
    The thermal or diffusive conductivity within a domain.
    <Continuity_stabilization_coefficient> real [0.0] </Continuity_stabilization_coefficient>
    The stabilizing coefficient associated to the continuity equation in the Variational Multiscale Method.
    <Density> real [0.5] </Density>
    The density property for a solid or fluid.
    <Dilational_penalty_model> string [none] </Dilational_penalty_mode>
    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

    <Elasticity_modulus> real [1e7] </Elasticity_modulus>
    The elastic modulus property for a solid or fluid.
    <Electrophysiology_model> string [none] </Electrophysiology_model>
    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] </Feedback_parameter_for_stretch_activated_currents>
    The constant for stretch-activated currents (SACs) used in cardiac electrophysiology simulations.
    <Fluid_density> real [0.5] </Fluid_density>
    The density property for a fluid.
    <Force_x> real [0.0] </Force_x>
    The x-component of a body force applied to a domain.
    <Force_y> real [0.0] </Force_y>
    The y-component of a body force applied to a domain.
    <Force_z> real [0.0] </Force_z>
    The z-component of a body force applied to a domain.
    <G_Na> real [14.838] </G_Na>
    The maximal sodium conductance used in the TenTusscher-Panfilov Ventricular Myocyte Model.
    <G_CaL> real [3.98E-5] </G_CaL>
    The maximal calcium conductance used in the TenTusscher-Panfilov Ventricular Myocyte Model.
    <G_Kr> real [0.153] </G_Kr>
    The maximal rapid delayed rectifier potassium current used in the TenTusscher-Panfilov Ventricular Myocyte Model.
    <G_Ks> real [0.392] </G_Ks>
    The maximal rapid delayed rectifier sodium current used in the TenTusscher-Panfilov Ventricular Myocyte Model.
    <G_to> real [0.294] </G_to>
    The maximal transient outward current conductance used in the TenTusscher-Panfilov Ventricular Myocyte Model.
    <Inverse_darcy_permeability> real [0.0] </Inverse_darcy_permeability>
    Sets the inverse darcy permeability constant.
    <Isotropic_conductivity> real [0.0] </Isotropic_conductivity>
    Sets the isotropic conductivity constant.
    <Mass_damping> real [0.0] </Mass_damping>
    The mass damping coefficient used to dissipate energy.
    <Momentum_stabilization_coefficient> real [0.0] </Momentum_stabilization_coefficient>
    The stabilizing coefficient associated to the momentum equation in the Variational Multiscale Method.
    <Myocardial_zone> string [none] </Myocardial_zone>
    The region of the heart the domain is representing. This is needed for some activation models. Permissible values are
    • endocardium
    • epicardium
    • myocardium

    <ODE_solver> string [euler] </ODE_solver>
    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] </Penalty_parameter>
    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.
    <Poisson_ratio> real [0.3] </Poisson_ratio>
    The Poisson ratio property for a solid. Permissible values are between 0.0 and 0.5.
    <Relative_tolerance> real [1e-4] </Relative_tolerance>
    Convergence is deemed to occur when the size of the residual becomes sufficiently small compared with the first error estimate calculated.
    <Shell_thickness> real [0.0] </Shell_thickness>
    The thickness of the shell elements in a domain.
    <Solid_density> real [0.5] </Solid_density>
    The density property for a solid.
    <Source_term> real [0.0] </Source_term>
    The volumetric source term (e.g. heat or dye) within a domain.
    <tau_fi> real [0.110] </tau_fi>
    The fast inward current time scale used in the Bueno-Orovio cell activation model.
    <tau_si> real [1.88750] </tau_si>
    The slow inward current time scale used in the Bueno-Orovio cell activation model.
    <Time_step_for_integration> real [0.0] </Time_step_for_integration>
    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
    <Amplitude> real [0.0] </Amplitude>
    The stimulus amplitude.
    <Cycle_length> real [0.0] </Cycle_length>
    The stimulus cycle length.
    <Duration> real [0.0] </Duration>
    The stimulus duration.
    <Start_time> real [0.0] </Start_time>
    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

    <LS type=linear_solver_type>

    <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

    <Max_iterations> integer [5] </Max_iterations>
    The maximum number of iterations to use for a linear solve.
    <Krylov_space_dimension> integer [200] </Krylov_space_dimension>
    The dimension of the Kyrlov matrix used to approximate a linear system.
    <Tolerance> real [0.4] </Tolerance>
    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] </Absolute_tolerance>
    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] </NS_GM_max_iterations>
    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] </NS_GM_tolerance>
    The tolerance used to control the convergence of the GMRES portion of the bi-partitioned iterative algorithm.
    <NS_CG_max_iterations> integer [5] </NS_CG_max_iterations>
    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] </NS_CG_tolerance>
    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
    <Preconditioner> string [none] </Preconditioner>
    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

    <svZeroDSolver_interface>
    [ svZeroDSolver Interface Parameters ]
    </svZeroDSolver_interface>

    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

    <Coupling_type> string [none] </Coupling_type>

    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 .
    <Configuration_file> string [none]
    </Configuration_file>

    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]
    </Shared_library>

    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
    On MacOS and Linux this library will be installed in the directory /usr/local/sv/svZeroDSolver/DATE/lib.

    <Initial_flows> real [0.0]
    </Initial_flows>

    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]
    </Initial_pressures>

    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

    <Output type=output_type>
    [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

    <Output type="Spatial"> <Displacement> true </Displacement> <Stress> true </Stress> <Output>

    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

    <Add_BC name=face_name>

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

    Some forms of boundary conditions only apply to certain equations.

    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 spatial distribution of values is used for only Neumann or Traction boundary condition types to specify a spatially varying load (pressure/traction). The values are stored as nodal values for a face in a VTK VTP format file.

    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

    <Bct_file_path> string [none] </Bct_file_path>

    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]
    </CST_shell_bc_type>

    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]
    </Impose_flux>

    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]
    </Impose_on_state_variable_integral>

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

    <Initial_displacements_file_path> string [none] </Initial_displacements_file_path>

    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]
    </Penalty_parameter>

    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.

    <Prestress_file_path> string [none] </Prestress_file_path>

    Use the VTK VTP-format file to initialize CMM using a prestressed wall under equilibrium with fluid traction.

    <Profile> string [flat] </Profile>

    Set the spatial distribution of a state variable on the face. Acceptable values: Flat, Parabolic or User_defined

    <Ramp_function> boolean [false]
    </Ramp_function>

    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.

    <Spatial_profile_file_path> string [none] </Spatial_profile_file_path>

    Use the given text file to set the spatial distribution of a state variable for a User_defined profile.

    <Spatial_values_file_path> string [none] </Spatial_values_file_path>

    The path to the VTK VTP-format file used to set the spatially varying body force to a face.

    <svZeroDSolver_block> string [none] </svZeroDSolver_block>

    The name of the svZeroDSolver block this boundary face is coupled to.

    <Temporal_and_spatial_values_file_path> string [none] </Temporal_and_spatial_values_file_path>

    The path to the text file containing temporal and spatial values.

    <Temporal_values_file_path> string [none] </Temporal_values_file_path>

    The path to the text file containing temporal values.

    <Time_dependence> string [steady]
    </Time_dependence>

    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]
    </Traction_values_file_path>

    The path to the VTK VTP-format file containing nodally varying traction values.

    <Traction_multiplier> real [1.0]
    </Traction_multiplier>

    The value used to scale the traction values read from from a file.

    <Type> string
    </Type>

    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
    </Value>

    The value of the state variable.

    RCR Boundary Condition Parameters

    <RCR_values>
    <Capacitance> real
    </Capacitance>
    Capacitance.
    <Distal_resistance> real </Distal_resistance>
    Distal resistance.
    <Distal_pressure> real </Distal_pressure>
    The distal pressure used to initialize an RCR boundary condition.
    <Initial_pressure> real </Initial_pressure>
    The initial pressure used to initialize an RCR boundary condition.
    <Proximal_resistance> real </Proximal_resistance>
    Proximal resistance.
    </RCR_values>

    Dirichlet Boundary Condition Parameters

    <Effective_direction> vector [none] </Effective_direction>
    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] </Penalty_parameter_normal>
    If the Dirichlet boundary condition is weakly applied then it is applied in a normal direction.
    <Penalty_parameter_tangential> real [0.0] </Penalty_parameter_tangential>
    If the Dirichlet boundary condition is weakly applied then it is applied in a tangential direction.
    <Weakly_applied> boolean [false] </Weakly_applied>
    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] </Zero_out_perimeter>
    If true then the solver will zero out the nodes shared by two adjacent faces.

    Neumann Boundary Condition Parameters

    <Distal_pressure> real [0.0] </Distal_pressure>
    The distal pressure used to initialize an RCR boundary condition.
    <Follower_pressure_load> boolean [false] </Follower_pressure_load>
    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.
    <Fourier_coefficients_file_path> string [none] </Fourier_coefficients_file_path>
    A text file containing the Fourier coefficients interpolating a time-dependent state variable.
    <Undeforming_neu_face> boolean [false] </Undeforming_neu_face>
    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

    <Apply_along_normal_direction> boolean [false] </Apply_along_normal_direction>
    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] </Damping>
    The damping constant for the spring-mass-damper force used to implement a Robin boundary condition.
    <Stiffness>> real [1.0] </Stiffness>
    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

    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 files

    See 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

    1. Number of processors used to run the simulation
    2. Number of equations
    3. Number of meshes
    4. Number of nodes
    5. Number of unknowns in the coupled 0D-3D problem
    6. Number of degrees of freedom
    7. An error flag

    These values are followed by an integer (4-byte) and two double (8-byte) values containing

    1. Simulation time step
    2. Simulation time
    3. 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..

    The restart file contains information about how the simulation data (e.g. finite element mesh) was partitioned to run on multiple processors during a parallel simulation. A simulation can therefore not be restarted using a different number of processors than was used for the original simulation.

    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 MultiPhysics Tool does not yet support writing svMultiPhysics solver XML files.

    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.

    See the Getting Started / Introduction documentation for a description of the SimVascular graphical user interface (GUI) controls and useful application Tools concepts.

    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.
    The Create MultiPhysics Simulation Job popup menu.

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

    Density TextBox - Fluid density
    Viscosity TextBox - Fluid viscosity



    Linear elasticity equation properties

    Density TextBox - Solid density
    Elastic modulus TextBox - Elastic modulus
    Poisson ratio TextBox - Poisson ratio. Value <= 0.5



    Structure equation properties

    Density TextBox - Solid density
    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

    Density TextBox - Fluid density
    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

    Add Button - Select to add a boundary condition. A Boundary Condition Setup popup window is displayed with a list of mesh boundary faces. Selecting a face name allows for setting the various parameters for a given type of boundary condition.

    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

    Time dependence

    Steady CheckBox - If selected

    Value TextBox -

    Unsteady CheckBox - If selected

    Temporal values file FileBrowser - Select to identify the time values

    Resistance - CheckBox - If selected

    Value TextBox -

    Coupled - CheckBox - If selected

    General - CheckBox - If selected

    Temporal and spatial values file FileBrowser - Select to identify the time values

    Projection - CheckBox - If selected

    Profiledependance

    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.