/[ascend]/trunk/solvers/dopri5/asc_dopri5.c
ViewVC logotype

Diff of /trunk/solvers/dopri5/asc_dopri5.c

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

revision 1756 by jpye, Wed Feb 6 04:12:44 2008 UTC revision 1757 by jpye, Wed Feb 27 05:40:40 2008 UTC
# Line 1  Line 1 
1  /*  ASCEND modelling environment  /*  ASCEND modelling environment
2      Copyright (C) 2007 Carnegie Mellon University      Copyright (C) 2007 Carnegie Mellon University
3    
4      This program is free software; you can redistribute it and/or modify      This program is free software; you can redistribute it and/or modify
5      it under the terms of the GNU General Public License as published by      it under the terms of the GNU General Public License as published by
6      the Free Software Foundation; either version 2, or (at your option)      the Free Software Foundation; either version 2, or (at your option)
7      any later version.      any later version.
8    
9      This program is distributed in the hope that it will be useful,      This program is distributed in the hope that it will be useful,
10      but WITHOUT ANY WARRANTY; without even the implied warranty of      but WITHOUT ANY WARRANTY; without even the implied warranty of
11      MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the      MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
12      GNU General Public License for more details.      GNU General Public License for more details.
13    
14      You should have received a copy of the GNU General Public License      You should have received a copy of the GNU General Public License
15      along with this program; if not, write to the Free Software      along with this program; if not, write to the Free Software
16      Foundation, Inc., 59 Temple Place - Suite 330,      Foundation, Inc., 59 Temple Place - Suite 330,
17      Boston, MA 02111-1307, USA.      Boston, MA 02111-1307, USA.
18  *//** @file  *//** @file
19      DOPRI5 Runge-Kutta integrator      DOPRI5 Runge-Kutta integrator
20    
21      Based on the implementation of LSODE integrator, but adapted to      Based on the implementation of LSODE integrator, but adapted to
22      an explicit one-step method.      an explicit one-step method.
23  *//*  *//*
24      by John Pye, May 2007.      by John Pye, May 2007.
# Line 291  static double *dopri5_get_artol(Integrat Line 291  static double *dopri5_get_artol(Integrat
291    
292      // single tol for all vars      // single tol for all vars
293      tolval = SLV_PARAM_REAL(&(blsys->params),DOPRI5_PARAM_RTOL);      tolval = SLV_PARAM_REAL(&(blsys->params),DOPRI5_PARAM_RTOL);
294      CONSOLE_DEBUG("Using RTOL = %f for all vars", tolval);      CONSOLE_DEBUG("Using RTOL = %g for all vars", tolval);
295    
296      tol = ASC_NEW(double);      tol = ASC_NEW(double);
297      if(tol == NULL){      if(tol == NULL){
# Line 336  static double *dopri5_get_artol(Integrat Line 336  static double *dopri5_get_artol(Integrat
336    STATS    STATS
337  */  */
338    
339  /*  /*
340      Several functions provide access to different values :      Several functions provide access to different values :
341    
342      xRead   x value for which the solution has been computed (x=xend after      xRead   x value for which the solution has been computed (x=xend after
343          successful return).          successful return).
344    
345      hRead   Predicted step size of the last accepted step (useful for a      hRead   Predicted step size of the last accepted step (useful for a
346          subsequent call to dopri5).          subsequent call to dopri5).
347    
348      nstepRead   Number of used steps.      nstepRead   Number of used steps.
349      naccptRead  Number of accepted steps.      naccptRead  Number of accepted steps.
350      nrejctRead  Number of rejected steps.      nrejctRead  Number of rejected steps.
# Line 408  static FcnEqDiff integrator_dopri5_fex; Line 408  static FcnEqDiff integrator_dopri5_fex;
408      @param t indep var value      @param t indep var value
409      @param y input vector of variable values      @param y input vector of variable values
410      @param ydot return vector of calculated derivatives      @param ydot return vector of calculated derivatives
411      @param user_data point to whatever we want, in this case the IntegratorSystem        @param user_data point to whatever we want, in this case the IntegratorSystem
412  */  */
413  static void integrator_dopri5_fex(  static void integrator_dopri5_fex(
414          unsigned n_eq, double t, double *y, double *ydot          unsigned n_eq, double t, double *y, double *ydot
# Line 454  static void integrator_dopri5_fex( Line 454  static void integrator_dopri5_fex(
454          ERROR_REPORTER_START_HERE(ASC_PROG_ERR);          ERROR_REPORTER_START_HERE(ASC_PROG_ERR);
455          FPRINTF(ASCERR,"Unable to compute the vector of derivatives with the following values for the state variables:\n");          FPRINTF(ASCERR,"Unable to compute the vector of derivatives with the following values for the state variables:\n");
456          for (i = 0; i< n_eq; i++) {          for (i = 0; i< n_eq; i++) {
457              FPRINTF(ASCERR,"y[%4d] = %f\n",i, y[i]);              FPRINTF(ASCERR,"y[%4d] = %g\n",i, y[i]);
458          }          }
459          error_reporter_end_flush();          error_reporter_end_flush();
460  #endif  #endif
# Line 536  static void integrator_dopri5_reporter( Line 536  static void integrator_dopri5_reporter(
536          //CONSOLE_DEBUG("t=%f > ts=%f (currentsample = %ld",t,ts,d->currentsample);          //CONSOLE_DEBUG("t=%f > ts=%f (currentsample = %ld",t,ts,d->currentsample);
537          integrator_output_write_obs(blsys);          integrator_output_write_obs(blsys);
538  #ifdef DOPRI5_DEBUG  #ifdef DOPRI5_DEBUG
539          CONSOLE_DEBUG("step = %ld", nr-1);            CONSOLE_DEBUG("step = %ld", nr-1);
540  #endif  #endif
541          while(t>ts){          while(t>ts){
542              d->currentsample++;              d->currentsample++;
# Line 596  int integrator_dopri5_solve(IntegratorSy Line 596  int integrator_dopri5_solve(IntegratorSy
596    
597      d->y_vars = ASC_NEW_ARRAY(struct var_variable *,d->n_eqns+1);      d->y_vars = ASC_NEW_ARRAY(struct var_variable *,d->n_eqns+1);
598      d->ydot_vars = ASC_NEW_ARRAY(struct var_variable *, d->n_eqns+1);      d->ydot_vars = ASC_NEW_ARRAY(struct var_variable *, d->n_eqns+1);
599        
600      d->yinter = ASC_NEW_ARRAY(double,d->n_eqns);      d->yinter = ASC_NEW_ARRAY(double,d->n_eqns);
601    
602      /* set up the NLA solver here */      /* set up the NLA solver here */
603    
604      /*      /*
605          DOPRI5 should be OK to deal with any linsol/linsolqr-based solver.          DOPRI5 should be OK to deal with any linsol/linsolqr-based solver.
606          But for the moment we restrict to just QRSlv :-(          But for the moment we restrict to just QRSlv :-(
607      */      */
# Line 615  int integrator_dopri5_solve(IntegratorSy Line 615  int integrator_dopri5_solve(IntegratorSy
615      slv_get_status(blsys->system, &status);      slv_get_status(blsys->system, &status);
616    
617      if(status.struct_singular){      if(status.struct_singular){
618          ERROR_REPORTER_HERE(ASC_USER_WARNING              ERROR_REPORTER_HERE(ASC_USER_WARNING
619              ,"The system (according to QRSlv) is structurally singular."              ,"The system (according to QRSlv) is structurally singular."
620              " The ODE system may also be singular, but not necessarily."              " The ODE system may also be singular, but not necessarily."
621          );          );
# Line 653  int integrator_dopri5_solve(IntegratorSy Line 653  int integrator_dopri5_solve(IntegratorSy
653      d->currentsample = 1;      d->currentsample = 1;
654    
655      y = integrator_get_y(blsys, NULL);      y = integrator_get_y(blsys, NULL);
656        
657      rtoler = dopri5_get_artol(blsys,1,tolvect);      rtoler = dopri5_get_artol(blsys,1,tolvect);
658      atoler = dopri5_get_artol(blsys,0,tolvect);      atoler = dopri5_get_artol(blsys,0,tolvect);
659    
Legend:
Removed from v.1756  
changed lines
  Added in v.1757
john.pye@anu.edu.au
 ViewVC Help
Powered by ViewVC 1.1.22