/[ascend]/trunk/base/generic/solver/diffvars.c
ViewVC logotype

Contents of /trunk/base/generic/solver/diffvars.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1315 - (show annotations) (download) (as text)
Mon Mar 5 06:04:15 2007 UTC (15 years, 3 months ago) by johnpye
File MIME type: text/x-csrc
File size: 8380 byte(s)
Adding flagging of differential variables (those which have derivatives) via var.h.
1 /* ASCEND modelling environment
2 Copyright (C) 2007 Carnegie Mellon University
3
4 This program is free software; you can redistribute it and/or modify
5 it under the terms of the GNU General Public License as published by
6 the Free Software Foundation; either version 2, or (at your option)
7 any later version.
8
9 This program is distributed in the hope that it will be useful,
10 but WITHOUT ANY WARRANTY; without even the implied warranty of
11 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12 GNU General Public License for more details.
13
14 You should have received a copy of the GNU General Public License
15 along with this program; if not, write to the Free Software
16 Foundation, Inc., 59 Temple Place - Suite 330,
17 Boston, MA 02111-1307, USA.
18 *//** @file
19 Routines for creating 'derivative chains' during system build.
20 *//*
21 By John Pye, Jan 2007
22 */
23 #include "diffvars.h"
24
25 #include <utilities/ascMalloc.h>
26 #include <utilities/ascPanic.h>
27 #include <utilities/error.h>
28
29 #include <compiler/instance_io.h>
30
31 #include "analyse_impl.h"
32 #include "system_impl.h"
33
34 /* #define DIFFVARS_DEBUG */
35
36 /*------------------------------------------------------------------------------
37 DERIVATIVES & DIFFERENTIAL VARIABLES
38 */
39
40 SolverDiffVarCollection *system_get_diffvars(slv_system_t sys){
41 return sys->diffvars;
42 }
43
44 static
45 int CmpDiffVars(const struct solver_ipdata *a, const struct solver_ipdata *b){
46 if(a->u.v.odeid < b->u.v.odeid)return -1;
47 if(a->u.v.odeid > b->u.v.odeid)return 1;
48 /* now, the ode_id is equal: next level of sorting... */
49 if(a->u.v.deriv < b->u.v.deriv)return -1;
50 if(a->u.v.deriv > b->u.v.deriv)return 1;
51 return 0;
52 }
53
54 /**
55 This function steals a lot of what was in integrator.c, but we try to make
56 it more general and capable of extracting info about high-order derivative
57 chains. One eye open to the possibility of doing Pantelides (or other?)
58 style index reduction.
59
60 Also, note that we try to only use data from the struct problem_t, NOT
61 to go back to the instance hierarchy. The idea is that we have already got
62 what we need in classify_instance, etc.
63
64 @return 0 on success
65 */
66 int system_generate_diffvars(slv_system_t sys, struct problem_t *prob){
67 SolverDiffVarCollection *diffvars = NULL;
68 struct solver_ipdata *vip, *vipnext;
69 SolverDiffVarSequence *seq;
70 struct gl_list_t *seqs;
71 long i, seqstart, nalg, ndiff;
72 short j;
73 char cont;
74 short maxorder = 0;
75
76 asc_assert(prob);
77
78 if(gl_length(prob->diffvars)==0){
79 CONSOLE_DEBUG("No differential variables were seen. Skipping generation of diffvars struct.");
80 return 0;
81 }
82
83 seqs = gl_create(prob->nr);
84
85 /* add the list of algebraic variables to the structure too */
86 asc_assert(prob->algebvars);
87 nalg = gl_length(prob->algebvars);
88 for(i=0; i<nalg; ++i){
89 vip = (struct solver_ipdata *)gl_fetch(prob->algebvars,i+1);
90 seq = ASC_NEW(SolverDiffVarSequence);
91 seq->n = 1;
92 seq->ode_id = 0;
93 seq->vars = ASC_NEW_ARRAY(struct var_variable*,1);
94 seq->vars[0] = vip->u.v.data;
95 gl_append_ptr(seqs,(void *)seq);
96 }
97 #ifdef DIFFVARS_DEBUG
98 CONSOLE_DEBUG("Added %ld algebraic vars to chains",nalg);
99 CONSOLE_DEBUG("Sorting %ld differential & derivative vars...",gl_length(prob->diffvars));
100 #endif
101
102 /* first sort the list of diffvars */
103 gl_sort(prob->diffvars, (CmpFunc)CmpDiffVars);
104
105 /* scan the list for derivs that are missing vars */
106 i = 1; cont = TRUE; seqstart =1; ndiff = 0;
107 vip = (struct solver_ipdata *)gl_fetch(prob->diffvars,i);
108 if(vip->u.v.deriv > 1){
109 ERROR_REPORTER_START_NOLINE(ASC_USER_ERROR);
110 FPRINTF(ASCERR,"Missing ode_type %d for ode_id %d (check var '",vip->u.v.deriv+1,vip->u.v.odeid);
111 WriteInstanceName(ASCERR,vip->i,prob->root);
112 FPRINTF(ASCERR,"')");
113 error_reporter_end_flush();
114 return 2;
115 }
116 while(cont){
117 /* CONSOLE_DEBUG("Working, seqstart=%ld",seqstart); */
118 if(i >= gl_length(prob->diffvars))cont = FALSE;
119 else vipnext = (struct solver_ipdata *)gl_fetch(prob->diffvars,i+1);
120
121 if(cont && vipnext->u.v.odeid == vip->u.v.odeid){
122 /* same sequence, check that it's the next derivative */
123 asc_assert(vip->u.v.deriv >= vip->u.v.deriv);
124
125 if(vipnext->u.v.deriv == vip->u.v.deriv){
126 ERROR_REPORTER_START_NOLINE(ASC_USER_ERROR);
127 FPRINTF(ASCERR,"Repeated ode_type %d for ode_id %d (var '%s')",vip->u.v.deriv,vip->u.v.odeid);
128 WriteInstanceName(ASCERR,vip->i,prob->root);
129 FPRINTF(ASCERR,"')");
130 error_reporter_end_flush();
131 return 1;
132 }else if(vipnext->u.v.deriv > vip->u.v.deriv + 1){
133 ERROR_REPORTER_START_NOLINE(ASC_USER_ERROR);
134 FPRINTF(ASCERR,"Missing ode_type %d for ode_id %d (check var '",vip->u.v.deriv+1,vip->u.v.odeid);
135 WriteInstanceName(ASCERR,vip->i,prob->root);
136 FPRINTF(ASCERR,"')");
137 error_reporter_end_flush();
138 return 2;
139 }
140
141 /* allow this var into the sequence, loop again */
142 vip = vipnext; ++i; continue;
143 }
144
145 /* close the current sequence */
146 if(i == seqstart){
147 ERROR_REPORTER_START_NOLINE(ASC_USER_ERROR);
148 FPRINTF(ASCERR,"Variable '");
149 WriteInstanceName(ASCWAR,vip->i,prob->root);
150 FPRINTF(ASCERR,"' declared differential without derivative being identified (ode_id=%d)",vip->u.v.odeid);
151 error_reporter_end_flush();
152 return 3;
153 }
154 seq = ASC_NEW(SolverDiffVarSequence);
155 seq->n = i - seqstart + 1;
156 seq->ode_id = vip->u.v.odeid;
157 seq->vars = ASC_NEW_ARRAY(struct var_variable*,seq->n);
158 #ifdef DIFFVARS_DEBUG
159 CONSOLE_DEBUG("Saving sequence ode_id = %ld, n = %d",seq->ode_id,seq->n);
160 #endif
161 for(j=0;j<seq->n;++j){
162 vip = (struct solver_ipdata *)gl_fetch(prob->diffvars,seqstart+j);
163 seq->vars[j]=vip->u.v.data;
164 /* set the VAR_ALGEB flag as req */
165 var_set_diff(seq->vars[j],j==0 && seq->n > 1);
166 #ifdef DIFFVARS_DEBUG
167 CONSOLE_DEBUG("seq, j=%d: is %s",j
168 ,(var_diff(seq->vars[j]) ? "DIFF" :
169 (var_deriv(seq->vars[j]) ? "deriv" : "alg")
170 )
171 );
172 #endif
173 }
174 gl_append_ptr(seqs,(void *)seq);
175 ndiff++;
176 #ifdef DIFFVARS_DEBUG
177 CONSOLE_DEBUG("Completed seq %ld",gl_length(seqs));
178 #endif
179
180 if(vip->u.v.deriv > maxorder){
181 maxorder = vip->u.v.deriv;
182 /* CONSOLE_DEBUG("maxorder increased to %d",maxorder); */
183 }
184
185 ++i; seqstart=i; vip=vipnext;
186 if(cont && vip->u.v.deriv != 1){
187 ERROR_REPORTER_START_NOLINE(ASC_USER_ERROR);
188 FPRINTF(ASCERR,"Missing undifferentiated form of var '");
189 WriteInstanceName(ASCWAR,vip->i,prob->root);
190 FPRINTF(ASCERR,"' (ode_id = %d, ode_type = %d)",vip->u.v.odeid,vip->u.v.deriv);
191 error_reporter_end_flush();
192 return 4;
193 }
194 continue;
195 }
196
197 #ifdef DIFFVARS_DEBUG
198 CONSOLE_DEBUG("Identified %ld derivative chains, maximum length %d...",gl_length(seqs),maxorder);
199 #endif
200
201 diffvars = ASC_NEW(SolverDiffVarCollection);
202 diffvars->maxorder = maxorder;
203 diffvars->nseqs = gl_length(seqs);
204 diffvars->seqs = ASC_NEW_ARRAY(SolverDiffVarSequence,diffvars->nseqs);
205 for(i=0;i<diffvars->nseqs;++i){
206 diffvars->seqs[i] = *(SolverDiffVarSequence *)gl_fetch(seqs,i+1);
207 }
208 gl_destroy(seqs);
209 diffvars->nalg = nalg;
210 diffvars->ndiff = ndiff;
211
212 diffvars->nindep = gl_length(prob->indepvars);
213 diffvars->indep = ASC_NEW_ARRAY(struct var_variable *,diffvars->nindep);
214 for(i=0;i<diffvars->nindep;++i){
215 vip = (struct solver_ipdata *)gl_fetch(prob->indepvars,i+1);
216 diffvars->indep[i] = vip->u.v.data;
217 }
218 #ifdef DIFFVARS_DEBUG
219 CONSOLE_DEBUG("Identified %ld indep vars",diffvars->nindep);
220 #endif
221
222 /* I know, I know, this might not be the best place for this... */
223 diffvars->nobs = gl_length(prob->obsvars);
224 diffvars->obs = ASC_NEW_ARRAY(struct var_variable *,diffvars->nobs);
225 for(i=0;i<diffvars->nobs;++i){
226 vip = (struct solver_ipdata *)gl_fetch(prob->obsvars,i+1);
227 diffvars->obs[i] = vip->u.v.data;
228 }
229 #ifdef DIFFVARS_DEBUG
230 CONSOLE_DEBUG("Identified %ld obs vers",diffvars->nobs);
231 CONSOLE_DEBUG("There were %ld rels",prob->nr);
232 #endif
233
234 slv_set_diffvars(sys,(void *)diffvars);
235
236 return 0;
237 }
238
239
240 int system_diffvars_debug(slv_system_t sys,FILE *fp){
241 int i, j;
242 char *varname;
243 const SolverDiffVarCollection *diffvars;
244 SolverDiffVarSequence seq;
245 diffvars = sys->diffvars;
246 if(diffvars==NULL){
247 fprintf(fp,"NO DIFFVARS (NULL)");
248 return 0;
249 }
250 for(i=0; i<diffvars->nseqs;++i){
251 seq = diffvars->seqs[i];
252 fprintf(fp,"%d: ",i);
253 for(j=0; j<seq.n; ++j){
254 if(j)fprintf(fp," <-- ");
255 fprintf(fp,"%d: (%p)",var_sindex(seq.vars[j]),seq.vars[j]);
256 varname = var_make_name(sys,seq.vars[j]);
257 fprintf(fp,"'%s'",varname);
258 ASC_FREE(varname);
259 }
260 fprintf(fp,"\n");
261 }
262 return 0;
263 }
264

john.pye@anu.edu.au
ViewVC Help
Powered by ViewVC 1.1.22