This page will look better in a graphical browser that supports web standards, but is accessible to any browser or internet device.

Served by Samwise.

Cardiac Physiome Society workshop: November 6-9, 2017 , Toronto

Using Partial-Differential Equations in MML

This page is for the current JSim version 2.0. Click here for the earlier JSim 1.6 version.

Introduction

This document describes how to use partial differential equations (PDEs) within MML models. JSim currently support only first degree PDEs in one spatial dimension.

Prerequisites:

Contents:

See also:

First Examples

JSim can solve a large class of 1st degree 1-dimensional partial differential equations (PDEs). The five examples illustrated here are the type of advection and diffusion equations used to model concentrations of materials in capillaries.

Example 1: The Advection Equation


(Java plugin required)

This is a pure advection (aka wave) equation in a tube. What comes in at the left hand boundary exits at the right hand boundary.

There are two independent variables, t in seconds, where t stands for time, and x in centimeters, where x is a spatial variable. The tube is parameterized by length L, volume V, and flow F. The velocity in the tube, U, is given by U=F*L/V. The transit time, Tmean (time for the input concentration to appear at the outflow) is given by Tmean = V/F.

The left hand Dirichlet boundary condition is given as -U*C = -U*Cin.

The right hand Neumann boundary condition, C:x=0, is, strictly speaking, extraneous, insofaras there is only one derivative in the x space. However, JSim requires two boundary conditions for PDE's.

The result is the advection of the wave form (Lagged-normal curve representing a bolus injection) in the downstream direction. See Figures 1 and 2 in the JSim applet.

Example 2: The Diffusion Equation

This models requires either the MacCormack or the Toms731 solver because it lacks flow.


(Java plugin required)

A simple mass balance check is performed by integrating C with respect at x at model end time, t.max.

The left and right Neumann boundary conditions can be written as either

D*C:x = 0;
or
C:x=0;

The result is an uncentered spike in the initial condition gradually becomes rounded, flattens out. An integral of C(x,t.max) with respect to x, shows that the amount of material is conserved.

Example 3: The Combined Advection/Diffusion Equation


(Java plugin required)

The left hand total flux boundary condition is written as -U*C + D*C:x = -U*Cin. The right hand boundary Neumann boundary condition is written as either D*C:x = 0 or C:x = 0.

Figure 1 in the JSim applet shows the outflowing concentration has a reduced peak height and is more spread out relative to the inflowing concentration. Figure 2, a contour plot of the concentration in time (vertical axis) and and space (horizontal axis) shows the concentration peak decaying and spreading out with time.

Example 4: A Multiple Variable Problem


(Java plugin required)

As with ODEs, PDE variables may interact with other PDE and ODE variables. This is an example of counter-current flow and exchange. Two tubes (capillaries) flowing in opposite direction exchange material with each other. The tube with flow from left to right (lr) is similar to the combined advection diffusion equation with exchange with the tube with flow from right to left (rl).

This model illustrates both left (lr) and right (rl) total flux boundary conditions and right (lr) and left (rl) Neumann boundary conditions.

Example 5: A Cautionary Note About Mixing PDEs with ODEs


(Java plugin required)

Mixing PDEs and ODEs can result in computational errors. This is an example of a two region model (BTEX20) which compares having the diffusion coefficient in the non-flowing region, Disf, set to zero with having no diffusion coefficient in the equation for the non-flowing region. The results are different, and omitting the diffusion coefficient in the PDE for the non-flowing region leads to erroneous results. In models where the non-flowing regions may represent transporters bound to membranes and no diffusion is possible, it is still necessary to include a diffusion term with the diffusion coefficient set to zero. It is permissible to have PDE and ODE models in series, but not in parallel.

Types of Boundary Conditions for advection diffusion problems

There are primarily three types of boundary conditions that JSim can use for advection-diffusion problems:

Initial and Boundary Consistency in PDEs

A PDE variable's initial and boundary condition on the inflow side of the problem must be consistent. Consider the following initial and boundary condition:

when(t=t.min) {C(x,t)=f(x);} //implies C(x.min,t.min)=f(x.min)
when(x=x.min) {C(x,t)=g(t);} //implies C(x.min,t.min)=g(t.min)

