[cherry]<Thyme>Scratch>spSolve.mesa!1

File: [Cherry]<Thyme>Cedar5.2>System>spSolve.mesa

Last Edited by: SChen, September 15, 1984 7:44:54 pm PDT

DIRECTORY
BasicTime USING [Now],
IO USING [char, int, PutF, PutFR, real, rope],
Labels USING [Set],
Real USING [Fix, FixC, FixI, Float, RoundLI, SmallestNormalizedNumber],
RealOps USING [SqRt],
RefText USING [Equal],
Rope USING [ROPE],
spGlobals USING [Aborted, argList, argNodes, branchLinkPtr, branchPtr, Canned, CheckPoint, dumpAll, error, findNodeOrBranch, getSignedNumber, Handle, initPlot, initPrint, maxCountsPerStep, maxLog, modBrPtr, modFuncPtr, next, NextStage, nodePtr, plotFromList, printFromList, PutMsgLine, RefFuncBody, RefModBody, RefC, RefI, RefL, RefR, RefV, ResetCheckPoint, ShowDetails];

spSolve: CEDAR PROGRAM
IMPORTS BasicTime, IO, Labels, Real, RealOps, RefText, spGlobals
EXPORTS spGlobals=

BEGIN OPEN spGlobals;

solveError: SIGNAL[s: Rope.ROPE]= CODE;
solveTest: SIGNAL[s: Rope.ROPE]= CODE;
Retreat: PUBLIC SIGNAL[cause: Rope.ROPE]= CODE;
Failure: PUBLIC SIGNAL[errorNum: NAT]= CODE;

nParms: NAT= 28;
parmNames: ARRAY[0..nParms) OF REF TEXT= ["maxiter", "printiter", "initstep", "squish", "tol", "dt", "tmax", "printstep", "rungekutta", "tmin", "printall", "amiter", "amup", "inttol", "amdown", "vmin", "vmax", "dvdtmax", "dvdtfac", "tbreak", "rkerr", "amerr", "amdelay", "floor", "minratio", "rkfac", "retratio", ""];

updateFunctions: PROC[handle: Handle, f: modFuncPtr]= {
mb: modBrPtr;
source: argNodes;
args, res: argList;
allInactive: BOOL;

handle.vars.modelsSkipped← 0;
UNTIL f=NIL DO
source← f.arguments;
args← f.argVector;
allInactive← TRUE;
IF source # NIL THEN
FOR i: NAT IN [0..source.size) DO
args[i]← source[i].nHist.y;
allInactive← (allInactive AND source[i].marked);
ENDLOOP;
WITH f.body SELECT FROM
fcn: RefFuncBody => {
fcn.branch.comVal← fcn.functionProc[handle.vars.t, args, f.parmVector];
};
mod: RefModBody =>
IF ~allInactive THEN {
-- since ~allInactive, from above we know that source is not Nil
oldArgs: argList← mod.oldArgVector;
FOR i: NAT IN [0..source.size) DO
oldArgs[i]← source[i].nHist.oldy;
ENDLOOP;
res← mod.modelResults;
mod.modelProc[args, oldArgs, f.parmVector, res];
FOR i: NAT IN [0..source.size) DO
source[i].nHist.y← args[i];
ENDLOOP;
mb← mod.modelBranches;
UNTIL mb=NIL DO
mb.b.comVal← res[mb.b.modelIndex];
mb← mb.nextBranch;
ENDLOOP;
}
ELSE handle.vars.modelsSkipped← handle.vars.modelsSkipped + 1;
ENDCASE;
f← f.nextMFPtr;
ENDLOOP;
}; -- updateFunctions

