/[ascend]/trunk/base/generic/solver/aww.c
ViewVC logotype

Contents of /trunk/base/generic/solver/aww.c

Parent Directory Parent Directory | Revision Log Revision Log


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

john.pye@anu.edu.au
ViewVC Help
Powered by ViewVC 1.1.22