generic/solver/ida.c

/*  ASCEND modelling environment
    Copyright (C) 2006 Carnegie Mellon University

    This program is free software; you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation; either version 2, or (at your option)
    any later version.

    This program is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details.

    You should have received a copy of the GNU General Public License
    along with this program; if not, write to the Free Software
    Foundation, Inc., 59 Temple Place - Suite 330,
    Boston, MA 02111-1307, USA.
*//**
    @file
    Access to the IDA integrator for ASCEND. IDA is a DAE solver that comes
    as part of the GPL-licensed SUNDIALS solver package from LLNL.
    @see http://www.llnl.gov/casc/sundials/
*//*
    by John Pye, May 2006
*/

/* 
    Be careful with the following. This file requires both the 'ida.h' from
    SUNDIALS as well as the 'ida.h' from ASCEND. Make sure that we're getting
    both of these; if you get problems check your build tool for the paths being
    passed to the C preprocessor.
*/

/* standard includes */
#include <signal.h>

/* ASCEND includes */
#include "ida.h"
#include <utilities/error.h>
#include <utilities/ascConfig.h>
#include <utilities/ascSignal.h>
#include <compiler/instance_enum.h>
#include "var.h"
#include "rel.h"
#include "discrete.h"
#include "conditional.h"
#include "logrel.h"
#include "bnd.h"
#include "linsol.h"
#include "linsolqr.h"
#include "slv_common.h"
#include "slv_client.h"
#include "relman.h"

/* SUNDIALS includes */
#ifdef ASC_WITH_IDA
# include <sundials/sundials_config.h>
# include <ida/ida.h>
# include <nvector/nvector_serial.h>
# include <ida/ida_spgmr.h>
# ifndef IDA_SUCCESS
#  error "Failed to include SUNDIALS IDA header file"
# endif
#endif

/*
    for the benefit of build tools that didn't sniff the SUNDIALS version, we
    assume version 2.2.x (and thence possible errors)
*/
#ifndef SUNDIALS_VERSION_MINOR
# warning "GUESSING SUNDIALS VERSION 2.2"
# define SUNDIALS_VERSION_MINOR 2
#endif
#ifndef SUNDIALS_VERSION_MAJOR
# define SUNDIALS_VERSION_MAJOR 2
#endif

/* check that we've got what we expect now */
#ifndef ASC_IDA_H
# error "Failed to include ASCEND IDA header file"
#endif

/**
    Struct containing any stuff that IDA needs that doesn't fit into the 
    common IntegratorSystem struct.
*/
typedef struct{
    struct rel_relation **rellist;   /**< NULL terminated list of rels */
    struct var_variable **varlist;   /**< NULL terminated list of vars. ONLY USED FOR DEBUGGING -- get rid of it! */
    int nrels;
    int safeeval;                    /**< whether to pass the 'safe' flag to relman_eval */
} IntegratorIdaData;

/*-------------------------------------------------------------
  FORWARD DECLS
*/
/* residual function forward declaration */
int integrator_ida_fex(realtype tt, N_Vector yy, N_Vector yp, N_Vector rr, void *res_data);

int integrator_ida_jvex(realtype tt, N_Vector yy, N_Vector yp, N_Vector rr
        , N_Vector v, N_Vector Jv, realtype c_j
        , void *jac_data, N_Vector tmp1, N_Vector tmp2
);

/* error handler forward declaration */
void integrator_ida_error(int error_code
        , const char *module, const char *function
        , char *msg, void *eh_data
);

/*-------------------------------------------------------------
  SETUP/TEARDOWN ROUTINES
*/
void integrator_ida_create(IntegratorSystem *blsys){
    CONSOLE_DEBUG("ALLOCATING IDA ENGINE DATA");
    IntegratorIdaData *enginedata;
    enginedata = ASC_NEW(IntegratorIdaData);
    enginedata->rellist = NULL;
    enginedata->varlist = NULL;
    enginedata->safeeval = 1;
    blsys->enginedata = (void *)enginedata;
}

void integrator_ida_free(void *enginedata){
    CONSOLE_DEBUG("DELETING IDA ENGINE DATA");
    IntegratorIdaData *d = (IntegratorIdaData *)enginedata;
    /* note, we don't own the rellist, so don't need to free it */
    ASC_FREE(d);
}

IntegratorIdaData *integrator_ida_enginedata(IntegratorSystem *blsys){
    IntegratorIdaData *d;
    assert(blsys!=NULL);
    assert(blsys->enginedata!=NULL);
    assert(blsys->engine==INTEG_IDA);
    d = ((IntegratorIdaData *)(blsys->enginedata));
    assert(d->safeeval = 1);
    return d;
}

/*-------------------------------------------------------------
  MAIN IDA SOLVER ROUTINE, see IDA manual, sec 5.4, p. 27 ff.
*/