GSiteration: PROC[handle: Handle, dT: REAL, print, firstTime, mark: BOOL← FALSE]
RETURNS[nIter: INT← 0]= {
E, newE, dvdt, I, newI, cond, csum, otherV, maxI, denom: REAL;
n: nodePtr;
b: branchPtr;
vCur: branchPtr;
vCurBody: RefV;
bLink: branchLinkPtr;
plus, nodeOK, OK: BOOL;
realMaxIter: INT = IF handle.vars.maxIter >= LAST[INT]/2 THEN LAST[INT] ELSE handle.vars.maxIter + handle.vars.maxIter;

updateFunctions[handle, handle.vars.intNodeModFunc];
IF print THEN printNodeEqn[handle];
UNTIL nIter >= realMaxIter DO
OK← TRUE;
IF print THEN
IO.PutF[handle.msgStream, "t= %10.4f, iteration %2d --\n",
IO.real[handle.vars.t], IO.int[nIter]];
IF handle.vars.otherModFunc # NIL THEN
updateFunctions[handle, handle.vars.otherModFunc];
IF firstTime THEN updateFunctions[handle, handle.vars.intNodeModFunc];
FOR n← handle.vars.nodeList, n.nextNode UNTIL n=NIL DO
IF n=handle.vars.gndNode THEN LOOP;
IF nIter > 1 AND n.converged THEN LOOP;
IF print THEN IO.PutF[handle.msgStream, "Node %g: ",
IO.rope[n.nodeName.name]];
bLink← n.branches;
IF n.integrate THEN {
maxI← csum← I← 0.0;
UNTIL bLink=NIL DO
b ← bLink.branch;
plus ← bLink.pos;
WITH b.body SELECT FROM
r: RefR =>
newI← (bLink.otherNode.nHist.y - n.nHist.y)/b.comVal;
c: RefC => {
csum← csum + b.comVal;
newI← b.comVal*bLink.otherNode.nHist.f0;
};
l: RefL =>
newI← (IF plus THEN -l.iHist.y ELSE l.iHist.y)/b.comVal;
v: RefV =>
newI← IF plus THEN -v.vsCurrent ELSE v.vsCurrent;
i: RefI =>
newI← IF plus THEN -b.comVal ELSE b.comVal
ENDCASE;
maxI← MAX[maxI, ABS[newI]];
I← I + newI;
bLink← bLink.nextLink;
ENDLOOP;
E← n.nHist.f0;
dvdt← n.nHist.f0← I/csum;
IF mark THEN n.marked← (ABS[I] <= handle.vars.dvdtTrunc*maxI);
denom← MAX[ABS[E], ABS[dvdt]];
nodeOK← IF denom <= Real.SmallestNormalizedNumber THEN TRUE
ELSE (handle.vars.tol >= ABS[(dvdt - E)/denom]);
IF print THEN
IO.PutF[handle.msgStream, "Old v'= %9.2f New v'= %9.2f C= %9.2f.\n",
IO.real[E], IO.real[dvdt], IO.real[csum]];
}
ELSE {
E← n.nHist.y;
IF n.curPtr # NIL THEN {
vCur← n.curPtr;
vCurBody← NARROW[vCur.body];
newE← IF vCur.posNode=n
THEN vCur.negNode.nHist.y + vCur.comVal
ELSE vCur.posNode.nHist.y - vCur.comVal;
I← vCurBody.vsCurrent;
newI← 0.0;
UNTIL bLink=NIL DO
b ← bLink.branch;
otherV← bLink.otherNode.nHist.y;
plus ← bLink.pos;
WITH b.body SELECT FROM
r: RefR =>
newI← newI + (otherV - newE)/b.comVal;
c: RefC =>
newI← newI + b.comVal*
(bLink.otherNode.nHist.f0 - n.nHist.f0);
l: RefL =>
newI← newI + (IF plus THEN -l.iHist.y
ELSE l.iHist.y)/b.comVal;
v: RefV =>
newI← IF plus THEN newI - v.vsCurrent
ELSE newI + v.vsCurrent;
i: RefI =>
newI← IF plus THEN newI - b.comVal
ELSE newI + b.comVal
ENDCASE;
bLink← bLink.nextLink
ENDLOOP;
newI← IF n=vCur.posNode THEN I + newI
ELSE I - newI;
denom← MAX[ABS[I], ABS[newI]];
nodeOK← IF denom <= Real.SmallestNormalizedNumber THEN TRUE
ELSE (handle.vars.tol >= ABS[(newI - I)/denom]);
n.nHist.f0← IF firstTime THEN 0.0 ELSE
handle.vars.squish*n.nHist.f1 +
(1.0-handle.vars.squish)*(newE-n.nHist.oldy)/dT +
(IF vCur.posNode=n THEN vCur.negNode.nHist.f0
ELSE vCur.posNode.nHist.f0);
vCurBody.vsCurrent← newI;
IF print THEN
IO.PutF[handle.msgStream, "Old I= %9.2f New I= %9.2f",
IO.real[I], IO.real[newI]];
}
ELSE {
cond← newE← 0.0;
UNTIL bLink=NIL DO
b ← bLink.branch;
plus ← bLink.pos;
WITH b.body SELECT FROM
r: RefR => {
cond← cond + 1.0/b.comVal;
newE← newE + bLink.otherNode.nHist.y/b.comVal;
};
l: RefL =>
newE← newE + (IF plus THEN -l.iHist.y ELSE l.iHist.y)/b.comVal;
v: RefV =>
newE← IF plus THEN newE - v.vsCurrent
ELSE newE + v.vsCurrent;
i: RefI =>
newE← IF plus THEN newE - b.comVal
ELSE newE + b.comVal
ENDCASE;
bLink← bLink.nextLink
ENDLOOP;
newE← newE/cond;
denom← MAX[ABS[E], ABS[newE]];
nodeOK← IF denom <= Real.SmallestNormalizedNumber THEN TRUE
ELSE (handle.vars.tol >= ABS[(newE - E)/denom]);
};
n.nHist.y← newE;
IF print THEN IO.PutF[handle.msgStream, "Old E= %9.2f New E= %9.2f.\n",
IO.real[E], IO.real[newE]];
};
OK← OK AND nodeOK;
n.converged← nodeOK;
IF ~nodeOK THEN
FOR bLink← n.branches, bLink.nextLink UNTIL bLink=NIL DO
bLink.otherNode.converged← FALSE;
ENDLOOP;
ENDLOOP;
nIter← nIter + 1;
IF handle.vars.canned AND NOT handle.vars.forcedDump THEN SIGNAL Aborted;
IF nIter=handle.vars.maxIter THEN
IO.PutF[handle.msgStream, "t= %9.3f, iteration count exceeded",
IO.real[handle.vars.t]];
IF OK THEN EXIT;
ENDLOOP;
IF nIter >= realMaxIter THEN SIGNAL solveError["GS did not converge"]
}; -- GSiteration

