westerberg/ivpDAE/ivpStepN.a4c

REQUIRE "atoms.a4l";
REQUIRE "plot.a4l";

(*
 *  ivpStep.a4c
 *  by Arthur W. Westerberg
 *  Part of the ASCEND Library
 *  $Date: 02/07/28 $
 *  $Revision: 1.7 $
 *  $Author:  mthomas $
 *  $Source: /afs/cs.cmu.edu/project/ascend/Repository/models/ivpsystem.a4l,v $
 *
 *  This file is part of the ASCEND Modeling Library.
 *
 *  Copyright (C) 1994 - 1998 Carnegie Mellon University
 *
 *  The ASCEND Modeling Library is free software; you can redistribute
 *  it and/or modify it under the terms of the GNU General Public
 *  License as published by the Free Software Foundation; either
 *  version 2 of the License, or (at your option) any later version.
 *
 *  The ASCEND Modeling Library is distributed in hope that it
 *  will be useful, but WITHOUT ANY WARRANTY; without even the implied
 *  warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
 *  See the GNU General Public License for more details.
 *
 *  You should have received a copy of the GNU General Public License
 *  along with the program; if not, write to the Free Software
 *  Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139 USA.  Check
 *  the file named COPYING.
 *)

(*============================================================================

    I V P S T E P  .  A 4 C
    -----------------------------

    AUTHOR:      Arthur W. Westerberg

    DATES:       07/2002 - Original code (AWW)
                 07/2003 - Tested code.  Works for simple test problem.
		 (AWW)
		 12/2003 - Added Tcl script that runs complete simulation
		 (AWW)

    CONTENTS:	 Models for the numerical integration equations for
                 the two multistep methods: a bdf for stiff problems
		 and the Adams Moulton for non-stiff problems.  There
		 is also a framework within which one creates the
		 model for the physical system.  These models are for
		 taking one step of the independent variable.

============================================================================*)


(* ---------------------------------------------------------- *)

