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

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

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

revision 2881 by jpye, Sun Mar 29 05:31:03 2015 UTC revision 2882 by jpye, Sun Mar 29 12:20:00 2015 UTC
# Line 739  int ida_setup_IC(IntegratorSystem *integ Line 739  int ida_setup_IC(IntegratorSystem *integ
739      char *varname;      char *varname;
740  #endif  #endif
741    
742    #ifdef IDA_BND_DEBUG
743        CONSOLE_DEBUG("Solving initial conditions...");
744    #endif
745    
746      icopt = 0;      icopt = 0;
747      if(strcmp(SLV_PARAM_CHAR(&integ->params,IDA_PARAM_CALCIC), "Y") == 0) {      if(strcmp(SLV_PARAM_CHAR(&integ->params,IDA_PARAM_CALCIC), "Y") == 0) {
748  #ifdef SOLVE_DEBUG  #ifdef SOLVE_DEBUG
# Line 746  int ida_setup_IC(IntegratorSystem *integ Line 750  int ida_setup_IC(IntegratorSystem *integ
750  #endif  #endif
751          icopt = IDA_Y_INIT;          icopt = IDA_Y_INIT;
752          asc_assert(icopt!=0);          asc_assert(icopt!=0);
753      }else if (strcmp(SLV_PARAM_CHAR(&integ->params,IDA_PARAM_CALCIC), "YA_YDP")      }else if(0==strcmp(SLV_PARAM_CHAR(&integ->params,IDA_PARAM_CALCIC), "YA_YDP")){
             == 0) {  
754  #ifdef SOLVE_DEBUG  #ifdef SOLVE_DEBUG
755          CONSOLE_DEBUG("Solving initial conditions using values of yd");          CONSOLE_DEBUG("Solving initial conditions using values of yd");
756  #endif  #endif
# Line 1016  static int integrator_ida_solve(Integrat Line 1019  static int integrator_ida_solve(Integrat
1019      integrator_ida_debug(integ, stderr);      integrator_ida_debug(integ, stderr);
1020  #endif  #endif
1021    
   
1022      /* Initialise boundary condition states if appropriate. Reconfigure if necessary */      /* Initialise boundary condition states if appropriate. Reconfigure if necessary */
1023      if(enginedata->nbnds){      if(enginedata->nbnds){
1024          CONSOLE_DEBUG("Initialising boundary states");          CONSOLE_DEBUG("Initialising boundary states");
# Line 1027  static int integrator_ida_solve(Integrat Line 1029  static int integrator_ida_solve(Integrat
1029  #endif  #endif
1030          bnd_cond_states = ASC_NEW_ARRAY_CLEAR(int,enginedata->nbnds);          bnd_cond_states = ASC_NEW_ARRAY_CLEAR(int,enginedata->nbnds);
1031    
1032            /* identify if we're exactly *on* any boundaries currently */
1033          for(i = 0; i < enginedata->nbnds; i++) {          for(i = 0; i < enginedata->nbnds; i++) {
1034    #ifdef IDA_BND_DEBUG
1035                relname = bnd_make_name(integ->system,enginedata->bndlist[i]);
1036    #endif
1037              bnd_cond_states[i] = bndman_calc_satisfied(enginedata->bndlist[i]);              bnd_cond_states[i] = bndman_calc_satisfied(enginedata->bndlist[i]);
1038                bnd_set_ida_first_cross(enginedata->bndlist[i],1);
1039              if(bndman_real_eval(enginedata->bndlist[i]) == 0) {              if(bndman_real_eval(enginedata->bndlist[i]) == 0) {
1040                    /* if the residual for the boundary is zero (ie looks like we are *on* the boundary?) JP */
1041    #ifdef IDA_BND_DEBUG
1042                    CONSOLE_DEBUG("Boundary '%s': not set",relname);
1043    #endif
1044                  bnd_not_set[i] = 1;                  bnd_not_set[i] = 1;
1045                  all_bnds_set = 0;                  all_bnds_set = 0;
1046              }else              }else{
1047                  bnd_not_set[i] = 0;                  bnd_not_set[i] = 0;
1048              bnd_set_ida_first_cross(enginedata->bndlist[i],1);              }
1049  #ifdef IDA_BND_DEBUG  #ifdef IDA_BND_DEBUG
1050                  char *n = bnd_make_name(integ->system,enginedata->bndlist[i]);              CONSOLE_DEBUG("Boundary '%s' is %d",relname,bnd_cond_states[i]);
1051                  CONSOLE_DEBUG("Conditional '%s' is %d",n,bnd_cond_states[i]);              ASC_FREE(relname);
                 ASC_FREE(n);  
1052  #endif  #endif
1053          }          }
1054          CONSOLE_DEBUG("Setting up LRSlv...");          CONSOLE_DEBUG("Setting up LRSlv...");
# Line 1057  static int integrator_ida_solve(Integrat Line 1067  static int integrator_ida_solve(Integrat
1067    
1068      /* Setup parameter inputs and initial conditions for IDA. */      /* Setup parameter inputs and initial conditions for IDA. */
1069      tout = samplelist_get(integ->samples, start_index + 1);      tout = samplelist_get(integ->samples, start_index + 1);
1070        /* solve the initial conditions, allocate memory, other stuff... */
1071      ida_prepare_integrator(integ, ida_mem, tout);      ida_prepare_integrator(integ, ida_mem, tout);
1072    
1073      tol = 0.0001*(samplelist_get(integ->samples, finish_index)      tol = 0.0001*(samplelist_get(integ->samples, finish_index)
# Line 1096  static int integrator_ida_solve(Integrat Line 1107  static int integrator_ida_solve(Integrat
1107                  rootdir[i] = 0;                  rootdir[i] = 0;
1108  #ifdef IDA_BND_DEBUG  #ifdef IDA_BND_DEBUG
1109                  char *n = bnd_make_name(integ->system,enginedata->bndlist[i]);                  char *n = bnd_make_name(integ->system,enginedata->bndlist[i]);
1110                  CONSOLE_DEBUG("Boundary '%s': value=%d; (trigger dirn=%s)"                  CONSOLE_DEBUG("Boundary '%s': bnd_cond_states[%d]=%d, bndman_calc_satisfied=%d; (trigger dirn=%s)"
1111                      ,n, bndman_calc_satisfied(enginedata->bndlist[i])                      ,n, i, bnd_cond_states[i], bndman_calc_satisfied(enginedata->bndlist[i])
1112                      ,rootdir[i]==1?"UP":(rootdir[i]==-1?"DOWN":"both")                      ,rootdir[i]==1?"UP":(rootdir[i]==-1?"DOWN":"both")
1113                  );                  );
1114                  ASC_FREE(n);                  ASC_FREE(n);
# Line 1187  static int integrator_ida_solve(Integrat Line 1198  static int integrator_ida_solve(Integrat
1198                      }                      }
1199  #endif  #endif
1200    
1201                        
1202                      if(all_bnds_set == 0){                      if(all_bnds_set == 0){
1203    #ifdef IDA_BND_DEBUG
1204                            CONSOLE_DEBUG("Unset bounds exist; evaluate them explicitly...");
1205    #endif
1206                          all_bnds_set = 1;                          all_bnds_set = 1;
1207                          for(i = 0; i < enginedata->nbnds; i++){                          for(i = 0; i < enginedata->nbnds; i++){
1208                              if(bnd_not_set[i]){                              if(bnd_not_set[i]){
1209                                  if(!rootsfound[i])                                  if(!rootsfound[i]){
1210                                      bnd_cond_states[i] = bndman_calc_satisfied(enginedata->bndlist[i]);                                      bnd_cond_states[i] = bndman_calc_satisfied(enginedata->bndlist[i]);
1211                                  else all_bnds_set = 0;  #ifdef IDA_BND_DEBUG
1212                                        relname = bnd_make_name(integ->system, enginedata->bndlist[i]);
1213                                        CONSOLE_DEBUG("Boundary '%s': bnd_cond_states[%d] = %d"
1214                                            ,relname,i,bnd_cond_states[i]
1215                                        );
1216                                        ASC_FREE(relname);
1217    #endif
1218                                    }else all_bnds_set = 0;
1219                              }                              }
1220                          }                          }
1221                      }                      }
1222    
1223                      if(!after_root){                      if(!after_root){
1224    #ifdef IDA_BND_DEBUG
1225                            CONSOLE_DEBUG("Just 'after_root'...");
1226                            for(i=0;i<enginedata->nbnds;++i){
1227                                relname = bnd_make_name(integ->system, enginedata->bndlist[i]);
1228                                CONSOLE_DEBUG("Boundary '%s': bnd_cond_states[%d] = %d"
1229                                    ,relname,i,bnd_cond_states[i]
1230                                );
1231                                ASC_FREE(relname);
1232                            }
1233    #endif
1234                          need_to_reconfigure = ida_cross_boundary(integ, rootsfound,                          need_to_reconfigure = ida_cross_boundary(integ, rootsfound,
1235                              bnd_cond_states, qrslv_ind, lrslv_ind);                              bnd_cond_states, qrslv_ind, lrslv_ind);
1236                      }                      }
# Line 1266  static int integrator_ida_solve(Integrat Line 1298  static int integrator_ida_solve(Integrat
1298              for(i = 0; i < enginedata->nbnds; i++) {              for(i = 0; i < enginedata->nbnds; i++) {
1299                  if(bnd_not_set[i]) {                  if(bnd_not_set[i]) {
1300                      bnd_cond_states[i] = bndman_calc_satisfied(enginedata->bndlist[i]);                      bnd_cond_states[i] = bndman_calc_satisfied(enginedata->bndlist[i]);
1301    #ifdef IDA_BND_DEBUG
1302                        char *n = bnd_make_name(integ->system,enginedata->bndlist[i]);
1303                        CONSOLE_DEBUG("Boundary '%s' not set; satisfied=%d",n,bnd_cond_states[i]);
1304                        ASC_FREE(n);
1305    #endif
1306                  }                  }
1307              }              }
1308          }          }

Legend:
Removed from v.2881  
changed lines
  Added in v.2882

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