printNodeEqn: PROC[handle: Handle]= {
n: nodePtr;
bLink: branchLinkPtr;
branch: branchPtr;
plus: BOOL;

FOR n← handle.vars.intNodeList, n.nextIntNode UNTIL n=NIL DO
IO.PutF[handle.msgStream, "\n%g -- %13.6f %13.6f\n",
IO.rope[n.nodeName.name],
IO.real[n.nHist.y],
IO.real[n.nHist.oldy]];
FOR bLink← n.branches, bLink.nextLink UNTIL bLink=NIL DO
branch← bLink.branch;
plus ← bLink.pos;
WITH branch.body SELECT FROM
r: RefR =>
IO.PutF[handle.msgStream, "(%13.6f - %13.6f)/%13.6f",
IO.real[bLink.otherNode.nHist.y],
IO.real[n.nHist.y],
IO.real[branch.comVal] ];
c: RefC =>
IO.PutF[handle.msgStream, "%13.6fx%13.6f",
IO.real[branch.comVal],
IO.real[bLink.otherNode.nHist.f0] ];
l: RefL =>
IF plus THEN
IO.PutF[handle.msgStream, "-1x%13.6f/%13.6f",
IO.real[l.iHist.y],
IO.real[branch.comVal]]
ELSE IO.PutF[handle.msgStream, "%13.6f/%13.6f",
IO.real[l.iHist.y],
IO.real[branch.comVal]];
v: RefV =>
IF plus THEN IO.PutF[handle.msgStream, "-1x%13.6f",
IO.real[v.vsCurrent]]
ELSE IO.PutF[handle.msgStream, "%13.6f", IO.real[v.vsCurrent]];
i: RefI =>
IF plus THEN IO.PutF[handle.msgStream, "-1x%13.6f",
IO.real[branch.comVal]]
ELSE IO.PutF[handle.msgStream, "%13.6f",IO.real[branch.comVal]];
ENDCASE;
IO.PutF[handle.msgStream, " %s(%13.6f, %13.6f)\n",
IO.rope[bLink.otherNode.nodeName.name],
IO.real[bLink.otherNode.nHist.y],
IO.real[bLink.otherNode.nHist.f0]];
ENDLOOP;
ENDLOOP;
}; -- printNodeEqn

