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

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

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

revision 950 by johnpye, Sun Nov 26 03:14:02 2006 UTC revision 951 by johnpye, Sun Nov 26 05:01:49 2006 UTC
# Line 85  Line 85 
85  # error "Failed to include ASCEND IDA header file"  # error "Failed to include ASCEND IDA header file"
86  #endif  #endif
87    
88  /* #define JEX_DEBUG */  #define FEX_DEBUG
89  /* #define FEX_DEBUG */  #define JEX_DEBUG
90    
91  /**  /**
92      Struct containing any stuff that IDA needs that doesn't fit into the      Struct containing any stuff that IDA needs that doesn't fit into the
93      common IntegratorSystem struct.      common IntegratorSystem struct.
# Line 162  enum ida_parameters{ Line 163  enum ida_parameters{
163      ,IDA_PARAM_RTOL      ,IDA_PARAM_RTOL
164      ,IDA_PARAM_ATOL      ,IDA_PARAM_ATOL
165      ,IDA_PARAM_ATOLVECT      ,IDA_PARAM_ATOLVECT
166      ,IDA_PARAMS_SIZE      ,IDA_PARAM_GSMODIFIED
167            ,IDA_PARAMS_SIZE
168  };  };
169    
170  /**  /**
# Line 230  int integrator_ida_params_default(Integr Line 232  int integrator_ida_params_default(Integr
232          }, "SPGMR"}, (char *[]){"DENSE","SPGMR",NULL}          }, "SPGMR"}, (char *[]){"DENSE","SPGMR",NULL}
233      );      );
234    
235        slv_param_bool(p,IDA_PARAM_GSMODIFIED
236                ,(SlvParameterInitBool){{"gsmodified"
237                ,"Use modified Gram-Schmidt orthogonalisation in SPGMR?",2
238                ,"TRUE = GS_MODIFIED, FALSE = GS_CLASSICAL. See IDA manual section 5.5.6.6"
239            }, TRUE}
240        );
241    
242      asc_assert(p->num_parms == IDA_PARAMS_SIZE);      asc_assert(p->num_parms == IDA_PARAMS_SIZE);
243    
244      CONSOLE_DEBUG("Created %d params", p->num_parms);      CONSOLE_DEBUG("Created %d params", p->num_parms);
# Line 333  int integrator_ida_solve( Line 342  int integrator_ida_solve(
342      IDASetInitStep(ida_mem, integrator_get_stepzero(blsys));      IDASetInitStep(ida_mem, integrator_get_stepzero(blsys));
343      IDASetMaxNumSteps(ida_mem, integrator_get_maxsubsteps(blsys));      IDASetMaxNumSteps(ida_mem, integrator_get_maxsubsteps(blsys));
344      if(integrator_get_minstep(blsys)>0){      if(integrator_get_minstep(blsys)>0){
345          ERROR_REPORTER_HERE(ASC_PROG_NOTE,"IDA does not support minstep (ignored)");          ERROR_REPORTER_HERE(ASC_PROG_NOTE,"IDA does not support minstep (ignored)\n");
346      }      }
347      /* there's no capability for setting *minimum* step size in IDA */      /* there's no capability for setting *minimum* step size in IDA */
348    
# Line 366  int integrator_ida_solve( Line 375  int integrator_ida_solve(
375          }else{          }else{
376              CONSOLE_DEBUG("USING NUMERICAL DIFF");              CONSOLE_DEBUG("USING NUMERICAL DIFF");
377          }          }
378    
379            /* select Gram-Schmidt orthogonalisation */
380            if(SLV_PARAM_BOOL(&(blsys->params),IDA_PARAM_GSMODIFIED)){
381                CONSOLE_DEBUG("USING MODIFIED GS");
382                flag = IDASpilsSetGSType(ida_mem,MODIFIED_GS);
383                if(flag!=IDASPILS_SUCCESS){
384                    ERROR_REPORTER_HERE(ASC_PROG_ERR,"Failed to set GS_MODIFIED");
385                    return 0;
386                }
387            }else{
388                CONSOLE_DEBUG("USING CLASSICAL GS");
389                flag = IDASpilsSetGSType(ida_mem,CLASSICAL_GS);
390                if(flag!=IDASPILS_SUCCESS){
391                    ERROR_REPORTER_HERE(ASC_PROG_ERR,"Failed to set GS_MODIFIED");
392                    return 0;
393                }
394            }
395      }else if(strcmp(linsolver,"DENSE")==0){      }else if(strcmp(linsolver,"DENSE")==0){
396          CONSOLE_DEBUG("DENSE DIRECT SOLVER, size = %d",size);          CONSOLE_DEBUG("DENSE DIRECT SOLVER, size = %d",size);
397          flag = IDADense(ida_mem, size);          flag = IDADense(ida_mem, size);
# Line 546  int integrator_ida_fex(realtype tt, N_Ve Line 572  int integrator_ida_fex(realtype tt, N_Ve
572              NV_Ith_S(rr,i) = resid;              NV_Ith_S(rr,i) = resid;
573              if(!calc_ok){              if(!calc_ok){
574                  relname = rel_make_name(blsys->system, *relptr);                  relname = rel_make_name(blsys->system, *relptr);
575                  ERROR_REPORTER_HERE(ASC_PROG_ERR,"Residual calculation error in rel '%s'",relname);                  ERROR_REPORTER_HERE(ASC_PROG_ERR,"Calculation error in rel '%s'",relname);
576                  ASC_FREE(relname);                  ASC_FREE(relname);
577                  /* presumable some output already made? */                  /* presumable some output already made? */
578                  is_error = 1;                  is_error = 1;
# Line 764  int integrator_ida_jvex(realtype tt, N_V Line 790  int integrator_ida_jvex(realtype tt, N_V
790      var_filter_t filter;      var_filter_t filter;
791      int count;      int count;
792    
793      /* fprintf(stderr,"\n--------------\n"); */  #ifdef JEX_DEBUG
794      /* CONSOLE_DEBUG("EVALUTING JACOBIAN..."); */      CONSOLE_DEBUG("EVALUATING JACOBIAN...");
795    #endif
796    
797      blsys = (IntegratorSystem *)jac_data;      blsys = (IntegratorSystem *)jac_data;
798      enginedata = integrator_ida_enginedata(blsys);      enginedata = integrator_ida_enginedata(blsys);
# Line 786  int integrator_ida_jvex(realtype tt, N_V Line 813  int integrator_ida_jvex(realtype tt, N_V
813      filter.matchbits = VAR_SVAR;      filter.matchbits = VAR_SVAR;
814      filter.matchvalue = VAR_SVAR;      filter.matchvalue = VAR_SVAR;
815    
     /* CONSOLE_DEBUG("PRINTING VALUES OF 'v' VECTOR (length %ld)",NV_LENGTH_S(v)); */  
     /* for(i=0; i<NV_LENGTH_S(v); ++i){  
         CONSOLE_DEBUG("v[%d] = %f",i,NV_Ith_S(v,i));  
     }*/  
   
816      Asc_SignalHandlerPush(SIGFPE,SIG_IGN);      Asc_SignalHandlerPush(SIGFPE,SIG_IGN);
817      if (setjmp(g_fpe_env)==0) {      if (setjmp(g_fpe_env)==0) {
818          for(i=0, relptr = enginedata->rellist;          for(i=0, relptr = enginedata->rellist;
819                  i< enginedata->nrels && relptr != NULL;                  i< enginedata->nrels && relptr != NULL;
820                  ++i, ++relptr                  ++i, ++relptr
821          ){          ){
             /* fprintf(stderr,"\n"); */  
             /* relname = rel_make_name(blsys->system, *relptr);  
             CONSOLE_DEBUG("RELATION %d '%s'",i,relname);  
             ASC_FREE(relname); */  
   
822              /* get derivatives for this particular relation */              /* get derivatives for this particular relation */
823              status = relman_diff2(*relptr, &filter, derivatives, variables, &count, enginedata->safeeval);              status = relman_diff2(*relptr, &filter, derivatives, variables, &count, enginedata->safeeval);
824              /* CONSOLE_DEBUG("Got derivatives against %d matching variables", count); */              /* CONSOLE_DEBUG("Got derivatives against %d matching variables", count); */
825    
826              for(j=0;j<count;++j){              if(status){
827                  /* varname = var_make_name(blsys->system, enginedata->varlist[variables[j]]);                  relname = rel_make_name(blsys->system, *relptr);
828                  CONSOLE_DEBUG("derivatives[%d] = %f (variable %d, '%s')",j,derivatives[j],variables[j],varname);                  ERROR_REPORTER_HERE(ASC_PROG_ERR,"Calculation error in rel '%s'",relname);
829                  ASC_FREE(varname); */                  ASC_FREE(relname);
830              }                  is_error = 1;
   
             if(!status){  
                 /* CONSOLE_DEBUG("Derivatives for relation %d OK",i); */  
             }else{  
                 CONSOLE_DEBUG("ERROR calculating derivatives for relation %d",i);  
831                  break;                  break;
832              }              }
833    
# Line 837  int integrator_ida_jvex(realtype tt, N_V Line 849  int integrator_ida_jvex(realtype tt, N_V
849                  }                  }
850                  */                  */
851                                    
852                  var_yindex = blsys->y_id[variables[j] - 1];                  var_yindex = blsys->y_id[variables[j]];
853                  /* CONSOLE_DEBUG("j = %d: variables[j] = %d, y_id = %d",j,variables[j],var_yindex); */                  /* CONSOLE_DEBUG("j = %d: variables[j] = %d, y_id = %d",j,variables[j],var_yindex); */
854    
855                  if(var_yindex >= 0){                  if(var_yindex >= 0){
856                      /* CONSOLE_DEBUG("j = %d: algebraic, deriv[j] = %f, v[%d] = %f",j,derivatives[j], var_yindex, NV_Ith_S(v,var_yindex)); */  #ifdef JEX_DEBUG
857                        asc_assert(blsys->y[var_yindex]==enginedata->varlist[variables[j]]);
858                        fprintf(stderr,"Jv[%d] += %f (dF[%d]/dy[%d] = %f, v[%d] = %f)\n", i
859                            , derivatives[j] * NV_Ith_S(v,var_yindex)
860                            , i, var_yindex, derivatives[j]
861                            , var_yindex, NV_Ith_S(v,var_yindex)
862                        );
863    #endif
864                      Jv_i += derivatives[j] * NV_Ith_S(v,var_yindex);                      Jv_i += derivatives[j] * NV_Ith_S(v,var_yindex);
865                  }else{                  }else{
866                      var_yindex = -var_yindex-1;  #ifdef JEX_DEBUG
867                      /* CONSOLE_DEBUG("j = %d: differential, deriv[j] = %f, v[%d] = %f",j,derivatives[j], var_yindex, NV_Ith_S(v,var_yindex)); */                      fprintf(stderr,"Jv[%d] += %f (dF[%d]/dydot[%d] = %f, v[%d] = %f)\n", i
868                      Jv_i += derivatives[j] * NV_Ith_S(v,var_yindex) / c_j;                          , derivatives[j] * NV_Ith_S(v,-var_yindex-1)
869                            , i, -var_yindex-1, derivatives[j]
870                            , -var_yindex-1, NV_Ith_S(v,-var_yindex-1)
871                        );
872    #endif
873                        asc_assert(blsys->ydot[-var_yindex-1]==enginedata->varlist[variables[j]]);
874                        Jv_i += derivatives[j] * NV_Ith_S(v,-var_yindex-1) * c_j;
875                  }                  }
876              }              }
877    
878              NV_Ith_S(Jv,i) = Jv_i;              NV_Ith_S(Jv,i) = Jv_i;
879              /* CONSOLE_DEBUG("(J*v)[%d] = %f", i, Jv_i); */  #ifdef JEX_DEBUG
880          }              relname = rel_make_name(blsys->system, *relptr);
881          if(status){              CONSOLE_DEBUG("'%s': Jv[%d] = %f", relname, i, NV_Ith_S(Jv,i));
882              CONSOLE_DEBUG("JVEX ERROR RESULT");              ASC_FREE(relname);
883              /* presumably some error_reporter will already have been made*/              return 1;          
884              is_error = 1;  #endif
885          }          }
886      }else{      }else{
887          CONSOLE_DEBUG("FLOATING POINT ERROR WITH i=%d",i);          relname = rel_make_name(blsys->system, *relptr);
888            ERROR_REPORTER_HERE(ASC_PROG_ERR,"Floating point error (SIGFPE) in rel '%s'",relname);
889            ASC_FREE(relname);
890            is_error = 1;
891      }      }
892      Asc_SignalHandlerPop(SIGFPE,SIG_IGN);      Asc_SignalHandlerPop(SIGFPE,SIG_IGN);
893    
894      if(is_error)CONSOLE_DEBUG("SOME ERRORS FOUND IN EVALUATION");      if(is_error){
895                CONSOLE_DEBUG("SOME ERRORS FOUND IN EVALUATION");
896      return is_error;          return 1;
897        }
898        return 0;
899  }  }
900    
901  /*----------------------------------------------  /*----------------------------------------------
Legend:
Removed from v.950  
changed lines
  Added in v.951
john.pye@anu.edu.au
 ViewVC Help
Powered by ViewVC 1.1.22