# Using Partial-Differential Equations in MML

*This page is for the current JSim version 2.0 and higher. 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:

- Introductory Survey of MML (required)
- Using ODEs in MML (recommended)
- Introduction to the JSim GUI (recommended)

Contents:

- First Examples
- Types of Boundary Conditions for Advection/Diffusion Problems
- IC/BC Consistency in PDEs
- Available Solvers
- The LSFEA Solver
- The MacCormack Solver
- The Toms 731 Solver
- Comments or Questions?

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

import nsrunit; unit conversion on; // SHORT DESCRIPTION // Model for the flow of a substance in a tube. math Advection { realDomain t sec; t.min=0; t.max=15; t.delta=0.1; real L = 0.1 cm, // Length of tube Ngridx = 51; // Number of grid points spatially realDomain x cm; x.min=0; x.max=L; x.ct=Ngridx; real C(x,t) mM; // Concentration of substance extern real Cin(t) mM; // Inflow Concentration extern real C0(x) mM; // Initial Distribution of C real Cout(t) mM, // OutflowConcentration F = 1 ml/min, // Flow rate V = 0.05 ml, // Volume of tube U cm/sec, // Flow Velocity Tmean sec; // Time for inflow to appear as outflow U = F*L/V; Tmean = V/F; // Initial Condition when(t=t.min) {C(x,t)=if(x=x.min) Cin else C0;} // Boundary Conditions when (x=x.min) { -U*C = -U*Cin;} // Left Hand BC when (x=x.max) { C:x = 0; Cout = C;} // Right Hand BC // Partial Differential Equation C:t = -U*C:x ; }

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

import nsrunit; unit conversion on; // SHORT DESCRIPTION // Axial diffusion in one dimension. Uses only MacCormack or Toms732 solver. math SimpleDiffusionEquation { realDomain t sec; t.min=0; t.max=15; t.delta=0.1; real L = 0.1 cm, // Length of tube Ngridx = 51; // Number of grid divisions realDomain x cm; x.min=0; x.max=L; x.ct=Ngridx; real C(x,t) mM; // Concentration of substance extern real C0(x) mM; // Initial Distribution of C real D = 0.00005 cm^2/sec; // Diffusion coefficient // Initial Condition when(t=t.min) {C(x,t) = C0;} // Boundary Conditions when (x=x.min) { D*C:x = 0;} // Left Hand BC when (x=x.max) { D*C:x = 0; } // Right Hand BC // Partial Differential Equation C:t = D*C:x:x; real integralCdx mM*cm; // Integral C(x,t.max)*dx from x.min to x.max integralCdx = integral(x=x.min to x.max, C(x,t.max)); }