/* return 1 on success */
int integrator_ida_solve(
        IntegratorSystem *blsys
        , unsigned long start_index
        , unsigned long finish_index
){
    void *ida_mem;
    int size, flag, t_index;
    realtype t0, reltol, t, tret, tout1;
    N_Vector y0, yp0, abstol, ypret, yret;
    IntegratorIdaData *enginedata;

    CONSOLE_DEBUG("STARTING IDA...");

    enginedata = integrator_ida_enginedata(blsys);
    CONSOLE_DEBUG("safeeval = %d",enginedata->safeeval);

    /* store reference to list of relations (in enginedata) */
    enginedata->nrels = slv_get_num_solvers_rels(blsys->system);
    enginedata->rellist = slv_get_solvers_rel_list(blsys->system);
    enginedata->varlist = slv_get_solvers_var_list(blsys->system);
    CONSOLE_DEBUG("Number of relations: %d",enginedata->nrels);
    CONSOLE_DEBUG("Number of dependent vars: %ld",blsys->n_y);
    size = blsys->n_y;

    if(enginedata->nrels!=size){
        ERROR_REPORTER_HERE(ASC_USER_ERROR,"Integration problem is not square (%d rels, %d vars)", enginedata->nrels, size);
        return 0; /* failure */
    }

    /* retrieve initial values from the system */
    t0 = samplelist_get(blsys->samples,start_index);

    y0 = N_VNew_Serial(size);
    integrator_get_y(blsys,NV_DATA_S(y0));

    yp0 = N_VNew_Serial(size);
    integrator_get_ydot(blsys,NV_DATA_S(yp0));

    N_VPrint_Serial(yp0);
    CONSOLE_DEBUG("yp0 is at %p",&yp0);

    /* create IDA object */
    ida_mem = IDACreate();

    /* retrieve the absolute tolerance values for each variable */
    abstol = N_VNew_Serial(size);
    N_VConst(0.1,abstol); /** @TODO fill in the abstol values from the model */
    reltol = 0.001;

    /* allocate internal memory */
    flag = IDAMalloc(ida_mem, &integrator_ida_fex, t0, y0, yp0, IDA_SV, reltol, abstol);
    if(flag==IDA_MEM_NULL){
        ERROR_REPORTER_HERE(ASC_PROG_ERR,"ida_mem is NULL");
        return 0;
    }else if(flag==IDA_MEM_FAIL){
        ERROR_REPORTER_HERE(ASC_PROG_ERR,"Unable to allocate memory (IDAMalloc)");
        return 0;
    }else if(flag==IDA_ILL_INPUT){
        ERROR_REPORTER_HERE(ASC_PROG_ERR,"Invalid input to IDAMalloc");
        return 0;
    }/* else success */

    /* set optional inputs... */
    IDASetErrHandlerFn(ida_mem, &integrator_ida_error, (void *)blsys);
    IDASetRdata(ida_mem, (void *)blsys);
    IDASetMaxStep(ida_mem, integrator_get_maxstep(blsys));
    IDASetInitStep(ida_mem, integrator_get_stepzero(blsys));
    IDASetMaxNumSteps(ida_mem, integrator_get_maxsubsteps(blsys));
    /* there's no capability for setting *minimum* step size in IDA */

    CONSOLE_DEBUG("ASSIGNING LINEAR SOLVER");

    /* attach linear solver module, using the default value of maxl */
    flag = IDASpgmr(ida_mem, 0);
    if(flag==IDASPILS_MEM_NULL){
        ERROR_REPORTER_HERE(ASC_PROG_ERR,"ida_mem is NULL");
        return 0;
    }else if(flag==IDASPILS_MEM_FAIL){
        ERROR_REPORTER_HERE(ASC_PROG_ERR,"Unable to allocate memory (IDASpgmr)");
        return 0;
    }/* else success */

    /* assign the J*v function */
    flag = IDASpilsSetJacTimesVecFn(ida_mem, &integrator_ida_jvex, (void *)blsys);
    if(flag==IDASPILS_MEM_NULL){
        ERROR_REPORTER_HERE(ASC_PROG_ERR,"ida_mem is NULL");
        return 0;
    }else if(flag==IDASPILS_LMEM_NULL){
        ERROR_REPORTER_HERE(ASC_PROG_ERR,"IDASPILS linear solver has not been initialized");
        return 0;
    }/* else success */

    /* set linear solver optional inputs...

        ...nothing here at the moment... 

    */

    /* correct initial values, given derivatives */
    blsys->currentstep=0;
    t_index=start_index+1;
    tout1 = samplelist_get(blsys->samples, t_index);

#if SUNDIALS_VERSION_MAJOR==2 && SUNDIALS_VERSION_MINOR==3
    flag = IDACalcIC(ida_mem, IDA_Y_INIT, tout1);
#else
    flag = IDACalcIC(ida_mem, t0, y0, yp0, IDA_Y_INIT, tout1);
#endif

    if(flag!=IDA_SUCCESS){
        ERROR_REPORTER_HERE(ASC_PROG_ERR,"Unable to solve initial values (IDACalcIC)");
        return 0;
    }/* else success */

    CONSOLE_DEBUG("INITIAL CONDITIONS SOLVED :-)");

    /* optionally, specify ROO-FINDING PROBLEM */

    /* -- set up the IntegratorReporter */
    integrator_output_init(blsys);

    /* -- store the initial values of all the stuff */
    integrator_output_write(blsys);
    integrator_output_write_obs(blsys);

    /* specify where the returned values should be stored */
    yret = y0;
    ypret = yp0;

    /* advance solution in time, return values as yret and derivatives as ypret */
    blsys->currentstep=1;
    for(t_index=start_index+1;t_index <= finish_index;++t_index, ++blsys->currentstep){
        t = samplelist_get(blsys->samples, t_index);

        CONSOLE_DEBUG("SOLVING UP TO t = %f", t);
    
        flag = IDASolve(ida_mem, t, &tret, yret, ypret, IDA_NORMAL);

        /* pass the values of everything back to the compiler */
        integrator_set_t(blsys, (double)tret);
        integrator_set_y(blsys, NV_DATA_S(yret));
        integrator_set_ydot(blsys, NV_DATA_S(ypret));

        if(flag<0){
            ERROR_REPORTER_HERE(ASC_PROG_ERR,"Failed to solve t = %f (IDASolve), error %d", t, flag);
            break;
        }

        /* -- do something so that blsys knows the values of tret, yret and ypret */

        /* -- store the current values of all the stuff */
        integrator_output_write(blsys);
        integrator_output_write_obs(blsys);
    }

    /* -- close the IntegratorReporter */
    integrator_output_close(blsys);

    if(flag < 0){
        ERROR_REPORTER_HERE(ASC_PROG_ERR,"Solving aborted while attempting t = %f", t);
        return 0;
    }

    /* get optional outputs */
    
    /* free solution memory */
    N_VDestroy_Serial(yret);
    N_VDestroy_Serial(ypret);
    
    /* free solver memory */
    IDAFree(ida_mem);

    /* all done */
    return 1;
}