runParms: PROC[handle: Handle]= {
OPEN handle.vars;

parmValue: PROC[booleanOrNot: {boolean, notBoolean}]
RETURNS[r: REAL, b: BOOL, n: NAT]= {
IF booleanOrNot=boolean THEN
IF item=name THEN {
IF RefText.Equal[newString, "FALSE", FALSE] THEN b← FALSE
ELSE
IF RefText.Equal[newString, "TRUE", FALSE] THEN b← TRUE
ELSE error[handle, 719];
next[handle]
}
ELSE error[handle, 704]
ELSE {
r← getSignedNumber[handle];
IF r >= 0.0 AND r < Real.Float[LAST[NAT]] THEN n← Real.FixC[r]; -- else?
};
}; -- parmValue

nameIndex: NAT;
tempR: REAL;
WHILE item=name DO
FOR i: NAT IN [0..nParms) DO
nameIndex← i;
IF RefText.Equal[newString, parmNames[nameIndex], FALSE] THEN EXIT;
ENDLOOP;
next[handle];
IF item # leftArrow THEN error[handle, 703] ELSE next[handle];
SELECT nameIndex FROM
0=> {
limit: INT ← LAST[INT];
[tempR,,]← parmValue[notBoolean];
IF tempR >= Real.Float[limit] THEN maxIter← limit
ELSE IF tempR > 0 THEN maxIter← Real.Fix[tempR];
};
1=> [, printIter,] ← parmValue[boolean];
2=> [initStep,,] ← parmValue[notBoolean];
3=> [initStep,,] ← parmValue[notBoolean];
4=> [tol,,] ← parmValue[notBoolean];
5=> [printdT,,] ← parmValue[notBoolean];
6=> [tMax,,] ← parmValue[notBoolean];
7=> [, printStep,] ← parmValue[boolean];
8=> [, doRungeKutta,]← parmValue[boolean];
9=> [tMin,,] ← parmValue[notBoolean];
10=> [, printAll,] ← parmValue[boolean];
11=> [,, AMiter] ← parmValue[notBoolean];
12=> [AMup,,] ← parmValue[notBoolean];
13=> [intTol,,] ← parmValue[notBoolean];
14=> [AMdown,,] ← parmValue[notBoolean];
15=> [Vmin,,] ← parmValue[notBoolean];
16=> [Vmax,,] ← parmValue[notBoolean];
17=> [dvdtmax,,] ← parmValue[notBoolean];
18=> [dvdtFac,,] ← parmValue[notBoolean];
19=> [tBreak,,] ← parmValue[notBoolean];
20=> [RKerr,,] ← parmValue[notBoolean];
21=> [AMerr,,] ← parmValue[notBoolean];
22=> [,, AMdelay] ← parmValue[notBoolean];
23=> [floor,,] ← parmValue[notBoolean];
24=> [minRatio,,] ← parmValue[notBoolean];
25=> [RKfac,,] ← parmValue[notBoolean];
26=> [retRatio,,] ← parmValue[notBoolean];
ENDCASE=> error[handle, 790];
IF item=comma THEN next[handle];
ENDLOOP;
}; -- runParms

predict: PROC[handle: Handle, h: REAL]= {
nodes: nodePtr ← handle.vars.intNodeList;
inds: branchPtr← handle.vars.inductorList;
incr: REAL;

h← h/24.0;

UNTIL nodes=NIL DO
IF ~nodes.marked THEN {
OPEN nodes.nHist;
incr← h*(55.0*f1 - 59.0*f2 + 37.0*f3 - 9.0*f4);
y← oldy + incr;
};
nodes← nodes.nextIntNode;
ENDLOOP;

UNTIL inds=NIL DO
indsBody: RefL← NARROW[inds.body];
{OPEN indsBody.iHist;
y← oldy + h*(55.0*f1 - 59.0*f2 + 37.0*f3 - 9.0*f4)};
inds← indsBody.nextInductor;
ENDLOOP;
}; -- predict

firstCorrector: PROC[handle: Handle, h: REAL]= {
nodes: nodePtr ← handle.vars.intNodeList;
inds: branchPtr← handle.vars.inductorList;
f: REAL;

h← h/24.0;

UNTIL nodes=NIL DO
IF ~nodes.marked THEN {
OPEN nodes.nHist;
y← oldy + h*(9.0*f0 + 19.0*f1 - 5.0*f2 + f3);
oldf← f0;
};
nodes← nodes.nextIntNode;
ENDLOOP;

UNTIL inds=NIL DO
indsBody: RefL← NARROW[inds.body];
{OPEN indsBody.iHist;
f← inds.posNode.nHist.y - inds.negNode.nHist.y;
y← oldy + h*(9.0*f + 19.0*f1 - 5.0*f2 + f3);
oldf← f};
inds← indsBody.nextInductor;
ENDLOOP;
}; -- firstCorrector