*(Java 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

import nsrunit; unit conversion on; // SHORT DESCRIPTION // Models flow (left to right) of a substance in tube with axial diffusion. math AdvectionDiffusion { realDomain t sec; t.min=0; t.max=15; t.delta=0.1; real L = 0.1 cm, // Length of tube Ngridx = 51; // Number of grid points realDomain x cm; x.min=0; x.max=L; x.ct=Ngridx; real C(x,t) mM; // Concentration of substance extern real Cin(t) mM; // Inflow Concentration extern real C0(x) mM; // Initial Distribution of C real Cout(t) mM, // Outflow Concentration F = 1 ml/min, // Flow rate in tube V = 0.05 ml, // Volume of tube D = 0.001 cm^2/sec, // Diffusion coefficient U cm/sec; // Velocity U = F*L/V; // Initial Condition when(t=t.min) {C(x,t)=if(x=x.min) Cin else C0;} // Boundary Conditions when (x=x.min) { -U*C + D*C:x = -U*Cin; } // Left Hand Total Flux BC when (x=x.max) { D*C:x=0; Cout=C;} // Right Hand Neumann BC // Partial Differential Equation C:t = -U*C:x + D*C:x:x; }

*(Java 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

import nsrunit; unit conversion on; // SHORT DESCRIPTION // Counter current exchanging flows in two tubes with axial diffusion. math AdvectionDiffusion { realDomain t sec; t.min=0; t.max=15; t.delta=0.1; real L = 0.1 cm, // Length of tubes Ngridx = 51; // Number of grid points realDomain x cm; x.min=0; x.max=L; x.ct=Ngridx; //Left to Right Concentration real Clr(x,t) mM; // Concentration of substance extern real Clrin(t) mM; // Inflow Concentration extern real Clr0(x) mM; // Initial Distribution of C real Clrout(t) mM; // Outflow Concentration //Right to Left Concentration real Crl(x,t) mM; // Concentration of substance extern real Crlin(t) mM; // Inflow Concentration extern real Crl0(x) mM; // Initial Distribution of C real Crlout(t) mM; // Outflow Concentration real PS = 3 ml/min, // Exchange rate between two tubes F = 1 ml/min, // Flow rate in tubes V = 0.05 ml, // Volume of tubes D = 0.0001 cm^2/sec, // Diffusion coefficient U cm/sec; // Velocity U = F*L/V; // Initial Conditions when(t=t.min) {Clr(x,t)=if(x=x.min) Clrin else Clr0; Crl(x,t)=if(x=x.max) Crlin else Crl0;} // Boundary Conditions when (x=x.min) { -U*Clr + D*Clr:x = -U*Clrin; } // Left Hand Total Flux BC when (x=x.max) { D*Clr:x=0; Clrout=Clr;} // Right Hand Neumann BC when (x=x.max) { U*Crl + D*Crl:x = U*Crlin; } // Right Hand Total Flux BC when (x=x.min) { D*Crl:x=0; Crlout=Crl;} // Left Hand Neumann BC // Partial Differential Equations Clr:t = -U*Clr:x + D*Clr:x:x +PS/V*(Crl-Clr); Crl:t = U*Crl:x + D*Crl:x:x +PS/V*(Clr-Crl); }

*(Java 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

/* DESCRIPTION: Shows the difference between setting Disf to zero vs. omitting Disf from the isf concentration equation. The omission causes an error. Run LOOPS with t.delta = 0.2, 0.1, 0.01; PDE solver equals LSFEA; and plot Cout and CoutNoDisf. Then repeat with the other PDE solvers. */ import nsrunit; unit conversion on; math BTEX20NoDisf { // INDEPENDENT VARIABLES realDomain t sec ; t.min=0; t.max=50; t.delta=0.2; realDomain x cm; real L=0.1 cm, Ngrid=61; x.min=0; x.max=L; x.ct=Ngrid; private x.min, x.max, x.ct; // PARAMETERS real Fp = 1 ml/(g*min), // plasma flow Vp = 0.05 ml/g, // plasma volume Visfp = 0.15 ml/g, // isf volume (Visf') PSg = 100 ml/(g*min), // exchange coefficient between plasma and isf Dp = 0.0 cm^2/sec, // plasma diffusion coefficient Disf = 0.0 cm^2/sec; // INFLOWING CONCENTRATION real Cin(t) mM; Cin = if (t>=1 and t<2) (1 mM) else (0 mM); // CONCENTRATION VARIABLES (Disf in equation for Cisf) real Cp(t,x) mM, // plasma Cisf(t,x) mM, // isf Cout(t) mM; // Outflow Concentration from plasma region // BOUNDARY CONDITIONS when (x=x.min) { (-Fp*L/Vp)*(Cp-Cin)+Dp*Cp:x = 0; Cisf:x = 0; } when (x=x.max) { Cp:x = 0; Cisf:x = 0; Cout = Cp;} // INITIAL CONDITIONS when (t=t.min) { Cp = 0; Cisf = 0; } // PARTIAL DIFFERENTIAL EQUATIONS Cp:t = -Fp*L/Vp*Cp:x + PSg/Vp*(Cisf-Cp) + Dp*Cp:x:x ; Cisf:t = PSg/Visfp*(Cp-Cisf) + Disf*Cisf:x:x ; // CONCENTRATION VARIABLES (Disf omitted from equation for CisfNoDisf) real CpNoDisf(t,x) mM, // plasma CisfNoDisf(t,x) mM, // isf CoutNoDisf(t) mM; // Outflow Concentration from plasma region // BOUNDARY CONDITIONS when (x=x.min) { (-Fp*L/Vp)*(CpNoDisf-Cin)+Dp*CpNoDisf:x = 0; } when (x=x.max) { CpNoDisf:x = 0; CoutNoDisf = CpNoDisf;} // INITIAL CONDITIONS when (t=t.min) { CpNoDisf = 0; CisfNoDisf = 0;} // PARTIAL DIFFERENTIAL EQUATIONS CpNoDisf:t = -Fp*L/Vp*CpNoDisf:x + PSg/Vp*(CisfNoDisf-CpNoDisf) + Dp*CpNoDisf:x:x ; CisfNoDisf:t = + PSg/Visfp*(CpNoDisf-CisfNoDisf) ; //<-No diffusion term }

*(Java 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:

- Dirichlet: C (at either x.min or x.max) = g(t).
- Neumann: C:x (at either x.min or x.max) = 0, primarily for the downstream boundary and required for LSFEA solver. C:x (at either x.min or x.max) = h(t), used with MacCormack and Toms731 solver.
- Total Flux: -U*C + D*C:x = -U*Cin (at either x.min or x.max), primarily for the upstream boundary.

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

- LSFEA - Lagrangian Sliding Fluid Element Algorithm;
- MacCormack - 2nd order finite difference method with optional flux-corrected transport for hyperbolic and parabolic equations;
- 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

- D*u:x:x is called the diffusion term, D the diffusion coefficient;
- B*u:x is called the advection term;
- S is called the source or sink term;
- Algebraic expressions B, D, S, f1, f2, f3, g1, g2 and g3 may not contain any derivative variables (e.g. u:x);
- B and D must be spatially invariant (no x dependency);
- D must be non-negative;
- f1, f2, f3, g1, g2 and g3 may not contain u;
- If B=0 then f1, f3, g1 and g3 must be zero;
- If B>0, then f1 must be non-zero;
- If B<0, then g1 must be non-zero;
- At least one of f1 and g1 must be zero;
- At least one of f2 and g2 must be non-zero;
- at least one PDE in the equation block must have a non-zero B value;

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

- B, D, S, f1, f2, f3, g1, g2, g3 are defined as in LSFEA and may no contain any derivative variables;
- D must be non-negative.
- f1, f2, f3, g1, g2 and g3 may not contain u;

For further information on the MacCormack solver see:

- "Computational Techniques for Fluid Dynamics 2" by C.A.J. Fletcher, Springer-Verlag, page 410.
- "Flux-Corrected Transport.I.Shasta, A Fluid Transport Algorithm That Works" by J.P. Boris and D.L. Book, J. Comp. Phys. 11, pages 38-69, 1973.

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

- D and P are expressions of u, u:x or other non-derivative model variables (note that this is a more general form than the one used for LSFEA and MacCormack);
- f1, f2, f3, g1, g2 and g3 may not contain u;
- D is symbolically differentiable (i.e. it may not contain extern variables or JSim F&P calls);

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 U01HL122199 Analyzing the Cardiac Power Grid, 09/15/2015 - 05/31/2020, 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.