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