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
# ifdef __GNUC__
#  warning "GUESSING SUNDIALS VERSION 2.2"
# endif
# 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 */
    }

    CONSOLE_DEBUG("RETRIEVING t0");

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

    CONSOLE_DEBUG("RETRIEVING y0");

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

    CONSOLE_DEBUG("RETRIEVING yp0");

    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 var_yindex;

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

    fprintf(stderr,"\n--------------\n");
    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
        ){
            fprintf(stderr,"\n");
            relname = rel_make_name(blsys->system, *relptr);
            CONSOLE_DEBUG("RELATION %d '%s'",i,relname);
            ASC_FREE(relname);

            /* 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);

            for(j=0;j<count;++j){
                varname = var_make_name(blsys->system, enginedata->varlist[variables[j]]);
                CONSOLE_DEBUG("derivatives[%d] = %f (variable %d, '%s')",j,derivatives[j],variables[j],varname);
                ASC_FREE(varname);
            }

            if(!status){
                CONSOLE_DEBUG("Derivatives for relation %d OK",i);
            }else{
                CONSOLE_DEBUG("ERROR calculating derivatives for relation %d",i);
                break;
            }

            /*
                Now we have the derivatives wrt each alg/diff variable in the
                present equation. variables[] points into the varlist. need
                a mapping from the varlist to the y and ydot lists.
            */

            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]);
                }
                
                /* this bit is still not right!!! */
                var_yindex = blsys->y_id[variables[j]];

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