/*--------------------------------------------------
  RESIDUALS AND JACOBIAN
*/
/**
    Function to evaluate system residuals, in the form required for IDA.

    Given tt, yy and yp, we need to evaluate and return rr.

    @param tt current value of indep variable (time)
    @param yy current values of dependent variable vector
    @param yp current values of derivatives of dependent variables
    @param rr the output residual vector (is we're returning data to)
    @param res_data pointer to our stuff (blsys in this case).

    @return 0 on success, positive on recoverable error, and
        negative on unrecoverable error.
*/
int integrator_ida_fex(realtype tt, N_Vector yy, N_Vector yp, N_Vector rr, void *res_data){
    IntegratorSystem *blsys;
    IntegratorIdaData *enginedata;
    int i, calc_ok, is_error;
    struct rel_relation** relptr;
    double resid;
    char *relname;

    blsys = (IntegratorSystem *)res_data;
    enginedata = integrator_ida_enginedata(blsys);

    /* fprintf(stderr,"\n\n"); */
    /* CONSOLE_DEBUG("ABOUT TO EVALUTE RESIDUALS..."); */

    /* pass the values of everything back to the compiler */
    integrator_set_t(blsys, (double)tt);
    integrator_set_y(blsys, NV_DATA_S(yy));
    integrator_set_ydot(blsys, NV_DATA_S(yp));

    /* revaluate the system residuals using the new data */
    is_error = 0; 
    relptr = enginedata->rellist;

    /* CONSOLE_DEBUG("IDA requests residuals of length %lu",NV_LENGTH_S(rr)); */
    if(NV_LENGTH_S(rr)!=enginedata->nrels){
        ERROR_REPORTER_HERE(ASC_PROG_ERR,"Invalid residuals nrels!=length(rr)");
        return -1; /* unrecoverable */
    }

    Asc_SignalHandlerPush(SIGFPE,SIG_IGN);
    if (setjmp(g_fpe_env)==0) {
        for(i=0, relptr = enginedata->rellist;
                i< enginedata->nrels && relptr != NULL;
                ++i, ++relptr
        ){
            resid = relman_eval(*relptr, &calc_ok, enginedata->safeeval);
            
            relname = rel_make_name(blsys->system, *relptr);
            /* if(calc_ok){
                CONSOLE_DEBUG("residual[%d:\"%s\"] = %f",i,relname,resid);
            }else{
                CONSOLE_DEBUG("residual[%d:\"%s\"] = %f (ERROR)",i,relname,resid);
            }*/
            ASC_FREE(relname);

            NV_Ith_S(rr,i) = resid;
            if(!calc_ok){
                /* presumable some output already made? */
                is_error = 1;
            }
        }
    }else{
        CONSOLE_DEBUG("FLOATING POINT ERROR WITH i=%d",i);
    }
    Asc_SignalHandlerPop(SIGFPE,SIG_IGN);

    if(is_error)CONSOLE_DEBUG("SOME ERRORS FOUND IN EVALUATION");
    return is_error;
}

