1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
|
REQUIRE "ivpNondimensional/ivpStepN.a4c";
(* The following is a simple tank model used to test ivpStep.a4c
This model contains the prototype of a set of models a user should
write when using ivpStep.a4c to solve initial value mixed DAE
systems. *)
(* ---------------------------------------------------------- *)
MODEL pastPoint REFINES integrationPoint;
(* This model is for saving all past points needed for
the fitting of the Taylor series *)
(* There is one state - the holdup in the tank. We will
need one predicted algebraic variable for a stop condition
based on too low a tank level. *)
nStates :== 1;
nPredVars :== 2;
METHODS
METHOD default_self;
RUN integrationPoint::default_self;
END default_self;
END pastPoint;
(* ---------------------------------------------------------- *)
MODEL currentPtDynTank REFINES pastPoint;
(* This model is for the current point. It requires the dynamic
and algebraic model equations for the model.
See the comments for MODEL integrationPoint.
Define states, state derivatives and dimensionless quantities
for the integration package. *)
t,tNom IS_A time;
xDefn: x*tNom = t;
dt,dtNom IS_A time;
dxDefn: dx*dtNom = dt;
holdup,holdupNom IS_A delta_mole;
y1Defn: y[1]*holdupNom = holdup;
dholdup_dt IS_A delta_molar_rate;
dydxDefn: dydx[1]*holdupNom = dholdup_dt*dtNom;
flowIn,flowOut,
flowStop,flowMin,
flowNom IS_A molar_rate;
Cv IS_A positive_factor; (* valve constant *)
y2Defn: y[2]*flowNom = flowStop;
eqnDynTankholdup: dholdup_dt = flowIn - flowOut;
eqnFlowIn: flowIn = 2.0{mol/s}*(1+sin(20.0{deg/s}*t));
eqnFlowOut: (flowOut/1{mol/s})^2 = Cv^2*(holdup/1.0{mol});
eqnStoppingCond: flowStop = flowOut - flowMin;
METHODS
(* ----- currentPtDynTank ----- *)
METHOD default_self;
RUN pastPoint::default_self;
END default_self;
(* ----- currentPtDynTank ----- *)
METHOD specifyForInitializing;
(* states, state derivatives and predicted algebraics *)
(* fix the time for the initial point *)
FIX t;
FIX tNom;
FREE x; (* computed using xDefn *)
(* fix the time step for the initial point *)
FIX dt;
FIX dtNom;
FREE dx; (* computed using dxDefn *)
(* fix the holdup for the initial point *)
FIX holdup;
FIX holdupNom;
FREE y[1]; (* computed using y1Defn *)
(* compute the time derivative for the holdup *)
FREE dholdup_dt; (* computed using eqnDynTankholdup *)
FREE dydx[1]; (* computed using dydxDefn *)
(* algebraic variables *)
FREE flowIn; (* computed by eqnFlowIn *)
FIX flowMin;
FREE flowStop; (* computed by eqnFlowStop *)
FIX flowNom;
FREE y[2]; (* computed by y2Defn *)
(* When initializing, we wish to compute the valve constant
that will give us a desired flowOut *)
FIX flowOut;
FREE Cv; (* computed by eqnFlowOut *)
END specifyForInitializing;
(* ----- currentPtDynTank ----- *)
METHOD valuesForInitializing;
(* provide values for all items whose values are fixed *)
t := 0.0{s};
tNom := 1.0 {s};
dt := 0.1{s}; (* value not used when initializing *)
dtNom := 1.0 {s};
holdup := 10.0{mol};
holdupNom := 1.0 {mol};
flowOut := 1.0{mol/s};
flowMin := 0.2{mol/s};
flowNom := 1.0 {mol/s};
END valuesForInitializing;
(* ----- currentPtDynTank ----- *)
METHOD specifyForStepping;
RUN specifyForInitializing;
(* to allow one to handle stopping conditions, the integration
package includes an equation to compute time for current point.
The dimensionless value for step size has to be fixed when
initializing for stepping. *)
FREE t; (* computed within integration package *)
FREE dt; (* computed by xDefn *)
FIX dx; (* value has to be fixed for stepping *)
(* taking a step is to compute the states, their derivatives
and the algebraic variables for the model - in general *)
FREE holdup; (* computed by the integration eqns *)
(* when initializing, we fixed flowOut and computed the
required valve constant, Cv. We need to reverse that for stepping *)
FREE flowOut;
FIX Cv;
END specifyForStepping;
(* ----- currentPtDynTank ----- *)
METHOD valuesForStepping;
(* there are no new values we need to set/estimate for stepping *)
END valuesForStepping;
METHOD testForIndexProblem;
(* this method should be run only for the instance of the
currentPt. *)
(* method testForIndexProblem will set all fixed flags for
the state variables y to TRUE and all the fixed flags
for the dydx and predicted algebraic variables to
FALSE. The user should assure that the fixed flags for
the remaining algebraic variables are set to make the
currentPt model square. *)
RUN integrationPoint::testForIndexProblem;
FIX Cv;
END testForIndexProblem;
END currentPtDynTank;
(* ---------------------------------------------------------- *)
MODEL ivpTest REFINES ivpStep;
currentPt IS_REFINED_TO currentPtDynTank;
iP[1..7] IS_REFINED_TO pastPoint;
deltaTime IS_A time;
stopTime IS_A time;
maxDeltaTime IS_A time;
deltaX*currentPt.dtNom = deltaTime;
maxDeltaX*currentPt.dtNom = maxDeltaTime;
stopX*currentPt.tNom = stopTime;
(* stopConds is a set containing the indices of the predicted
variables (y[stopConds]) whose sign change will cause the
integration to stop.*)
stopConds :== [2];
METHODS
(* ----- ivpTest ----- *)
METHOD default_self;
(* run first *)
RUN ivpStep::default_self;
RUN currentPt.default_self;
RUN iP[1..7].default_self;
END default_self;
(* ----- ivpTest ----- *)
METHOD valuesForInitializing;
(* run after default_self. *)
(* set values for initial point *)
RUN currentPt.valuesForInitializing;
(* after running this method,
run specifyForInitializing *)
END valuesForInitializing;
(* ----- ivpTest ----- *)
METHOD specifyForInitializing;
(* run after valuesForInitializing *)
(* set fixed flags to initialize currentPt *)
RUN currentPt.specifyForInitializing;
END specifyForInitializing;
(* after running specifyForInitializing,
solve the currentPt to set the initial conditions
for the problem*)
(* ----- ivpTest ----- *)
METHOD valuesForStepping;
(* run after solving for the initial conditions. *)
RUN ivpStep::values;
eps := 1.0e-6;
(* The values following are default values. The
script running the integration will likely
overwrite these using values the user will supply
when invoking the script. *)
deltaX := 0.01;
stopX := 2.0;
maxNominalSteppingError := 0.001;
END valuesForStepping;
(* ----- ivpTest ----- *)
METHOD specifyForStepping;
(* run after running the method valuesForStepping *)
RUN ivpStep::specify;
RUN currentPt.specifyForStepping;
FIX maxDeltaX;
FIX stopX;
END specifyForStepping;
(* at this point select method for solving (Adams-Moulton
or BDF and then run stepX.
Then to march forward one step at a time, carry out the
following
1. solve the next step
2. adjust polynomial order and method if needed
3. repeat from 1 to integrate. *)
END ivpTest;
(* ---------------------------------------------------------- *)
|