MODEL ivpBase() REFINES cmumodel();

    (* ----- ivpBase -----
       The base model for solving DAE initial value problems using the
       multistep methods: (a) a bdf method for stiff problems and (b)
       the Adams Moulton method for non-stiff problems.  A DAE ivp
       model has two types of sets of equations:

       (1) the model of the physical system at the k-th step

	      f(dydx(k), y(k), z(k), u(k), x(k)) = 0
	      h(y(k), z(k), u(k), x(k))          = 0
	      u(k) = u(x(k))

	  where y are the state, z the algebraic and u the forcing
	  variable vectors and x the scalar independent variable, and
	  where f is a vector of functions equal to the number of
	  states, h is a vector of algebraic equations equal in number
	  to the algebraic variables, and where one typically
	  specifies the forcing variables as functions only of the
	  independent variable.

      (2) the numerical multistep integration equations
      
	      y(k) = y(k-1) + s(x[k]-x[k-1];
                                y[k-1],y[k-2],...;
                                dydx[k],dydx[k-1],...)
	      
	  These equations relate the states at the current point to
	  the time derivatives at current and past points and states
	  at past points.  They are not related to the model of the
	  physical system but only to the numerical method being used.

       A common solving approach would be to solve the integration
       equations in an outside loop, and, for each iteration of that
       outside loop, use a Newton-based method to solve the physical
       system model for the state derivatives and the algebraic
       variables.  Often this approach solves the outer loop using a
       simple substitution method - which can be very slow to
       converge.  The following figure illustrates.
 

                       ------------------     
            u(x) ->---|   (inner loop)   |---- z  -->
                      | solve model eqns |
            ------>---|                  |-- dydx-->-  
           |           ------------------            |
           |                                         |
           |           ------------------            |
           |          |   (outer loop)   |           |
            --- y,x =-|  iterate  inte-  |---<-------
                      |   gration eqns   |
                      |                  |
                       ------------------


       By creating both the model and integration equations through
       ASCEND models, ASCEND can simultaneously converge the model
       and integration equations using its sparse Newton-based solver
       - thus making the solving much more efficient.  We wrote this
       model precisely to gain this efficiency.
       

       The integration module - degrees of freedom considerations
       
       For the q-th order bdf and Adams-Moulton (am) methods used
       here, we write a q-th order polynomial in the independent
       variable x for each state variable y[k][i].  This polynomial is
       written in the form of a Taylor Series expanded around the
       current point x[k].

       A q-th order polynomial has the q+1 polynomial coefficients in
       it a[k][i]/i!.  The a[k][i] are estimate the q-th derivative of
       the function at the current point, with the q-th term being an
       estimate of the error in the fit.  We determine these
       coefficients by fitting this polynomial and/or the first
       derivative of this polynomial to past y[j][i] and dy/dt[j][i],
       j<k, data, respectively, and to the computed value for
       dy/dt[k][i] at the current point.  We would need to write
       precisely q+1 such equations for this fitting.  The bdf and am
       methods differ in which equations we write to do this fitting.
       See comments in later in the code for the actual polynomials we
       write.

       We then use this polynomial evaluated at the current point to
       estimate the state y[k][i] at the current point x[k].

       Looking back at the two boxes in the previous figure, the lower
       integration box will have these q+1 equations in it, tabulated
       values for previous points, and the means to use past data and
       the current estimate for dy/dt[k][i] (the input we show above
       for this box) to fit the coefficients for each state variable.
       It will also write the polynomial one additional time to
       provide the estimate for y[k][i] at the current point, which is
       the output from this box.  It thus contains q+2 polynomial
       equations in it, in terms of q+1 unknown coefficients for each
       state.  The model equations introduce the states y[k][i] and their
       derivatives dy/dx[k][i] for the current point and the equation to
       compute the latter in terms of the former, often in the form

           dy/dt[k][i] = f(all y[k],t).  

       So it introduces, in principle two variables and one equation
       for each state.

       To summarize, the integration box introduces q+1 equations and
       q+2 variables per state while the model box introduces 1
       equation and 2 new variables.  The model is in square and
       solvable.


       Predictor behavior

       The polynomial we fit for the current point can easily be used
       to predict y[k+1][i], the values for the states i at the next
       time step.  We can also create polynomials to predict any of
       the non state variables in the problem.  We can include these
       non state variables in the integration model and use the
       current and q past values to determine the coefficients for a
       Taylor Series polynomial that passes through these points.
       This polynomial would only be useful for predicting as it play
       no role in solving the problem.  But often predicting one of
       the algebraic variables might improve the chance to converge
       the equations.


       Error control

       The package checks periodically to see discover the size step
       it should take to keep the error at that the user has
       prescribed.  It can also change the order of the polynomial up
       or down by one if such a move will allow for a larger step.

       We write the polynomial in the form of a Taylor Series around
       the current point.  So the error in estimating the current
       point, as a function of the step size we take, is the highest
       order term in that series.  We back compute the step size
       needed to obtain the error we allow and use it as the step size
       for the next integration step.  We keep this step size for the
       next q steps before checking again.  Changing the order of the
       polynomial can also affect the error.  This package will change
       the order up or down by one (i.e., to q-1 or q+1 (keeping in
       the limits of 1 and 6) if that change allows it to take a
       larger step.  The next to last term in the Taylor Series
       estimates the error when using a polynomial of order q-1.  We
       estimate the derivative for order q+1 as the difference in the
       q-th derivative between the current point and the previous
       point divided by the change in x between these two points -
       i.e., a divided difference which is essentially a free
       calculation.


       Stopping conditions

       We have created this package so one can stop whenever any of
       the variables we are predicting (all the states and whichever
       of the non states we select) crosses through zero.  In taking a
       normal step in this package, we use the error control criterion
       to fix the step size in the independent variable, x, prior to
       taking the next step.  If after we solve the model for this
       step, any one of the variables we have selected to use as
       stopping conditions changes sign, we fix its value to zero and
       resolve the entire model to compute the step size.  As more
       than one such variable may change sign in this final interval,
       we test after each solution to see if another variable changed
       sign in the reduced interval we just computed.  If so, we
       repeat by fixing this new variable to zero and resolving.  In
       this manner we will find the variable that changed sign at the
       earliest point in the final interval.

     *)

     (* ----- ivpBase -----
       TO USE THIS PACKAGE, the user must provide three models that
       are refinements of models we provide here.  The first two
       will be refinements of integrationPoint, while the third must be
       a refinement of ivpStep.  See the comments for those models.
       Also see the model ivpStep.testModel.a4c, which is a prototype
       for what the user must write.
       
       NOTE: THE MODELS IN THIS LIBRARY ARE WRITTEN USING DIMENSIONLESS
       QUANTITIES.  THE USER MUST CONVERT ANY DIMENSIONAL STATES
       AND PREDICTED ALGEBRAIC VARIABLES, AND ALL STOPPING CONDITIONS,
       ERRORS, ETC., TO DIMENSIONLESS VARIABLES FOR TYING INTO THESE
       MODELS.
       
     *)

    
END ivpBase;

(* ---------------------------------------------------------- *)

MODEL factorial(
    N         IS_A integer_constant;
    )
WHERE(
    N>=0;
    N<=12;  (* Note that 12! = 0.48E9, 2^31=2.1E9 *)
) REFINES ivpBase();
    
    (* create a vector of factorials, where fac[i] contains
      i!, where i goes from 0 to N *)

    fac[0..N] IS_A factor;

METHODS

    (* ----- factorial ----- *)

    METHOD default_self;
	fac[0] := 1;

	(* if N=0, the following loop produces no code *)
	FOR i IN [1..N] DO
	    fac[i] := i*fac[i-1];
	END FOR;
    END default_self;

    (* ----- factorial ----- *)

    METHOD specify;
	FIX fac[0..N];
    END specify;

END factorial;

(* ---------------------------------------------------------- *)

MODEL pTrow(
    set IS_A set OF integer_constant;
) REFINES ivpBase();

    (* This model is needed and instanced in
      "MODEL pascalTriangle."  No user should
      instance it by itself. *)

    col[set] IS_A factor;
    
END pTrow;

(* ---------------------------------------------------------- *)

MODEL pascalTriangle(
    dim IS_A integer_constant;
) REFINES ivpBase();
    
    (* Creates a Pascal Triangle of dimension "dim."
      A Pascal triangle gives the coefficients for the
      polynomials resulting from expanding (x-y)^q,
      for q = 0..dim.  For dim=5, the triangle we form
      here is
      
      row(1): 1  1  1  1  1  1
      row(2): 1  2  3  4  5
      row(3): 1  3  6 10
      row(4): 1  4 10
      row(5): 1  5
      row(6): 1
      
      The coefficients for q=4 for (x-y)^4 are then in
      the next to last place in each row, namely
      1, 4, 6, 4, 1
      as in
      (1)x^4 - (4)x^3y + (6)x^2y^2 - (4)xy^3 + (1)y^4.
      *)
 
    FOR i IN [0..dim] CREATE
	row[i] IS_A pTrow([0..dim-i]);
    END FOR;
    
    METHODS
    
    (* ----- pascalTriangle ----- *)

    METHOD default_self;
	FOR i IN [0..dim] DO
	    row[0].col[i] := 1;
	    row[i].col[0] := 1;
	END FOR;
	FOR i IN [1..dim] DO
	    FOR j IN [1..dim-i] DO
		row[i].col[j] := row[i-1].col[j]+row[i].col[j-1];
	    END FOR;
	END FOR;
    END default_self;

    (* ----- pascalTriangle ----- *)

    METHOD specify;
	FOR i IN [0..dim] DO
	    FOR j IN [0..dim-i] DO
		FIX row[i].col[j];
	    END FOR;
	END FOR;
    END specify;
    
END pascalTriangle;

(* ---------------------------------------------------------- *)
(* ---------------------------------------------------------- *)
(* ---------------------------------------------------------- *)
(* ---------------------------------------------------------- *)
(* ---------------------------------------------------------- *)
(* ---------------------------------------------------------- *)

MODEL yPolysDef(
    y         WILL_BE factor;
    dx        WILL_BE factor;
    dim       WILL_BE integer_constant;
    a[0..dim] WILL_BE factor;
) WHERE(
    (dim >= 3) == TRUE;
);

(*
  create a polynomial equation for y of order dim at the
  independent variable dx.
*)

    f                  IS_A factorial(dim);
    yTerm[0..dim]      IS_A factor;
    
    yTerm[dim] = a[dim]/f.fac[dim];
    FOR i IN [2..dim] CREATE
	yTermEqn[i-1]: yTerm[i-1] =
	    a[i-1]/f.fac[i-1] + dx*yTerm[i];
    END FOR;
    yTermEqn0: y = a[0] + dx*yTerm[1];
    
END yPolysDef;

(* ---------------------------------------------------------- *)

MODEL dydxPolysDef(
    dydx      WILL_BE factor;
    dx        WILL_BE factor;
    dim       WILL_BE integer_constant;
    a[0..dim] WILL_BE factor;
) WHERE(
    (dim >= 3) == TRUE;
);

(*
  create a polynomial equation for dydx of order dim at the
  independent variable dx.
*)

    f                  IS_A factorial(dim);
    dydxTerm[0..dim-1] IS_A factor;
    
    dydxTerm[dim] = a[dim]/f.fac[dim-1];
    FOR i IN [3..dim] CREATE
	dydxTermEqn[i-1]: dydxTerm[i-1] =
	    a[i-1]/f.fac[i-2] + dx*dydxTerm[i];
    END FOR;
    dydxTermEqn1: dydx = a[1] + dx*dydxTerm[2];

END dydxPolysDef;

(* ---------------------------------------------------------- *)

MODEL state(
    y             WILL_BE solver_var;
    dydx          WILL_BE solver_var;
    dim           WILL_BE integer_constant;
) WHERE(
    (dim >= 3) == TRUE;
);

(*
  Model to provide storage for the state variable y, its derivative
  dydt, and the  independent variable dx, written at current and
  past points.  Note that dx = 0 at the current point.  All variables
  are stored in dimensionless form (of type factor), by  dividing
  values by nominal values when storing or retrieving.
  
  It also provides polynomials of order dim at current and past points
  for the integration equations written "around" the current point.
 *)
  

    (* storage *)
    yStore[0..dim]     IS_A factor;
    dydxStore[0..dim]  IS_A factor;
    dxStore[0..dim]    IS_A factor;
    
    (* polynomials *)
    a[0..dim]          IS_A factor;
    FOR i IN [1..dim] CREATE
	yPolys[i]      IS_A yPolysDef(yStore[i],    dxStore[i], dim, a);
	dydxPolys[i-1] IS_A yPolysDef(dydxStore[i-1], dxStore[i-1], dim, a);
    END FOR;
    
    (* Note that a[0] and a[1] predict y and dydx at the current point *)
    yPredEqn:    y/y.nominal       = a[0];
    dydxPredEqn: dydx/dydx.nominal = a[1];
  
(*
  
  The solver can now choose which order and which method to use as
  follows for an order q method 
  
  Set coefficients a[q+1..dim] to zero and fix.
  
  To choose the AM method, set only the following included flags to TRUE
    polys.dydxTermEqn[0..q-1]
    polys.yTermEqn[1]

  To choose the BDF method, set only the following included flags to TRUE
    polys.dydxTermEqn[0]
    polys.yTermEqn[1..q]
    
  One will have written q+1 eqns in either case in the q+1 coefficients
  a[0..q].
  
  *)
    
    METHODS
    
    METHOD default_self;
    END default_self;

    METHOD specify;
    END specify;
    
    METHOD values;
    END values;

    METHOD setUseMethodAM;
    END setUseMethodAM;

    METHOD setUseMethodBDF;
    END setUseMethodBDF;

    METHOD incrementPolyOrder;
    END incrementPolyOrder;

    METHOD decrementPolyOrder;
    END decrementPolyOrder;
    
    METHOD stepX;
    END stepX;
    
END state;

(* ---------------------------------------------------------- *)
    
MODEL predictedAlgebraic(
    z WILL_BE solver_var;
    dim WILL_BE integer_constant;
) WHERE(
    (dim >= 3) == TRUE;
);

    (* storage *)
    zStore[0..dim]    IS_A factor;
    dxStore[0..dim]   IS_A factor;
    
    (* polynomials *)
    a[0..dim]         IS_A factor;
    FOR i IN [0..dim] CREATE
	zPolys[i]    IS_A yPolysDef(zStore[i],    dxStore[i], dim, a);
    END FOR;

(*    
  Note that algebraic prediction is only to improve initial guesses
  for z at the current point so z will not necessarily equal its
  converged current value.  Do not include a prediction equation.
  
  The solver can now choose which order to use as follows for an order
  q method
  
  Set coefficients a[q+1..dim] to zero and fix.

  Set only the following include flags to TRUE
    polys.yTermEqn[0..q]
  
  One will have written q+1 eqns in the q+1 coefficients a[0..q].
*)

END predictedAlgebraic;

(* ---------------------------------------------------------- *)

MODEL myDynamicModel;
    
    maxPolyOrder          IS_A integer_constant;
    maxPolyOrder          :== 6;

    t,dt                  IS_A time;
    holdup                IS_A mole;
    dholdup_dt            IS_A delta_molar_rate;
    flowIn, flowOut,
    flowStop, flowMin     IS_A molar_rate;
    Cv                    IS_A positive_factor; (* valve constant *)
    
    holdupState           IS_A state(holdup, dholdup_dt, maxPolyOrder);
    predFlowStop          IS_A predictedAlgebraic(flowStop, maxPolyOrder);

    eqnDynTankholdup:     dholdup_dt = flowIn - flowOut;
    eqnFlowOut:           (flowOut/1{mol/s})^2 = Cv^2*(holdup/1.0{mol});
    eqnStoppingCond:      flowStop = flowOut - flowMin;

    METHODS

    METHOD default_self;
    END default_self;
	
    METHOD specifyForInitializing;
    END specifyForInitializing;
    
    METHOD valuesForInitializing;
    END valuesForInitializing;
    
    METHOD specifyForStepping;
    END specifyForStepping;
    
    METHOD valuesForStepping;
    END valuesForStepping;

    METHOD testForIndexProblem;
    END testForIndexProblem;

END myDynamicModel;

(* ---------------------------------------------------------- *)
1	REQUIRE "atoms.a4l";
2	REQUIRE "plot.a4l";
3
4	(*
5	* ivpStep.a4c
6	* by Arthur W. Westerberg
7	* Part of the ASCEND Library
8	* $Date: 02/07/28 $
9	* $Revision: 1.7 $
10	* $Author: mthomas $
11	* $Source: /afs/cs.cmu.edu/project/ascend/Repository/models/ivpsystem.a4l,v $
12	*
13	* This file is part of the ASCEND Modeling Library.
14	*
15	* Copyright (C) 1994 - 1998 Carnegie Mellon University
16	*
17	* The ASCEND Modeling Library is free software; you can redistribute
18	* it and/or modify it under the terms of the GNU General Public
19	* License as published by the Free Software Foundation; either
20	* version 2 of the License, or (at your option) any later version.
21	*
22	* The ASCEND Modeling Library is distributed in hope that it
23	* will be useful, but WITHOUT ANY WARRANTY; without even the implied
24	* warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
25	* See the GNU General Public License for more details.
26	*
27	* You should have received a copy of the GNU General Public License
28	* along with the program; if not, write to the Free Software
29	* Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139 USA. Check
30	* the file named COPYING.
31	*)
32
33	(*============================================================================
34
35	I V P S T E P . A 4 C
36	-----------------------------
37
38	AUTHOR: Arthur W. Westerberg
39
40	DATES: 07/2002 - Original code (AWW)
41	07/2003 - Tested code. Works for simple test problem.
42	(AWW)
43	12/2003 - Added Tcl script that runs complete simulation
44	(AWW)
45
46	CONTENTS: Models for the numerical integration equations for
47	the two multistep methods: a bdf for stiff problems
48	and the Adams Moulton for non-stiff problems. There
49	is also a framework within which one creates the
50	model for the physical system. These models are for
51	taking one step of the independent variable.
52
53	============================================================================*)
54
55
56	(* ---------------------------------------------------------- *)
57
58	MODEL ivpBase() REFINES cmumodel();
59
60	(* ----- ivpBase -----
61	The base model for solving DAE initial value problems using the
62	multistep methods: (a) a bdf method for stiff problems and (b)
63	the Adams Moulton method for non-stiff problems. A DAE ivp
64	model has two types of sets of equations:
65
66	(1) the model of the physical system at the k-th step
67
68	f(dydx(k), y(k), z(k), u(k), x(k)) = 0
69	h(y(k), z(k), u(k), x(k)) = 0
70	u(k) = u(x(k))
71
72	where y are the state, z the algebraic and u the forcing
73	variable vectors and x the scalar independent variable, and
74	where f is a vector of functions equal to the number of
75	states, h is a vector of algebraic equations equal in number
76	to the algebraic variables, and where one typically
77	specifies the forcing variables as functions only of the
78	independent variable.
79
80	(2) the numerical multistep integration equations
81
82	y(k) = y(k-1) + s(x[k]-x[k-1];
83	y[k-1],y[k-2],...;
84	dydx[k],dydx[k-1],...)
85
86	These equations relate the states at the current point to
87	the time derivatives at current and past points and states
88	at past points. They are not related to the model of the
89	physical system but only to the numerical method being used.
90
91	A common solving approach would be to solve the integration
92	equations in an outside loop, and, for each iteration of that
93	outside loop, use a Newton-based method to solve the physical
94	system model for the state derivatives and the algebraic
95	variables. Often this approach solves the outer loop using a
96	simple substitution method - which can be very slow to
97	converge. The following figure illustrates.
98
99
100	------------------
101	u(x) ->---\| (inner loop) \|---- z -->
102	\| solve model eqns \|
103	------>---\| \|-- dydx-->-
104	\| ------------------ \|
105	\| \|
106	\| ------------------ \|
107	\| \| (outer loop) \| \|
108	--- y,x =-\| iterate inte- \|---<-------
109	\| gration eqns \|
110	\| \|
111	------------------
112
113
114	By creating both the model and integration equations through
115	ASCEND models, ASCEND can simultaneously converge the model
116	and integration equations using its sparse Newton-based solver
117	- thus making the solving much more efficient. We wrote this
118	model precisely to gain this efficiency.
119
120
121	The integration module - degrees of freedom considerations
122
123	For the q-th order bdf and Adams-Moulton (am) methods used
124	here, we write a q-th order polynomial in the independent
125	variable x for each state variable y[k][i]. This polynomial is
126	written in the form of a Taylor Series expanded around the
127	current point x[k].
128
129	A q-th order polynomial has the q+1 polynomial coefficients in
130	it a[k][i]/i!. The a[k][i] are estimate the q-th derivative of
131	the function at the current point, with the q-th term being an
132	estimate of the error in the fit. We determine these
133	coefficients by fitting this polynomial and/or the first
134	derivative of this polynomial to past y[j][i] and dy/dt[j][i],
135	j<k, data, respectively, and to the computed value for
136	dy/dt[k][i] at the current point. We would need to write
137	precisely q+1 such equations for this fitting. The bdf and am
138	methods differ in which equations we write to do this fitting.
139	See comments in later in the code for the actual polynomials we
140	write.
141
142	We then use this polynomial evaluated at the current point to
143	estimate the state y[k][i] at the current point x[k].
144
145	Looking back at the two boxes in the previous figure, the lower
146	integration box will have these q+1 equations in it, tabulated
147	values for previous points, and the means to use past data and
148	the current estimate for dy/dt[k][i] (the input we show above
149	for this box) to fit the coefficients for each state variable.
150	It will also write the polynomial one additional time to
151	provide the estimate for y[k][i] at the current point, which is
152	the output from this box. It thus contains q+2 polynomial
153	equations in it, in terms of q+1 unknown coefficients for each
154	state. The model equations introduce the states y[k][i] and their
155	derivatives dy/dx[k][i] for the current point and the equation to
156	compute the latter in terms of the former, often in the form
157
158	dy/dt[k][i] = f(all y[k],t).
159
160	So it introduces, in principle two variables and one equation
161	for each state.
162
163	To summarize, the integration box introduces q+1 equations and
164	q+2 variables per state while the model box introduces 1
165	equation and 2 new variables. The model is in square and
166	solvable.
167
168
169	Predictor behavior
170
171	The polynomial we fit for the current point can easily be used
172	to predict y[k+1][i], the values for the states i at the next
173	time step. We can also create polynomials to predict any of
174	the non state variables in the problem. We can include these
175	non state variables in the integration model and use the
176	current and q past values to determine the coefficients for a
177	Taylor Series polynomial that passes through these points.
178	This polynomial would only be useful for predicting as it play
179	no role in solving the problem. But often predicting one of
180	the algebraic variables might improve the chance to converge
181	the equations.
182
183
184	Error control
185
186	The package checks periodically to see discover the size step
187	it should take to keep the error at that the user has
188	prescribed. It can also change the order of the polynomial up
189	or down by one if such a move will allow for a larger step.
190
191	We write the polynomial in the form of a Taylor Series around
192	the current point. So the error in estimating the current
193	point, as a function of the step size we take, is the highest
194	order term in that series. We back compute the step size
195	needed to obtain the error we allow and use it as the step size
196	for the next integration step. We keep this step size for the
197	next q steps before checking again. Changing the order of the
198	polynomial can also affect the error. This package will change
199	the order up or down by one (i.e., to q-1 or q+1 (keeping in
200	the limits of 1 and 6) if that change allows it to take a
201	larger step. The next to last term in the Taylor Series
202	estimates the error when using a polynomial of order q-1. We
203	estimate the derivative for order q+1 as the difference in the
204	q-th derivative between the current point and the previous
205	point divided by the change in x between these two points -
206	i.e., a divided difference which is essentially a free
207	calculation.
208
209
210	Stopping conditions
211
212	We have created this package so one can stop whenever any of
213	the variables we are predicting (all the states and whichever
214	of the non states we select) crosses through zero. In taking a
215	normal step in this package, we use the error control criterion
216	to fix the step size in the independent variable, x, prior to
217	taking the next step. If after we solve the model for this
218	step, any one of the variables we have selected to use as
219	stopping conditions changes sign, we fix its value to zero and
220	resolve the entire model to compute the step size. As more
221	than one such variable may change sign in this final interval,
222	we test after each solution to see if another variable changed
223	sign in the reduced interval we just computed. If so, we
224	repeat by fixing this new variable to zero and resolving. In
225	this manner we will find the variable that changed sign at the
226	earliest point in the final interval.
227
228	*)
229
230	(* ----- ivpBase -----
231	TO USE THIS PACKAGE, the user must provide three models that
232	are refinements of models we provide here. The first two
233	will be refinements of integrationPoint, while the third must be
234	a refinement of ivpStep. See the comments for those models.
235	Also see the model ivpStep.testModel.a4c, which is a prototype
236	for what the user must write.
237
238	NOTE: THE MODELS IN THIS LIBRARY ARE WRITTEN USING DIMENSIONLESS
239	QUANTITIES. THE USER MUST CONVERT ANY DIMENSIONAL STATES
240	AND PREDICTED ALGEBRAIC VARIABLES, AND ALL STOPPING CONDITIONS,
241	ERRORS, ETC., TO DIMENSIONLESS VARIABLES FOR TYING INTO THESE
242	MODELS.
243
244	*)
245
246
247	END ivpBase;
248
249	(* ---------------------------------------------------------- *)
250
251	MODEL factorial(
252	N IS_A integer_constant;
253	)
254	WHERE(
255	N>=0;
256	N<=12; (* Note that 12! = 0.48E9, 2^31=2.1E9 *)
257	) REFINES ivpBase();
258
259	(* create a vector of factorials, where fac[i] contains
260	i!, where i goes from 0 to N *)
261
262	fac[0..N] IS_A factor;
263
264	METHODS
265
266	(* ----- factorial ----- *)
267
268	METHOD default_self;
269	fac[0] := 1;
270
271	(* if N=0, the following loop produces no code *)
272	FOR i IN [1..N] DO
273	fac[i] := i*fac[i-1];
274	END FOR;
275	END default_self;
276
277	(* ----- factorial ----- *)
278
279	METHOD specify;
280	FIX fac[0..N];
281	END specify;
282
283	END factorial;
284
285	(* ---------------------------------------------------------- *)
286
287	MODEL pTrow(
288	set IS_A set OF integer_constant;
289	) REFINES ivpBase();
290
291	(* This model is needed and instanced in
292	"MODEL pascalTriangle." No user should
293	instance it by itself. *)
294
295	col[set] IS_A factor;
296
297	END pTrow;
298
299	(* ---------------------------------------------------------- *)
300
301	MODEL pascalTriangle(
302	dim IS_A integer_constant;
303	) REFINES ivpBase();
304
305	(* Creates a Pascal Triangle of dimension "dim."
306	A Pascal triangle gives the coefficients for the
307	polynomials resulting from expanding (x-y)^q,
308	for q = 0..dim. For dim=5, the triangle we form
309	here is
310
311	row(1): 1 1 1 1 1 1
312	row(2): 1 2 3 4 5
313	row(3): 1 3 6 10
314	row(4): 1 4 10
315	row(5): 1 5
316	row(6): 1
317
318	The coefficients for q=4 for (x-y)^4 are then in
319	the next to last place in each row, namely
320	1, 4, 6, 4, 1
321	as in
322	(1)x^4 - (4)x^3y + (6)x^2y^2 - (4)xy^3 + (1)y^4.
323	*)
324
325	FOR i IN [0..dim] CREATE
326	row[i] IS_A pTrow([0..dim-i]);
327	END FOR;
328
329	METHODS
330
331	(* ----- pascalTriangle ----- *)
332
333	METHOD default_self;
334	FOR i IN [0..dim] DO
335	row[0].col[i] := 1;
336	row[i].col[0] := 1;
337	END FOR;
338	FOR i IN [1..dim] DO
339	FOR j IN [1..dim-i] DO
340	row[i].col[j] := row[i-1].col[j]+row[i].col[j-1];
341	END FOR;
342	END FOR;
343	END default_self;
344
345	(* ----- pascalTriangle ----- *)
346
347	METHOD specify;
348	FOR i IN [0..dim] DO
349	FOR j IN [0..dim-i] DO
350	FIX row[i].col[j];
351	END FOR;
352	END FOR;
353	END specify;
354
355	END pascalTriangle;
356
357	(* ---------------------------------------------------------- *)
358	(* ---------------------------------------------------------- *)
359	(* ---------------------------------------------------------- *)
360	(* ---------------------------------------------------------- *)
361	(* ---------------------------------------------------------- *)
362	(* ---------------------------------------------------------- *)
363
364	MODEL yPolysDef(
365	y WILL_BE factor;
366	dx WILL_BE factor;
367	dim WILL_BE integer_constant;
368	a[0..dim] WILL_BE factor;
369	) WHERE(
370	(dim >= 3) == TRUE;
371	);
372
373	(*
374	create a polynomial equation for y of order dim at the
375	independent variable dx.
376	*)
377
378	f IS_A factorial(dim);
379	yTerm[0..dim] IS_A factor;
380
381	yTerm[dim] = a[dim]/f.fac[dim];
382	FOR i IN [2..dim] CREATE
383	yTermEqn[i-1]: yTerm[i-1] =
384	a[i-1]/f.fac[i-1] + dx*yTerm[i];
385	END FOR;
386	yTermEqn0: y = a[0] + dx*yTerm[1];
387
388	END yPolysDef;
389
390	(* ---------------------------------------------------------- *)
391
392	MODEL dydxPolysDef(
393	dydx WILL_BE factor;
394	dx WILL_BE factor;
395	dim WILL_BE integer_constant;
396	a[0..dim] WILL_BE factor;
397	) WHERE(
398	(dim >= 3) == TRUE;
399	);
400
401	(*
402	create a polynomial equation for dydx of order dim at the
403	independent variable dx.
404	*)
405
406	f IS_A factorial(dim);
407	dydxTerm[0..dim-1] IS_A factor;
408
409	dydxTerm[dim] = a[dim]/f.fac[dim-1];
410	FOR i IN [3..dim] CREATE
411	dydxTermEqn[i-1]: dydxTerm[i-1] =
412	a[i-1]/f.fac[i-2] + dx*dydxTerm[i];
413	END FOR;
414	dydxTermEqn1: dydx = a[1] + dx*dydxTerm[2];
415
416	END dydxPolysDef;
417
418	(* ---------------------------------------------------------- *)
419
420	MODEL state(
421	y WILL_BE solver_var;
422	dydx WILL_BE solver_var;
423	dim WILL_BE integer_constant;
424	) WHERE(
425	(dim >= 3) == TRUE;
426	);
427
428	(*
429	Model to provide storage for the state variable y, its derivative
430	dydt, and the independent variable dx, written at current and
431	past points. Note that dx = 0 at the current point. All variables
432	are stored in dimensionless form (of type factor), by dividing
433	values by nominal values when storing or retrieving.
434
435	It also provides polynomials of order dim at current and past points
436	for the integration equations written "around" the current point.
437	*)
438
439
440	(* storage *)
441	yStore[0..dim] IS_A factor;
442	dydxStore[0..dim] IS_A factor;
443	dxStore[0..dim] IS_A factor;
444
445	(* polynomials *)
446	a[0..dim] IS_A factor;
447	FOR i IN [1..dim] CREATE
448	yPolys[i] IS_A yPolysDef(yStore[i], dxStore[i], dim, a);
449	dydxPolys[i-1] IS_A yPolysDef(dydxStore[i-1], dxStore[i-1], dim, a);
450	END FOR;
451
452	(* Note that a[0] and a[1] predict y and dydx at the current point *)
453	yPredEqn: y/y.nominal = a[0];
454	dydxPredEqn: dydx/dydx.nominal = a[1];
455
456	(*
457
458	The solver can now choose which order and which method to use as
459	follows for an order q method
460
461	Set coefficients a[q+1..dim] to zero and fix.
462
463	To choose the AM method, set only the following included flags to TRUE
464	polys.dydxTermEqn[0..q-1]
465	polys.yTermEqn[1]
466
467	To choose the BDF method, set only the following included flags to TRUE
468	polys.dydxTermEqn[0]
469	polys.yTermEqn[1..q]
470
471	One will have written q+1 eqns in either case in the q+1 coefficients
472	a[0..q].
473
474	*)
475
476	METHODS
477
478	METHOD default_self;
479	END default_self;
480
481	METHOD specify;
482	END specify;
483
484	METHOD values;
485	END values;
486
487	METHOD setUseMethodAM;
488	END setUseMethodAM;
489
490	METHOD setUseMethodBDF;
491	END setUseMethodBDF;
492
493	METHOD incrementPolyOrder;
494	END incrementPolyOrder;
495
496	METHOD decrementPolyOrder;
497	END decrementPolyOrder;
498
499	METHOD stepX;
500	END stepX;
501
502	END state;
503
504	(* ---------------------------------------------------------- *)
505
506	MODEL predictedAlgebraic(
507	z WILL_BE solver_var;
508	dim WILL_BE integer_constant;
509	) WHERE(
510	(dim >= 3) == TRUE;
511	);
512
513	(* storage *)
514	zStore[0..dim] IS_A factor;
515	dxStore[0..dim] IS_A factor;
516
517	(* polynomials *)
518	a[0..dim] IS_A factor;
519	FOR i IN [0..dim] CREATE
520	zPolys[i] IS_A yPolysDef(zStore[i], dxStore[i], dim, a);
521	END FOR;
522
523	(*
524	Note that algebraic prediction is only to improve initial guesses
525	for z at the current point so z will not necessarily equal its
526	converged current value. Do not include a prediction equation.
527
528	The solver can now choose which order to use as follows for an order
529	q method
530
531	Set coefficients a[q+1..dim] to zero and fix.
532
533	Set only the following include flags to TRUE
534	polys.yTermEqn[0..q]
535
536	One will have written q+1 eqns in the q+1 coefficients a[0..q].
537	*)
538
539	END predictedAlgebraic;
540
541	(* ---------------------------------------------------------- *)
542
543	MODEL myDynamicModel;
544
545	maxPolyOrder IS_A integer_constant;
546	maxPolyOrder :== 6;
547
548	t,dt IS_A time;
549	holdup IS_A mole;
550	dholdup_dt IS_A delta_molar_rate;
551	flowIn, flowOut,
552	flowStop, flowMin IS_A molar_rate;
553	Cv IS_A positive_factor; (* valve constant *)
554
555	holdupState IS_A state(holdup, dholdup_dt, maxPolyOrder);
556	predFlowStop IS_A predictedAlgebraic(flowStop, maxPolyOrder);
557
558	eqnDynTankholdup: dholdup_dt = flowIn - flowOut;
559	eqnFlowOut: (flowOut/1{mol/s})^2 = Cv^2*(holdup/1.0{mol});
560	eqnStoppingCond: flowStop = flowOut - flowMin;
561
562	METHODS
563
564	METHOD default_self;
565	END default_self;
566
567	METHOD specifyForInitializing;
568	END specifyForInitializing;
569
570	METHOD valuesForInitializing;
571	END valuesForInitializing;
572
573	METHOD specifyForStepping;
574	END specifyForStepping;
575
576	METHOD valuesForStepping;
577	END valuesForStepping;
578
579	METHOD testForIndexProblem;
580	END testForIndexProblem;
581
582	END myDynamicModel;
583
584	(* ---------------------------------------------------------- *)