/**
    Function to evaluate the product J*v, in the form required for IDA (see IDASpilsSetJacTimesVecFn)

    Given tt, yy, yp, rr and v, we need to evaluate and return Jv.

    @param tt current value of the independent variable (time, t)
    @param yy current value of the dependent variable vector, y(t).
    @param yp current value of y'(t).
    @param rr current value of the residual vector F(t, y, y').
    @param v  the vector by which the Jacobian must be multiplied to the right.
    @param Jv the output vector computed
    @param c_j the scalar in the system Jacobian, proportional to the inverse of the step size ($ \alpha$ in Eq. (3.5) ).
    @param jac_data pointer to our stuff (blsys in this case, passed into IDA via IDASp*SetJacTimesVecFn.)
    @param tmp1 @see tmp2
    @param tmp2 (as well as tmp1) pointers to memory allocated for variables of type N_Vector for use here as temporary storage or work space.
    @return 0 on success
*/
int integrator_ida_jvex(realtype tt, N_Vector yy, N_Vector yp, N_Vector rr
        , N_Vector v, N_Vector Jv, realtype c_j
        , void *jac_data, N_Vector tmp1, N_Vector tmp2
){
    IntegratorSystem *blsys;
    IntegratorIdaData *enginedata;
    int i, j, is_error=0;
    struct rel_relation** relptr;
    char *relname, *varname;
    int status;
    double Jv_i;

    int *variables;
    double *derivatives;
    var_filter_t filter;
    int count;

    CONSOLE_DEBUG("EVALUTING JACOBIAN...");

    blsys = (IntegratorSystem *)jac_data;
    enginedata = integrator_ida_enginedata(blsys);

    /* pass the values of everything back to the compiler */
    integrator_set_t(blsys, (double)tt);
    integrator_set_y(blsys, NV_DATA_S(yy));
    integrator_set_ydot(blsys, NV_DATA_S(yp));
    /* no real use for residuals (rr) here, I don't think? */

    /* allocate space for returns from relman_diff2: we *should* be able to use 'tmp1' and 'tmp2' here... */
    variables = ASC_NEW_ARRAY(int, NV_LENGTH_S(yy) * 2);
    derivatives = ASC_NEW_ARRAY(double, NV_LENGTH_S(yy) * 2);

    /* evaluate the derivatives... */
    /* J = dG_dy = dF_dy + alpha * dF_dyp */

    filter.matchbits = VAR_SVAR;
    filter.matchvalue = VAR_SVAR;

    CONSOLE_DEBUG("PRINTING VALUES OF 'v' VECTOR (length %ld)",NV_LENGTH_S(v));
    for(i=0; i<NV_LENGTH_S(v); ++i){
        varname = var_make_name(blsys->system, enginedata->varlist[i]);
        if(varname==NULL)varname="UNKNOWN";
        CONSOLE_DEBUG("v[%d] = '%s' = %f",i,varname,NV_Ith_S(v,i));
        ASC_FREE(varname);
    }

    Asc_SignalHandlerPush(SIGFPE,SIG_IGN);
    if (setjmp(g_fpe_env)==0) {
        for(i=0, relptr = enginedata->rellist;
                i< enginedata->nrels && relptr != NULL;
                ++i, ++relptr
        ){
            /* get derivatives for this particular relation */
            status = relman_diff2(*relptr, &filter, derivatives, variables, &count, enginedata->safeeval);

            CONSOLE_DEBUG("Got derivatives against %d matching variables", count);

            relname = rel_make_name(blsys->system, *relptr);
            if(!status){
                CONSOLE_DEBUG("Derivatives for relation %d '%s' OK",i,relname);
            }else{
                CONSOLE_DEBUG("ERROR calculating derivatives for relation %d '%s'",i,relname);
                break;
            }
            ASC_FREE(relname);

            Jv_i = 0;
            for(j=0; j < count; ++j){
                /* CONSOLE_DEBUG("j = %d, variables[j] = %d, n_y = %ld", j, variables[j], blsys->n_y); */
                varname = var_make_name(blsys->system, enginedata->varlist[variables[j]]);
                if(varname){
                    CONSOLE_DEBUG("Variable %d '%s' derivative = %f", variables[j],varname,derivatives[j]);
                    ASC_FREE(varname);
                }else{
                    CONSOLE_DEBUG("Variable %d (UNKNOWN!): derivative = %f",variables[j],derivatives[j]);
                }
                
                Jv_i += derivatives[j] * NV_Ith_S(v,variables[j]);
            }

            NV_Ith_S(Jv,i) = Jv_i;
            CONSOLE_DEBUG("(J*v)[%d] = %f", i, Jv_i);

            if(status){
                /* presumably some error_reporter will already have been made*/
                is_error = 1;
            }
        }
    }else{
        CONSOLE_DEBUG("FLOATING POINT ERROR WITH i=%d",i);
    }
    Asc_SignalHandlerPop(SIGFPE,SIG_IGN);

    if(is_error)CONSOLE_DEBUG("SOME ERRORS FOUND IN EVALUATION");

    
    
    return is_error;
}

