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
- The Advection Equation
- The Diffusion Equation
- The Combined Advection/Diffusion Equation
- A Multiple Variable Problem
- A Cautionary Note About Mixing PDEs with ODEs
- 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 ;
}
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));
}
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;
}
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);
}
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
}
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?
Model development and archiving support at https://www.imagwiki.nibib.nih.gov/physiome 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.