/[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 1249 - (show annotations) (download) (as text)
Sat Jan 27 04:14:58 2007 UTC (15 years, 5 months ago) by johnpye
File MIME type: text/x-csrc
File size: 8086 byte(s)
Replacement integrator_ida_check_diffindex function. Fixed some errors in the old 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 #include "system_impl.h"
32
33 /*------------------------------------------------------------------------------
34 DERIVATIVES & DIFFERENTIAL VARIABLES
35 */
36
37 SolverDiffVarCollection *system_get_diffvars(slv_system_t sys){
38 return sys->diffvars;
39 }
40
41 static
42 int CmpDiffVars(const struct solver_ipdata *a, const struct solver_ipdata *b){
43 if(a->u.v.odeid < b->u.v.odeid)return -1;
44 if(a->u.v.odeid > b->u.v.odeid)return 1;
45 /* now, the ode_id is equal: next level of sorting... */
46 if(a->u.v.deriv < b->u.v.deriv)return -1;
47 if(a->u.v.deriv > b->u.v.deriv)return 1;
48 return 0;
49 }
50
51 /**
52 This function steals a lot of what was in integrator.c, but we try to make
53 it more general and capable of extracting info about high-order derivative
54 chains. One eye open to the possibility of doing Pantelides (or other?)
55 style index reduction.
56
57 Also, note that we try to only use data from the struct problem_t, NOT
58 to go back to the instance hierarchy. The idea is that we have already got
59 what we need in classify_instance, etc.
60
61 @return 0 on success
62 */
63 int system_generate_diffvars(slv_system_t sys, struct problem_t *prob){
64 SolverDiffVarCollection *diffvars = NULL;
65 struct solver_ipdata *vip, *vipnext;
66 SolverDiffVarSequence *seq;
67 struct gl_list_t *seqs;
68 long i, seqstart, nalg, ndiff;
69 short j;
70 char cont/*, *varname*/;
71 short maxorder = 0;
72
73 asc_assert(prob);
74
75 if(gl_length(prob->diffvars)==0){
76 CONSOLE_DEBUG("No differential variables were seen. Skipping generation of diffvars struct.");
77 return 0;
78 }
79
80 seqs = gl_create(prob->nr);
81
82 /* add the list of algebraic variables to the structure too */
83 asc_assert(prob->algebvars);
84 nalg = gl_length(prob->algebvars);
85 for(i=0; i<nalg; ++i){
86 vip = (struct solver_ipdata *)gl_fetch(prob->algebvars,i+1);
87 seq = ASC_NEW(SolverDiffVarSequence);
88 seq->n = 1;
89 seq->ode_id = 0;
90 seq->vars = ASC_NEW_ARRAY(struct var_variable*,1);
91 seq->vars[0] = vip->u.v.data;
92 gl_append_ptr(seqs,(void *)seq);
93 }
94 CONSOLE_DEBUG("Added %ld algebraic vars to chains",nalg);
95
96 CONSOLE_DEBUG("Sorting %ld differential & derivative vars...",gl_length(prob->diffvars));
97
98 /* first sort the list of diffvars */
99 gl_sort(prob->diffvars, (CmpFunc)CmpDiffVars);
100
101 /* scan the list for derivs that are missing vars */
102 i = 1; cont = TRUE; seqstart =1; ndiff = 0;
103 vip = (struct solver_ipdata *)gl_fetch(prob->diffvars,i);
104 if(vip->u.v.deriv > 1){
105 ERROR_REPORTER_START_NOLINE(ASC_USER_ERROR);
106 FPRINTF(ASCERR,"Missing ode_type %d for ode_id %d (check var '",vip->u.v.deriv+1,vip->u.v.odeid);
107 WriteInstanceName(ASCWAR,vip->i,prob->root);
108 FPRINTF(ASCERR,"')");
109 error_reporter_end_flush();
110 return 2;
111 }
112 while(cont){
113 /* CONSOLE_DEBUG("Working, seqstart=%ld",seqstart); */
114 if(i >= gl_length(prob->diffvars))cont = FALSE;
115 else vipnext = (struct solver_ipdata *)gl_fetch(prob->diffvars,i+1);
116
117 if(cont && vipnext->u.v.odeid == vip->u.v.odeid){
118 /* same sequence, check that it's the next derivative */
119 asc_assert(vipnext->u.v.odeid == vip->u.v.odeid);
120 asc_assert(vip->u.v.deriv >= vip->u.v.deriv);
121
122 if(vipnext->u.v.deriv == vip->u.v.deriv){
123 ERROR_REPORTER_START_NOLINE(ASC_USER_ERROR);
124 FPRINTF(ASCERR,"Repeated ode_type %d for ode_id %d (var '",vip->u.v.deriv,vip->u.v.odeid);
125 WriteInstanceName(ASCWAR,vip->i,prob->root);
126 FPRINTF(ASCERR,"')");
127 error_reporter_end_flush();
128 return 1;
129 }else if(vipnext->u.v.deriv > vip->u.v.deriv + 1){
130 ERROR_REPORTER_START_NOLINE(ASC_USER_ERROR);
131 FPRINTF(ASCERR,"Missing ode_type %d for ode_id %d (check var '",vip->u.v.deriv+1,vip->u.v.odeid);
132 WriteInstanceName(ASCWAR,vip->i,prob->root);
133 FPRINTF(ASCERR,"')");
134 error_reporter_end_flush();
135 return 2;
136 }
137
138 /* allow this var into the sequence, loop again */
139 vip = vipnext; ++i; continue;
140 }
141
142 /* close the current sequence */
143 if(i == seqstart){
144 ERROR_REPORTER_START_NOLINE(ASC_USER_ERROR);
145 FPRINTF(ASCERR,"Variable '");
146 WriteInstanceName(ASCWAR,vip->i,prob->root);
147 FPRINTF(ASCERR,"' declared differential without derivative being identified (ode_id=%d)",vip->u.v.odeid);
148 error_reporter_end_flush();
149 return 3;
150 }
151 seq = ASC_NEW(SolverDiffVarSequence);
152 seq->n = i - seqstart + 1;
153 seq->ode_id = vip->u.v.odeid;
154 /* CONSOLE_DEBUG("Saving sequence ode_id = %ld, n = %d",seq->ode_id,seq->n); */
155 seq->vars = ASC_NEW_ARRAY(struct var_variable*,seq->n);
156 for(j=0;j<seq->n;++j){
157 vip = (struct solver_ipdata *)gl_fetch(prob->diffvars,seqstart+j);
158 seq->vars[j]=vip->u.v.data;
159 /* varname = var_make_name(sys,seq->vars[j]);
160 CONSOLE_DEBUG("seq, j=%d: '%s'",j,varname);
161 ASC_FREE(varname); */
162 }
163 gl_append_ptr(seqs,(void *)seq);
164 ndiff++;
165 /* CONSOLE_DEBUG("Completed seq %ld",gl_length(seqs)); */
166
167 if(vip->u.v.deriv > maxorder){
168 maxorder = vip->u.v.deriv;
169 /* CONSOLE_DEBUG("maxorder increased to %d",maxorder); */
170 }
171
172 ++i; seqstart=i; vip=vipnext;
173 if(cont && vip->u.v.deriv != 1){
174 ERROR_REPORTER_START_NOLINE(ASC_USER_ERROR);
175 FPRINTF(ASCERR,"Missing undifferentiated form of var '");
176 WriteInstanceName(ASCWAR,vip->i,prob->root);
177 FPRINTF(ASCERR,"' (ode_id = %d, ode_type = %d)",vip->u.v.odeid,vip->u.v.deriv);
178 error_reporter_end_flush();
179 return 4;
180 }
181 continue;
182 }
183
184 CONSOLE_DEBUG("Identified %ld derivative chains, maximum length %d...",gl_length(seqs),maxorder);
185
186 diffvars = ASC_NEW(SolverDiffVarCollection);
187 diffvars->maxorder = maxorder;
188 diffvars->nseqs = gl_length(seqs);
189 diffvars->seqs = ASC_NEW_ARRAY(SolverDiffVarSequence,diffvars->nseqs);
190 for(i=0;i<diffvars->nseqs;++i){
191 diffvars->seqs[i] = *(SolverDiffVarSequence *)gl_fetch(seqs,i+1);
192 }
193 gl_destroy(seqs);
194 diffvars->nalg = nalg;
195 diffvars->ndiff = ndiff;
196
197 diffvars->nindep = gl_length(prob->indepvars);
198 diffvars->indep = ASC_NEW_ARRAY(struct var_variable *,diffvars->nindep);
199 for(i=0;i<diffvars->nindep;++i){
200 vip = (struct solver_ipdata *)gl_fetch(prob->indepvars,i+1);
201 diffvars->indep[i] = vip->u.v.data;
202 }
203 CONSOLE_DEBUG("Identified %ld indep vars",diffvars->nindep);
204
205 /* I know, I know, this might not be the best place for this... */
206 diffvars->nobs = gl_length(prob->obsvars);
207 diffvars->obs = ASC_NEW_ARRAY(struct var_variable *,diffvars->nobs);
208 for(i=0;i<diffvars->nobs;++i){
209 vip = (struct solver_ipdata *)gl_fetch(prob->obsvars,i+1);
210 diffvars->obs[i] = vip->u.v.data;
211 }
212 CONSOLE_DEBUG("Identified %ld obs vers",diffvars->nobs);
213
214 CONSOLE_DEBUG("There were %ld rels",prob->nr);
215
216 slv_set_diffvars(sys,(void *)diffvars);
217
218 return 0;
219 }
220
221
222 int system_diffvars_debug(slv_system_t sys,FILE *fp){
223 int i, j;
224 char *varname;
225 const SolverDiffVarCollection *diffvars;
226 SolverDiffVarSequence seq;
227 diffvars = sys->diffvars;
228 if(diffvars==NULL){
229 fprintf(fp,"NO DIFFVARS (NULL)");
230 return 0;
231 }
232 for(i=0; i<diffvars->nseqs;++i){
233 seq = diffvars->seqs[i];
234 fprintf(fp,"%d: ",i);
235 for(j=0; j<seq.n; ++j){
236 if(j)fprintf(fp," <-- ");
237 fprintf(fp,"%d: (%p)",var_sindex(seq.vars[j]),seq.vars[j]);
238 varname = var_make_name(sys,seq.vars[j]);
239 fprintf(fp,"'%s'",varname);
240 ASC_FREE(varname);
241 }
242 fprintf(fp,"\n");
243 }
244 return 0;
245 }
246

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