/[ascend]/branches/ksenija2/solvers/ida/idaboundary.c
ViewVC logotype

Contents of /branches/ksenija2/solvers/ida/idaboundary.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2881 - (show annotations) (download) (as text)
Sun Mar 29 05:31:03 2015 UTC (7 years, 10 months ago) by jpye
File MIME type: text/x-csrc
File size: 13813 byte(s)
added lots of debugging output to try to figure out issue with test/test integrator_idaevent.solardynamics.


1 #include "ida.h"
2 #include "idalinear.h"
3 #include "idaanalyse.h"
4 #include "idatypes.h"
5 #include "idaprec.h"
6 #include "idacalc.h"
7 #include "idaio.h"
8 #include "idaboundary.h"
9 #include <stdio.h>
10
11 #include <ascend/general/platform.h>
12 #include <ascend/general/list.h>
13 #include <ascend/general/panic.h>
14
15 #include <ascend/solver/solver.h>
16
17 #include <ascend/system/slv_client.h>
18 #include <ascend/system/cond_config.h>
19 #include <ascend/system/discrete.h>
20 #include <ascend/system/slv_types.h>
21 #include <ascend/system/slv_common.h>
22 #include <ascend/system/logrel.h>
23 #include <ascend/system/rel.h>
24 #include <ascend/system/discrete.h>
25 #include <ascend/system/var.h>
26 #include <ascend/system/block.h>
27 #include <ascend/system/bndman.h>
28
29 #define IDA_BND_DEBUG
30
31 /**
32 * Check to see if any of the system discrete variables have changed.
33 *
34 * @return 1 if any of the values have changed (does not necessarily mean
35 * the system needs to be reconfigured)
36 */
37 int some_dis_vars_changed(slv_system_t sys) {
38 struct dis_discrete **dvlist, *cur_dis;
39 int numDVs, i;
40 #ifdef IDA_BND_DEBUG
41 char *dis_name;
42 #endif
43
44 dvlist = slv_get_solvers_dvar_list(sys);
45 numDVs = slv_get_num_solvers_dvars(sys);
46
47 for (i = 0; i < numDVs; i++) {
48 cur_dis = dvlist[i];
49 if((dis_kind(cur_dis) == e_dis_boolean_t)
50 && (dis_inwhen(cur_dis) || dis_inevent(cur_dis))
51 ){
52 if(dis_value(cur_dis) != dis_previous_value(cur_dis)){
53 #ifdef IDA_BND_DEBUG
54 dis_name = dis_make_name(sys, cur_dis);
55 CONSOLE_DEBUG("Boundary %s (i=%d) has changed (current=%d, prev=%d)", dis_name,
56 i, dis_value(cur_dis), dis_previous_value(cur_dis));
57 ASC_FREE(dis_name);
58 #endif
59 return 1;
60 }
61 }
62 }
63 CONSOLE_DEBUG("No boundary vars have changed");
64 return 0;
65
66 }
67
68 static void ida_write_values(IntegratorSystem *integ) {
69 struct var_variable **vl;
70 struct dis_discrete **dvl;
71 int32 c;
72 vl = slv_get_solvers_var_list(integ->system);
73 dvl = slv_get_solvers_dvar_list(integ->system);
74 for (c = 0; vl[c] != NULL; c++)
75 CONSOLE_DEBUG("Value of %s is %f", var_make_name(integ->system,vl[c]),var_value(vl[c]));
76 for (c = 0; dvl[c] != NULL; c++)
77 CONSOLE_DEBUG("Value of %s is %d", dis_make_name(integ->system,dvl[c]),dis_value(dvl[c]));
78 }
79
80 int ida_setup_lrslv(IntegratorSystem *integ, int qrslv_ind, int lrslv_ind) {
81 CONSOLE_DEBUG("Running logical solver...");
82 ida_log_solve(integ, lrslv_ind);
83 if(some_dis_vars_changed(integ->system)){
84 CONSOLE_DEBUG("Some discrete vars changed; reanalysing");
85 return ida_bnd_reanalyse_cont(integ);
86 }
87 return 0;
88 }
89
90 int ida_bnd_reanalyse(IntegratorSystem *integ){
91
92 if (integ->y_id != NULL) {
93 ASC_FREE(integ->y_id);
94 integ->y_id = NULL;
95 }
96
97 if (integ->obs_id != NULL){
98 ASC_FREE(integ->obs_id);
99 integ->obs_id = NULL;
100 }
101 if (integ->y != NULL) {
102 ASC_FREE(integ->y);
103 integ->y = NULL;
104 }
105 if (integ->ydot != NULL) {
106 ASC_FREE(integ->ydot);
107 integ->ydot = NULL;
108 }
109 if (integ->obs != NULL) {
110 ASC_FREE(integ->obs);
111 integ->obs = NULL;
112 }
113
114 integ->n_y = 0;
115
116
117 return reanalyze_solver_lists(integ->system);
118 }
119
120 int ida_bnd_reanalyse_cont(IntegratorSystem *integ){
121 CONSOLE_DEBUG("Clearing memory from previous calls...");
122 if (integ->y_id != NULL) {
123 ASC_FREE(integ->y_id);
124 integ->y_id = NULL;
125 }
126
127 if (integ->obs_id != NULL){
128 ASC_FREE(integ->obs_id);
129 integ->obs_id = NULL;
130 }
131 if (integ->y != NULL) {
132 ASC_FREE(integ->y);
133 integ->y = NULL;
134 }
135 if (integ->ydot != NULL) {
136 ASC_FREE(integ->ydot);
137 integ->ydot = NULL;
138 }
139 if (integ->obs != NULL) {
140 ASC_FREE(integ->obs);
141 integ->obs = NULL;
142 }
143
144 integ->n_y = 0;
145
146 return integrator_ida_analyse(integ);
147 }
148
149 int ida_bnd_update_relist(IntegratorSystem *integ){
150 IntegratorIdaData *enginedata;
151 struct rel_relation **rels;
152 int i,j,n_solverrels,n_active_rels;
153
154 enginedata = integrator_ida_enginedata(integ);
155
156 n_solverrels = slv_get_num_solvers_rels(integ->system);
157 n_active_rels = slv_count_solvers_rels(integ->system, &integrator_ida_rel);
158 rels = slv_get_solvers_rel_list(integ->system);
159
160 if(enginedata->rellist != NULL){
161 ASC_FREE(enginedata->rellist);
162 enginedata->rellist = NULL;
163 enginedata->rellist = ASC_NEW_ARRAY(struct rel_relation *, n_active_rels);
164 }
165
166 j = 0;
167 for (i = 0; i < n_solverrels; ++i) {
168 if (rel_apply_filter(rels[i], &integrator_ida_rel)) {
169 #ifdef IDA_BND_DEBUG
170 char *relname;
171 relname = rel_make_name(integ->system, rels[i]);
172 CONSOLE_DEBUG("rel '%s': 0x%x", relname, rel_flags(rels[i]));
173 ASC_FREE(relname);
174 #endif
175 enginedata->rellist[j++] = rels[i];
176 }
177 }
178 asc_assert(j == n_active_rels);
179 enginedata->nrels = n_active_rels;
180
181 if (enginedata->nrels != integ->n_y) {
182 ERROR_REPORTER_HERE(ASC_USER_ERROR
183 ,"Integration problem is not square (%d active rels, %d vars)"
184 ,n_active_rels, integ->n_y
185 );
186 return 1; /* failure */
187 }
188
189 return 0;
190 }
191
192 N_Vector ida_bnd_new_zero_NV(long int vec_length){
193 int i;
194 N_Vector nv = N_VNew_Serial(vec_length);
195 for(i= 0; i< vec_length; i++) {
196 NV_Ith_S(nv,i) = 0.0;
197 }
198
199 return nv;
200 }
201
202
203 void ida_bnd_update_IC(IntegratorSystem *integ, realtype t0, N_Vector y0, N_Vector yp0) {
204 /* First destroy since n_y may have changed */
205 N_VDestroy_Serial(y0);
206 N_VDestroy_Serial(yp0);
207 /* retrieve new initial values from the system */
208 t0 = integrator_get_t(integ);
209 y0 = ida_bnd_new_zero_NV(integ->n_y);
210 integrator_get_y(integ, NV_DATA_S(y0));
211
212 yp0 = ida_bnd_new_zero_NV(integ->n_y);
213 integrator_get_ydot(integ, NV_DATA_S(yp0));
214
215 #ifdef IDA_BND_DEBUG
216 CONSOLE_DEBUG("BEFORE IC SOLVING:");
217 CONSOLE_DEBUG("TIME: %f", t0);
218 CONSOLE_DEBUG("Y");
219 N_VPrint_Serial(y0);
220 CONSOLE_DEBUG("Yp");
221 N_VPrint_Serial(yp0);
222 CONSOLE_DEBUG("Press any to continue...");
223 getchar();
224 #endif
225
226 }
227
228 int ida_log_solve(IntegratorSystem *integ, int lrslv_ind) {
229 slv_parameters_t parameters;
230 int num_params;
231 char *pname;
232 slv_status_t status;
233 int i;
234 if (slv_select_solver(integ->system, lrslv_ind) == -1) {
235 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Error attempting to load LRSlv");
236 }
237 #ifdef IDA_BND_DEBUG
238 CONSOLE_DEBUG("Solver selected is '%s'"
239 ,slv_solver_name(slv_get_selected_solver(integ->system))
240 );
241 #endif
242
243 /* Flip LRSlv into ida mode */
244 slv_get_parameters(integ->system, &parameters);
245 num_params = parameters.num_parms;
246 for (i = 0; i<num_params; i++) {
247 pname = parameters.parms[i].name;
248 if(strcmp(pname, "withida") == 0) {
249 parameters.parms[i].info.b.value = 1;
250 }
251 }
252 /* solve the logical relations in the model, if possible */
253 slv_presolve(integ->system);
254 slv_solve(integ->system);
255
256 /* Check for convergence */
257 slv_get_status(integ->system, &status);
258 if (!status.converged) {
259 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Non-convergence in logical solver.");
260 return 0;
261 #ifdef IDA_BND_DEBUG
262 }else{
263 CONSOLE_DEBUG("Converged");
264 #endif
265 }
266 return 1;
267 }
268
269 /*
270 * Uses LRSlv to check if any of the logical conditions in the model need
271 * updating after a boundary crossing. If so, the model is reconfigured and
272 * new initial conditions are determined at the point of the crossing.
273 *
274 * @return 1 if the crossing causes a change in the system
275 */
276
277 int ida_cross_boundary(IntegratorSystem *integ, int *rootsfound,
278 int *bnd_cond_states, int qrslv_ind, int lrslv_ind) {
279
280 IntegratorIdaData *enginedata;
281 slv_status_t status;
282
283 struct bnd_boundary *bnd = NULL;
284 int i, num_bnds, num_dvars;
285
286 struct dis_discrete **dvl;
287 int32 c, *prevals = NULL;
288
289 double *bnd_prev_eval;
290
291 CONSOLE_DEBUG("Crossing boundary...");
292
293 integrator_output_write(integ);
294 integrator_output_write_obs(integ);
295 integrator_output_write_event(integ);
296 /*^^^ FIXME should we write the above event? It may not have crossed in correct direction...? */
297
298 /* Flag the crossed boundary and update bnd_cond_states */
299 enginedata = integ->enginedata;
300 num_bnds = enginedata->nbnds;
301 bnd_prev_eval = ASC_NEW_ARRAY(double,num_bnds);
302 for(i = 0; i < num_bnds; i++) {
303 /* get the value of the boundary condition before crossing */
304 bnd_prev_eval[i] = bndman_real_eval(enginedata->bndlist[i]);
305 if(rootsfound[i]){
306 /* reminder: 'rootsfound[i] == 1 for UP or -1 for DOWN. */
307 /* this boundary was one of the boundaries triggered */
308 bnd = enginedata->bndlist[i];
309 bnd_set_ida_crossed(bnd, 1);
310 if(bnd_ida_first_cross(bnd) /* not crossed before */
311 || !((rootsfound[i] == 1 && bnd_ida_incr(bnd)) /* not crossing upwards twice in a row */
312 || (rootsfound[i] == -1 && !bnd_ida_incr(bnd))) /* not crossing downwards twice in a row */
313 ){
314 if(bnd_cond_states[i] == 0){
315 bnd_set_ida_value(bnd, 1);
316 bnd_cond_states[i] = 1;
317 }else{
318 bnd_set_ida_value(bnd, 0);
319 bnd_cond_states[i] = 0;
320 }
321 #ifdef IDA_BND_DEBUG
322 char *n = bnd_make_name(integ->system,bnd);
323 CONSOLE_DEBUG("Set boundary '%s' to %d",n,bnd_ida_value(bnd)?1:0);
324 ASC_FREE(n);
325 #endif
326 }else{
327 /* Boundary crossed twice in one direction. Very unlikey! */
328
329 /* The aim of the following two lines is to set the value of
330 the boolean variable to such a value as if the boundary
331 was crossed in the opposite direction before */
332 CONSOLE_DEBUG("DOUBLE CROSS");
333 if(!prevals){
334 num_dvars = slv_get_num_solvers_dvars(integ->system);
335 dvl = slv_get_solvers_dvar_list(integ->system);
336 prevals = ASC_NEW_ARRAY(int32,num_dvars);
337 for(c = 0; dvl[c] != NULL; c++)
338 prevals[c] = dis_value(dvl[c]);
339 }
340 bnd_set_ida_value(bnd, !bnd_cond_states[i]);
341 if(!ida_log_solve(integ,lrslv_ind)){
342 CONSOLE_DEBUG("Error in logic solve in double-cross");
343 return -1;
344 }
345 bnd_set_ida_value(bnd,bnd_cond_states[i]);
346 #ifdef IDA_BND_DEBUG
347 char *n = bnd_make_name(integ->system,bnd);
348 CONSOLE_DEBUG("After double-cross, set boundary '%s' to %d",n,bnd_ida_value(bnd)?1:0);
349 ASC_FREE(n);
350 #endif
351 }
352 /* flag this boundary as having previously been crossed */
353 bnd_set_ida_first_cross(bnd,0);
354 /* store this recent crossing-direction for future reference */
355 bnd_set_ida_incr(bnd,(rootsfound[i] > 0));
356 }
357 }
358 if(!ida_log_solve(integ,lrslv_ind)) return -1;
359
360 /* If there was a double crossing, because of ida_log_solve the previous values
361 of discrete variables may be equal to their current values, which would mean that
362 the discrete vars haven't changed their values, but actually they did. So here we
363 restore the previous values of those variables, which are not connected with the
364 boundary which was double-crossed. */
365 if(prevals){
366 for(c = 0; dvl[c] != NULL; c++)
367 if(!(dis_value(dvl[c]) == prevals[c] && dis_value(dvl[c]) != dis_previous_value(dvl[c])))
368 dis_set_previous_value(dvl[c],prevals[c]);
369 ASC_FREE(prevals);
370 }
371
372 integrator_output_write(integ);
373 integrator_output_write_obs(integ);
374 integrator_output_write_event(integ);
375 if(!some_dis_vars_changed(integ->system)){
376 CONSOLE_DEBUG("Crossed boundary but no effect on system");
377 /* Boundary crossing that has no effect on system */
378 ASC_FREE(bnd_prev_eval);
379
380 /* Reset the boundary flag */
381 for(i = 0; i < num_bnds; i++)
382 bnd_set_ida_crossed(enginedata->bndlist[i], 0);
383 return 0;
384 }
385 /* update the main system if required */
386 int events_triggered = 1;
387 if(slv_get_num_solvers_events(integ->system)!=0) {
388 while(some_dis_vars_changed(integ->system) && events_triggered) {
389 if(ida_bnd_reanalyse(integ)) {
390 /* select QRSlv solver, and solve the system */
391 if(slv_select_solver(integ->system, qrslv_ind) == -1) {
392 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Error attempting to load QRSlv");
393 }
394 #ifdef IDA_BND_DEBUG
395 {
396 struct var_variable **vl = slv_get_solvers_var_list(integ->system);
397 int i, n=slv_get_num_solvers_vars(integ->system), first=1;
398 CONSOLE_DEBUG("In boundary problem: variables (active, incident, free):");
399 for(i=0;i<n;++i){
400 if(var_incident(vl[i])&&!var_fixed(vl[i])&&var_active(vl[i])){
401 char *name = var_make_name(integ->system,vl[i]);
402 fprintf(stderr,"%s%s",(first?"\t":", "),name);
403 first=0; ASC_FREE(name);
404 }
405 }
406 fprintf(stderr,"\n");
407 }
408 {
409 struct rel_relation **rl = slv_get_solvers_rel_list(integ->system);
410 int i, n=slv_get_num_solvers_rels(integ->system), first=1;
411 CONSOLE_DEBUG("...relations (equality, included, active):");
412 for(i=0;i<n;++i){
413 if(rel_equality(rl[i])&&rel_active(rl[i])&&rel_included(rl[i])){
414 char *name = rel_make_name(integ->system,rl[i]);
415 fprintf(stderr,"%s%s",(first?"\t":", "),name);
416 first=0; ASC_FREE(name);
417 }
418 }
419 fprintf(stderr,"\n");
420 }
421 #endif
422 slv_presolve(integ->system);
423 slv_solve(integ->system);
424
425 slv_get_status(integ->system, &status);
426 if (!status.converged) {
427 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Non-convergence in "
428 "non-linear solver at boundary");
429 #ifdef IDA_BND_DEBUG
430 }else{
431 CONSOLE_DEBUG("Converged");
432 #endif
433 }
434
435 for (i = 0; i < num_bnds; i++) {
436 if (!rootsfound[i]) {
437 if ((bndman_real_eval(enginedata->bndlist[i]) > 0 && bnd_prev_eval[i] < 0) ||
438 (bndman_real_eval(enginedata->bndlist[i]) < 0 && bnd_prev_eval[i] > 0)) {
439 bnd_cond_states[i] = !bnd_cond_states[i];
440 if (bnd_prev_eval[i] > 0) {
441 bnd_set_ida_incr(enginedata->bndlist[i],0);
442 rootsfound[i] = -1;
443 }
444 else {
445 bnd_set_ida_incr(enginedata->bndlist[i],1);
446 rootsfound[i] = 1;
447 }
448 bnd_set_ida_value(enginedata->bndlist[i],bnd_cond_states[i]);
449 }
450 }
451 }
452 if (!ida_log_solve(integ,lrslv_ind)) return -1;
453
454 }else events_triggered = 0;
455 if (ida_bnd_reanalyse_cont(integ)) return 2;
456 if (events_triggered) {
457 integrator_output_write(integ);
458 integrator_output_write_obs(integ);
459 integrator_output_write_event(integ);
460 }
461 }
462 }
463
464 integrator_output_write(integ);
465 integrator_output_write_obs(integ);
466 integrator_output_write_event(integ);
467
468 ASC_FREE(bnd_prev_eval);
469
470 /* Reset the boundary flag */
471 for (i = 0; i < num_bnds; i++)
472 bnd_set_ida_crossed(enginedata->bndlist[i], 0);
473 return 1;
474 }
475

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