/[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 951 by johnpye, Sun Nov 26 01:36:49 2006 UTC revision 952 by johnpye, Tue Nov 28 23:01:50 2006 UTC
# Line 40  Line 40 
40      @TODO this needs to go away.      @TODO this needs to go away.
41  */  */
42  #define INTEG_DEBUG TRUE  #define INTEG_DEBUG TRUE
43    /* #define ANALYSE_DEBUG */
44    
45  /**  /**
46      Print a debug message with value if INTEG_DEBUG is TRUE.      Print a debug message with value if INTEG_DEBUG is TRUE.
# Line 398  int integrator_analyse_dae(IntegratorSys Line 399  int integrator_analyse_dae(IntegratorSys
399    
400      integrator_visit_system_vars(sys,&integrator_dae_classify_var);      integrator_visit_system_vars(sys,&integrator_dae_classify_var);
401    
402    #ifdef ANALYSE_DEBUG
403      CONSOLE_DEBUG("Found %lu observation variables:",gl_length(sys->obslist));      CONSOLE_DEBUG("Found %lu observation variables:",gl_length(sys->obslist));
404      for(i=1; i<=gl_length(sys->obslist); ++i){      for(i=1; i<=gl_length(sys->obslist); ++i){
405          info = (struct Integ_var_t *)gl_fetch(sys->obslist, i);          info = (struct Integ_var_t *)gl_fetch(sys->obslist, i);
# Line 405  int integrator_analyse_dae(IntegratorSys Line 407  int integrator_analyse_dae(IntegratorSys
407          CONSOLE_DEBUG("observation[%d] = \"%s\"",i,varname);          CONSOLE_DEBUG("observation[%d] = \"%s\"",i,varname);
408          ASC_FREE(varname);          ASC_FREE(varname);
409      }      }
410    #endif
411    
412      /* CONSOLE_DEBUG("Checking found vars..."); */      /* CONSOLE_DEBUG("Checking found vars..."); */
413      if(gl_length(sys->dynvars)==0){      if(gl_length(sys->dynvars)==0){
# Line 438  int integrator_analyse_dae(IntegratorSys Line 441  int integrator_analyse_dae(IntegratorSys
441      gl_sort(sys->dynvars,(CmpFunc)Integ_CmpDynVars);      gl_sort(sys->dynvars,(CmpFunc)Integ_CmpDynVars);
442      /* !!! NOTE THAT dynvars IS NOW REARRANGED !!! *sigh* bug hunting(!) */      /* !!! NOTE THAT dynvars IS NOW REARRANGED !!! *sigh* bug hunting(!) */
443    
444    #ifdef ANALYSE_DEBUG
445      CONSOLE_DEBUG("Variables rearranged to increasing ODE ID & TYPE (varindx = matrix order)");      CONSOLE_DEBUG("Variables rearranged to increasing ODE ID & TYPE (varindx = matrix order)");
446      for(i=1; i<=gl_length(sys->dynvars); ++i){      for(i=1; i<=gl_length(sys->dynvars); ++i){
447          info = (struct Integ_var_t *)gl_fetch(sys->dynvars, i);          info = (struct Integ_var_t *)gl_fetch(sys->dynvars, i);
# Line 445  int integrator_analyse_dae(IntegratorSys Line 449  int integrator_analyse_dae(IntegratorSys
449          CONSOLE_DEBUG("var[%d] = \"%s\": ode_type = %ld (varindx = %d)",i-1,varname,info->type,info->varindx);          CONSOLE_DEBUG("var[%d] = \"%s\": ode_type = %ld (varindx = %d)",i-1,varname,info->type,info->varindx);
450          ASC_FREE(varname);          ASC_FREE(varname);
451      }      }
452    #endif
453    
454      /* link up variables with their derivatives */      /* link up variables with their derivatives */
455      prev = NULL;      prev = NULL;
# Line 452  int integrator_analyse_dae(IntegratorSys Line 457  int integrator_analyse_dae(IntegratorSys
457          info = (struct Integ_var_t *)gl_fetch(sys->dynvars, i);          info = (struct Integ_var_t *)gl_fetch(sys->dynvars, i);
458                    
459          if(info->type == INTEG_STATE_VAR || info->type == INTEG_ALGEBRAIC_VAR){          if(info->type == INTEG_STATE_VAR || info->type == INTEG_ALGEBRAIC_VAR){
460    #ifdef ANALYSE_DEBUG
461              varname = var_make_name(sys->system,info->i);              varname = var_make_name(sys->system,info->i);
462              CONSOLE_DEBUG("Var \"%s\" is an algebraic variable",varname);              CONSOLE_DEBUG("Var \"%s\" is an algebraic variable",varname);
463              ASC_FREE(varname);              ASC_FREE(varname);
464    #endif
465              info->type = INTEG_STATE_VAR;              info->type = INTEG_STATE_VAR;
466              info->derivative_of = NULL;              info->derivative_of = NULL;
467          }else{          }else{
468              if(prev==NULL || info->index != prev->index){              if(prev==NULL || info->index != prev->index){
469                  /* derivative, but without undifferentiated var present in model */                  /* derivative, but without undifferentiated var present in model */
470    #ifdef ANALYSE_DEBUG
471                  varname = var_make_name(sys->system,info->i);                  varname = var_make_name(sys->system,info->i);
472                  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"
473                      , info->type-1                      , info->type-1
474                      , varname                      , varname
475                  );                  );
476                  ASC_FREE(varname);                  ASC_FREE(varname);
477    #endif
478                  return 0;                  return 0;
479              }else if(info->type != prev->type + 1){              }else if(info->type != prev->type + 1){
480                  /* derivative, but missing the next-lower-order derivative */                  /* derivative, but missing the next-lower-order derivative */
481    #ifdef ANALYSE_DEBUG
482                  derivname = var_make_name(sys->system,info->i);                  derivname = var_make_name(sys->system,info->i);
483                  varname = var_make_name(sys->system,prev->i);                  varname = var_make_name(sys->system,prev->i);
484                  ERROR_REPORTER_HERE(ASC_USER_ERROR                  ERROR_REPORTER_HERE(ASC_USER_ERROR
# Line 479  int integrator_analyse_dae(IntegratorSys Line 489  int integrator_analyse_dae(IntegratorSys
489                  );                  );
490                  ASC_FREE(varname);                  ASC_FREE(varname);
491                  ASC_FREE(derivname);                  ASC_FREE(derivname);
492    #endif
493                  return 0;                  return 0;
494              }else{              }else{
495                  /* variable with derivative */                  /* variable with derivative */
496    #ifdef ANALYSE_DEBUG
497                  varname = var_make_name(sys->system,prev->i);                  varname = var_make_name(sys->system,prev->i);
498                  derivname = var_make_name(sys->system,info->i);                  derivname = var_make_name(sys->system,info->i);
499                  CONSOLE_DEBUG("Var \"%s\" is the derivative of \"%s\"",derivname,varname);                  CONSOLE_DEBUG("Var \"%s\" is the derivative of \"%s\"",derivname,varname);
500                  ASC_FREE(varname);                  ASC_FREE(varname);
501                  ASC_FREE(derivname);                  ASC_FREE(derivname);
502    #endif
503                  info->derivative_of = prev;                  info->derivative_of = prev;
504              }              }
505          }          }
# Line 516  int integrator_analyse_dae(IntegratorSys Line 529  int integrator_analyse_dae(IntegratorSys
529      /* now add variables and their derivatives to 'ydot' and 'y' */      /* now add variables and their derivatives to 'ydot' and 'y' */
530      yindex = 0;      yindex = 0;
531            
532    #ifdef ANALYSE_DEBUG
533      CONSOLE_DEBUG("VARS IN THEIR MATRIX ORDER");      CONSOLE_DEBUG("VARS IN THEIR MATRIX ORDER");
534    #endif
535      for(i=1; i<=gl_length(sys->dynvars); ++i){      for(i=1; i<=gl_length(sys->dynvars); ++i){
536          info = (struct Integ_var_t *)gl_fetch(sys->dynvars, i);          info = (struct Integ_var_t *)gl_fetch(sys->dynvars, i);
537          if(info->derivative_of)continue;          if(info->derivative_of)continue;
# Line 526  int integrator_analyse_dae(IntegratorSys Line 541  int integrator_analyse_dae(IntegratorSys
541              sys->ydot[yindex] = info->derivative->i;              sys->ydot[yindex] = info->derivative->i;
542              if(info->varindx >= 0){              if(info->varindx >= 0){
543                  sys->y_id[info->varindx] = yindex;                  sys->y_id[info->varindx] = yindex;
544    #ifdef ANALYSE_DEBUG
545                  CONSOLE_DEBUG("y_id[%d] = %d",info->varindx,yindex);                  CONSOLE_DEBUG("y_id[%d] = %d",info->varindx,yindex);
546    #endif
547              }              }
548              if(info->derivative->varindx >= 0){              if(info->derivative->varindx >= 0){
549                  sys->y_id[info->derivative->varindx] = -1-yindex;                  sys->y_id[info->derivative->varindx] = -1-yindex;
550    #ifdef ANALYSE_DEBUG
551                  CONSOLE_DEBUG("y_id[%d] = %d",info->derivative->varindx,-1-yindex);                  CONSOLE_DEBUG("y_id[%d] = %d",info->derivative->varindx,-1-yindex);
552    #endif
553              }              }
554          }else{          }else{
555              sys->y[yindex] = info ->i;              sys->y[yindex] = info ->i;
556              sys->ydot[yindex] = NULL;              sys->ydot[yindex] = NULL;
557              if(info->varindx >= 0){              if(info->varindx >= 0){
558                  sys->y_id[info->varindx] = yindex;                  sys->y_id[info->varindx] = yindex;
559    #ifdef ANALYSE_DEBUG
560                  CONSOLE_DEBUG("y_id[%d] = %d",info->varindx,yindex);                  CONSOLE_DEBUG("y_id[%d] = %d",info->varindx,yindex);
561    #endif
562              }              }
563          }          }
564          yindex++;          yindex++;
# Line 558  int integrator_analyse_dae(IntegratorSys Line 579  int integrator_analyse_dae(IntegratorSys
579    
580      if(!integrator_sort_obs_vars(sys))return 0;      if(!integrator_sort_obs_vars(sys))return 0;
581    
582    #ifdef ANALYSE_DEBUG
583      CONSOLE_DEBUG("RESULTS OF ANALYSIS");      CONSOLE_DEBUG("RESULTS OF ANALYSIS");
584      fprintf(stderr,"index\ty\tydot\n");      fprintf(stderr,"index\ty\tydot\n");
585      fprintf(stderr,"-----\t-----\t-----\n");      fprintf(stderr,"-----\t-----\t-----\n");
# Line 579  int integrator_analyse_dae(IntegratorSys Line 601  int integrator_analyse_dae(IntegratorSys
601      fprintf(stderr,"index\tname\ty_id\ty\tydot\n");      fprintf(stderr,"index\tname\ty_id\ty\tydot\n");
602      fprintf(stderr,"-----\t-----\t-----\t-----\t-----\n");      fprintf(stderr,"-----\t-----\t-----\t-----\t-----\n");
603      integrator_visit_system_vars(sys,integrator_dae_show_var);      integrator_visit_system_vars(sys,integrator_dae_show_var);
604    #endif
605    
606      return 1;      return 1;
607  }  }

Legend:
Removed from v.951  
changed lines
  Added in v.952

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