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

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

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 2868 by jpye, Tue Feb 17 07:34:59 2015 UTC revision 2869 by jpye, Tue Mar 24 13:25:23 2015 UTC
# Line 26  Line 26 
26  #include <ascend/system/block.h>  #include <ascend/system/block.h>
27  #include <ascend/system/bndman.h>  #include <ascend/system/bndman.h>
28    
29  /* #define IDA_BND_DEBUG */  #define IDA_BND_DEBUG
30    
31  /*  /*
32   *   *
# Line 81  static void ida_write_values(IntegratorS Line 81  static void ida_write_values(IntegratorS
81  }  }
82    
83  int ida_setup_lrslv(IntegratorSystem *integ, int qrslv_ind, int lrslv_ind) {  int ida_setup_lrslv(IntegratorSystem *integ, int qrslv_ind, int lrslv_ind) {
     int32 c;  
     struct dis_discrete **dvl;  
     slv_status_t status;  
   
84      ida_log_solve(integ, lrslv_ind);      ida_log_solve(integ, lrslv_ind);
85        if(some_dis_vars_changed(integ->system)){
86      dvl = slv_get_solvers_dvar_list(integ->system);          return ida_bnd_reanalyse_cont(integ);
87        }
     if (some_dis_vars_changed(integ->system)) return ida_bnd_reanalyse_cont(integ);  
88      return 0;      return 0;
89  }  }
90    
# Line 240  int ida_log_solve(IntegratorSystem *inte Line 235  int ida_log_solve(IntegratorSystem *inte
235          ERROR_REPORTER_HERE(ASC_PROG_ERR,"Error attempting to load LRSlv");          ERROR_REPORTER_HERE(ASC_PROG_ERR,"Error attempting to load LRSlv");
236      }      }
237  #ifdef IDA_BND_DEBUG  #ifdef IDA_BND_DEBUG
238      CONSOLE_DEBUG("Solver selected is '%s'",slv_solver_name      CONSOLE_DEBUG("Solver selected is '%s'"
239          (slv_get_selected_solver(integ->system)));          ,slv_solver_name(slv_get_selected_solver(integ->system))
240        );
241  #endif  #endif
242    
243      /* Flip LRSlv into ida mode */      /* Flip LRSlv into ida mode */
# Line 262  int ida_log_solve(IntegratorSystem *inte Line 258  int ida_log_solve(IntegratorSystem *inte
258      if (!status.converged) {      if (!status.converged) {
259          ERROR_REPORTER_HERE(ASC_PROG_ERR,"Non-convergence in logical solver.");          ERROR_REPORTER_HERE(ASC_PROG_ERR,"Non-convergence in logical solver.");
260          return 0;          return 0;
261    #ifdef IDA_BND_DEBUG
262        }else{
263            CONSOLE_DEBUG("Converged");
264    #endif
265      }      }
266      return 1;      return 1;
267  }  }
# Line 296  int ida_cross_boundary(IntegratorSystem Line 296  int ida_cross_boundary(IntegratorSystem
296      enginedata = integ->enginedata;      enginedata = integ->enginedata;
297      num_bnds = enginedata->nbnds;      num_bnds = enginedata->nbnds;
298      bnd_prev_eval = ASC_NEW_ARRAY(double,num_bnds);      bnd_prev_eval = ASC_NEW_ARRAY(double,num_bnds);
299      for (i = 0; i < num_bnds; i++) {      for(i = 0; i < num_bnds; i++) {
300          bnd_prev_eval[i] = bndman_real_eval(enginedata->bndlist[i]);          bnd_prev_eval[i] = bndman_real_eval(enginedata->bndlist[i]);
301          if (rootsfound[i]) {          if(rootsfound[i]) {
302              bnd = enginedata->bndlist[i];              bnd = enginedata->bndlist[i];
303              bnd_set_ida_crossed(bnd, 1);              bnd_set_ida_crossed(bnd, 1);
304              if (bnd_ida_first_cross(bnd) || !((rootsfound[i] == 1 && bnd_ida_incr(bnd)) || (rootsfound[i] == -1 && !bnd_ida_incr(bnd))) ) {              if(bnd_ida_first_cross(bnd) || !((rootsfound[i] == 1 && bnd_ida_incr(bnd)) || (rootsfound[i] == -1 && !bnd_ida_incr(bnd))) ) {
305                  if (bnd_cond_states[i] == 0) {                  if(bnd_cond_states[i] == 0) {
306                      bnd_set_ida_value(bnd, 1);                      bnd_set_ida_value(bnd, 1);
307                      bnd_cond_states[i] = 1;                      bnd_cond_states[i] = 1;
308                  } else {                  }else{
309                      bnd_set_ida_value(bnd, 0);                      bnd_set_ida_value(bnd, 0);
310                      bnd_cond_states[i] = 0;                      bnd_cond_states[i] = 0;
311                  }                  }
312              }else{ /* Boundary crossed twice in one direction              }else{ /* Boundary crossed twice in one direction
313                    This is very unlikely to happen */                    This is very unlikely to happen */
314                  /* The aim of the following two lines is to set                  /* The aim of the following two lines is to set
315                                     the value of the boolean variable to such a                      the value of the boolean variable to such a
316                                     value as if the boundary was crossed in the                      value as if the boundary was crossed in the
317                                     opposite direction before */                      opposite direction before */
318                  if (!prevals) {                  if(!prevals){
319                      num_dvars = slv_get_num_solvers_dvars(integ->system);                      num_dvars = slv_get_num_solvers_dvars(integ->system);
320                      dvl = slv_get_solvers_dvar_list(integ->system);                      dvl = slv_get_solvers_dvar_list(integ->system);
321                      prevals = ASC_NEW_ARRAY(int32,num_dvars);                      prevals = ASC_NEW_ARRAY(int32,num_dvars);
322                      for (c = 0; dvl[c] != NULL; c++)                      for(c = 0; dvl[c] != NULL; c++)
323                          prevals[c] = dis_value(dvl[c]);                          prevals[c] = dis_value(dvl[c]);
324                  }                  }
325                  bnd_set_ida_value(bnd, !bnd_cond_states[i]);                  bnd_set_ida_value(bnd, !bnd_cond_states[i]);
326                  if (!ida_log_solve(integ,lrslv_ind)) return -1;                  if(!ida_log_solve(integ,lrslv_ind))
327                        return -1;
328                  bnd_set_ida_value(bnd,bnd_cond_states[i]);                  bnd_set_ida_value(bnd,bnd_cond_states[i]);
329              }              }
330              bnd_set_ida_first_cross(bnd,0);              bnd_set_ida_first_cross(bnd,0);
331              if (rootsfound[i]>0) bnd_set_ida_incr(bnd,1);              if
332              else bnd_set_ida_incr(bnd,0);                  (rootsfound[i]>0) bnd_set_ida_incr(bnd,1);
333                else
334                    bnd_set_ida_incr(bnd,0);
335          }          }
336      }      }
337      if (!ida_log_solve(integ,lrslv_ind)) return -1;      if(!ida_log_solve(integ,lrslv_ind)) return -1;
338    
339      /* If there was a double crossing, because of ida_log_solve the previous values      /* If there was a double crossing, because of ida_log_solve the previous values
340      of discrete variables may be equal to their current values, which would mean that      of discrete variables may be equal to their current values, which would mean that
341      the discrete vars haven't changed their values, but actually they did. So here we      the discrete vars haven't changed their values, but actually they did. So here we
342      restore the previous values of those variables, which are not connected with the      restore the previous values of those variables, which are not connected with the
343      boundary which was double-crossed. */      boundary which was double-crossed. */
344      if (prevals) {      if(prevals){
345          for (c = 0; dvl[c] != NULL; c++)          for(c = 0; dvl[c] != NULL; c++)
346              if (!(dis_value(dvl[c]) == prevals[c] && dis_value(dvl[c]) != dis_previous_value(dvl[c])))              if(!(dis_value(dvl[c]) == prevals[c] && dis_value(dvl[c]) != dis_previous_value(dvl[c])))
347                  dis_set_previous_value(dvl[c],prevals[c]);                  dis_set_previous_value(dvl[c],prevals[c]);
348          ASC_FREE(prevals);          ASC_FREE(prevals);
349      }      }
# Line 348  int ida_cross_boundary(IntegratorSystem Line 351  int ida_cross_boundary(IntegratorSystem
351      integrator_output_write(integ);      integrator_output_write(integ);
352      integrator_output_write_obs(integ);      integrator_output_write_obs(integ);
353      integrator_output_write_event(integ);      integrator_output_write_event(integ);
354      if (!some_dis_vars_changed(integ->system)){      if(!some_dis_vars_changed(integ->system)){
355          /* Boundary crossing that has no effect on system */          /* Boundary crossing that has no effect on system */
   
356          ASC_FREE(bnd_prev_eval);          ASC_FREE(bnd_prev_eval);
357    
358          /* Reset the boundary flag */          /* Reset the boundary flag */
359          for (i = 0; i < num_bnds; i++)          for(i = 0; i < num_bnds; i++)
360              bnd_set_ida_crossed(enginedata->bndlist[i], 0);              bnd_set_ida_crossed(enginedata->bndlist[i], 0);
361          return 0;          return 0;
362      }      }
363      /* update the main system if required */      /* update the main system if required */
364      int events_triggered = 1;      int events_triggered = 1;
365      if (slv_get_num_solvers_events(integ->system)!=0) {      if(slv_get_num_solvers_events(integ->system)!=0) {
366          while (some_dis_vars_changed(integ->system) && events_triggered)  {          while(some_dis_vars_changed(integ->system) && events_triggered)  {
367                if(ida_bnd_reanalyse(integ)) {
368              if (ida_bnd_reanalyse(integ)) {                  /* select QRSlv solver, and solve the system */
369                    if(slv_select_solver(integ->system, qrslv_ind) == -1) {
                 if (slv_select_solver(integ->system, qrslv_ind) == -1) {  
370                      ERROR_REPORTER_HERE(ASC_PROG_ERR,"Error attempting to load QRSlv");                      ERROR_REPORTER_HERE(ASC_PROG_ERR,"Error attempting to load QRSlv");
371                  }                  }
372  #ifdef IDA_BND_DEBUG  #ifdef IDA_BND_DEBUG
373                  CONSOLE_DEBUG("Solver selected is '%s'",slv_solver_name  # if 0
374                      (slv_get_selected_solver(integ->system)));                  CONSOLE_DEBUG("Solver selected is '%s'"
375                        ,slv_solver_name(slv_get_selected_solver(integ->system))
376                    );
377    #endif
378                    {
379                        struct var_variable **vl = slv_get_solvers_var_list(integ->system);
380                        int i, n=slv_get_num_solvers_vars(integ->system), first=1;
381                        CONSOLE_DEBUG("In boundary problem: variables (active, incident, free):");
382                        for(i=0;i<n;++i){
383                            if(var_incident(vl[i])&&!var_fixed(vl[i])&&var_active(vl[i])){
384                                char *name = var_make_name(integ->system,vl[i]);
385                                fprintf(stderr,"%s%s",(first?"\t":", "),name);
386                                first=0; ASC_FREE(name);
387                            }
388                        }
389                        fprintf(stderr,"\n");
390                    }
391                    {
392                        struct rel_relation **rl = slv_get_solvers_rel_list(integ->system);
393                        int i, n=slv_get_num_solvers_rels(integ->system), first=1;
394                        CONSOLE_DEBUG("...relations (equality, included, active):");
395                        for(i=0;i<n;++i){
396                            if(rel_equality(rl[i])&&rel_active(rl[i])&&rel_included(rl[i])){
397                                char *name = rel_make_name(integ->system,rl[i]);
398                                fprintf(stderr,"%s%s",(first?"\t":", "),name);
399                                first=0; ASC_FREE(name);
400                            }
401                        }
402                        fprintf(stderr,"\n");
403                    }
404  #endif  #endif
405                  slv_presolve(integ->system);                  slv_presolve(integ->system);
406                  slv_solve(integ->system);                  slv_solve(integ->system);
407    
408                  slv_get_status(integ->system, &status);                  slv_get_status(integ->system, &status);
409                  if (!status.converged) {                  if (!status.converged) {
410                      ERROR_REPORTER_HERE(ASC_PROG_ERR,"Non-convergence in non-linear solver at "                      ERROR_REPORTER_HERE(ASC_PROG_ERR,"Non-convergence in "
411                          "boundary");                          "non-linear solver at boundary");
412    #ifdef IDA_BND_DEBUG
413                    }else{
414                        CONSOLE_DEBUG("Converged");
415    #endif
416                  }                  }
417    
418                  for (i = 0; i < num_bnds; i++) {                  for (i = 0; i < num_bnds; i++) {

Legend:
Removed from v.2868  
changed lines
  Added in v.2869

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