/*----------------------------------------------
  ERROR REPORTING
*/
/**
    Error message reporter function to be passed to IDA. All error messages
    will trigger a call to this function, so we should find everything
    appearing on the console (in the case of Tcl/Tk) or in the errors/warnings
    panel (in the case of PyGTK).
*/
void integrator_ida_error(int error_code
        , const char *module, const char *function
        , char *msg, void *eh_data
){
    IntegratorSystem *blsys;
    error_severity_t sev;

    /* cast back the IntegratorSystem, just in case we need it */
    blsys = (IntegratorSystem *)eh_data;

    /* severity depends on the sign of the error_code value */
    if(error_code <= 0){
        sev = ASC_PROG_ERR;
    }else{
        sev = ASC_PROG_WARNING;
    }

    /* use our all-purpose error reporting to get stuff back to the GUI */
    error_reporter(sev,module,0,function,"%s (error %d)",msg,error_code);
}
1	johnpye	669	/* ASCEND modelling environment
2			Copyright (C) 2006 Carnegie Mellon University
3
4			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
6			the Free Software Foundation; either version 2, or (at your option)
7			any later version.
8
9			This program is distributed in the hope that it will be useful,
10			but WITHOUT ANY WARRANTY; without even the implied warranty of
11			MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12			GNU General Public License for more details.
13
14			You should have received a copy of the GNU General Public License
15			along with this program; if not, write to the Free Software
16			Foundation, Inc., 59 Temple Place - Suite 330,
17			Boston, MA 02111-1307, USA.
18			//*
19			@file
20			Access to the IDA integrator for ASCEND. IDA is a DAE solver that comes
21			as part of the GPL-licensed SUNDIALS solver package from LLNL.
22			@see http://www.llnl.gov/casc/sundials/
23			//
24			by John Pye, May 2006
25			*/
26
27			/*
28	johnpye	912	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
30			both of these; if you get problems check your build tool for the paths being
31			passed to the C preprocessor.
32	johnpye	669	*/
33
34			/* standard includes */
35			#include <signal.h>
36
37			/* ASCEND includes */
38			#include "ida.h"
39			#include <utilities/error.h>
40			#include <utilities/ascConfig.h>
41			#include <utilities/ascSignal.h>
42			#include <compiler/instance_enum.h>
43			#include "var.h"
44			#include "rel.h"
45			#include "discrete.h"
46			#include "conditional.h"
47			#include "logrel.h"
48			#include "bnd.h"
49			#include "linsol.h"
50			#include "linsolqr.h"
51			#include "slv_common.h"
52			#include "slv_client.h"
53	johnpye	719	#include "relman.h"
54	johnpye	669
55			/* SUNDIALS includes */
56			#ifdef ASC_WITH_IDA
57	johnpye	913	# include <sundials/sundials_config.h>
58	johnpye	782	# include <ida/ida.h>
59			# include <nvector/nvector_serial.h>
60			# include <ida/ida_spgmr.h>
61	johnpye	669	# ifndef IDA_SUCCESS
62			# error "Failed to include SUNDIALS IDA header file"
63			# endif
64			#endif
65
66	johnpye	913	/*
67			for the benefit of build tools that didn't sniff the SUNDIALS version, we
68			assume version 2.2.x (and thence possible errors)
69			*/
70			#ifndef SUNDIALS_VERSION_MINOR
71			# warning "GUESSING SUNDIALS VERSION 2.2"
72			# define SUNDIALS_VERSION_MINOR 2
73			#endif
74			#ifndef SUNDIALS_VERSION_MAJOR
75			# define SUNDIALS_VERSION_MAJOR 2
76			#endif
77
78	johnpye	669	/* check that we've got what we expect now */
79			#ifndef ASC_IDA_H
80			# error "Failed to include ASCEND IDA header file"
81			#endif
82
83			/**
84			Struct containing any stuff that IDA needs that doesn't fit into the
85			common IntegratorSystem struct.
86			*/
87			typedef struct{
88			struct rel_relation rellist; /< NULL terminated list of rels */
89	johnpye	912	struct var_variable varlist; /< NULL terminated list of vars. ONLY USED FOR DEBUGGING -- get rid of it! */
90	johnpye	669	int nrels;
91			int safeeval; /*< whether to pass the 'safe' flag to relman_eval /
92			} IntegratorIdaData;
93
94			/*-------------------------------------------------------------
95			FORWARD DECLS
96			*/
97			/* residual function forward declaration */
98	johnpye	903	int integrator_ida_fex(realtype tt, N_Vector yy, N_Vector yp, N_Vector rr, void *res_data);
99	johnpye	669
100	johnpye	909	int integrator_ida_jvex(realtype tt, N_Vector yy, N_Vector yp, N_Vector rr
101			, N_Vector v, N_Vector Jv, realtype c_j
102			, void *jac_data, N_Vector tmp1, N_Vector tmp2
103			);
104
105	johnpye	669	/* error handler forward declaration */
106			void integrator_ida_error(int error_code
107			, const char module, const char function
108			, char msg, void eh_data
109			);
110
111			/*-------------------------------------------------------------
112			SETUP/TEARDOWN ROUTINES
113			*/
114			void integrator_ida_create(IntegratorSystem *blsys){
115			CONSOLE_DEBUG("ALLOCATING IDA ENGINE DATA");
116			IntegratorIdaData *enginedata;
117			enginedata = ASC_NEW(IntegratorIdaData);
118			enginedata->rellist = NULL;
119	johnpye	909	enginedata->varlist = NULL;
120	johnpye	669	enginedata->safeeval = 1;
121	ben.allan	704	blsys->enginedata = (void *)enginedata;
122			}
123	johnpye	669
124			void integrator_ida_free(void *enginedata){
125			CONSOLE_DEBUG("DELETING IDA ENGINE DATA");
126			IntegratorIdaData d = (IntegratorIdaData )enginedata;
127			/* note, we don't own the rellist, so don't need to free it */
128			ASC_FREE(d);
129			}
130
131			IntegratorIdaData integrator_ida_enginedata(IntegratorSystem blsys){
132			IntegratorIdaData *d;
133			assert(blsys!=NULL);
134			assert(blsys->enginedata!=NULL);
135			assert(blsys->engine==INTEG_IDA);
136			d = ((IntegratorIdaData *)(blsys->enginedata));
137			assert(d->safeeval = 1);
138			return d;
139			}
140
141			/*-------------------------------------------------------------
142			MAIN IDA SOLVER ROUTINE, see IDA manual, sec 5.4, p. 27 ff.
143			*/
144
145			/* return 1 on success */
146			int integrator_ida_solve(
147			IntegratorSystem *blsys
148			, unsigned long start_index
149			, unsigned long finish_index
150			){
151			void *ida_mem;
152	johnpye	719	int size, flag, t_index;
153	johnpye	669	realtype t0, reltol, t, tret, tout1;
154			N_Vector y0, yp0, abstol, ypret, yret;
155			IntegratorIdaData *enginedata;
156
157			CONSOLE_DEBUG("STARTING IDA...");
158
159			enginedata = integrator_ida_enginedata(blsys);
160			CONSOLE_DEBUG("safeeval = %d",enginedata->safeeval);
161
162			/* store reference to list of relations (in enginedata) */
163			enginedata->nrels = slv_get_num_solvers_rels(blsys->system);
164			enginedata->rellist = slv_get_solvers_rel_list(blsys->system);
165	johnpye	909	enginedata->varlist = slv_get_solvers_var_list(blsys->system);
166	johnpye	669	CONSOLE_DEBUG("Number of relations: %d",enginedata->nrels);
167	johnpye	719	CONSOLE_DEBUG("Number of dependent vars: %ld",blsys->n_y);
168	johnpye	669	size = blsys->n_y;
169
170			if(enginedata->nrels!=size){
171			ERROR_REPORTER_HERE(ASC_USER_ERROR,"Integration problem is not square (%d rels, %d vars)", enginedata->nrels, size);
172			return 0; /* failure */
173			}
174
175			/* retrieve initial values from the system */
176			t0 = samplelist_get(blsys->samples,start_index);
177
178			y0 = N_VNew_Serial(size);
179			integrator_get_y(blsys,NV_DATA_S(y0));
180
181			yp0 = N_VNew_Serial(size);
182			integrator_get_ydot(blsys,NV_DATA_S(yp0));
183
184			N_VPrint_Serial(yp0);
185			CONSOLE_DEBUG("yp0 is at %p",&yp0);
186
187			/* create IDA object */
188			ida_mem = IDACreate();
189
190			/* retrieve the absolute tolerance values for each variable */
191			abstol = N_VNew_Serial(size);
192			N_VConst(0.1,abstol); /** @TODO fill in the abstol values from the model */
193			reltol = 0.001;
194
195			/* allocate internal memory */
196	johnpye	903	flag = IDAMalloc(ida_mem, &integrator_ida_fex, t0, y0, yp0, IDA_SV, reltol, abstol);
197	johnpye	669	if(flag==IDA_MEM_NULL){
198			ERROR_REPORTER_HERE(ASC_PROG_ERR,"ida_mem is NULL");
199			return 0;
200			}else if(flag==IDA_MEM_FAIL){
201			ERROR_REPORTER_HERE(ASC_PROG_ERR,"Unable to allocate memory (IDAMalloc)");
202			return 0;
203			}else if(flag==IDA_ILL_INPUT){
204			ERROR_REPORTER_HERE(ASC_PROG_ERR,"Invalid input to IDAMalloc");
205			return 0;
206			}/* else success */
207
208			/* set optional inputs... */
209			IDASetErrHandlerFn(ida_mem, &integrator_ida_error, (void *)blsys);
210			IDASetRdata(ida_mem, (void *)blsys);
211			IDASetMaxStep(ida_mem, integrator_get_maxstep(blsys));
212			IDASetInitStep(ida_mem, integrator_get_stepzero(blsys));
213			IDASetMaxNumSteps(ida_mem, integrator_get_maxsubsteps(blsys));
214			/* there's no capability for setting minimum step size in IDA */
215
216	johnpye	912	CONSOLE_DEBUG("ASSIGNING LINEAR SOLVER");
217
218	johnpye	669	/* attach linear solver module, using the default value of maxl */
219			flag = IDASpgmr(ida_mem, 0);
220			if(flag==IDASPILS_MEM_NULL){
221			ERROR_REPORTER_HERE(ASC_PROG_ERR,"ida_mem is NULL");
222			return 0;
223			}else if(flag==IDASPILS_MEM_FAIL){
224			ERROR_REPORTER_HERE(ASC_PROG_ERR,"Unable to allocate memory (IDASpgmr)");
225			return 0;
226			}/* else success */
227
228	johnpye	909	/* assign the Jv function /
229			flag = IDASpilsSetJacTimesVecFn(ida_mem, &integrator_ida_jvex, (void *)blsys);
230			if(flag==IDASPILS_MEM_NULL){
231			ERROR_REPORTER_HERE(ASC_PROG_ERR,"ida_mem is NULL");
232			return 0;
233			}else if(flag==IDASPILS_LMEM_NULL){
234			ERROR_REPORTER_HERE(ASC_PROG_ERR,"IDASPILS linear solver has not been initialized");
235			return 0;
236			}/* else success */
237	johnpye	669
238	johnpye	909	/* set linear solver optional inputs...
239	johnpye	669
240	johnpye	909	...nothing here at the moment...
241
242			*/
243
244	johnpye	669	/* correct initial values, given derivatives */
245			blsys->currentstep=0;
246			t_index=start_index+1;
247			tout1 = samplelist_get(blsys->samples, t_index);
248	johnpye	913
249			#if SUNDIALS_VERSION_MAJOR==2 && SUNDIALS_VERSION_MINOR==3
250			flag = IDACalcIC(ida_mem, IDA_Y_INIT, tout1);
251			#else
252	johnpye	669	flag = IDACalcIC(ida_mem, t0, y0, yp0, IDA_Y_INIT, tout1);
253	johnpye	913	#endif
254
255	johnpye	669	if(flag!=IDA_SUCCESS){
256			ERROR_REPORTER_HERE(ASC_PROG_ERR,"Unable to solve initial values (IDACalcIC)");
257			return 0;
258			}/* else success */
259
260			CONSOLE_DEBUG("INITIAL CONDITIONS SOLVED :-)");
261
262			/* optionally, specify ROO-FINDING PROBLEM */
263
264			/* -- set up the IntegratorReporter */
265			integrator_output_init(blsys);
266
267			/* -- store the initial values of all the stuff */
268			integrator_output_write(blsys);
269			integrator_output_write_obs(blsys);
270
271			/* specify where the returned values should be stored */
272			yret = y0;
273			ypret = yp0;
274
275			/* advance solution in time, return values as yret and derivatives as ypret */
276			blsys->currentstep=1;
277			for(t_index=start_index+1;t_index <= finish_index;++t_index, ++blsys->currentstep){
278			t = samplelist_get(blsys->samples, t_index);
279
280			CONSOLE_DEBUG("SOLVING UP TO t = %f", t);
281
282			flag = IDASolve(ida_mem, t, &tret, yret, ypret, IDA_NORMAL);
283
284			/* pass the values of everything back to the compiler */
285			integrator_set_t(blsys, (double)tret);
286			integrator_set_y(blsys, NV_DATA_S(yret));
287			integrator_set_ydot(blsys, NV_DATA_S(ypret));
288
289			if(flag<0){
290			ERROR_REPORTER_HERE(ASC_PROG_ERR,"Failed to solve t = %f (IDASolve), error %d", t, flag);
291			break;
292			}
293
294			/* -- do something so that blsys knows the values of tret, yret and ypret */
295
296			/* -- store the current values of all the stuff */
297			integrator_output_write(blsys);
298			integrator_output_write_obs(blsys);
299			}
300
301			/* -- close the IntegratorReporter */
302			integrator_output_close(blsys);
303
304			if(flag < 0){
305			ERROR_REPORTER_HERE(ASC_PROG_ERR,"Solving aborted while attempting t = %f", t);
306			return 0;
307			}
308
309			/* get optional outputs */
310
311			/* free solution memory */
312			N_VDestroy_Serial(yret);
313			N_VDestroy_Serial(ypret);
314
315			/* free solver memory */
316			IDAFree(ida_mem);
317
318			/* all done */
319			return 1;
320			}
321
322			/*--------------------------------------------------
323			RESIDUALS AND JACOBIAN
324			*/
325			/**
326			Function to evaluate system residuals, in the form required for IDA.
327
328			Given tt, yy and yp, we need to evaluate and return rr.
329
330			@param tt current value of indep variable (time)
331			@param yy current values of dependent variable vector
332			@param yp current values of derivatives of dependent variables
333			@param rr the output residual vector (is we're returning data to)
334			@param res_data pointer to our stuff (blsys in this case).
335
336			@return 0 on success, positive on recoverable error, and
337			negative on unrecoverable error.
338			*/
339	johnpye	903	int integrator_ida_fex(realtype tt, N_Vector yy, N_Vector yp, N_Vector rr, void *res_data){
340	johnpye	669	IntegratorSystem *blsys;
341			IntegratorIdaData *enginedata;
342			int i, calc_ok, is_error;
343			struct rel_relation** relptr;
344			double resid;
345			char *relname;
346
347			blsys = (IntegratorSystem *)res_data;
348			enginedata = integrator_ida_enginedata(blsys);
349
350	johnpye	903	/* fprintf(stderr,"\n\n"); */
351			/* CONSOLE_DEBUG("ABOUT TO EVALUTE RESIDUALS..."); */
352	johnpye	669
353			/* pass the values of everything back to the compiler */
354			integrator_set_t(blsys, (double)tt);
355			integrator_set_y(blsys, NV_DATA_S(yy));
356			integrator_set_ydot(blsys, NV_DATA_S(yp));
357
358			/* revaluate the system residuals using the new data */
359			is_error = 0;
360			relptr = enginedata->rellist;
361
362	johnpye	903	/* CONSOLE_DEBUG("IDA requests residuals of length %lu",NV_LENGTH_S(rr)); */
363	johnpye	669	if(NV_LENGTH_S(rr)!=enginedata->nrels){
364			ERROR_REPORTER_HERE(ASC_PROG_ERR,"Invalid residuals nrels!=length(rr)");
365			return -1; /* unrecoverable */
366			}
367
368			Asc_SignalHandlerPush(SIGFPE,SIG_IGN);
369			if (setjmp(g_fpe_env)==0) {
370			for(i=0, relptr = enginedata->rellist;
371			i< enginedata->nrels && relptr != NULL;
372			++i, ++relptr
373			){
374			resid = relman_eval(*relptr, &calc_ok, enginedata->safeeval);
375
376			relname = rel_make_name(blsys->system, *relptr);
377	johnpye	903	/* if(calc_ok){
378	johnpye	669	CONSOLE_DEBUG("residual[%d:\"%s\"] = %f",i,relname,resid);
379			}else{
380			CONSOLE_DEBUG("residual[%d:\"%s\"] = %f (ERROR)",i,relname,resid);
381	johnpye	903	}*/
382	johnpye	669	ASC_FREE(relname);
383
384			NV_Ith_S(rr,i) = resid;
385			if(!calc_ok){
386			/* presumable some output already made? */
387			is_error = 1;
388			}
389			}
390			}else{
391			CONSOLE_DEBUG("FLOATING POINT ERROR WITH i=%d",i);
392			}
393			Asc_SignalHandlerPop(SIGFPE,SIG_IGN);
394
395			if(is_error)CONSOLE_DEBUG("SOME ERRORS FOUND IN EVALUATION");
396			return is_error;
397			}
398
399	johnpye	903	/**
400	johnpye	912	Function to evaluate the product J*v, in the form required for IDA (see IDASpilsSetJacTimesVecFn)
401	johnpye	909
402			Given tt, yy, yp, rr and v, we need to evaluate and return Jv.
403
404			@param tt current value of the independent variable (time, t)
405			@param yy current value of the dependent variable vector, y(t).
406			@param yp current value of y'(t).
407			@param rr current value of the residual vector F(t, y, y').
408			@param v the vector by which the Jacobian must be multiplied to the right.
409			@param Jv the output vector computed
410			@param c_j the scalar in the system Jacobian, proportional to the inverse of the step size ($ \alpha$ in Eq. (3.5) ).
411			@param jac_data pointer to our stuff (blsys in this case, passed into IDA via IDASp*SetJacTimesVecFn.)
412			@param tmp1 @see tmp2
413			@param tmp2 (as well as tmp1) pointers to memory allocated for variables of type N_Vector for use here as temporary storage or work space.
414			@return 0 on success
415	johnpye	903	*/
416	johnpye	909	int integrator_ida_jvex(realtype tt, N_Vector yy, N_Vector yp, N_Vector rr
417			, N_Vector v, N_Vector Jv, realtype c_j
418			, void *jac_data, N_Vector tmp1, N_Vector tmp2
419			){
420			IntegratorSystem *blsys;
421			IntegratorIdaData *enginedata;
422			int i, j, is_error=0;
423			struct rel_relation** relptr;
424			char relname, varname;
425			int status;
426			double Jv_i;
427	johnpye	903
428	johnpye	909	int *variables;
429			double *derivatives;
430			var_filter_t filter;
431			int count;
432
433			CONSOLE_DEBUG("EVALUTING JACOBIAN...");
434
435			blsys = (IntegratorSystem *)jac_data;
436			enginedata = integrator_ida_enginedata(blsys);
437
438			/* pass the values of everything back to the compiler */
439			integrator_set_t(blsys, (double)tt);
440			integrator_set_y(blsys, NV_DATA_S(yy));
441			integrator_set_ydot(blsys, NV_DATA_S(yp));
442			/* no real use for residuals (rr) here, I don't think? */
443
444	johnpye	912	/* allocate space for returns from relman_diff2: we should be able to use 'tmp1' and 'tmp2' here... */
445	johnpye	909	variables = ASC_NEW_ARRAY(int, NV_LENGTH_S(yy) * 2);
446			derivatives = ASC_NEW_ARRAY(double, NV_LENGTH_S(yy) * 2);
447
448			/* evaluate the derivatives... */
449			/* J = dG_dy = dF_dy + alpha * dF_dyp */
450
451			filter.matchbits = VAR_SVAR;
452			filter.matchvalue = VAR_SVAR;
453
454	johnpye	912	CONSOLE_DEBUG("PRINTING VALUES OF 'v' VECTOR (length %ld)",NV_LENGTH_S(v));
455			for(i=0; i<NV_LENGTH_S(v); ++i){
456			varname = var_make_name(blsys->system, enginedata->varlist[i]);
457			if(varname==NULL)varname="UNKNOWN";
458			CONSOLE_DEBUG("v[%d] = '%s' = %f",i,varname,NV_Ith_S(v,i));
459			ASC_FREE(varname);
460			}
461
462	johnpye	909	Asc_SignalHandlerPush(SIGFPE,SIG_IGN);
463			if (setjmp(g_fpe_env)==0) {
464			for(i=0, relptr = enginedata->rellist;
465			i< enginedata->nrels && relptr != NULL;
466			++i, ++relptr
467			){
468			/* get derivatives for this particular relation */
469			status = relman_diff2(*relptr, &filter, derivatives, variables, &count, enginedata->safeeval);
470
471			CONSOLE_DEBUG("Got derivatives against %d matching variables", count);
472
473			relname = rel_make_name(blsys->system, *relptr);
474			if(!status){
475			CONSOLE_DEBUG("Derivatives for relation %d '%s' OK",i,relname);
476			}else{
477			CONSOLE_DEBUG("ERROR calculating derivatives for relation %d '%s'",i,relname);
478			break;
479			}
480			ASC_FREE(relname);
481
482			Jv_i = 0;
483			for(j=0; j < count; ++j){
484	johnpye	910	/* CONSOLE_DEBUG("j = %d, variables[j] = %d, n_y = %ld", j, variables[j], blsys->n_y); */
485	johnpye	909	varname = var_make_name(blsys->system, enginedata->varlist[variables[j]]);
486			if(varname){
487			CONSOLE_DEBUG("Variable %d '%s' derivative = %f", variables[j],varname,derivatives[j]);
488			ASC_FREE(varname);
489			}else{
490	johnpye	912	CONSOLE_DEBUG("Variable %d (UNKNOWN!): derivative = %f",variables[j],derivatives[j]);
491	johnpye	909	}
492	johnpye	912
493	johnpye	909	Jv_i += derivatives[j] * NV_Ith_S(v,variables[j]);
494			}
495
496			NV_Ith_S(Jv,i) = Jv_i;
497	johnpye	910	CONSOLE_DEBUG("(J*v)[%d] = %f", i, Jv_i);
498
499	johnpye	909	if(status){
500			/* presumably some error_reporter will already have been made*/
501			is_error = 1;
502			}
503			}
504			}else{
505			CONSOLE_DEBUG("FLOATING POINT ERROR WITH i=%d",i);
506			}
507			Asc_SignalHandlerPop(SIGFPE,SIG_IGN);
508
509			if(is_error)CONSOLE_DEBUG("SOME ERRORS FOUND IN EVALUATION");
510	johnpye	912
511
512
513	johnpye	909	return is_error;
514			}
515
516	johnpye	669	/*----------------------------------------------
517			ERROR REPORTING
518	johnpye	782	*/
519	johnpye	669	/**
520			Error message reporter function to be passed to IDA. All error messages
521			will trigger a call to this function, so we should find everything
522			appearing on the console (in the case of Tcl/Tk) or in the errors/warnings
523			panel (in the case of PyGTK).
524			*/
525			void integrator_ida_error(int error_code
526			, const char module, const char function
527			, char msg, void eh_data
528			){
529			IntegratorSystem *blsys;
530			error_severity_t sev;
531
532			/* cast back the IntegratorSystem, just in case we need it */
533			blsys = (IntegratorSystem *)eh_data;
534
535			/* severity depends on the sign of the error_code value */
536			if(error_code <= 0){
537			sev = ASC_PROG_ERR;
538			}else{
539			sev = ASC_PROG_WARNING;
540			}
541
542			/* use our all-purpose error reporting to get stuff back to the GUI */
543			error_reporter(sev,module,0,function,"%s (error %d)",msg,error_code);
544			}