/[ascend]/trunk/base/generic/solver/integrator.c
ViewVC logotype

Diff of /trunk/base/generic/solver/integrator.c

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

revision 902 by johnpye, Mon Oct 23 01:07:58 2006 UTC revision 903 by johnpye, Wed Oct 25 13:07:12 2006 UTC
# Line 330  int integrator_analyse_dae(IntegratorSys Line 330  int integrator_analyse_dae(IntegratorSys
330      assert(blsys->indepvars==NULL);      assert(blsys->indepvars==NULL);
331    
332      blsys->indepvars = gl_create(10L);  /* t var info */      blsys->indepvars = gl_create(10L);  /* t var info */
333      blsys->dynvars = gl_create(200L);  /* y ydot var info */      blsys->dynvars = gl_create(200L);  /* y and ydot var info */
334      blsys->obslist = gl_create(100L);  /* obs info */      blsys->obslist = gl_create(100L);  /* obs info */
335    
336      if(blsys->indepvars==NULL      if(blsys->indepvars==NULL
# Line 363  int integrator_analyse_dae(IntegratorSys Line 363  int integrator_analyse_dae(IntegratorSys
363      numstates = 0;      numstates = 0;
364      for(i=1; i<=gl_length(blsys->dynvars); ++i){      for(i=1; i<=gl_length(blsys->dynvars); ++i){
365          info = (struct Integ_var_t *)gl_fetch(blsys->dynvars, i);          info = (struct Integ_var_t *)gl_fetch(blsys->dynvars, i);
366          if(info->index==1 || info->index==0)numstates++;          if(info->type==1 || info->type==0)numstates++;
367          if(maxderiv < info->type - 1)maxderiv = info->type - 1;          if(maxderiv < info->type - 1)maxderiv = info->type - 1;
368          varname = var_make_name(blsys->system,info->i);          varname = var_make_name(blsys->system,info->i);
369          CONSOLE_DEBUG("var[%d] = \"%s\": order = %ld",i,varname,info->type-1);          /* CONSOLE_DEBUG("var[%d] = \"%s\": ode_index = %ld",i,varname,info->type); */
370          ASC_FREE(varname);          ASC_FREE(varname);
371      }      }
372      if(maxderiv == 0){      if(maxderiv == 0){
373          ERROR_REPORTER_HERE(ASC_USER_ERROR,"No derivatives found (check 'ode_type' values in your vars).");          ERROR_REPORTER_HERE(ASC_USER_ERROR,"No derivatives found (check 'ode_type' values for your vars).");
374          return 0;          return 0;
375      }      }
376      if(numstates == 0){      if(numstates == 0){
377          ERROR_REPORTER_HERE(ASC_USER_ERROR,"No states found (check 'odetype' values in your vars).");          ERROR_REPORTER_HERE(ASC_USER_ERROR,"No states found (check 'odetype' values for your vars).");
378          return 0;          return 0;
379      }      }
380    
# Line 388  int integrator_analyse_dae(IntegratorSys Line 388  int integrator_analyse_dae(IntegratorSys
388      for(i=1; i<=gl_length(blsys->dynvars); ++i){      for(i=1; i<=gl_length(blsys->dynvars); ++i){
389          info = (struct Integ_var_t *)gl_fetch(blsys->dynvars, i);          info = (struct Integ_var_t *)gl_fetch(blsys->dynvars, i);
390          varname = var_make_name(blsys->system,info->i);          varname = var_make_name(blsys->system,info->i);
391          CONSOLE_DEBUG("var[%d] = \"%s\": order = %ld",i,varname,info->type-1);          /* CONSOLE_DEBUG("var[%d] = \"%s\": ode_type = %ld",i,varname,info->type); */
392          ASC_FREE(varname);          ASC_FREE(varname);
393      }      }
394    
# Line 405  int integrator_analyse_dae(IntegratorSys Line 405  int integrator_analyse_dae(IntegratorSys
405          }else{          }else{
406              varname = NULL;              varname = NULL;
407          }          }
408          CONSOLE_DEBUG("current = \"%s\", previous = \"%s\"",derivname,varname);          /* CONSOLE_DEBUG("current = \"%s\", previous = \"%s\"",derivname,varname); */
409          ASC_FREE(derivname);          ASC_FREE(derivname);
410          if(varname)ASC_FREE(varname);          if(varname)ASC_FREE(varname);
411    
412          if(info->type == INTEG_STATE_VAR || info->type == INTEG_ALGEBRAIC_VAR){          if(info->type == INTEG_STATE_VAR || info->type == INTEG_ALGEBRAIC_VAR){
413              varname = var_make_name(blsys->system,info->i);              varname = var_make_name(blsys->system,info->i);
414              CONSOLE_DEBUG("Var \"%s\" is not a derivative",varname);              /* CONSOLE_DEBUG("Var \"%s\" is not a derivative",varname); */
415              ASC_FREE(varname);              ASC_FREE(varname);
416              info->derivative_of = NULL;              info->derivative_of = NULL;
417              info->type = INTEG_STATE_VAR;              info->type = INTEG_STATE_VAR;
418          }else{          }else{
419              if(prev==NULL || info->index != prev->index){              if(prev==NULL || info->index != prev->index){
420                  CONSOLE_DEBUG("current current type = %ld",info->type);                  /* CONSOLE_DEBUG("current current type = %ld",info->type); */
421                  if(prev!=NULL){                  /* if(prev!=NULL){
422                      CONSOLE_DEBUG("current index = %ld, previous = %ld",info->index,prev->index);                      CONSOLE_DEBUG("current index = %ld, previous = %ld",info->index,prev->index);
423                  }else{                  }else{
424                      CONSOLE_DEBUG("current index = %ld, current type = %ld",info->index,info->type);                      CONSOLE_DEBUG("current index = %ld, current type = %ld",info->index,info->type);
425                  }                  } */
426                  varname = var_make_name(blsys->system,info->i);                  varname = var_make_name(blsys->system,info->i);
427                  ERROR_REPORTER_HERE(ASC_USER_ERROR,"Derivative %d of \"%s\" is present without its un-differentiated equivalent"                  ERROR_REPORTER_HERE(ASC_USER_ERROR,"Derivative %d of \"%s\" is present without its un-differentiated equivalent"
428                      , info->type-1                      , info->type-1
# Line 460  int integrator_analyse_dae(IntegratorSys Line 460  int integrator_analyse_dae(IntegratorSys
460          info = (struct Integ_var_t *)gl_fetch(blsys->dynvars, i);          info = (struct Integ_var_t *)gl_fetch(blsys->dynvars, i);
461          if(info->derivative_of){          if(info->derivative_of){
462              info->derivative_of->derivative = info->i;              info->derivative_of->derivative = info->i;
463                /* varname = var_make_name(blsys->system,info->derivative_of->i);
464                derivname = var_make_name(blsys->system,info->derivative_of->derivative);
465                CONSOLE_DEBUG("Var \"%s\" is the derivative of \"%s\"",derivname,varname);
466                ASC_FREE(varname);
467                ASC_FREE(derivname); */
468          }          }
469      }      }
470    
# Line 471  int integrator_analyse_dae(IntegratorSys Line 476  int integrator_analyse_dae(IntegratorSys
476          info = (struct Integ_var_t *)gl_fetch(blsys->dynvars, i);          info = (struct Integ_var_t *)gl_fetch(blsys->dynvars, i);
477          if(info->type == INTEG_STATE_VAR || info->type == INTEG_ALGEBRAIC_VAR || info->derivative != NULL){          if(info->type == INTEG_STATE_VAR || info->type == INTEG_ALGEBRAIC_VAR || info->derivative != NULL){
478              varname = var_make_name(blsys->system,info->i);              varname = var_make_name(blsys->system,info->i);
479              CONSOLE_DEBUG("Var \"%s\" is a state variable",varname);              if(var_fixed(info->i)){
480                    CONSOLE_DEBUG("Var \"%s\" is a FIXED state variable",varname);
481                }else{
482                    CONSOLE_DEBUG("Var \"%s\" is a state variable",varname);
483                }          
484              ASC_FREE(varname);              ASC_FREE(varname);
485              info->isstate = 1;              info->isstate = 1;
486              numy++;              numy++;
487          }else{          }else{
488                varname = var_make_name(blsys->system,info->i);
489                CONSOLE_DEBUG("Var \"%s\" is a NON-STATE derivative",varname);
490                ASC_FREE(varname);
491              info->isstate = 0;              info->isstate = 0;
492          }          }
493      }      }
# Line 816  void integrator_dae_classify_var(Integra Line 828  void integrator_dae_classify_var(Integra
828                  /* any other type of var is in the DAE system, at least for now */                  /* any other type of var is in the DAE system, at least for now */
829                  INTEG_ADD_TO_LIST(info,type,index,var,blsys->dynvars);                  INTEG_ADD_TO_LIST(info,type,index,var,blsys->dynvars);
830              }              }
831            }else{
832                /* fixed variable, only include it if ode_type == 1 */
833                type = DynamicVarInfo(var,&index);
834                if(type==INTEG_STATE_VAR){
835                    INTEG_ADD_TO_LIST(info,type,index,var,blsys->dynvars);
836                }
837          }          }
838    
839          /* if the var's obs_id > 0, add it to the observation list */          /* if the var's obs_id > 0, add it to the observation list */
# Line 1164  void integrator_set_ydot(IntegratorSyste Line 1182  void integrator_set_ydot(IntegratorSyste
1182          if(blsys->ydot[i]!=NULL){          if(blsys->ydot[i]!=NULL){
1183              var_set_value(blsys->ydot[i],dydx[i]);              var_set_value(blsys->ydot[i],dydx[i]);
1184  #ifndef NDEBUG  #ifndef NDEBUG
1185              varname = var_make_name(blsys->system, blsys->ydot[i]);              /* varname = var_make_name(blsys->system, blsys->ydot[i]);
1186              CONSOLE_DEBUG("ydot[%ld] = \"%s\" = %g --> ASCEND", i+1, varname, dydx[i]);              CONSOLE_DEBUG("ydot[%ld] = \"%s\" = %g --> ASCEND", i+1, varname, dydx[i]);
1187              ASC_FREE(varname);              ASC_FREE(varname); */
1188  #endif  #endif
1189          }          }
1190  #ifndef NDEBUG  #ifndef NDEBUG
1191          else{          /*else{
1192              CONSOLE_DEBUG("ydot[%ld] = %g (internal)", i+1, dydx[i]);              CONSOLE_DEBUG("ydot[%ld] = %g (internal)", i+1, dydx[i]);
1193          }          }*/
1194  #endif  #endif
1195      }      }
1196  }  }
# Line 1339  int integrator_checkstatus(slv_status_t Line 1357  int integrator_checkstatus(slv_status_t
1357      return 1;      return 1;
1358    }    }
1359    if (status.diverged) {    if (status.diverged) {
1360      FPRINTF(stderr, "The derivative system did not converge.\nIntegration "      FPRINTF(stderr, "The derivative system did not converge. Integration will terminate.");
         "will be terminated at the end of the current step.");  
1361      return 0;      return 0;
1362    }    }
1363    if (status.inconsistent) {    if (status.inconsistent) {
1364      FPRINTF(stderr, "A numerical inconsistency was discovered while "      FPRINTF(stderr, "A numerically inconsistent state was discovered while "
1365          "calculating derivatives. Integration will be terminated at "          "calculating derivatives. Integration will terminate.");
         "the end of the current step.");  
1366      return 0;      return 0;
1367    }    }
1368    if (status.time_limit_exceeded) {    if (status.time_limit_exceeded) {
1369      FPRINTF(stderr, "The time limit was exceeded while calculating "      FPRINTF(stderr, "The time limit was exceeded while calculating "
1370          "derivatives.\nIntegration will be terminated at the end of "          "derivatives. Integration will terminate.");
         "the current step.");  
1371      return 0;      return 0;
1372    }    }
1373    if (status.iteration_limit_exceeded) {    if (status.iteration_limit_exceeded) {
1374      FPRINTF(stderr, "The iteration limit was exceeded while calculating "      FPRINTF(stderr, "The iteration limit was exceeded while calculating "
1375          "derivatives.\nIntegration will be terminated at "          "derivatives. Integration will terminate.");
         "the end of the current step.");  
1376      return 0;      return 0;
1377    }    }
1378    if (status.panic) {    if (status.panic) {
1379      FPRINTF(stderr, "The user patience limit was exceeded while "      FPRINTF(stderr, "The user patience limit was exceeded while "
1380          "calculating derivatives.\nIntegration will be terminated "          "calculating derivatives. Integration will terminate.");
         "at the end of the current step.");  
1381      return 0;      return 0;
1382    }    }
1383    return 0;    return 0;
Legend:
Removed from v.902  
changed lines
  Added in v.903
john.pye@anu.edu.au
 ViewVC Help
Powered by ViewVC 1.1.22