File: [Cherry]<Thyme>Cedar5.2>System>spSolve.mesa
Last Edited by: SChen, May 17, 1985 9:09:04 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],
RealFns USING [SqRt],
RefText USING [Equal],
Rope USING [Equal, Length, ROPE],
spGlobals USING [Aborted, argList, argNodes, branchLinkPtr, branchPtr, Canned, CheckPoint, dumpAll, error, findNodeOrBranch, getSignedNumber, Handle, initPlot, initPrint, makeStringNB, 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, RealFns, 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= 29;
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", "ignorDvDtBelow", ""];
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, tempDvDt: 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 handle.msgStream.PutF["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 handle.msgStream.PutF["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;
tempDvDt← I/csum;
dvdt← n.nHist.f0← IF ABS[tempDvDt] < handle.vars.ignorDvDtBelow THEN 0.0 ELSE tempDvDt;
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 handle.msgStream.PutF["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);
n.nHist.f0← 0.0;
IF NOT firstTime THEN {
tempDvDt ← 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);
IF ABS[tempDvDt] > handle.vars.ignorDvDtBelow THEN n.nHist.f0← tempDvDt;
};
vCurBody.vsCurrent← newI;
IF print THEN handle.msgStream.PutF["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 handle.msgStream.PutF["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
IF nIter = realMaxIter THEN
SIGNAL solveError[IO.PutFR["Node %g GS did not converge.",
IO.rope[spGlobals.makeStringNB[handle, n, NIL]]]];
FOR bLink← n.branches, bLink.nextLink UNTIL bLink=NIL DO
bLink.otherNode.converged← FALSE;
ENDLOOP;
ENDLOOP; -- n
nIter← nIter + 1;
IF handle.vars.canned AND NOT handle.vars.forcedDump THEN SIGNAL Aborted;
IF nIter=handle.vars.maxIter THEN
handle.msgStream.PutF["t= %9.3f, iteration count exceeded",
IO.real[handle.vars.t]];
IF OK THEN EXIT;
ENDLOOP; -- nIter
}; -- GSiteration
printNodeEqn: PROC[handle: Handle]= {
n: nodePtr;
bLink: branchLinkPtr;
branch: branchPtr;
plus: BOOL;
FOR n← handle.vars.intNodeList, n.nextIntNode UNTIL n=NIL DO
handle.msgStream.PutF["\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 => handle.msgStream.PutF["(%13.6f - %13.6f)/%13.6f",
IO.real[bLink.otherNode.nHist.y],
IO.real[n.nHist.y],
IO.real[branch.comVal] ];
c: RefC => handle.msgStream.PutF["%13.6fx%13.6f",
IO.real[branch.comVal],
IO.real[bLink.otherNode.nHist.f0] ];
l: RefL =>
IF plus THEN handle.msgStream.PutF["-1x%13.6f/%13.6f",
IO.real[l.iHist.y],
IO.real[branch.comVal]]
ELSE handle.msgStream.PutF["%13.6f/%13.6f",
IO.real[l.iHist.y],
IO.real[branch.comVal]];
v: RefV =>
IF plus THEN handle.msgStream.PutF["-1x%13.6f",
IO.real[v.vsCurrent]]
ELSE handle.msgStream.PutF["%13.6f", IO.real[v.vsCurrent]];
i: RefI =>
IF plus THEN handle.msgStream.PutF["-1x%13.6f",
IO.real[branch.comVal]]
ELSE handle.msgStream.PutF["%13.6f",IO.real[branch.comVal]];
ENDCASE;
handle.msgStream.PutF[" %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];
27=> [ignorDvDtBelow, , ]← 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 AND deltay > 1.0E-30 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 > handle.vars.floor AND deltay > 1.0E-30 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;
IF handle.vars.countsPerStep=1 THEN
SIGNAL solveError[IO.PutFR["\nNode %g convergence problems.",
IO.rope[spGlobals.makeStringNB[handle, nodes, NIL]]]];
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
[visited: FALSE,
node: handle.vars.worstNode,
t: handle.vars.t,
v: handle.vars.worstNode.nHist.y,
dvdt: handle.vars.worstNode.nHist.f0]
ELSE [FALSE, NIL, 0.0, 0.0, 0.0];
handle.vars.curLog← (handle.vars.curLog + 1) MOD maxLog;
IF handle.vars.countsPerStep=1 THEN SIGNAL solveError["Step too small"];
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]]];
};
};
}; -- changeStepSize
saveState: PROC[handle: Handle]= {
inds: branchPtr← handle.vars.inductorList;
volts: branchPtr← handle.vars.vSourceList;
indsBody: RefL;
voltsBody: RefV;
FOR nodes: nodePtr← handle.vars.nodeList, nodes.nextNode UNTIL nodes=NIL DO
OPEN nodes.nHist;
SolveTrouble: PROC [explanation: Rope.ROPE] = {
SIGNAL solveTest[IO.PutFR["Node %g %g", IO.rope[spGlobals.makeStringNB[handle, nodes, NIL]], IO.rope[explanation]]];
};
oldy← y;
f4← f3; f3← f2; f2← f1; f1← f0;
SELECT TRUE FROM
y < handle.vars.Vmin => SolveTrouble["below Vmin"];
y > handle.vars.Vmax => SolveTrouble["above Vmax"];
ABS[f0] > handle.vars.dvdtmax => SolveTrouble["change rate too large"];
ENDCASE => NULL;
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;
IF ABS[nodes.nHist.f0] < handle.vars.ignorDvDtBelow THEN nodes.nHist.f0← 0.0;
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]= {
lastRetreatCause: Rope.ROPE← NIL;
nGS: INT;
out, ok, showRetreat: BOOL← FALSE;
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=> {
If you really want to see what the cause is, use the debugger to set showRetreat to TRUE.
IF ~Rope.Equal[cause, lastRetreatCause] THEN {
IF Rope.Length[cause] > 0 AND showRetreat THEN
handle.msgStream.PutF["%g at t= %12.6e.\n",
IO.rope[cause], IO.real[handle.vars.t]];
lastRetreatCause← cause;
IF handle.vars.countsPerStep = 1 THEN REJECT;
};
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← RealFns.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← RealFns.SqRt[RealFns.SqRt[ratio]];
changeStepSize[handle, ratio, FALSE];
handle.vars.AMcount← 0;
};
}
ELSE ratio← RealFns.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;
handle.msgStream.PutF["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.
Barth, February 8, 1985 3:17:34 pm PST
McCreight, April 4, 1985 1:57:23 pm PST
Chen, May 9, 1985 6:40:28 pm PDT:
1. RealOps.SqRt => RealFns.SqRt;
2. changed lastRetreatCause from global to local.
3. ignor dvdt when its magnitude is smaller than ignorDvDtBelow.