westerberg/ivpNondimensional/ivpStepN.dynTank.a4c

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;

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