/[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 979 by johnpye, Wed Dec 20 14:34:16 2006 UTC revision 980 by johnpye, Wed Dec 20 23:23:50 2006 UTC
# Line 24  Line 24 
24      by John Pye, May 2006      by John Pye, May 2006
25  */  */
26    
27  /*  /*
28      Be careful with the following. This file requires both the 'ida.h' from      Be careful with the following. This file requires both the 'ida.h' from
29      SUNDIALS as well as the 'ida.h' from ASCEND. Make sure that we're getting      SUNDIALS as well as the 'ida.h' from ASCEND. Make sure that we're getting
30      both of these; if you get problems check your build tool for the paths being      both of these; if you get problems check your build tool for the paths being
# Line 101  const IntegratorInternals integrator_ida Line 101  const IntegratorInternals integrator_ida
101  };  };
102    
103  /**  /**
104      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
105      common IntegratorSystem struct.      common IntegratorSystem struct.
106  */  */
107  typedef struct{  typedef struct{
# Line 331  int integrator_ida_solve( Line 331  int integrator_ida_solve(
331      }      }
332    
333      /* retrieve initial values from the system */      /* retrieve initial values from the system */
334        
335      /** @TODO fix this, the starting time != first sample */      /** @TODO fix this, the starting time != first sample */
336      t0 = integrator_get_t(blsys);      t0 = integrator_get_t(blsys);
337      CONSOLE_DEBUG("RETRIEVED t0 = %f",t0);      CONSOLE_DEBUG("RETRIEVED t0 = %f",t0);
# Line 352  int integrator_ida_solve( Line 352  int integrator_ida_solve(
352      /* create IDA object */      /* create IDA object */
353      ida_mem = IDACreate();      ida_mem = IDACreate();
354    
355      /* relative error tolerance */        /* relative error tolerance */
356      reltol = SLV_PARAM_REAL(&(blsys->params),IDA_PARAM_RTOL);      reltol = SLV_PARAM_REAL(&(blsys->params),IDA_PARAM_RTOL);
357      CONSOLE_DEBUG("rtol = %8.2e",reltol);      CONSOLE_DEBUG("rtol = %8.2e",reltol);
358    
# Line 387  int integrator_ida_solve( Line 387  int integrator_ida_solve(
387      /* set optional inputs... */      /* set optional inputs... */
388      IDASetErrHandlerFn(ida_mem, &integrator_ida_error, (void *)blsys);      IDASetErrHandlerFn(ida_mem, &integrator_ida_error, (void *)blsys);
389      IDASetRdata(ida_mem, (void *)blsys);      IDASetRdata(ida_mem, (void *)blsys);
390      IDASetMaxStep(ida_mem, integrator_get_maxstep(blsys));        IDASetMaxStep(ida_mem, integrator_get_maxstep(blsys));
391      IDASetInitStep(ida_mem, integrator_get_stepzero(blsys));      IDASetInitStep(ida_mem, integrator_get_stepzero(blsys));
392      IDASetMaxNumSteps(ida_mem, integrator_get_maxsubsteps(blsys));      IDASetMaxNumSteps(ida_mem, integrator_get_maxsubsteps(blsys));
393      if(integrator_get_minstep(blsys)>0){      if(integrator_get_minstep(blsys)>0){
# Line 409  int integrator_ida_solve( Line 409  int integrator_ida_solve(
409              case IDADENSE_MEM_FAIL: ERROR_REPORTER_HERE(ASC_PROG_ERR,"Memory allocation failed for IDADENSE"); return 0;              case IDADENSE_MEM_FAIL: ERROR_REPORTER_HERE(ASC_PROG_ERR,"Memory allocation failed for IDADENSE"); return 0;
410              default: ERROR_REPORTER_HERE(ASC_PROG_ERR,"bad return"); return 0;              default: ERROR_REPORTER_HERE(ASC_PROG_ERR,"bad return"); return 0;
411          }          }
412            
413          if(SLV_PARAM_BOOL(&(blsys->params),IDA_PARAM_AUTODIFF)){          if(SLV_PARAM_BOOL(&(blsys->params),IDA_PARAM_AUTODIFF)){
414              CONSOLE_DEBUG("USING AUTODIFF");              CONSOLE_DEBUG("USING AUTODIFF");
415              flag = IDADenseSetJacFn(ida_mem, &integrator_ida_djex, (void *)blsys);              flag = IDADenseSetJacFn(ida_mem, &integrator_ida_djex, (void *)blsys);
# Line 439  int integrator_ida_solve( Line 439  int integrator_ida_solve(
439          }else{          }else{
440              ERROR_REPORTER_HERE(ASC_PROG_ERR,"Unknown IDA linear solver choice '%s'",linsolver);              ERROR_REPORTER_HERE(ASC_PROG_ERR,"Unknown IDA linear solver choice '%s'",linsolver);
441              return 0;              return 0;
442          }                        }
443                
444          if(flag==IDASPILS_MEM_NULL){          if(flag==IDASPILS_MEM_NULL){
445              ERROR_REPORTER_HERE(ASC_PROG_ERR,"ida_mem is NULL");              ERROR_REPORTER_HERE(ASC_PROG_ERR,"ida_mem is NULL");
446              return 0;              return 0;
# Line 485  int integrator_ida_solve( Line 485  int integrator_ida_solve(
485      }      }
486    
487      /* set linear solver optional inputs...      /* set linear solver optional inputs...
488          ...nothing here at the moment...          ...nothing here at the moment...
489      */      */
490    
491      /* calculate initial conditions */      /* calculate initial conditions */
# Line 535  int integrator_ida_solve( Line 535  int integrator_ida_solve(
535          t = samplelist_get(blsys->samples, t_index);          t = samplelist_get(blsys->samples, t_index);
536          t0 = integrator_get_t(blsys);          t0 = integrator_get_t(blsys);
537          asc_assert(t > t0);          asc_assert(t > t0);
538            
539  #ifdef SOLVE_DEBUG  #ifdef SOLVE_DEBUG
540          CONSOLE_DEBUG("Integratoring from t0 = %f to t = %f", t0, t);          CONSOLE_DEBUG("Integratoring from t0 = %f to t = %f", t0, t);
541  #endif  #endif
542        
543          flag = IDASolve(ida_mem, t, &tret, yret, ypret, IDA_NORMAL);          flag = IDASolve(ida_mem, t, &tret, yret, ypret, IDA_NORMAL);
544    
545          /* pass the values of everything back to the compiler */          /* pass the values of everything back to the compiler */
# Line 568  int integrator_ida_solve( Line 568  int integrator_ida_solve(
568      }      }
569    
570      /* get optional outputs */      /* get optional outputs */
571        
572      /* free solution memory */      /* free solution memory */
573      N_VDestroy_Serial(yret);      N_VDestroy_Serial(yret);
574      N_VDestroy_Serial(ypret);      N_VDestroy_Serial(ypret);
575        
576      /* free solver memory */      /* free solver memory */
577      IDAFree(ida_mem);      IDAFree(ida_mem);
578    
# Line 631  int integrator_ida_fex(realtype tt, N_Ve Line 631  int integrator_ida_fex(realtype tt, N_Ve
631      */      */
632    
633      /* evaluate each residual in the rellist */      /* evaluate each residual in the rellist */
634      is_error = 0;      is_error = 0;
635      relptr = enginedata->rellist;      relptr = enginedata->rellist;
636    
637      if(enginedata->safeeval){      if(enginedata->safeeval){
# Line 697  int integrator_ida_fex(realtype tt, N_Ve Line 697  int integrator_ida_fex(realtype tt, N_Ve
697          return 1;          return 1;
698      }      }
699    
700  #ifdef FEX_DEBUG      #ifdef FEX_DEBUG
701      CONSOLE_DEBUG("RESIDUAL OK");      CONSOLE_DEBUG("RESIDUAL OK");
702  #endif  #endif
703      return 0;      return 0;
# Line 762  int integrator_ida_djex(long int Neq, re Line 762  int integrator_ida_djex(long int Neq, re
762      /* print step size */      /* print step size */
763      CONSOLE_DEBUG("<c_j> = %f",c_j);      CONSOLE_DEBUG("<c_j> = %f",c_j);
764  #endif  #endif
765        
766      /* build up the dense jacobian matrix... */      /* build up the dense jacobian matrix... */
767      status = 0;      status = 0;
768      for(i=0, relptr = enginedata->rellist;      for(i=0, relptr = enginedata->rellist;
# Line 796  int integrator_ida_djex(long int Neq, re Line 796  int integrator_ida_djex(long int Neq, re
796              ASC_FREE(varname);              ASC_FREE(varname);
797          }          }
798          fprintf(stderr,"\n");          fprintf(stderr,"\n");
799  #endif    #endif
800          /* insert values into the Jacobian row in appropriate spots (can assume Jac starts with zeros -- IDA manual) */          /* insert values into the Jacobian row in appropriate spots (can assume Jac starts with zeros -- IDA manual) */
801          for(j=0; j < count; ++j){                    for(j=0; j < count; ++j){
802              var_yindex = blsys->y_id[variables[j]];              var_yindex = blsys->y_id[variables[j]];
803    #ifndef __WIN32__
804                /* the SUNDIALS headers seem not to store 'N' on Windows */
805              ASC_ASSERT_RANGE(var_yindex, -Jac->N, Jac->N);              ASC_ASSERT_RANGE(var_yindex, -Jac->N, Jac->N);
806    #endif
807              if(var_yindex >= 0){              if(var_yindex >= 0){
808                  asc_assert(blsys->y[var_yindex]==enginedata->varlist[variables[j]]);                  asc_assert(blsys->y[var_yindex]==enginedata->varlist[variables[j]]);
809                  DENSE_ELEM(Jac,i,var_yindex) += derivatives[j];                  DENSE_ELEM(Jac,i,var_yindex) += derivatives[j];
# Line 831  int integrator_ida_djex(long int Neq, re Line 834  int integrator_ida_djex(long int Neq, re
834              fprintf(stderr,"%11.2e",DENSE_ELEM(Jac,i,j));              fprintf(stderr,"%11.2e",DENSE_ELEM(Jac,i,j));
835          }          }
836          fprintf(stderr,"\n");          fprintf(stderr,"\n");
837      }        }
838  #endif  #endif
839    
840      if(status){      if(status){
# Line 893  int integrator_ida_jvex(realtype tt, N_V Line 896  int integrator_ida_jvex(realtype tt, N_V
896      /* no real use for residuals (rr) here, I don't think? */      /* no real use for residuals (rr) here, I don't think? */
897    
898      /* allocate space for returns from relman_diff2: we *should* be able to use 'tmp1' and 'tmp2' here... */      /* allocate space for returns from relman_diff2: we *should* be able to use 'tmp1' and 'tmp2' here... */
899        
900      i = NV_LENGTH_S(yy) * 2;      i = NV_LENGTH_S(yy) * 2;
901  #ifdef JEX_DEBUG  #ifdef JEX_DEBUG
902      CONSOLE_DEBUG("Allocating 'variables' with length %d",i);      CONSOLE_DEBUG("Allocating 'variables' with length %d",i);
# Line 975  int integrator_ida_jvex(realtype tt, N_V Line 978  int integrator_ida_jvex(realtype tt, N_V
978                      );                      );
979  #endif  #endif
980                      asc_assert(blsys->ydot[-var_yindex-1]==enginedata->varlist[variables[j]]);                      asc_assert(blsys->ydot[-var_yindex-1]==enginedata->varlist[variables[j]]);
981                      Jv_i += derivatives[j] * NV_Ith_S(v,-var_yindex-1) * c_j;                      Jv_i += derivatives[j] * NV_Ith_S(v,-var_yindex-1) * c_j;
982                  }                  }
983              }              }
984    
# Line 985  int integrator_ida_jvex(realtype tt, N_V Line 988  int integrator_ida_jvex(realtype tt, N_V
988              relname = rel_make_name(blsys->system, *relptr);              relname = rel_make_name(blsys->system, *relptr);
989              CONSOLE_DEBUG("'%s': Jv[%d] = %f", relname, i, NV_Ith_S(Jv,i));              CONSOLE_DEBUG("'%s': Jv[%d] = %f", relname, i, NV_Ith_S(Jv,i));
990              //ASC_FREE(relname);              //ASC_FREE(relname);
991              return 1;                        return 1;
992  #endif  #endif
993          }          }
994      }else{      }else{

Legend:
Removed from v.979  
changed lines
  Added in v.980

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