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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2882 - (hide annotations) (download) (as text)
Sun Mar 29 12:20:00 2015 UTC (7 years, 10 months ago) by jpye
File MIME type: text/x-csrc
File size: 14074 byte(s)
more debugging info

1 jpye 2559 #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 jpye 2368 #include "idaboundary.h"
9 jpye 2559 #include <stdio.h>
10 jpye 2368
11 jpye 2559 #include <ascend/general/platform.h>
12     #include <ascend/general/list.h>
13     #include <ascend/general/panic.h>
14 jpye 2368
15 jpye 2559 #include <ascend/solver/solver.h>
16 jpye 2368
17 jpye 2559 #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 jpye 2836 #include <ascend/system/discrete.h>
25     #include <ascend/system/var.h>
26     #include <ascend/system/block.h>
27     #include <ascend/system/bndman.h>
28 jpye 2368
29 jpye 2869 #define IDA_BND_DEBUG
30 jpye 2836
31 jpye 2881 /**
32     * Check to see if any of the system discrete variables have changed.
33 jpye 2559 *
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 jpye 2836 int numDVs, i;
40     #ifdef IDA_BND_DEBUG
41 jpye 2559 char *dis_name;
42 jpye 2836 #endif
43 jpye 2559
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 jpye 2881 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 jpye 2559 #ifdef IDA_BND_DEBUG
54 jpye 2881 dis_name = dis_make_name(sys, cur_dis);
55 jpye 2882 CONSOLE_DEBUG("Boolean %s (i=%d) has changed (current=%d, prev=%d)", dis_name,
56 jpye 2881 i, dis_value(cur_dis), dis_previous_value(cur_dis));
57     ASC_FREE(dis_name);
58 jpye 2559 #endif
59 jpye 2836 return 1;
60 jpye 2559 }
61     }
62 jpye 2368 }
63 jpye 2881 CONSOLE_DEBUG("No boundary vars have changed");
64 jpye 2836 return 0;
65 jpye 2368
66 jpye 2559 }
67    
68 jpye 2836 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 jpye 2881 CONSOLE_DEBUG("Running logical solver...");
82 jpye 2836 ida_log_solve(integ, lrslv_ind);
83 jpye 2869 if(some_dis_vars_changed(integ->system)){
84 jpye 2881 CONSOLE_DEBUG("Some discrete vars changed; reanalysing");
85 jpye 2869 return ida_bnd_reanalyse_cont(integ);
86     }
87 jpye 2836 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 jpye 2368 }
96    
97 jpye 2836 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 jpye 2368
114 jpye 2836 integ->n_y = 0;
115 jpye 2559
116 jpye 2836
117     return reanalyze_solver_lists(integ->system);
118 jpye 2368 }
119 jpye 2559
120 jpye 2836 int ida_bnd_reanalyse_cont(IntegratorSystem *integ){
121 jpye 2881 CONSOLE_DEBUG("Clearing memory from previous calls...");
122 jpye 2559 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 jpye 2836 return integrator_ida_analyse(integ);
147 jpye 2559 }
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 jpye 2836 char *relname;
171 jpye 2559 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 jpye 2836 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 jpye 2869 CONSOLE_DEBUG("Solver selected is '%s'"
239     ,slv_solver_name(slv_get_selected_solver(integ->system))
240     );
241 jpye 2836 #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 jpye 2869 #ifdef IDA_BND_DEBUG
262     }else{
263     CONSOLE_DEBUG("Converged");
264     #endif
265 jpye 2836 }
266     return 1;
267     }
268    
269 jpye 2559 /*
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 jpye 2836 int *bnd_cond_states, int qrslv_ind, int lrslv_ind) {
279 jpye 2559
280     IntegratorIdaData *enginedata;
281     slv_status_t status;
282    
283 jpye 2836 struct bnd_boundary *bnd = NULL;
284     int i, num_bnds, num_dvars;
285 jpye 2559
286 jpye 2836 struct dis_discrete **dvl;
287     int32 c, *prevals = NULL;
288    
289     double *bnd_prev_eval;
290    
291 jpye 2881 CONSOLE_DEBUG("Crossing boundary...");
292    
293 jpye 2836 integrator_output_write(integ);
294     integrator_output_write_obs(integ);
295     integrator_output_write_event(integ);
296 jpye 2881 /*^^^ FIXME should we write the above event? It may not have crossed in correct direction...? */
297 jpye 2836
298 jpye 2559 /* Flag the crossed boundary and update bnd_cond_states */
299     enginedata = integ->enginedata;
300     num_bnds = enginedata->nbnds;
301 jpye 2836 bnd_prev_eval = ASC_NEW_ARRAY(double,num_bnds);
302 jpye 2869 for(i = 0; i < num_bnds; i++) {
303 jpye 2881 /* get the value of the boundary condition before crossing */
304 jpye 2836 bnd_prev_eval[i] = bndman_real_eval(enginedata->bndlist[i]);
305 jpye 2881 if(rootsfound[i]){
306     /* reminder: 'rootsfound[i] == 1 for UP or -1 for DOWN. */
307     /* this boundary was one of the boundaries triggered */
308 jpye 2559 bnd = enginedata->bndlist[i];
309 jpye 2882 #ifdef IDA_BND_DEBUG
310     char *name = bnd_make_name(integ->system,enginedata->bndlist[i]);
311     CONSOLE_DEBUG("Boundary '%s': ida_incr=%d, dirn=%d,ida_value=%d,ida_first_cross=%d,bnd_cond_states[i]=%d"
312     ,name,bnd_ida_incr(bnd),rootsfound[i],bnd_ida_value(bnd)
313     ,bnd_ida_first_cross(bnd)!=0,bnd_cond_states[i]
314     );
315     #endif
316 jpye 2559 bnd_set_ida_crossed(bnd, 1);
317 jpye 2881 if(bnd_ida_first_cross(bnd) /* not crossed before */
318     || !((rootsfound[i] == 1 && bnd_ida_incr(bnd)) /* not crossing upwards twice in a row */
319     || (rootsfound[i] == -1 && !bnd_ida_incr(bnd))) /* not crossing downwards twice in a row */
320     ){
321     if(bnd_cond_states[i] == 0){
322 jpye 2836 bnd_set_ida_value(bnd, 1);
323     bnd_cond_states[i] = 1;
324 jpye 2869 }else{
325 jpye 2836 bnd_set_ida_value(bnd, 0);
326     bnd_cond_states[i] = 0;
327     }
328 jpye 2881 #ifdef IDA_BND_DEBUG
329 jpye 2882 CONSOLE_DEBUG("Set boundary '%s' to %d (single cross)",name,bnd_ida_value(bnd)!=0);
330 jpye 2881 #endif
331     }else{
332     /* Boundary crossed twice in one direction. Very unlikey! */
333    
334     /* The aim of the following two lines is to set the value of
335     the boolean variable to such a value as if the boundary
336     was crossed in the opposite direction before */
337     CONSOLE_DEBUG("DOUBLE CROSS");
338 jpye 2869 if(!prevals){
339 jpye 2836 num_dvars = slv_get_num_solvers_dvars(integ->system);
340     dvl = slv_get_solvers_dvar_list(integ->system);
341     prevals = ASC_NEW_ARRAY(int32,num_dvars);
342 jpye 2869 for(c = 0; dvl[c] != NULL; c++)
343 jpye 2881 prevals[c] = dis_value(dvl[c]);
344 jpye 2836 }
345     bnd_set_ida_value(bnd, !bnd_cond_states[i]);
346 jpye 2881 if(!ida_log_solve(integ,lrslv_ind)){
347     CONSOLE_DEBUG("Error in logic solve in double-cross");
348 jpye 2869 return -1;
349 jpye 2881 }
350 jpye 2836 bnd_set_ida_value(bnd,bnd_cond_states[i]);
351 jpye 2881 #ifdef IDA_BND_DEBUG
352 jpye 2882 CONSOLE_DEBUG("After double-cross, set boundary '%s' to %d",name,bnd_ida_value(bnd)?1:0);
353 jpye 2881 #endif
354 jpye 2559 }
355 jpye 2881 /* flag this boundary as having previously been crossed */
356 jpye 2836 bnd_set_ida_first_cross(bnd,0);
357 jpye 2881 /* store this recent crossing-direction for future reference */
358     bnd_set_ida_incr(bnd,(rootsfound[i] > 0));
359 jpye 2882 #ifdef IDA_BND_DEBUG
360     ASC_FREE(name);
361     #endif
362 jpye 2559 }
363     }
364 jpye 2869 if(!ida_log_solve(integ,lrslv_ind)) return -1;
365 jpye 2559
366 jpye 2836 /* If there was a double crossing, because of ida_log_solve the previous values
367     of discrete variables may be equal to their current values, which would mean that
368     the discrete vars haven't changed their values, but actually they did. So here we
369     restore the previous values of those variables, which are not connected with the
370     boundary which was double-crossed. */
371 jpye 2869 if(prevals){
372     for(c = 0; dvl[c] != NULL; c++)
373     if(!(dis_value(dvl[c]) == prevals[c] && dis_value(dvl[c]) != dis_previous_value(dvl[c])))
374 jpye 2836 dis_set_previous_value(dvl[c],prevals[c]);
375     ASC_FREE(prevals);
376 jpye 2559 }
377    
378 jpye 2836 integrator_output_write(integ);
379     integrator_output_write_obs(integ);
380     integrator_output_write_event(integ);
381 jpye 2869 if(!some_dis_vars_changed(integ->system)){
382 jpye 2881 CONSOLE_DEBUG("Crossed boundary but no effect on system");
383 jpye 2836 /* Boundary crossing that has no effect on system */
384     ASC_FREE(bnd_prev_eval);
385 jpye 2559
386 jpye 2836 /* Reset the boundary flag */
387 jpye 2869 for(i = 0; i < num_bnds; i++)
388 jpye 2836 bnd_set_ida_crossed(enginedata->bndlist[i], 0);
389 jpye 2559 return 0;
390     }
391 jpye 2836 /* update the main system if required */
392     int events_triggered = 1;
393 jpye 2869 if(slv_get_num_solvers_events(integ->system)!=0) {
394     while(some_dis_vars_changed(integ->system) && events_triggered) {
395     if(ida_bnd_reanalyse(integ)) {
396     /* select QRSlv solver, and solve the system */
397     if(slv_select_solver(integ->system, qrslv_ind) == -1) {
398 jpye 2836 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Error attempting to load QRSlv");
399     }
400     #ifdef IDA_BND_DEBUG
401 jpye 2869 {
402     struct var_variable **vl = slv_get_solvers_var_list(integ->system);
403     int i, n=slv_get_num_solvers_vars(integ->system), first=1;
404     CONSOLE_DEBUG("In boundary problem: variables (active, incident, free):");
405     for(i=0;i<n;++i){
406     if(var_incident(vl[i])&&!var_fixed(vl[i])&&var_active(vl[i])){
407     char *name = var_make_name(integ->system,vl[i]);
408     fprintf(stderr,"%s%s",(first?"\t":", "),name);
409     first=0; ASC_FREE(name);
410     }
411     }
412     fprintf(stderr,"\n");
413     }
414     {
415     struct rel_relation **rl = slv_get_solvers_rel_list(integ->system);
416     int i, n=slv_get_num_solvers_rels(integ->system), first=1;
417     CONSOLE_DEBUG("...relations (equality, included, active):");
418     for(i=0;i<n;++i){
419     if(rel_equality(rl[i])&&rel_active(rl[i])&&rel_included(rl[i])){
420     char *name = rel_make_name(integ->system,rl[i]);
421     fprintf(stderr,"%s%s",(first?"\t":", "),name);
422     first=0; ASC_FREE(name);
423     }
424     }
425     fprintf(stderr,"\n");
426     }
427     #endif
428 jpye 2836 slv_presolve(integ->system);
429     slv_solve(integ->system);
430    
431     slv_get_status(integ->system, &status);
432     if (!status.converged) {
433 jpye 2869 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Non-convergence in "
434     "non-linear solver at boundary");
435     #ifdef IDA_BND_DEBUG
436     }else{
437     CONSOLE_DEBUG("Converged");
438     #endif
439 jpye 2836 }
440    
441     for (i = 0; i < num_bnds; i++) {
442     if (!rootsfound[i]) {
443     if ((bndman_real_eval(enginedata->bndlist[i]) > 0 && bnd_prev_eval[i] < 0) ||
444     (bndman_real_eval(enginedata->bndlist[i]) < 0 && bnd_prev_eval[i] > 0)) {
445     bnd_cond_states[i] = !bnd_cond_states[i];
446     if (bnd_prev_eval[i] > 0) {
447     bnd_set_ida_incr(enginedata->bndlist[i],0);
448     rootsfound[i] = -1;
449     }
450     else {
451     bnd_set_ida_incr(enginedata->bndlist[i],1);
452     rootsfound[i] = 1;
453     }
454     bnd_set_ida_value(enginedata->bndlist[i],bnd_cond_states[i]);
455     }
456     }
457     }
458     if (!ida_log_solve(integ,lrslv_ind)) return -1;
459    
460     }else events_triggered = 0;
461     if (ida_bnd_reanalyse_cont(integ)) return 2;
462     if (events_triggered) {
463     integrator_output_write(integ);
464     integrator_output_write_obs(integ);
465     integrator_output_write_event(integ);
466     }
467     }
468     }
469    
470     integrator_output_write(integ);
471     integrator_output_write_obs(integ);
472     integrator_output_write_event(integ);
473    
474     ASC_FREE(bnd_prev_eval);
475    
476     /* Reset the boundary flag */
477     for (i = 0; i < num_bnds; i++)
478     bnd_set_ida_crossed(enginedata->bndlist[i], 0);
479     return 1;
480 jpye 2559 }
481 jpye 2836

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