ascend/integrator/aww.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
    C-code implementation of Art Westerberg's DAE integrator
*//*
    by John Pye, Dec 2006.
*/

#include <utilities/config.h>
#include <utilities/ascConfig.h>
#include <utilities/ascPanic.h>
#include "aww.h"

const IntegratorInternals integrator_aww_internals = {
    integrator_aww_create
    ,integrator_aww_params_default
    ,integrator_aww_analyse
    ,integrator_aww_solve
    ,NULL /* writematrixfn */
    ,NULL /* debugfn */
    ,integrator_aww_free
    ,INTEG_AWW
    ,"AWW"
};

typedef struct{
    unsigned nothing;
}  IntegratorAWWData;

/*-------------------------------------------------------------
  PARAMETERS FOR AWW INTEGRATOR
*/

enum aww_parameters{
    AWW_PARAM_METHOD
    ,AWW_PARAM_RTOL
    ,AWW_PARAMS_SIZE
};

enum aww_methods{
    AWW_BDF
    ,AWW_AM
};
    
/**
    Here the full set of parameters is defined, along with upper/lower bounds,
    etc. The values are stuck into the blsys->params structure.

    @return 0 on success
*/
int integrator_aww_params_default(IntegratorSystem *blsys){
    asc_assert(blsys!=NULL);
    asc_assert(blsys->engine==INTEG_AWW);
    slv_parameters_t *p;
    p = &(blsys->params);

    slv_destroy_parms(p);

    if(p->parms==NULL){
        CONSOLE_DEBUG("params NULL");
        p->parms = ASC_NEW_ARRAY(struct slv_parameter, AWW_PARAMS_SIZE);
        if(p->parms==NULL)return -1;
        p->dynamic_parms = 1;
    }else{
        CONSOLE_DEBUG("params not NULL");
    }

    /* reset the number of parameters to zero so that we can check it at the end */
    p->num_parms = 0;

    slv_param_char(p,AWW_PARAM_METHOD
            ,(SlvParameterInitChar){{"method"
            ,"Stepping method",1
            ,"See Art's notes, sec. 15.2"
        }, "BDF"}, (char *[]){"BDF","AM",NULL}
    );

    slv_param_real(p,AWW_PARAM_RTOL
            ,(SlvParameterInitReal){{"rtol"
            ,"Scalar relative error tolerance",1
            ,"Value of the scalar relative error tolerance."
        }, 1e-4, 0, DBL_MAX }
    );

    asc_assert(p->num_parms == AWW_PARAMS_SIZE);

    CONSOLE_DEBUG("Created %d params", p->num_parms);

    return 0;
}

void integrator_aww_create(IntegratorSystem *blsys){
    CONSOLE_DEBUG("ALLOCATING AWW ENGINE DATA");
    IntegratorAWWData *enginedata;
    enginedata = ASC_NEW(IntegratorAWWData);
    enginedata->nothing = 0;
    blsys->enginedata = (void *)enginedata;
    integrator_aww_params_default(blsys);
}

void integrator_aww_free(void *enginedata){
    CONSOLE_DEBUG("DELETING AWW ENGINE DATA");
    IntegratorAWWData *d = (IntegratorAWWData *)enginedata;
    ASC_FREE(d);
}

/** return 1 on success */
int integrator_aww_analyse(IntegratorSystem *blsys){
    return 0;
}

/** return 0 on success */
int integrator_aww_solve(IntegratorSystem *blsys
        , unsigned long start_index, unsigned long finish_index
){
    double initstep, maxstep, rtol;
    unsigned maxsubsteps;
    const char *methodname;

    // solve initial point

    // run valuesForInit
    // run specifyForInit
    // solve with QRSlv
    // if not converged throw error

    // 'DELETE SYSTEM'
    // run valuesForStep
    // run specifyForStep
    
    initstep = integrator_get_stepzero(blsys);
    maxstep = integrator_get_maxstep(blsys);
    maxsubsteps = integrator_get_maxsubsteps(blsys);
    rtol = SLV_PARAM_REAL(&(blsys->params),AWW_PARAM_RTOL);

    methodname = SLV_PARAM_CHAR(&(blsys->params),AWW_PARAM_METHOD);
    CONSOLE_DEBUG("ASSIGNING STEPPING METHOD '%s'",methodname);
    
    enum aww_methods method;
    if(strcmp(methodname,"BDF")==0)method = AWW_BDF;
    else if(strcmp(methodname,"AM")==0)method = AWW_AM;
    else{
        ERROR_REPORTER_HERE(ASC_USER_ERROR,"Invalid method name '%s'",methodname);
    }

    /* initialise the integrator reporter */
    integrator_output_init(blsys);

    // while not too many solver calls {
    //    capture plotting points
    //    do something with polynomials
    //    for(i=0; i<= poly_oder; ++i){
    //      if(!in the final step){
    //          if(!have hit stopping configtion) UN 'stepX'
    //          SOLVE with QRSlv
    //          if(!converged)abort 'failed to convert
    //      }
    //      if(in the final step)abort 'STOP reached'
    //      RUN setStopConditions
    //      retrieve status of stop confitions
    //      retureve status of 'last step'
    //    }
    //    compute maximum steps for each variable
    // }

    ERROR_REPORTER_HERE(ASC_PROG_ERR,"Not implemented");
    return 1;
}

