/[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 1241 - (show annotations) (download) (as text)
Fri Jan 26 11:59:45 2007 UTC (15 years, 5 months ago) by johnpye
File MIME type: text/x-csrc
File size: 9293 byte(s)
Fixed a bug with idaanalyse, hunting another one.
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
28 #include <compiler/instance_io.h>
29
30 #include "analyse_impl.h"
31
32 /*------------------------------------------------------------------------------
33 DERIVATIVES & DIFFERENTIAL VARIABLES
34 */
35
36 #ifdef ASC_IDA_NEW_ANALYSE
37
38 SolverDiffVarCollection *analyse_get_diffvars(slv_system_t sys){
39 return (SolverDiffVarCollection *)slv_get_diffvars(sys);
40 }
41
42 static
43 int CmpDiffVars(const struct solver_ipdata *a, const struct solver_ipdata *b){
44 if(a->u.v.odeid < b->u.v.odeid)return -1;
45 if(a->u.v.odeid > b->u.v.odeid)return 1;
46 /* now, the ode_id is equal: next level of sorting... */
47 if(a->u.v.deriv < b->u.v.deriv)return -1;
48 if(a->u.v.deriv > b->u.v.deriv)return 1;
49 return 0;
50 }
51
52 /**
53 This function steals a lot of what was in integrator.c, but we try to make
54 it more general and capable of extracting info about high-order derivative
55 chains. One eye open to the possibility of doing Pantelides (or other?)
56 style index reduction.
57
58 Also, note that we try to only use data from the struct problem_t, NOT
59 to go back to the instance hierarchy. The idea is that we have already got
60 what we need in classify_instance, etc.
61
62 @return 0 on success
63 */
64 int analyse_generate_diffvars(slv_system_t sys, struct problem_t *prob){
65 SolverDiffVarCollection *diffvars = NULL;
66 struct solver_ipdata *vip, *vipnext;
67 SolverDiffVarSequence *seq;
68 struct gl_list_t *seqs;
69 long i, seqstart, nalg, ndiff;
70 short j;
71 char cont/*, *varname*/;
72 short maxorder = 0;
73
74 asc_assert(prob);
75
76 if(gl_length(prob->diffvars)==0){
77 CONSOLE_DEBUG("No differential variables were seen. Skipping generation of diffvars struct.");
78 return 0;
79 }
80
81 seqs = gl_create(prob->nr);
82
83 /* add the list of algebraic variables to the structure too */
84 asc_assert(prob->algebvars);
85 nalg = gl_length(prob->algebvars);
86 for(i=0; i<nalg; ++i){
87 vip = (struct solver_ipdata *)gl_fetch(prob->algebvars,i+1);
88 seq = ASC_NEW(SolverDiffVarSequence);
89 seq->n = 1;
90 seq->ode_id = 0;
91 seq->vars = ASC_NEW_ARRAY(struct var_variable*,1);
92 seq->vars[0] = vip->u.v.data;
93 gl_append_ptr(seqs,(void *)seq);
94 }
95 CONSOLE_DEBUG("Added %ld algebraic vars to chains",nalg);
96
97 CONSOLE_DEBUG("Sorting %ld differential & derivative vars...",gl_length(prob->diffvars));
98
99 /* first sort the list of diffvars */
100 gl_sort(prob->diffvars, (CmpFunc)CmpDiffVars);
101
102 /* scan the list for derivs that are missing vars */
103 i = 1; cont = TRUE; seqstart =1; ndiff = 0;
104 vip = (struct solver_ipdata *)gl_fetch(prob->diffvars,i);
105 if(vip->u.v.deriv > 1){
106 ERROR_REPORTER_START_NOLINE(ASC_USER_ERROR);
107 FPRINTF(ASCERR,"Missing ode_type %d for ode_id %d (check var '",vip->u.v.deriv+1,vip->u.v.odeid);
108 WriteInstanceName(ASCWAR,vip->i,prob->root);
109 FPRINTF(ASCERR,"')");
110 error_reporter_end_flush();
111 return 2;
112 }
113 while(cont){
114 /* CONSOLE_DEBUG("Working, seqstart=%ld",seqstart); */
115 if(i >= gl_length(prob->diffvars))cont = FALSE;
116 else vipnext = (struct solver_ipdata *)gl_fetch(prob->diffvars,i+1);
117
118 if(cont && vipnext->u.v.odeid == vip->u.v.odeid){
119 /* same sequence, check that it's the next derivative */
120 asc_assert(vipnext->u.v.odeid == vip->u.v.odeid);
121 asc_assert(vip->u.v.deriv >= vip->u.v.deriv);
122
123 if(vipnext->u.v.deriv == vip->u.v.deriv){
124 ERROR_REPORTER_START_NOLINE(ASC_USER_ERROR);
125 FPRINTF(ASCERR,"Repeated ode_type %d for ode_id %d (var '",vip->u.v.deriv,vip->u.v.odeid);
126 WriteInstanceName(ASCWAR,vip->i,prob->root);
127 FPRINTF(ASCERR,"')");
128 error_reporter_end_flush();
129 return 1;
130 }else if(vipnext->u.v.deriv > vip->u.v.deriv + 1){
131 ERROR_REPORTER_START_NOLINE(ASC_USER_ERROR);
132 FPRINTF(ASCERR,"Missing ode_type %d for ode_id %d (check var '",vip->u.v.deriv+1,vip->u.v.odeid);
133 WriteInstanceName(ASCWAR,vip->i,prob->root);
134 FPRINTF(ASCERR,"')");
135 error_reporter_end_flush();
136 return 2;
137 }
138
139 /* allow this var into the sequence, loop again */
140 vip = vipnext; ++i; continue;
141 }
142
143 /* close the current sequence */
144 if(i == seqstart){
145 ERROR_REPORTER_START_NOLINE(ASC_USER_ERROR);
146 FPRINTF(ASCERR,"Variable '");
147 WriteInstanceName(ASCWAR,vip->i,prob->root);
148 FPRINTF(ASCERR,"' declared differential without derivative being identified (ode_id=%d)",vip->u.v.odeid);
149 error_reporter_end_flush();
150 return 3;
151 }
152 seq = ASC_NEW(SolverDiffVarSequence);
153 seq->n = i - seqstart + 1;
154 seq->ode_id = vip->u.v.odeid;
155 /* CONSOLE_DEBUG("Saving sequence ode_id = %ld, n = %d",seq->ode_id,seq->n); */
156 seq->vars = ASC_NEW_ARRAY(struct var_variable*,seq->n);
157 for(j=0;j<seq->n;++j){
158 vip = (struct solver_ipdata *)gl_fetch(prob->diffvars,seqstart+j);
159 seq->vars[j]=vip->u.v.data;
160 /* varname = var_make_name(sys,seq->vars[j]);
161 CONSOLE_DEBUG("seq, j=%d: '%s'",j,varname);
162 ASC_FREE(varname); */
163 }
164 gl_append_ptr(seqs,(void *)seq);
165 ndiff++;
166 /* CONSOLE_DEBUG("Completed seq %ld",gl_length(seqs)); */
167
168 if(vip->u.v.deriv > maxorder){
169 maxorder = vip->u.v.deriv;
170 /* CONSOLE_DEBUG("maxorder increased to %d",maxorder); */
171 }
172
173 ++i; seqstart=i; vip=vipnext;
174 if(cont && vip->u.v.deriv != 1){
175 ERROR_REPORTER_START_NOLINE(ASC_USER_ERROR);
176 FPRINTF(ASCERR,"Missing undifferentiated form of var '");
177 WriteInstanceName(ASCWAR,vip->i,prob->root);
178 FPRINTF(ASCERR,"' (ode_id = %d, ode_type = %d)",vip->u.v.odeid,vip->u.v.deriv);
179 error_reporter_end_flush();
180 return 4;
181 }
182 continue;
183 }
184
185 CONSOLE_DEBUG("Identified %ld derivative chains, maximum length %d...",gl_length(seqs),maxorder);
186
187 diffvars = ASC_NEW(SolverDiffVarCollection);
188 diffvars->maxorder = maxorder;
189 diffvars->nseqs = gl_length(seqs);
190 diffvars->seqs = ASC_NEW_ARRAY(SolverDiffVarSequence,diffvars->nseqs);
191 for(i=0;i<diffvars->nseqs;++i){
192 diffvars->seqs[i] = *(SolverDiffVarSequence *)gl_fetch(seqs,i+1);
193 }
194 gl_destroy(seqs);
195 diffvars->nalg = nalg;
196 diffvars->ndiff = ndiff;
197
198 diffvars->nindep = gl_length(prob->indepvars);
199 diffvars->indep = ASC_NEW_ARRAY(struct var_variable *,diffvars->nindep);
200 for(i=0;i<diffvars->nindep;++i){
201 vip = (struct solver_ipdata *)gl_fetch(prob->indepvars,i+1);
202 diffvars->indep[i] = vip->u.v.data;
203 }
204 CONSOLE_DEBUG("Identified %ld indep vars",diffvars->nindep);
205
206 /* I know, I know, this might not be the best place for this... */
207 diffvars->nobs = gl_length(prob->obsvars);
208 diffvars->obs = ASC_NEW_ARRAY(struct var_variable *,diffvars->nobs);
209 for(i=0;i<diffvars->nobs;++i){
210 vip = (struct solver_ipdata *)gl_fetch(prob->obsvars,i+1);
211 diffvars->obs[i] = vip->u.v.data;
212 }
213 CONSOLE_DEBUG("Identified %ld obs vers",diffvars->nobs);
214
215 CONSOLE_DEBUG("There were %ld rels",prob->nr);
216
217 slv_set_diffvars(sys,(void *)diffvars);
218
219 return 0;
220 }
221
222
223 int analyse_diffvars_debug(slv_system_t sys,FILE *fp){
224 int i, j;
225 char *varname;
226 const SolverDiffVarCollection *diffvars;
227 SolverDiffVarSequence seq;
228 diffvars = analyse_get_diffvars(sys);
229
230 for(i=0; i<diffvars->nseqs;++i){
231 seq = diffvars->seqs[i];
232 fprintf(fp,"%d: ",i);
233 for(j=0; j<seq.n; ++j){
234 if(j)fprintf(fp," <-- ");
235 fprintf(fp,"%d: (%p)",var_sindex(seq.vars[j]),seq.vars[j]);
236 varname = var_make_name(sys,seq.vars[j]);
237 fprintf(fp,"'%s'",varname);
238 ASC_FREE(varname);
239 }
240 fprintf(fp,"\n");
241 }
242 return 0;
243 }
244
245 typedef int CmpFn(const void *,const void *);
246 /**
247 Comparison function for use by analyse_diffvars_sort.
248 All we're trying to do is to reproduce the ordering in var_sindex, so
249 no need to actually look at the flags.
250 */
251 static int analyse_diffvars_cmp(const SolverDiffVarSequence *a, const SolverDiffVarSequence *b){
252 const struct var_variable *va;
253 const struct var_variable *vb;
254 int ia;
255 int ib;
256 #if 0
257 int aok;
258
259 var_filter_t vf = {VAR_FIXED|VAR_INCIDENT|VAR_ACTIVE|VAR_SVAR|VAR_DERIV
260 ,0 |VAR_INCIDENT|VAR_ACTIVE|VAR_SVAR|0 };
261 #endif
262
263 asc_assert(a->n >= 1);
264 asc_assert(b->n >= 1);
265
266 va = a->vars[0];
267 vb = b->vars[0];
268
269 #if 0
270 aok = var_apply_filter(va,&vf);
271 if(var_apply_filter(vb,&vf)){
272 if(!aok){
273 /* 'b' better */
274 return 1;
275 }
276 }else if(aok){
277 /* 'a' better */
278 return -1;
279 }
280 #endif
281
282 ia = var_sindex(va);
283 ib = var_sindex(vb);
284
285 if(ia<ib)return -1;
286 if(ia==ib)return 0;
287
288 return 1;
289 }
290
291 int analyse_diffvars_sort(slv_system_t sys){
292 int i, j;
293 char *varname;
294 const SolverDiffVarCollection *diffvars;
295 SolverDiffVarSequence seq;
296 diffvars = analyse_get_diffvars(sys);
297
298 qsort(diffvars->seqs,diffvars->nseqs,sizeof(SolverDiffVarSequence),(CmpFn*)(&analyse_diffvars_cmp));
299 }
300
301
302 #endif

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