/[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 927 by johnpye, Mon Nov 6 07:49:06 2006 UTC revision 928 by johnpye, Wed Nov 22 10:32:18 2006 UTC
# Line 83  static symchar *g_symbols[4]; Line 83  static symchar *g_symbols[4];
83  struct Integ_var_t {  struct Integ_var_t {
84    long index;    long index;
85    long type;    long type;
86    struct var_variable *i;    struct Integ_var_t *derivative;
87    struct Integ_var_t *derivative_of;    struct Integ_var_t *derivative_of;
88    struct var_variable *derivative;    struct var_variable *i;
89    int varindx; /**< index into slv_get_master_vars_list, or -1 if not there */    int varindx; /**< index into slv_get_master_vars_list, or -1 if not there */
90    int isstate;    int isstate;
91  };  };
# Line 334  int integrator_analyse_dae(IntegratorSys Line 334  int integrator_analyse_dae(IntegratorSys
334      int i, j;      int i, j;
335      int numstates;      int numstates;
336      int numy, nrels;      int numy, nrels;
337        int yindex;
338      int maxderiv;      int maxderiv;
339    
340      CONSOLE_DEBUG("Starting DAE analysis");      CONSOLE_DEBUG("Starting DAE analysis");
# Line 377  int integrator_analyse_dae(IntegratorSys Line 378  int integrator_analyse_dae(IntegratorSys
378          info = (struct Integ_var_t *)gl_fetch(sys->dynvars, i);          info = (struct Integ_var_t *)gl_fetch(sys->dynvars, i);
379          if(info->type==1 || info->type==0)numstates++;          if(info->type==1 || info->type==0)numstates++;
380          if(maxderiv < info->type - 1)maxderiv = info->type - 1;          if(maxderiv < info->type - 1)maxderiv = info->type - 1;
381          varname = var_make_name(sys->system,info->i);          /* varname = var_make_name(sys->system,info->i);
382          /* CONSOLE_DEBUG("var[%d] = \"%s\": ode_index = %ld",i,varname,info->type); */          CONSOLE_DEBUG("var[%d] = \"%s\": ode_index = %ld",i,varname,info->type);
383          ASC_FREE(varname);          ASC_FREE(varname); */
384      }      }
385      if(maxderiv == 0){      if(maxderiv == 0){
386          ERROR_REPORTER_HERE(ASC_USER_ERROR,"No derivatives found (check 'ode_type' values for your vars).");          ERROR_REPORTER_HERE(ASC_USER_ERROR,"No derivatives found (check 'ode_type' values for your vars).");
387          return 0;          return 0;
388      }      }
389      if(numstates == 0){      if(maxderiv > 1){
390          ERROR_REPORTER_HERE(ASC_USER_ERROR,"No states found (check 'odetype' values for your vars).");          ERROR_REPORTER_HERE(ASC_USER_ERROR,"Higher-order derivatives found. You must provide a reduced order formulation for your model.");
391          return 0;          return 0;
392      }      }
393    
   
394      if(!integrator_check_indep_var(sys))return 0;      if(!integrator_check_indep_var(sys))return 0;
395    
396      gl_sort(sys->dynvars,(CmpFunc)Integ_CmpDynVars);      gl_sort(sys->dynvars,(CmpFunc)Integ_CmpDynVars);
397    
398      /*      fprintf(stderr,"\n\n\nSORTED VARS\n");
399      for(i=1; i<=gl_length(sys->dynvars); ++i){      for(i=1; i<=gl_length(sys->dynvars); ++i){
400          info = (struct Integ_var_t *)gl_fetch(sys->dynvars, i);          info = (struct Integ_var_t *)gl_fetch(sys->dynvars, i);
401          varname = var_make_name(sys->system,info->i);          varname = var_make_name(sys->system,info->i);
402          // CONSOLE_DEBUG("var[%d] = \"%s\": ode_type = %ld",i,varname,info->type);          CONSOLE_DEBUG("var[%d] = \"%s\": ode_type = %ld",i,varname,info->type);
403          ASC_FREE(varname);          ASC_FREE(varname);
404      }*/      }
   
     /* link up derivative chains */  