iteratingCorrector: PROC[handle: Handle, h: REAL] RETURNS[newRatio: REAL← 4.0,
OK: BOOL← TRUE]= {
nodes: nodePtr ← handle.vars.intNodeList;
inds: branchPtr← handle.vars.inductorList;
f, deltay, newy, err, denom: REAL;

h← 3.0*h/8.0;
UNTIL nodes=NIL DO
IF ~nodes.marked THEN {
OPEN nodes.nHist;
deltay← h*(f0 - oldf);
newy ← y + deltay;
denom ← MAX[ABS[newy], ABS[y]];
IF denom > handle.vars.floor THEN {
err← handle.vars.AMerr*ABS[deltay/denom];
IF handle.vars.intTol < err THEN {
OK← FALSE;
newRatio← MIN[newRatio, handle.vars.intTol/err];
handle.vars.worstNode← nodes;
};
};
y ← newy;
oldf← f0;
};
nodes← nodes.nextIntNode;
ENDLOOP;

UNTIL inds=NIL DO
indsBody: RefL← NARROW[inds.body];
{ OPEN indsBody.iHist;
f ← inds.posNode.nHist.y - inds.negNode.nHist.y;
deltay← h*(f - oldf);
newy ← y + deltay;
denom ← MAX[ABS[newy], ABS[y]];
IF denom # 0.0 THEN {
err← handle.vars.AMerr*ABS[deltay/denom];
IF handle.vars.intTol < err THEN {
OK← FALSE;
newRatio← MIN[newRatio, handle.vars.intTol/err];
};
};
y ← newy;
oldf← f;
};
inds← indsBody.nextInductor;
ENDLOOP;
}; -- iteratingCorrector

AdamsMoulton: PROC[handle: Handle, dT: REAL]
RETURNS[newRatio: REAL, nGS: INT, ok: BOOL]= {
nAM: NAT;

predict[handle, dT];
nGS← GSiteration[handle, dT, handle.vars.printAll];
firstCorrector[handle, dT];
nGS← GSiteration[handle, dT, handle.vars.printAll] + nGS;
FOR nAM← 2, nAM + 1 UNTIL nAM=handle.vars.AMiter DO
[]← iteratingCorrector[handle, dT];
nGS← GSiteration[handle, dT, handle.vars.printAll] + nGS
ENDLOOP;
[newRatio, ok]← iteratingCorrector[handle, dT];
}; -- AdamsMoulton

RKupdate: PROC[handle: Handle, c, h: REAL, last: BOOL]
RETURNS[newRatio: REAL← 1.0, ok: BOOL← TRUE]= {
nodes: nodePtr ← handle.vars.intNodeList;
inds: branchPtr← handle.vars.inductorList;
k, incr, y, err, denom: REAL;

UNTIL nodes=NIL DO
IF ~nodes.marked THEN {
nodes.nHist.sumk← c*nodes.nHist.f0 + nodes.nHist.sumk;
IF last THEN {
incr ← h*nodes.nHist.sumk;
y ← nodes.nHist.oldy + incr;
denom← MAX[ABS[y], ABS[nodes.nHist.y]];
IF denom > handle.vars.floor THEN {
err← handle.vars.RKerr*ABS[(y - nodes.nHist.y)/denom];
IF handle.vars.intTol < err THEN {
ok← FALSE;
newRatio← MIN[newRatio, handle.vars.intTol/err]
};
};
nodes.nHist.y← y;
nodes.nHist.sumk← 0.0;
}
ELSE {
incr← h*nodes.nHist.f0;
nodes.nHist.y← nodes.nHist.oldy + incr;
};
};
nodes← nodes.nextIntNode;
ENDLOOP;

UNTIL inds=NIL DO
indsBody: RefL← NARROW[inds.body];
k← inds.posNode.nHist.y - inds.negNode.nHist.y;
indsBody.iHist.sumk← c*k + indsBody.iHist.sumk;
IF last THEN {
indsBody.iHist.y← indsBody.iHist.oldy + h*indsBody.iHist.sumk;
indsBody.iHist.sumk← 0.0
}
ELSE indsBody.iHist.y← indsBody.iHist.oldy + h*k;
inds← indsBody.nextInductor;
ENDLOOP;
}; -- RKupdate

RungeKutta: PROC[handle: Handle, dT: REAL]
RETURNS[newRatio: REAL, nGS: NAT, ok: BOOL]= {
time: REAL;

time← handle.vars.t;
handle.vars.t← handle.vars.t - 0.5*dT;
[] ← RKupdate[handle, 1.0, 0.5*dT, FALSE];
nGS← GSiteration[handle, 0.5*dT, handle.vars.printAll];
[] ← RKupdate[handle, 2.0, 0.5*dT, FALSE];
nGS← GSiteration[handle, 0.5*dT, handle.vars.printAll] + nGS;
handle.vars.t← time;
[] ← RKupdate[handle, 2.0, dT, FALSE];
nGS← GSiteration[handle, dT, handle.vars.printAll] + nGS;
[newRatio, ok]← RKupdate[handle, 1.0, 0.1666667*dT, TRUE]
}; -- RungeKutta