1	johnpye	970	/* 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			//* @file
19			C-code implementation of Art Westerberg's DAE integrator
20			//
21			by John Pye, Dec 2006.
22			*/
23
24			#include <utilities/config.h>
25			#include <utilities/ascConfig.h>
26			#include <utilities/ascPanic.h>
27	johnpye	977	#include "aww.h"
28	johnpye	970
29	johnpye	977	const IntegratorInternals integrator_aww_internals = {
30			integrator_aww_create
31			,integrator_aww_params_default
32			,integrator_aww_analyse
33			,integrator_aww_solve
34	johnpye	1129	,NULL /* writematrixfn */
35	johnpye	1196	,NULL /* debugfn */
36	johnpye	977	,integrator_aww_free
37			,INTEG_AWW
38			,"AWW"
39			};
40
41	johnpye	970	typedef struct{
42			unsigned nothing;
43			} IntegratorAWWData;
44
45			/*-------------------------------------------------------------
46			PARAMETERS FOR AWW INTEGRATOR
47			*/
48
49			enum aww_parameters{
50			AWW_PARAM_METHOD
51			,AWW_PARAM_RTOL
52			,AWW_PARAMS_SIZE
53			};
54
55			enum aww_methods{
56			AWW_BDF
57			,AWW_AM
58			};
59
60			/**
61			Here the full set of parameters is defined, along with upper/lower bounds,
62			etc. The values are stuck into the blsys->params structure.
63
64			@return 0 on success
65			*/
66			int integrator_aww_params_default(IntegratorSystem *blsys){
67			asc_assert(blsys!=NULL);
68			asc_assert(blsys->engine==INTEG_AWW);
69			slv_parameters_t *p;
70			p = &(blsys->params);
71
72			slv_destroy_parms(p);
73
74			if(p->parms==NULL){
75			CONSOLE_DEBUG("params NULL");
76			p->parms = ASC_NEW_ARRAY(struct slv_parameter, AWW_PARAMS_SIZE);
77			if(p->parms==NULL)return -1;
78			p->dynamic_parms = 1;
79			}else{
80			CONSOLE_DEBUG("params not NULL");
81			}
82
83			/* reset the number of parameters to zero so that we can check it at the end */
84			p->num_parms = 0;
85
86			slv_param_char(p,AWW_PARAM_METHOD
87			,(SlvParameterInitChar){{"method"
88			,"Stepping method",1
89			,"See Art's notes, sec. 15.2"
90			}, "BDF"}, (char *[]){"BDF","AM",NULL}
91			);
92
93			slv_param_real(p,AWW_PARAM_RTOL
94			,(SlvParameterInitReal){{"rtol"
95			,"Scalar relative error tolerance",1
96			,"Value of the scalar relative error tolerance."
97			}, 1e-4, 0, DBL_MAX }
98			);
99
100			asc_assert(p->num_parms == AWW_PARAMS_SIZE);
101
102			CONSOLE_DEBUG("Created %d params", p->num_parms);
103
104			return 0;
105			}
106
107			void integrator_aww_create(IntegratorSystem *blsys){
108			CONSOLE_DEBUG("ALLOCATING AWW ENGINE DATA");
109			IntegratorAWWData *enginedata;
110			enginedata = ASC_NEW(IntegratorAWWData);
111			enginedata->nothing = 0;
112			blsys->enginedata = (void *)enginedata;
113			integrator_aww_params_default(blsys);
114			}
115
116			void integrator_aww_free(void *enginedata){
117			CONSOLE_DEBUG("DELETING AWW ENGINE DATA");
118			IntegratorAWWData d = (IntegratorAWWData )enginedata;
119			ASC_FREE(d);
120			}
121
122			/** return 1 on success */
123			int integrator_aww_analyse(IntegratorSystem *blsys){
124			return 0;
125			}
126
127	johnpye	1049	/** return 0 on success */
128	johnpye	970	int integrator_aww_solve(IntegratorSystem *blsys
129			, unsigned long start_index, unsigned long finish_index
130			){
131			double initstep, maxstep, rtol;
132			unsigned maxsubsteps;
133			const char *methodname;
134
135			// solve initial point
136
137			// run valuesForInit
138			// run specifyForInit
139			// solve with QRSlv
140			// if not converged throw error
141
142			// 'DELETE SYSTEM'
143			// run valuesForStep
144			// run specifyForStep
145
146			initstep = integrator_get_stepzero(blsys);
147			maxstep = integrator_get_maxstep(blsys);
148			maxsubsteps = integrator_get_maxsubsteps(blsys);
149			rtol = SLV_PARAM_REAL(&(blsys->params),AWW_PARAM_RTOL);
150
151			methodname = SLV_PARAM_CHAR(&(blsys->params),AWW_PARAM_METHOD);
152			CONSOLE_DEBUG("ASSIGNING STEPPING METHOD '%s'",methodname);
153
154			enum aww_methods method;
155			if(strcmp(methodname,"BDF")==0)method = AWW_BDF;
156			else if(strcmp(methodname,"AM")==0)method = AWW_AM;
157			else{
158			ERROR_REPORTER_HERE(ASC_USER_ERROR,"Invalid method name '%s'",methodname);
159			}
160
161			/* initialise the integrator reporter */
162			integrator_output_init(blsys);
163
164			// while not too many solver calls {
165			// capture plotting points
166			// do something with polynomials
167			// for(i=0; i<= poly_oder; ++i){
168			// if(!in the final step){
169			// if(!have hit stopping configtion) UN 'stepX'
170			// SOLVE with QRSlv
171			// if(!converged)abort 'failed to convert
172			// }
173			// if(in the final step)abort 'STOP reached'
174			// RUN setStopConditions
175			// retrieve status of stop confitions
176			// retureve status of 'last step'
177			// }
178			// compute maximum steps for each variable
179			// }
180
181			ERROR_REPORTER_HERE(ASC_PROG_ERR,"Not implemented");
182	johnpye	1049	return 1;
183	johnpye	970	}
184