/[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 976 by johnpye, Tue Dec 19 00:57:26 2006 UTC revision 977 by johnpye, Wed Dec 20 00:39:52 2006 UTC
# Line 41  Line 41 
41      @TODO this needs to go away.      @TODO this needs to go away.
42  */  */
43  #define INTEG_DEBUG TRUE  #define INTEG_DEBUG TRUE
44  /* #define ANALYSE_DEBUG */  #define ANALYSE_DEBUG
45    
46  /**  /**
47      Print a debug message with value if INTEG_DEBUG is TRUE.      Print a debug message with value if INTEG_DEBUG is TRUE.
# Line 79  static symchar *g_symbols[3]; Line 79  static symchar *g_symbols[3];
79   * that maintains this info by backpointers, but oh well until then.   * that maintains this info by backpointers, but oh well until then.
80   */   */
81  #define OBSINDEX g_symbols[2]  #define OBSINDEX g_symbols[2]
82  /* Integer child. All variables with OBSINDEX !=0 will be recorded in  /* Integer child. All variables with OBSINDEX !=0 will be sent to the
83   * the blsode output file. Tis someone else's job to grok this output.      IntegratorOutputWriteObsFn allowing output to a file, graph, console, etc.
84   */   */
85    
86    
# Line 103  struct Integ_var_t { Line 103  struct Integ_var_t {
103  void integrator_create_engine(IntegratorSystem *sys);  void integrator_create_engine(IntegratorSystem *sys);
104  void integrator_free_engine(IntegratorSystem *sys);  void integrator_free_engine(IntegratorSystem *sys);
105    
106  static int integrator_analyse_ode(IntegratorSystem *sys);  IntegratorAnalyseFn integrator_analyse_ode;
107  static int integrator_analyse_dae(IntegratorSystem *sys);  IntegratorAnalyseFn integrator_analyse_dae;
108    
109  typedef void (IntegratorVarVisitorFn)(IntegratorSystem *sys, struct var_variable *var, const int *varindx);  typedef void (IntegratorVarVisitorFn)(IntegratorSystem *sys, struct var_variable *var, const int *varindx);
110  static void integrator_visit_system_vars(IntegratorSystem *sys,IntegratorVarVisitorFn *visitor);  static void integrator_visit_system_vars(IntegratorSystem *sys,IntegratorVarVisitorFn *visitor);
# Line 144  IntegratorSystem *integrator_new(slv_sys Line 144  IntegratorSystem *integrator_new(slv_sys
144      sys->system = slvsys;      sys->system = slvsys;
145      sys->instance = inst;      sys->instance = inst;
146    
147        sys->engine = INTEG_UNKNOWN;
148        sys->internals = NULL;
149    
150      sys->states = NULL; sys->derivs = NULL;      sys->states = NULL; sys->derivs = NULL;
151      sys->dynvars = NULL; sys->obslist = NULL; sys->indepvars = NULL;      sys->dynvars = NULL; sys->obslist = NULL; sys->indepvars = NULL;
152    
# Line 207  static void IntegInitSymbols(void){ Line 210  static void IntegInitSymbols(void){
210      At present, integrators are all present at compile-time so this is a static      At present, integrators are all present at compile-time so this is a static
211      list; we just return the list.      list; we just return the list.
212  */  */
213  void integrator_get_engines(IntegratorLookup **listptr){  const IntegratorLookup *integrator_get_engines(){
214  #define S ,  #define S ,
215  #define I(N) {INTEG_##N, #N}  #define I(N,P) {INTEG_##N, P.name}
216      static const IntegratorLookup lookup[] = {      static const IntegratorLookup list[] = {
217          INTEG_LIST          INTEG_LIST
218          ,{INTEG_UNKNOWN,NULL}          ,{INTEG_UNKNOWN,NULL}
219      };      };
     *listptr = lookup;  
220  #undef S  #undef S
221  #undef I  #undef I
222        return list;
223  }  }
224    
225  /* return 0 on success */  /* return 0 on success */
226  int integrator_set_engine(IntegratorSystem *sys, IntegratorEngine engine){  int integrator_set_engine(IntegratorSystem *sys, IntegratorEngine engine){
227    
228        CONSOLE_DEBUG("Setting engine...");
229    
230      /* verify integrator type ok. always passes for nonNULL inst. */      /* verify integrator type ok. always passes for nonNULL inst. */
231      if(engine==INTEG_UNKNOWN){      if(engine==INTEG_UNKNOWN){
232          ERROR_REPORTER_NOLINE(ASC_USER_ERROR          ERROR_REPORTER_NOLINE(ASC_USER_ERROR
# Line 239  int integrator_set_engine(IntegratorSyst Line 244  int integrator_set_engine(IntegratorSyst
244          integrator_free_engine(sys);          integrator_free_engine(sys);
245      }      }
246      sys->engine = engine;      sys->engine = engine;
247        switch(sys->engine){
248    #ifdef ASC_WITH_IDA
249            case INTEG_IDA: sys->internals = &integrator_ida_internals; break;
250    #endif
251            case INTEG_LSODE: sys->internals = &integrator_lsode_internals; break;
252            case INTEG_AWW: sys->internals = &integrator_aww_internals; break;
253            default:
254                ERROR_REPORTER_HERE(ASC_PROG_ERR,"Unknown integrator engine");
255                sys->internals = NULL;
256                return 1;
257        };
258        
259        asc_assert(sys->internals);
260      integrator_create_engine(sys);      integrator_create_engine(sys);
   
261      return 0;      return 0;
262  }  }
263    
# Line 253  IntegratorEngine integrator_get_engine(c Line 270  IntegratorEngine integrator_get_engine(c
270      this system. Note that this data is pointed to by sys->enginedata.      this system. Note that this data is pointed to by sys->enginedata.
271  */  */
272  void integrator_free_engine(IntegratorSystem *sys){  void integrator_free_engine(IntegratorSystem *sys){
273      switch(sys->engine){      if(sys->engine==INTEG_UNKNOWN)return;
274          case INTEG_LSODE: integrator_lsode_free(sys->enginedata); break;      if(sys->enginedata){
275  #ifdef ASC_WITH_IDA          if(sys->internals){
276          case INTEG_IDA: integrator_ida_free(sys->enginedata); break;              (sys->internals->freefn)(sys->enginedata);
277  #endif              sys->enginedata=NULL;
278          case INTEG_AWW: integrator_aww_free(sys->enginedata); break;          }else{
279          default: break;              ERROR_REPORTER_HERE(ASC_PROG_ERR,"Unable to free engine data: no sys->internals");
280            }
281      }      }
     sys->enginedata=NULL;  
282  }  }
283    
284  /**  /**
# Line 272  void integrator_free_engine(IntegratorSy Line 289  void integrator_free_engine(IntegratorSy
289      during the solve stage (and freed inside integrator_free_engine)      during the solve stage (and freed inside integrator_free_engine)
290  */  */
291  void integrator_create_engine(IntegratorSystem *sys){  void integrator_create_engine(IntegratorSystem *sys){
292      if(sys->enginedata!=NULL)return;      asc_assert(sys->engine!=INTEG_UNKNOWN);
293      switch(sys->engine){      asc_assert(sys->internals);
294          case INTEG_LSODE: integrator_lsode_create(sys); break;      asc_assert(sys->enginedata==NULL);
295  #ifdef ASC_WITH_IDA      (sys->internals->createfn)(sys);
         case INTEG_IDA: integrator_ida_create(sys); break;  
 #endif  
         case INTEG_AWW: integrator_aww_create(sys); break;  
         default: break;  
     }  
296  }  }
297    
298  /*------------------------------------------------------------------------------  /*------------------------------------------------------------------------------
# Line 293  void integrator_create_engine(Integrator Line 305  void integrator_create_engine(Integrator
305    
306      @return 0 on success, 1 on error      @return 0 on success, 1 on error
307  */  */
308  static int integrator_params_default(IntegratorSystem *sys){  int integrator_params_default(IntegratorSystem *sys){
309      switch(sys->engine){      asc_assert(sys->engine!=INTEG_UNKNOWN);
310          case INTEG_LSODE: return integrator_lsode_params_default(sys);      asc_assert(sys->internals);
311  #ifdef ASC_WITH_IDA      return (sys->internals->paramsdefaultfn)(sys);
         case INTEG_IDA: return integrator_ida_params_default(sys);  
 #endif  
         case INTEG_AWW: return integrator_aww_params_default(sys);  
         default: return 0;  
     }  
312  }  }
313    
314  int integrator_params_get(const IntegratorSystem *sys, slv_parameters_t *parameters){  int integrator_params_get(const IntegratorSystem *sys, slv_parameters_t *parameters){
# Line 364  int integrator_find_indep_var(Integrator Line 371  int integrator_find_indep_var(Integrator
371      @return 1 on success      @return 1 on success
372  */  */
373  int integrator_analyse(IntegratorSystem *sys){  int integrator_analyse(IntegratorSystem *sys){
374      switch(sys->engine){      CONSOLE_DEBUG("Analysing integration system...");
375          case INTEG_LSODE: return integrator_analyse_ode(sys);      asc_assert(sys);
376  #ifdef ASC_WITH_IDA      if(sys->engine==INTEG_UNKNOWN){
377          case INTEG_IDA: return integrator_analyse_dae(sys);          ERROR_REPORTER_HERE(ASC_PROG_ERR,"No engine selected: can't analyse");
378  #endif          return 0;
         case INTEG_AWW: return integrator_aww_analyse(sys);  
         case INTEG_UNKNOWN:  
             ERROR_REPORTER_HERE(ASC_PROG_ERR,"No engine selected: can't analyse");  
         default:  
             ERROR_REPORTER_HERE(ASC_PROG_ERR  
                 ,"The selected integration engine (%d) is not available"  
                 ,sys->engine  
             );  
379      }      }
380      return 0;      asc_assert(sys->engine!=INTEG_UNKNOWN);
381        asc_assert(sys->internals);
382        return (sys->internals->analysefn)(sys);
383  }  }
384    
385  /**  /**
# Line 546  int integrator_analyse_dae(IntegratorSys Line 547  int integrator_analyse_dae(IntegratorSys
547      for(i=0; i<numy; ++i){      for(i=0; i<numy; ++i){
548          asc_assert(sys->ydot[i]==NULL);          asc_assert(sys->ydot[i]==NULL);
549      }      }
550        
551        for(i=0; i<numy; ++i){
552            sys->y_id[i] = 9999999L;
553        }
554    
555      /* now add variables and their derivatives to 'ydot' and 'y' */      /* now add variables and their derivatives to 'ydot' and 'y' */
556      yindex = 0;      yindex = 0;
# Line 561  int integrator_analyse_dae(IntegratorSys Line 566  int integrator_analyse_dae(IntegratorSys
566              assert(info->derivative);              assert(info->derivative);
567              sys->ydot[yindex] = info->derivative->i;              sys->ydot[yindex] = info->derivative->i;
568              if(info->varindx >= 0){              if(info->varindx >= 0){
569                    ASC_ASSERT_RANGE(yindex, -1e7L, 1e7L);
570                  sys->y_id[info->varindx] = yindex;                  sys->y_id[info->varindx] = yindex;
571  #ifdef ANALYSE_DEBUG  #ifdef ANALYSE_DEBUG
572                  CONSOLE_DEBUG("y_id[%d] = %d",info->varindx,yindex);                  CONSOLE_DEBUG("y_id[%d] = %d",info->varindx,yindex);
# Line 568  int integrator_analyse_dae(IntegratorSys Line 574  int integrator_analyse_dae(IntegratorSys
574              }              }
575              if(info->derivative->varindx >= 0){              if(info->derivative->varindx >= 0){
576                  sys->y_id[info->derivative->varindx] = -1-yindex;                  sys->y_id[info->derivative->varindx] = -1-yindex;
577                    ASC_ASSERT_RANGE(-1-yindex,-numy,0);
578  #ifdef ANALYSE_DEBUG  #ifdef ANALYSE_DEBUG
579                  CONSOLE_DEBUG("y_id[%d] = %d",info->derivative->varindx,-1-yindex);                  CONSOLE_DEBUG("y_id[%d] = %d",info->derivative->varindx,-1-yindex);
580  #endif  #endif
# Line 577  int integrator_analyse_dae(IntegratorSys Line 584  int integrator_analyse_dae(IntegratorSys
584              sys->ydot[yindex] = NULL;              sys->ydot[yindex] = NULL;
585              if(info->varindx >= 0){              if(info->varindx >= 0){
586                  sys->y_id[info->varindx] = yindex;                  sys->y_id[info->varindx] = yindex;
587                    ASC_ASSERT_RANGE(yindex,0,numy);
588  #ifdef ANALYSE_DEBUG  #ifdef ANALYSE_DEBUG
589                  CONSOLE_DEBUG("y_id[%d] = %d",info->varindx,yindex);                  CONSOLE_DEBUG("y_id[%d] = %d",info->varindx,yindex);
590  #endif  #endif
# Line 631  void integrator_dae_show_var(IntegratorS Line 639  void integrator_dae_show_var(IntegratorS
639          , struct var_variable *var, const int *varindx          , struct var_variable *var, const int *varindx
640  ){  ){
641      char *varname;      char *varname;
642      int y_id;      long y_id;
643      varname = var_make_name(sys->system, var);      varname = var_make_name(sys->system, var);
644      if(varindx==NULL){      if(varindx==NULL){
645          fprintf(stderr,".\t%s\n",varname);          fprintf(stderr,".\t%s\n",varname);
# Line 639  void integrator_dae_show_var(IntegratorS Line 647  void integrator_dae_show_var(IntegratorS
647          return;          return;
648      }      }
649      y_id = sys->y_id[*varindx];      y_id = sys->y_id[*varindx];
650      fprintf(stderr,"%d\t%s\t%d", *varindx, varname,y_id);  
651        fprintf(stderr,"%d\t%s\t%ld", *varindx, varname,y_id);
652      ASC_FREE(varname);      ASC_FREE(varname);
653      if(y_id >= 0){      if(y_id >= 0){
654          fprintf(stderr,"\ty[%d]\t.\n",y_id);          fprintf(stderr,"\ty[%ld]\t.\n",y_id);
655      }else{      }else{
656          fprintf(stderr,"\t.\tydot[%d]\n",-y_id-1);          fprintf(stderr,"\t.\tydot[%ld]\n",-y_id-1);
657      }      }
658    
659        ASC_ASSERT_LT(*varindx,1e7L);
660        ASC_ASSERT_LT(y_id, 9999999L);
661        ASC_ASSERT_LT(-9999999L, y_id);
662  }  }
663    
664  void integrator_visit_system_vars(IntegratorSystem *sys,IntegratorVarVisitorFn *visitfn){  void integrator_visit_system_vars(IntegratorSystem *sys,IntegratorVarVisitorFn *visitfn){
# Line 1121  int integrator_solve(IntegratorSystem *s Line 1134  int integrator_solve(IntegratorSystem *s
1134      unsigned long start_index=0, finish_index=0;      unsigned long start_index=0, finish_index=0;
1135      assert(sys!=NULL);      assert(sys!=NULL);
1136    
1137        assert(sys->internals);
1138        assert(sys->engine!=INTEG_UNKNOWN);
1139    
1140      nstep = integrator_getnsamples(sys)-1;      nstep = integrator_getnsamples(sys)-1;
1141      /* check for at least 2 steps and dimensionality of x vs steps here */      /* check for at least 2 steps and dimensionality of x vs steps here */
1142    
# Line 1154  int integrator_solve(IntegratorSystem *s Line 1170  int integrator_solve(IntegratorSystem *s
1170    
1171      CONSOLE_DEBUG("RUNNING INTEGRATION...");      CONSOLE_DEBUG("RUNNING INTEGRATION...");
1172    
1173      /* now go and run the integrator */      return (sys->internals->solvefn)(sys,start_index,finish_index);
     switch (sys->engine) {  
         case INTEG_LSODE:  
             return integrator_lsode_solve(sys, start_index, finish_index); break;  
 #ifdef ASC_WITH_IDA  
         case INTEG_IDA:  
             return integrator_ida_solve(sys,start_index, finish_index); break;  
 #endif  
         case INTEG_AWW:  
             return integrator_aww_solve(sys,start_index, finish_index); break;  
         default:  
             ERROR_REPORTER_HERE(ASC_PROG_ERR,"Unknown integrator (invalid, or not implemented yet)");  
             return 0;  
     }  
1174  }  }
1175    
1176  /*---------------------------------------------------------------  /*---------------------------------------------------------------
Legend:
Removed from v.976  
changed lines
  Added in v.977
john.pye@anu.edu.au
 ViewVC Help
Powered by ViewVC 1.1.22