405    
406        /* link up variables with their derivatives */
407      prev = NULL;      prev = NULL;
408      for(i=1; i<=gl_length(sys->dynvars); ++i){ /* why does gl_list index with base 1??? */      for(i=1; i<=gl_length(sys->dynvars); ++i){ /* why does gl_list index with base 1??? */
409          info = (struct Integ_var_t *)gl_fetch(sys->dynvars, i);          info = (struct Integ_var_t *)gl_fetch(sys->dynvars, i);
410          info->derivative = NULL;          
   
         derivname = var_make_name(sys->system,info->i);  
         if(prev!=NULL){  
             varname = var_make_name(sys->system,prev->i);  
         }else{  
             varname = NULL;  
         }  
         /* CONSOLE_DEBUG("current = \"%s\", previous = \"%s\"",derivname,varname); */  
         ASC_FREE(derivname);  
         if(varname)ASC_FREE(varname);  
   
411          if(info->type == INTEG_STATE_VAR || info->type == INTEG_ALGEBRAIC_VAR){          if(info->type == INTEG_STATE_VAR || info->type == INTEG_ALGEBRAIC_VAR){
412              varname = var_make_name(sys->system,info->i);              varname = var_make_name(sys->system,info->i);
413              /* CONSOLE_DEBUG("Var \"%s\" is not a derivative",varname); */              CONSOLE_DEBUG("Var \"%s\" is an algebraic variable",varname);
414              ASC_FREE(varname);              ASC_FREE(varname);
             info->derivative_of = NULL;  
415              info->type = INTEG_STATE_VAR;              info->type = INTEG_STATE_VAR;
416                info->derivative_of = NULL;
417          }else{          }else{
418              if(prev==NULL || info->index != prev->index){              if(prev==NULL || info->index != prev->index){
419                  /* CONSOLE_DEBUG("current current type = %ld",info->type); */                  /* derivative, but without undifferentiated var present in model */
                 /* if(prev!=NULL){  
                     CONSOLE_DEBUG("current index = %ld, previous = %ld",info->index,prev->index);  
                 }else{  
                     CONSOLE_DEBUG("current index = %ld, current type = %ld",info->index,info->type);  
                 } */  
420                  varname = var_make_name(sys->system,info->i);                  varname = var_make_name(sys->system,info->i);
421                  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"
422                      , info->type-1                      , info->type-1
# Line 442  int integrator_analyse_dae(IntegratorSys Line 425  int integrator_analyse_dae(IntegratorSys
425                  ASC_FREE(varname);                  ASC_FREE(varname);
426                  return 0;                  return 0;
427              }else if(info->type != prev->type + 1){              }else if(info->type != prev->type + 1){
428                    /* derivative, but missing the next-lower-order derivative */
429                  derivname = var_make_name(sys->system,info->i);                  derivname = var_make_name(sys->system,info->i);
430                  varname = var_make_name(sys->system,prev->i);                  varname = var_make_name(sys->system,prev->i);
431                  ERROR_REPORTER_HERE(ASC_USER_ERROR                  ERROR_REPORTER_HERE(ASC_USER_ERROR
# Line 454  int integrator_analyse_dae(IntegratorSys Line 438  int integrator_analyse_dae(IntegratorSys
438                  ASC_FREE(derivname);                  ASC_FREE(derivname);
439                  return 0;                  return 0;
440              }else{              }else{
441                    /* variable with derivative */
442                  varname = var_make_name(sys->system,prev->i);                  varname = var_make_name(sys->system,prev->i);
443                  derivname = var_make_name(sys->system,info->i);                  derivname = var_make_name(sys->system,info->i);
444                  CONSOLE_DEBUG("Var \"%s\" is the derivative of \"%s\"",derivname,varname);                  CONSOLE_DEBUG("Var \"%s\" is the derivative of \"%s\"",derivname,varname);
445                  ASC_FREE(varname);                  ASC_FREE(varname);
446                  ASC_FREE(derivname);                  ASC_FREE(derivname);
447                  info->derivative_of = prev;                  info->derivative_of = prev;
                 numy++;  
448              }              }
449          }          }
450          prev = info;          prev = info;
451      }      }
452    
453      /* record which vars have derivatives and which don't */      /* record which vars have derivatives and which don't, and count 'states' */
     for(i=1; i<=gl_length(sys->dynvars); ++i){  
         info = (struct Integ_var_t *)gl_fetch(sys->dynvars, i);  
         if(info->derivative_of){  
             info->derivative_of->derivative = info->i;  
             /* varname = var_make_name(sys->system,info->derivative_of->i);  
             derivname = var_make_name(sys->system,info->derivative_of->derivative);  
             CONSOLE_DEBUG("Var \"%s\" is the derivative of \"%s\"",derivname,varname);  
             ASC_FREE(varname);  
             ASC_FREE(derivname); */  
         }  
     }  
   
     CONSOLE_DEBUG("Indentifying states...");  
   
     /* count numy: either it's a state, or it has a higher-order derivative */  
454      numy = 0;      numy = 0;
455      for(i=1; i<=gl_length(sys->dynvars); ++i){      for(i=1; i<=gl_length(sys->dynvars); ++i){
456          info = (struct Integ_var_t *)gl_fetch(sys->dynvars, i);          info = (struct Integ_var_t *)gl_fetch(sys->dynvars, i);
457          if(info->type == INTEG_STATE_VAR || info->type == INTEG_ALGEBRAIC_VAR || info->derivative != NULL){          if(info->derivative_of){
458              varname = var_make_name(sys->system,info->i);              info->derivative_of->derivative = info;
             if(var_fixed(info->i)){  
                 CONSOLE_DEBUG("Var \"%s\" is a FIXED state variable",varname);  
             }else{  
                 CONSOLE_DEBUG("Var \"%s\" is a state variable",varname);  
             }  
             ASC_FREE(varname);  
             info->isstate = 1;  
             numy++;  
459          }else{          }else{
460              varname = var_make_name(sys->system,info->i);              numy++;
             CONSOLE_DEBUG("Var \"%s\" is a NON-STATE derivative",varname);  
             ASC_FREE(varname);  
             info->isstate = 0;  
461          }          }
462      }      }
463    
464      /*      /* allocate storage for the 'y' and 'ydot' arrays */
         create lists 'y' and 'ydot'. some elements of ydot don't correspond  
         to variables in our model: these are the algebraic vars.  
     */  
   
     CONSOLE_DEBUG("Identified %d state variables", numy);  
   
465      sys->y = ASC_NEW_ARRAY(struct var_variable *,numy);      sys->y = ASC_NEW_ARRAY(struct var_variable *,numy);
466      sys->ydot = ASC_NEW_ARRAY(struct var_variable *,numy);      sys->ydot = ASC_NEW_ARRAY(struct var_variable *,numy);
467        sys->y_id = ASC_NEW_ARRAY(int, slv_get_num_master_vars(sys->system));
468    
469      /*      /* now add variables and their derivatives to 'ydot' and 'y' */
470          at this point we know there are no missing derivatives etc, so we      yindex = 0;
471          can use (i-1) as the index into y and ydot. any variable with      
472          'derivative_of' set to null is a state variable... but it might already      for(i=1; i<=gl_length(sys->dynvars); ++i){
         be getting added  
     */  
     for(j=0, i=1; i<=gl_length(sys->dynvars); ++i){  
         info = (struct Integ_var_t *)gl_fetch(sys->dynvars, i);  
         if(!info->isstate)continue;  
         varname = var_make_name(sys->system,info->i);  
         if(info->derivative == NULL){  
             CONSOLE_DEBUG("Pure algebraic: %s",varname);  
             sys->y[j] = info->i; /* a pure algebraic variable */  
             sys->ydot[j] = NULL;  
         }else{  
             CONSOLE_DEBUG("Differential variable: %s",varname);  
             sys->y[j] = info->i; /* a variable whose derivative is present in the model */  
             sys->ydot[j] = info->derivative;  
         }  
         ASC_FREE(varname);  
         ++j;  
     }  
   
     /*  
         set up the y_id table so that given a 'variable number' from relman_diff2,  
         we can work out where that fits in our y and ydot vectors.  
   
         There's something a bit fishy about fetching the varlist at this late stage...  
     */  
   
     varlist = slv_get_solvers_var_list(sys->system);  
     nvarlist = slv_get_num_solvers_vars(sys->system);  
   
     CONSOLE_DEBUG("WORKING THROUGH THE SOLVER'S VAR LIST %d",nvarlist);  
   
     sys->y_id = ASC_NEW_ARRAY(long, nvarlist);  
     for(i=0; i< nvarlist; ++i){  
         sys->y_id[i] = -2;  
     }  
   
     CONSOLE_DEBUG("WORKING THROUGH %ld DYNVARS",gl_length(sys->dynvars));  
   
     for(i=1; i <= gl_length(sys->dynvars); ++i){  
473          info = (struct Integ_var_t *)gl_fetch(sys->dynvars, i);          info = (struct Integ_var_t *)gl_fetch(sys->dynvars, i);
474          sys->y_id[info->varindx] = i;          if(info->derivative_of)continue;
475          varname = var_make_name(sys->system,info->i);                    if(info->derivative){
476          if(info->varindx > 0){              sys->y[yindex] = info->i;
477              CONSOLE_DEBUG("Variable dynvars[%d] = y_id[%d] = '%s'",i,info->varindx,varname);              sys->ydot[yindex] = info->derivative->i;
478                if(info->varindx >= 0){
479                    sys->y_id[info->varindx] = yindex;
480                    CONSOLE_DEBUG("y_id[%d] = %d",info->varindx,yindex);
481                }
482                if(info->derivative->varindx >= 0){
483                    sys->y_id[info->derivative->varindx] = -1-yindex;
484                    CONSOLE_DEBUG("y_id[%d] = %d",info->derivative->varindx,-1-yindex);
485                }
486          }else{          }else{
487              CONSOLE_DEBUG("VARIABLE dynvars[%d] = y_id[%d] = '%s' NOT FOUND IN VARLIST",i,info->varindx,varname);              sys->y[yindex] = info ->i;
488          }              sys->ydot[yindex] = NULL;
489          ASC_FREE(varname);              if(info->varindx >= 0){
490      }                  sys->y_id[info->varindx] = yindex;
491      for(i=0; i< nvarlist; ++i){                  CONSOLE_DEBUG("y_id[%d] = %d",info->varindx,yindex);
492          if(sys->y_id[i] < 0){              }
             varname = var_make_name(sys->system,varlist[i]);  
             CONSOLE_DEBUG("UNCONNECTED SOLVER VAR varlist[%d] = '%s' (%s)",i,varname,var_fixed(varlist[i])?"fixed":"not fixed");  
             ASC_FREE(varname);  
493          }          }
494            yindex++;
495      }      }
496    
497      nrels = slv_get_num_solvers_rels(sys->system);      nrels = slv_get_num_solvers_rels(sys->system);
# Line 583  int integrator_analyse_dae(IntegratorSys Line 503  int integrator_analyse_dae(IntegratorSys
503          return 0;          return 0;
504      }      }
505    
506        CONSOLE_DEBUG("THERE ARE %d VARIABLES IN THE INTEGRATION SYSTEM",numy);
507        
508      sys->n_y = numy;      sys->n_y = numy;
509    
510      if(!integrator_sort_obs_vars(sys))return 0;      if(!integrator_sort_obs_vars(sys))return 0;
# Line 795  static void integrator_print_var_stats(I Line 717  static void integrator_print_var_stats(I
717  }  }
718    
719  /**  /**
720        Check sanity of the independent variable.
721    
722      @return 1 on success      @return 1 on success
723  */  */
724  static int integrator_check_indep_var(IntegratorSystem *sys){  static int integrator_check_indep_var(IntegratorSystem *sys){
# Line 879  void integrator_dae_classify_var(Integra Line 803  void integrator_dae_classify_var(Integra
803              CONSOLE_DEBUG("VARIABLE IS NOT ACTIVE");              CONSOLE_DEBUG("VARIABLE IS NOT ACTIVE");
804              return;              return;
805          }          }
806            
807            /* only non-fixed variables are accepted */
808          if(!var_fixed(var)){          if(!var_fixed(var)){
809              /* get the ode_type and ode_id of this solver_var */              /* get the ode_type and ode_id of this solver_var */
810              type = DynamicVarInfo(var,&index);              type = DynamicVarInfo(var,&index);
# Line 893  void integrator_dae_classify_var(Integra Line 818  void integrator_dae_classify_var(Integra
818                  INTEG_ADD_TO_LIST(info,type,index,var,varindx,sys->dynvars);                  INTEG_ADD_TO_LIST(info,type,index,var,varindx,sys->dynvars);
819              }              }
820          }          }
821  #if 1  #if 0
822  else{          else{
823              /* fixed variable, only include it if ode_type == 1 */              /* fixed variable, only include it if ode_type == 1 */
824              type = DynamicVarInfo(var,&index);              type = DynamicVarInfo(var,&index);
825              if(type==INTEG_STATE_VAR){              if(type==INTEG_STATE_VAR){
Legend:
Removed from v.927  
changed lines
  Added in v.928
john.pye@anu.edu.au
 ViewVC Help
Powered by ViewVC 1.1.22