changeStepSize: PROC[handle: Handle, ratio: REAL, RK: BOOL]= {
int: CHAR;

IF ratio < 1.0 OR handle.vars.countsPerStep < maxCountsPerStep THEN {
handle.vars.numGoodSteps← 1;

handle.vars.worstNodeLog[handle.vars.curLog]←
IF handle.vars.worstNode # NIL THEN
[node: handle.vars.worstNode,
t: handle.vars.t,
v: handle.vars.worstNode.nHist.y,
dvdt: handle.vars.worstNode.nHist.f0]
ELSE [NIL, 0.0, 0.0, 0.0];

handle.vars.curLog← (handle.vars.curLog + 1) MOD maxLog;

handle.vars.countsPerStep←
IF ratio > Real.Float[maxCountsPerStep/handle.vars.countsPerStep] THEN
maxCountsPerStep
ELSE MAX[1, Real.RoundLI[handle.vars.countsPerStep*ratio]];

int← IF RK THEN 'R ELSE 'A;
IF ShowDetails[handle] THEN {
Labels.Set[handle.time,
IO.PutFR["%12.6e", IO.real[handle.vars.t]]];
Labels.Set[handle.step,
IO.PutFR["%10d %g", IO.int[handle.vars.countsPerStep], IO.char[int]]];
};
IF handle.vars.countsPerStep=1 THEN SIGNAL solveError["Step too small"];
};
}; -- changeStepSize

saveState: PROC[handle: Handle]= {
nodes: nodePtr← handle.vars.nodeList;
inds: branchPtr← handle.vars.inductorList;
volts: branchPtr← handle.vars.vSourceList;
indsBody: RefL;
voltsBody: RefV;

UNTIL nodes=NIL DO
OPEN nodes.nHist;
oldy← y;
f4← f3; f3← f2; f2← f1; f1← f0;
IF y < handle.vars.Vmin THEN SIGNAL solveTest["Node below Vmin"]
ELSE
IF y > handle.vars.Vmax THEN SIGNAL solveTest["Node above Vmax"]
ELSE
IF ABS[f0] > handle.vars.dvdtmax THEN
SIGNAL solveTest["Rate above max"];
nodes← nodes.nextNode;
ENDLOOP;

UNTIL inds=NIL DO
indsBody← NARROW[inds.body];
{ OPEN indsBody.iHist;
oldy← y;
f4← f3; f3← f2; f2← f1;
f1← inds.posNode.nHist.y - inds.negNode.nHist.y;
};
inds← indsBody.nextInductor;
ENDLOOP;

UNTIL volts=NIL DO
voltsBody← NARROW[volts.body];
voltsBody.oldCurrent← voltsBody.vsCurrent;
volts← voltsBody.nextvSource;
ENDLOOP;
}; -- saveState

restoreState: PROC[handle: Handle]= {
nodes: nodePtr← handle.vars.nodeList;
inds: branchPtr← handle.vars.inductorList;
volts: branchPtr← handle.vars.vSourceList;
indsBody: RefL;
voltsBody: RefV;

UNTIL nodes=NIL DO
nodes.nHist.y← nodes.nHist.oldy;
nodes.nHist.f0← nodes.nHist.f1;
nodes← nodes.nextNode;
ENDLOOP;

UNTIL inds=NIL DO
indsBody← NARROW[inds.body];
indsBody.iHist.y← indsBody.iHist.oldy;
inds← indsBody.nextInductor;
ENDLOOP;

UNTIL volts=NIL DO
voltsBody← NARROW[volts.body];
voltsBody.vsCurrent← voltsBody.oldCurrent;
volts← voltsBody.nextvSource;
ENDLOOP;
}; -- restoreState

