This document explains variable functions (VFs), which are MML's syntax for treating a model variable v as a function v(). For example, values of a dynamic scalar variable v(t) may queried via v(expression), such as v(t-delay) or v(w(t)). This document discusses VFs in their full generality. Those whose interest is limited to creating delay lines, should instead see the document Using Delay Lines in MML.
Prerequisites:
- Introductory Survey of MML (required)
- Introduction to the JSim GUI (recommended)
Contents:
- VF Basics
- Quick Summary of VF Basics
- Creating Delay Lines using VFs
- Data Driven Functional Relationships and VFs
- VFs and PDE Boundaries
- Comments or Questions?
VF Basics
In MML, declaring a variable 1-dimensional v implies a function v(expr), where "expr" is an algebraic expression of other model variables, that can be used in model equations. Similarly, an N-dimensional variable v implies a function v(expr1, expr2, ... exprN). In the following example, w is calculated using the VF call v(t^2):
// Variable Function introductory example
math vfb1 {
realDomain t;
t.min=0; t.max=1; t.delta=0.1;
real u(t) = sin(10*t);
real v(t) = u(t^2);
}
When the model runs, each value for v(t) will be calculated by linearly interpolating u values at the t grid points closest to t^2. Since sin(t) is not linear, the interpolation introduces some error. v(t) is a better approximation of sin(t^2) when t.delta is small than when t.delta is large. In general, VFs of non-linear expressions always introduce some interpolation error, which may be minimized by choosing an appropriate t.delta.
Using the default time grid, the required interpolation values are always well-defined because 0 < t < 1 implies t^2 is also. However, if we set t.max=2, running the model will generate an error because the required interpolation value u(4) is out of bounds. We can prevent this run-time error by using an "if" clause:
// Variable Function with grid-range protection
math vfb1 {
realDomain t;
t.min=0; t.max=1; t.delta=0.1;
real u(t) = sin(10*t);
real v(t) = if (t>1) u(1) else u(t^2);
}
Another type of problem occurs if you set t.min=-1, all the values of v(t) for t < 0 are NaNs. This is because the calculations of u and v proceed in tandem from t.min to t.max. At t=t.min=-1, calculating v(-1) requires u(1), which hasn't been calculated yet, thus resulting in a NaN. In an ideal world, the JSim compiler would recognize this calculation order dependency and reorder the run-time calculations to accomodate it. In reality, unfortunately, this sort of order recognition is a very hard problem in complex models, and is presently without satisfactory solution. For now, modeler must be aware of calculation order dependency when writing VFs. One simple way to do this is by following the examples in this document.
Quick Summary of VF Basics
- Compensate for VF non-linearity by choosing an appropriately fine grid spacing.
- Protect VF calls from domain value overflow (or underflow) by using "if" clauses;
- Be aware of calculation order dependencies, following examples in this document whenever possible.
Creating Delay Lines using VFs
VFs are one of several mechanism for creating delay lines. All are discussed in the document Using Delay Lines in MML.
Data Driven Functional Relationships and VFs
VFs also provide a mechanism for incorporating external (often experimental) data into a model. Consider an initially charged capacitor(C), draining its charge(Q) over a coupled resistor(R). A simple model for this system, assuming constant R, is as follows:
// RC circuit, R constant
import nsrunit;
unit conversion on;
math rc1 {
realDomain t sec; // time
t.min=0; t.max=5; t.delta=0.1;
real R = 1 ohm; // resistance
real C = 1 farad; // capacitance
real Q0 = 1 coulomb; // initial charge on capacitor
real Q(t) coulomb; // charge on capacitor at time t
real V(t) = Q/C; // voltage drop over capacitor
when (t=t.min) Q=Q0;
Q:t = -V/R; // charge dissipated over resistor
}
Now suppose R is a non-linear function of the applied voltage that has been measured experimentally and stored in a JSim-compatible data set (see Data Files and Project Data Sets). MML for this system is below:
// RC circuit, non-linear R driven by external data
import nsrunit;
unit conversion on;
math rc1 {
realDomain t sec; // time
t.min=0; t.max=5; t.delta=0.1; // time grid definition
real C = 1 farad; // capacitance
real Q0 = 1 coulomb; // initial charge on capacitor
realDomain vr volt; // voltage over resistor
vr.min=0; vr.max=Q0/C; vr.ct=101; // vr grid defintion
extern real R(vr) ohm;// resistance
real Q(t) coulomb; // charge on capacitor at time t
real V(t) = Q/C; // voltage drop over capacitor
when (t=t.min) Q=Q0;
Q:t = -V/R(V); // charge dissipated over resistor
}
This model introduces a new realDomain vr representing the voltage over R and that R itself is now declared "extern" with domain vr. This construct allows R's values to be drawn from the external data set via a JSim function generator. Furthermore, the ODE for Q uses the VF R(V) to calculate the resistance appropriate to the current voltage. Also note that the current voltage V(t) is a "real", not a "realDomain", and so R cannot declare V as its domain.
VFs and PDE Boundaries
Models containing PDEs in series can exhibit a VF-related problem with NaNs if boundary conditions are not properly handled. Consider flow through a stirred tank C2(t) and an spatial distributed pipe C1(t,x). The pipe's input (C1in) is drawn from the tank, and the pipe's output (C1out) empties back into the tank. A first attempt at coding might look as follows:
// tank/pipe system with pipe output calculated using VF
import nsrunit;
unit conversion on;
math vfbc1 {
realDomain t sec; // time grid
t.min=0; t.max=10; t.delta=.1;
realDomain x cm; // pipe spatial grid
real L=1 cm; x.min=-L; x.max=L; x.ct=51;
real C1(t,x) mM, C2(t) mM; // pipe (C1) & tank(C2) conc.
real C1in(t) mM, C1out(t) mM; // pipe input & output conc.
real F = 1 ml/sec; // flow rate
real V1 = 1 ml, V2 = 1 ml; // pipe & tank volumes
real D = 1 cm^2/sec; // diffusion constant
when (t=t.min) { C1=C2*(x/x.max)^2; C2=1; }// initial pipe & tank conc.
when (x=x.min) C1=C1in; // pipe left-hand BC
when (x=x.max) C1:x = 0; // pipe right-hand BC
C1:t = D*C1:x:x - 2*F*L/V1*C1:x;// pipe PDE
C2:t = F/V2*(C1out-C2); // tank ODE
C1in = C2; // pipe input calculate from C2
C1out = C1(t, x.max); // pipe output calculated via VF
}
Here we calculate C1out at the pipe's right-hand boundary using a VF C1out=C1(t,x.max). Unfortunately, this code will fail at run-time with a "NaNs detected" message because the VF doesn't provide enough information to the compiler to sort out the circular dependencies of C1 and C2. If, however, we use a "when" clause to calculate C1out, the compiler can recognize the C1/C2 dependencies, and order the calculations properly:
// tank/pipe system with pipe output calculated using "when" clause
import nsrunit;
unit conversion on;
math vfbc2 {
realDomain t sec; // time grid
t.min=0; t.max=10; t.delta=.1;
realDomain x cm; // pipe spatial grid
real L=1 cm; x.min=-L; x.max=L; x.ct=51;
real C1(t,x) mM, C2(t) mM; // pipe (C1) & tank(C2) conc.
real C1in(t) mM, C1out(t) mM; // pipe input & output conc.
real F = 1 ml/sec; // flow rate
real V1 = 1 ml, V2 = 1 ml; // pipe & tank volumes
real D = 1 cm^2/sec; // diffusion constant
when (t=t.min) { C1=C2*(x/x.max)^2; C2=1; }// initial pipe & tank conc.
when (x=x.min) C1=C1in; // pipe left-hand BC
when (x=x.max) C1:x = 0; // pipe right-hand BC
C1:t = D*C1:x:x - 2*F*L/V1*C1:x;// pipe PDE
C2:t = F/V2*(C1out-C2); // tank ODE
C1in = C2; // pipe input calculate from C2
when (x=x.max) C1out = C1; // pipe output calculated via "when" clause
}
This model will run without the NaN error. Here, the proper solution is to not use a VF, but to use a "when" clause instead.
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.