Unless f(x.min)=g(t.min), the specification is inconsistent. This problem is fixed by either rewriting the initial condition as

when(t=t.min) {C(x,t)=if(x=x.min) g(t) else f(x); } 
// implies C(x.min,t.min) = g(t.min) not f(x.min),

or rewriting the boundary condition as

when(x=x.min) {C(x,t) = if(t=t.min) f(x) else g(t); }
// implies C(x.min,t.min) = f(x.min) not g(t.min).

Jsim rejects model runs with inconsistent boundary specifications.

Available PDE Solvers

JSim currently has three PDE solvers available:

  1. LSFEA - Lagrangian Sliding Fluid Element Algorithm;
  2. MacCormack - 2nd order finite difference method with optional flux-corrected transport for hyperbolic and parabolic equations;
  3. Toms731 - adaptive moving grid method for univariate partial differential equations;

LSFEA is the fastest solver, but is applicable to the narrowest range of equations. For example, LSFEA does not handle pure diffusion equations. Toms731 somewhat slower than LSFEA, but handles some equations (e.g. pure diffusion) LSFEA can't. Toms731 is slower still, but handles some non-linear problems the other solvers can't. An exact mathematical description of what equations each solver can handle follows in the remaining sections of this document. However, the user need not be overly concerned with the exact solver requirements, since the JSim compiler determines which solvers are appropriate for a set of PDEs, the fastest of which becomes the default solver. For more advanced users, the JSim user environment allows the user to select amoung available solvers.

In order to describe JSim's PDE solvers further, we must first describe a bit about how JSim processes PDEs. When the JSim compiler encounters a PDE, it attempts to factor the equation and boundary conditions into the standard equation form for each PDE solver, using simple algebraic manipulations. The standard form varies betweens solvers. If the factorization for a given solver succeeds, JSim will allow that solver to be used at run-time. Bear in mind that successful numeric calculation depends both equations stability and the solver used. If the factorization fails for all available solvers, then the compiler will abort with an appropriate error message.

We define a PDE equation block as a set of PDEs whose state equation or boundary conditions share state variables. PDE solvers operate on blocks of equations, rather than on individual equations, and some properties of the block are relevant to solver performance. In the following example, u and v share an equation block because their state equations depend upon each other:

      u:t = u:x:x - u:x - u + v;
      v:t = v - u;

The LSFEA Solver

The LSFEA solver supports PDE state and boundary equations that can be factored as follows:

      u:t = D*u:x:x - B*u:x + S;
      when (x=x.min) f1*u + f2*u:x = f3;
      when (x=x.max) g1*u + g2*u:x = g3;

where

The MacCormack Solver

The MacCormack solver requires the PDE state and boundary equations that can be factored as follows:

      u:t = D*u:x:x - B*u:x + S;
      when (x=x.min) f1*u + f2*u:x = f3;
      when (x=x.max) g1*u + g2*u:x = g3;

where

For further information on the MacCormack solver see:

The Toms731 Solver

The Toms731 solver requires the PDE state and boundary equations that can be factored as follows:

      u:t = D*u:x:x + P;
      when (x=x.min) f1*u + f2*u:x = f3;
      when (x=x.max) g1*u + g2*u:x = g3;

where

For further information on the Toms731 solver see:

Comments or Questions?

[This page was last modified 06Jul12, 3:13 pm.]

Model development and archiving support at physiome.org provided by the following grants: NIH/NIBIB BE08407 Software Integration, JSim and SBW 6/1/09-5/31/13; NIH/NHLBI T15 HL88516-01 Modeling for Heart, Lung and Blood: From Cell to Organ, 4/1/07-3/31/11; NSF BES-0506477 Adaptive Multi-Scale Model Simulation, 8/15/05-7/31/08; NIH/NHLBI R01 HL073598 Core 3: 3D Imaging and Computer Modeling of the Respiratory Tract, 9/1/04-8/31/09; as well as prior support from NIH/NCRR P41 RR01243 Simulation Resource in Circulatory Mass Transport and Exchange, 12/1/1980-11/30/01 and NIH/NIBIB R01 EB001973 JSim: A Simulation Analysis Platform, 3/1/02-2/28/07.