integration: PROC[handle: Handle, RK: BOOL]= {
nGS: INT;
out, ok: BOOL;
ratio, oldTime, DT: REAL;
oldLowStepcount, oldHighStepCount: INT;

saveState[handle];
oldTime← handle.vars.t;
oldLowStepcount← handle.vars.lowStepCount;
oldHighStepCount← handle.vars.highStepCount;
DO
handle.vars.lowStepCount← handle.vars.lowStepCount + handle.vars.countsPerStep;
IF handle.vars.lowStepCount > maxCountsPerStep THEN {
handle.vars.highStepCount← handle.vars.highStepCount + 1;
handle.vars.lowStepCount← handle.vars.lowStepCount - maxCountsPerStep;
out← TRUE;
}
ELSE out← handle.vars.printStep;

handle.vars.t← handle.vars.tMin + handle.vars.dT*handle.vars.lowStepCount +
(handle.vars.dT*maxCountsPerStep)*handle.vars.highStepCount;
DT← handle.vars.dT*handle.vars.countsPerStep;

{ ENABLE Retreat=> {
-- handle.msgStream.PutF["%g at t= %12.6e.\n",
-- IO.rope[cause], IO.real[handle.vars.t]];
ok← FALSE;
ratio← handle.vars.retRatio;
GOTO FlushIt;
};

IF handle.vars.numGoodSteps < 4 THEN {
[ratio, nGS, ok]← RungeKutta[handle, DT];
IF ok THEN nGS← nGS +
GSiteration[handle, DT, handle.vars.printAll,, TRUE];
ratio← RealOps.SqRt[ratio]*handle.vars.RKfac;
IF ok AND ~RK THEN
handle.vars.numGoodSteps← 1 + handle.vars.numGoodSteps;
}
ELSE {
handle.vars.AMcount← 1 + handle.vars.AMcount;
[ratio, nGS, ok]← AdamsMoulton[handle, DT];
IF ok THEN {
nGS← nGS + GSiteration[handle, DT, handle.vars.printAll,, TRUE];
IF ratio > handle.vars.AMup AND
handle.vars.AMcount >= handle.vars.AMdelay THEN {
ratio← RealOps.SqRt[RealOps.SqRt[ratio]];
changeStepSize[handle, ratio, FALSE];
handle.vars.AMcount← 0;
};
}
ELSE ratio← RealOps.SqRt[ratio];
};
EXITS
FlushIt=> NULL;
};
IF ok THEN EXIT;
restoreState[handle];
handle.vars.t← oldTime;
handle.vars.highStepCount← oldHighStepCount;
handle.vars.lowStepCount← oldLowStepcount;
changeStepSize[handle, MIN[ratio, handle.vars.minRatio],
(handle.vars.numGoodSteps < 4)];
handle.vars.AMcount← 0
ENDLOOP;

NextStage[handle];
IF out THEN printFromList[handle, nGS, handle.vars.t, handle.vars.printIter];
IF out THEN plotFromList[handle, handle.vars.t];
WHILE handle.vars.tBreak > 0.0 AND handle.vars.t >= handle.vars.tBreak DO
oldTime← handle.vars.t ENDLOOP;
}; -- integration

zeroVars: PROC[handle: Handle]= {
n: nodePtr← handle.vars.nodeList;
inds: branchPtr← handle.vars.inductorList;
vSs: branchPtr← handle.vars.vSourceList;
indsBody: RefL;
vSsBody: RefV;

handle.vars.gndNode.marked← TRUE;

UNTIL n=NIL DO
n.nHist← [];
n.marked← FALSE;
n← n.nextNode;
ENDLOOP;

UNTIL inds=NIL DO
indsBody← NARROW[inds.body];
indsBody.iHist← [];
inds← indsBody.nextInductor;
ENDLOOP;

UNTIL vSs=NIL DO
vSsBody← NARROW[vSs.body];
vSsBody.vsCurrent← vSsBody.oldCurrent← 0.0;
IF vSs.controller=NIL THEN
IF vSs.posNode=handle.vars.gndNode THEN vSs.negNode.marked← TRUE
ELSE
IF vSs.negNode=handle.vars.gndNode THEN vSs.posNode.marked← TRUE;
vSs← vSsBody.nextvSource;
ENDLOOP;
}; -- zeroVars

getStepSize: PROC[handle: Handle] RETURNS[stepOK: BOOL]= {
vMax, dvdtMax, frac: REAL← 0.0;
inds: branchPtr← handle.vars.inductorList;
nodes: nodePtr← handle.vars.nodeList;

handle.vars.lowStepCount← 0;
handle.vars.highStepCount← 0;
handle.vars.numberOfSteps← Real.FixI[handle.vars.tMax/handle.vars.printdT];
handle.vars.worstNode← NIL;

UNTIL nodes=NIL DO
vMax← MAX[ABS[nodes.nHist.y], vMax];
IF ABS[nodes.nHist.f0] > dvdtMax THEN {
dvdtMax← ABS[nodes.nHist.f0];
handle.vars.worstNode← nodes;
};
nodes← nodes.nextNode;
ENDLOOP;

frac← IF dvdtMax > 0.0 AND vMax > 0.0 THEN
handle.vars.intTol*vMax/(dvdtMax*handle.vars.printdT)
ELSE handle.vars.initStep;
handle.vars.dT← handle.vars.printdT/maxCountsPerStep;
handle.vars.dvdtTrunc← handle.vars.dvdtFac*handle.vars.intTol;
handle.vars.countsPerStep←
Real.RoundLI[maxCountsPerStep*MIN[frac, handle.vars.initStep]];
handle.vars.countsPerStep←
MAX[handle.vars.countsPerStep, 1];
IF handle.vars.countsPerStep=1 THEN {
stepOK← FALSE;
PutMsgLine[handle, "Initial step is too small!"];
}
ELSE {
stepOK← TRUE;
-- IO.PutF[handle.msgStream, "Initial step is %d.\n",
-- IO.int[handle.vars.countsPerStep]];
};
}; -- getStepSize

setICs: PUBLIC PROC[handle: Handle]= {
n: nodePtr;
b: branchPtr;
v: REAL;
allNodes: BOOL;

zeroVars[handle];
IF handle.vars.item # leftB THEN error[handle, 700] ELSE next[handle];
IF handle.vars.item=number THEN {
handle.vars.tMin← handle.vars.value;
next[handle];
IF handle.vars.item=comma THEN next[handle] ELSE error[handle, 702];
};
DO
allNodes← (handle.vars.item=star);
IF allNodes THEN next[handle] ELSE [n, b]← findNodeOrBranch[handle];
IF handle.vars.item=leftArrow THEN next[handle] ELSE error[handle, 703];
v← getSignedNumber[handle];
IF allNodes THEN
FOR n← handle.vars.nodeList, n.nextNode UNTIL n=NIL DO
n.nHist.y← v;
ENDLOOP
ELSE
IF n # NIL THEN n.nHist.y← v
ELSE
IF b # NIL THEN
WITH b.body SELECT FROM
l: RefL => l.iHist.y← v;
ENDCASE=> error[handle, 721]
ELSE error[handle, 720];
IF handle.vars.item=comma THEN next[handle] ELSE EXIT;
ENDLOOP;
IF handle.vars.item # rightB THEN error[handle, 701, TRUE] ELSE next[handle];
handle.vars.icsSet← TRUE;
}; -- setICs

runIt: PUBLIC PROC[handle: Handle] RETURNS[REAL]= {
n: INT;

IF handle.vars.item # leftB THEN error[handle, 700] ELSE next[handle];
runParms[handle];
IF handle.vars.item # rightB THEN error[handle, 701] ELSE next[handle];
IF handle.vars.tMax <= handle.vars.tMin THEN error[handle, 731];
IF handle.vars.errCount <= 0 THEN {
ENABLE Failure => {error[handle, errorNum]; GOTO quit};
IF ~handle.vars.icsSet THEN zeroVars[handle];
handle.vars.simTime← BasicTime.Now[];
initPrint[handle, handle.vars.simTime];
initPlot[handle, handle.vars.tMin, handle.vars.tMin + handle.vars.tMax];
handle.vars.t← handle.vars.tMin;
IF handle.vars.printdT <= 0.0 THEN
handle.vars.printdT← handle.vars.tMax/handle.vars.dfltInts;
handle.vars.dT← handle.vars.printdT;
n← GSiteration[handle, 0.0, handle.vars.printAll, TRUE, TRUE];
printFromList[handle, n, handle.vars.t, handle.vars.printIter];
plotFromList[handle, handle.vars.t];
IF ~getStepSize[handle] THEN RETURN[handle.vars.tMin];
handle.vars.numGoodSteps← 1;
WHILE handle.vars.highStepCount <= handle.vars.numberOfSteps DO
integration[handle, handle.vars.doRungeKutta];
IF Canned[handle] THEN EXIT;
IF CheckPoint[handle] THEN {
ResetCheckPoint[handle];
dumpAll[handle, handle.vars.t];
};
ENDLOOP;
handle.vars.icsSet← FALSE;
RETURN[handle.vars.t];
EXITS
quit=> RETURN[handle.vars.t];
}
ELSE RETURN[handle.vars.tMin]
}; -- runIt

END.

CHANGE LOG

Wilhelm, April 27, 1982 4:04 PM

Barth, 7-May-82 10:32:52 PDT

Chen, April 22, 1983 10:56 AM. Changed dvdtmax from 1E13 to 1E15.

Barth, July 11, 1983 1:42 PM

Chen, February 12, 1984 8:13 PM modified to support oldArgVector.

Chen, June 12, 1984 6:10:21 pm PDT, cedarized.