/[ascend]/trunk/pygtk/integrator.cpp
ViewVC logotype

Annotation of /trunk/pygtk/integrator.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1306 - (hide annotations) (download) (as text)
Sat Mar 3 11:50:47 2007 UTC (15 years, 5 months ago) by johnpye
File MIME type: text/x-c++src
File size: 7799 byte(s)
Added and performed basic tests of integrator_ida_write_matrix. Generalised the write_matrix
routine so that *any* requested output can be retrieved from the integrator (for the case
of IDA this is y and y', but it could equally be more complicated stuff.)
1 johnpye 669 #include "integrator.h"
2     #include "integratorreporter.h"
3 johnpye 942 #include "solverparameters.h"
4 johnpye 669 #include <stdexcept>
5 johnpye 941 #include <sstream>
6 johnpye 942 #include <cmath>
7 johnpye 669 using namespace std;
8    
9 johnpye 1303 // #define DESTROY_DEBUG
10    
11 johnpye 669 /**
12     'creating' an integrator in the context of the GUI just means an object
13     we can store the parameters that will be later sent to the underlying
14     C-code API.
15 johnpye 1199
16     @TODO at present the integrator requires a slv_system_t. This is the wrong
17     way around.
18 johnpye 669 */
19     Integrator::Integrator(Simulation &simulation)
20     : simulation(simulation)
21 johnpye 1087 {
22 johnpye 669 // create the C-level object
23     this->blsys = integrator_new(simulation.getSystem(),simulation.getModel().getInternalType());
24    
25     samplelist = NULL;
26    
27     // set default steps
28     setMinSubStep(0);
29     setMaxSubStep(0);
30     setMaxSubSteps(0);
31     setInitialSubStep(0);
32     }
33    
34     Integrator::~Integrator(){
35 johnpye 1303 #ifdef DESTROY_DEBUG
36 johnpye 921 CONSOLE_DEBUG("DESTROYING Integrator (C++) at %p",this);
37     CONSOLE_DEBUG("DESTROYING IntegratorSystem at %p",blsys);
38 johnpye 1303 #endif
39 johnpye 669 integrator_free(blsys);
40 johnpye 1303 #ifdef DESTROY_DEBUG
41 johnpye 1221 CONSOLE_DEBUG("done");
42 johnpye 921 CONSOLE_DEBUG("DESTROYING samplelist at %p",samplelist);
43 johnpye 1303 #endif
44 johnpye 1201 if(samplelist)samplelist_free(samplelist);
45 johnpye 669 }
46    
47 johnpye 1087 SolverParameters
48 johnpye 942 Integrator::getParameters() const{
49     SolverParameters params;
50     int res = integrator_params_get(blsys,&(params.getInternalType() ) );
51     if(res)throw runtime_error("Failed to get integrator parameters");
52     return params;
53     }
54 johnpye 669
55     void
56 johnpye 942 Integrator::setParameters(const SolverParameters &params){
57     int res = integrator_params_set(blsys,&(params.getInternalTypeConst() ) );
58     if(res)throw runtime_error("Failed to set integrator parameters");
59     }
60    
61     void
62 johnpye 669 Integrator::setReporter(IntegratorReporterCxx *reporter){
63     this->blsys->clientdata = reporter;
64     integrator_set_reporter(blsys,reporter->getInternalType());
65 johnpye 1087 //CONSOLE_DEBUG("REPORTER HAS BEEN SET");
66 johnpye 669 (*(this->blsys->reporter->init))(blsys);
67 johnpye 966 //CONSOLE_DEBUG("DONE TESTING OUTPUT_INIT");
68 johnpye 669 }
69    
70     double
71     Integrator::getCurrentTime(){
72     return integrator_get_t(blsys);
73     }
74    
75     long
76     Integrator::getCurrentStep(){
77     return integrator_getcurrentstep(blsys);
78     }
79    
80     long
81     Integrator::getNumSteps(){
82     return integrator_getnsamples(blsys);
83     }
84    
85 johnpye 1087 /**
86     Find the independent variable in the system, or throw an exception if not found.
87     */
88     void
89 johnpye 669 Integrator::findIndependentVar(){
90 johnpye 1087 int res = integrator_find_indep_var(blsys);
91    
92     if(res){
93     stringstream ss;
94     ss << "Independent variable not found (" << res << ")";
95     throw runtime_error(ss.str());
96     }
97 johnpye 669 }
98    
99 johnpye 979 void
100 johnpye 669 Integrator::analyse(){
101    
102     int res;
103 johnpye 1069 /*
104     Note, we never need to call analyze_make_system in any of the Integrator
105     code, as it gets called by Simulation::build.
106     */
107 johnpye 669 res = integrator_analyse(blsys);
108    
109 johnpye 1069 if(res){
110 johnpye 1097 stringstream ss; ss << "Failed system analysis (error " << res << ")";
111     throw runtime_error(ss.str());
112 johnpye 669 }
113     }
114    
115     /**
116     @TODO what about root detection?
117 johnpye 942
118     Integrate the function for the timesteps specified.
119    
120 johnpye 1049 Method will throw a runtime_error if integrator_solve returns error (non zero)
121 johnpye 953
122     @TODO does simulation.processVarStatus work for integrators like IDA???
123 johnpye 669 */
124 johnpye 942 void
125 johnpye 669 Integrator::solve(){
126    
127     // check the integration limits
128     // trigger of the solution process
129     // report errors?
130 johnpye 1087
131 johnpye 669 assert(samplelist!=NULL);
132     assert(samplelist->ns>0);
133     assert(blsys->reporter!=NULL);
134     assert(blsys->clientdata!=NULL);
135    
136     int res;
137     res = integrator_solve(blsys, 0, samplelist_length(samplelist)-1);
138 johnpye 915
139 johnpye 1049 if(res){
140     stringstream ss;
141     ss << "Failed integration (integrator_solve returned " << res << ")";
142     throw runtime_error(ss.str());
143 johnpye 669 }
144 johnpye 953
145     // communicate solver variable status back to the instance tree via 'interface_ptr'
146     simulation.processVarStatus();
147 johnpye 669 }
148    
149 johnpye 941 void
150 johnpye 1306 Integrator::writeMatrix(FILE *fp,const char *type) const{
151     if(integrator_write_matrix(this->blsys, fp, type)){
152 johnpye 1129 throw runtime_error("Failed to write matrix");
153     }
154     }
155    
156     void
157 johnpye 1196 Integrator::writeDebug(FILE *fp) const{
158     if(integrator_debug(this->blsys, fp)){
159     throw runtime_error("Failed to write debug output");
160     }
161     }
162    
163     void
164 johnpye 669 Integrator::setEngine(IntegratorEngine engine){
165 johnpye 941 int res = integrator_set_engine(this->blsys, engine);
166     if(!res)return;
167     if(res==1)throw range_error("Unknown integrator");
168     if(res==2)throw range_error("Invalid integrator");
169     stringstream ss;
170     ss << "Unknown error in setEngine (res = " << res << ")";
171     throw runtime_error(ss.str());
172 johnpye 669 }
173    
174 johnpye 941 void
175 johnpye 669 Integrator::setEngine(int engine){
176 johnpye 941 setEngine((IntegratorEngine)engine);
177 johnpye 669 }
178    
179 johnpye 941 void
180 johnpye 938 Integrator::setEngine(const string &name){
181 johnpye 941 CONSOLE_DEBUG("Setting integration engine to '%s'",name.c_str());
182 johnpye 938 IntegratorEngine engine = INTEG_UNKNOWN;
183     #ifdef ASC_WITH_LSODE
184     if(name=="LSODE")engine = INTEG_LSODE;
185     #endif
186     #ifdef ASC_WITH_IDA
187     if(name=="IDA")engine = INTEG_IDA;
188     #endif
189 johnpye 952 if(engine==INTEG_UNKNOWN){
190     throw runtime_error("Unkown integrator name");
191     }
192 johnpye 941 setEngine(engine);
193 johnpye 938 }
194    
195 johnpye 669 /**
196     Ideally this list would be dynamically generated based on what solvers
197     are available or are in memory.
198     */
199     map<int,string>
200 johnpye 941 Integrator::getEngines(){
201 johnpye 669 map<int,string> m;
202 johnpye 977 const IntegratorLookup *list = integrator_get_engines();
203 johnpye 969 while(list->id != INTEG_UNKNOWN){
204 johnpye 1071 if(list->name==NULL)throw runtime_error("list->name is NULL");
205 johnpye 969 m.insert(pair<int,string>(list->id,list->name));
206     ++list;
207     }
208 johnpye 669 return m;
209     }
210    
211     string
212 johnpye 1140 Integrator::getName() const{
213 johnpye 669 map<int,string> m=getEngines();
214     map<int,string>::iterator f = m.find(integrator_get_engine(blsys));
215     if(f==m.end()){
216     throw runtime_error("No engine selected");
217     }
218     return f->second;
219 johnpye 1087 }
220 johnpye 669
221 johnpye 942 /**
222     @TODO what about conversion factors? Is an allowance being made?
223     */
224 johnpye 669 void
225     Integrator::setLinearTimesteps(UnitsM units, double start, double end, unsigned long num){
226     if(samplelist!=NULL){
227     ASC_FREE(samplelist);
228     }
229     const dim_type *d = units.getDimensions().getInternalType();
230     samplelist = samplelist_new(num+1, d);
231     double val = start;
232     double inc = (end-start)/(num);
233 johnpye 746 for(unsigned long i=0;i<=num;++i){
234 johnpye 669 samplelist_set(samplelist,i,val);
235     val += inc;
236     }
237     integrator_set_samples(blsys,samplelist);
238     }
239    
240 johnpye 942 /**
241     @TODO what about conversion factors? Is an allowance being made?
242     */
243     void
244     Integrator::setLogTimesteps(UnitsM units, double start, double end, unsigned long num){
245     if(samplelist!=NULL){
246     ASC_FREE(samplelist);
247     }
248     const dim_type *d = units.getDimensions().getInternalType();
249    
250     if(start<=0)throw runtime_error("starting timestep needs to be > 0");
251     if(end<=0)throw runtime_error("end timestep needs to be > 0");
252     if(end <= start)throw runtime_error("end timestep needs to be > starting timestep");
253    
254     samplelist = samplelist_new(num+1, d);
255     double val = start;
256     double inc = exp((log(end)-log(start))/num);
257     for(unsigned long i=0;i<=num;++i){
258     samplelist_set(samplelist,i,val);
259     CONSOLE_DEBUG("samplelist[%lu] = %f",i,val);
260     val *= inc;
261     }
262     integrator_set_samples(blsys,samplelist);
263     }
264    
265 johnpye 669 vector<double>
266     Integrator::getCurrentObservations(){
267     double *d = ASC_NEW_ARRAY(double,getNumObservedVars());
268     integrator_get_observations(blsys,d);
269     vector<double> v=vector<double>(d,d+getNumObservedVars());
270     // do I need to free d?
271     // can I do this in such a way as I avoid all this memory-copying?
272     return v;
273     }
274    
275     Variable
276     Integrator::getObservedVariable(const long &i){
277     var_variable *v = integrator_get_observed_var(blsys,i);
278 johnpye 709 return Variable(&simulation,v);
279 johnpye 669 }
280    
281 johnpye 854 Variable
282     Integrator::getIndependentVariable(){
283     var_variable *v = integrator_get_independent_var(blsys);
284     if(v==NULL){
285     throw runtime_error("independent variable is null");
286     }
287     return Variable(&simulation,v);
288     }
289    
290 johnpye 669 int
291     Integrator::getNumVars(){
292     return blsys->n_y;
293     }
294    
295 johnpye 1087 int
296 johnpye 669 Integrator::getNumObservedVars(){
297     return blsys->n_obs;
298     }
299    
300     void
301     Integrator::setMinSubStep(double n){
302     integrator_set_minstep(blsys,n);
303     }
304    
305     void
306     Integrator::setMaxSubStep(double n){
307     integrator_set_maxstep(blsys,n);
308     }
309    
310     void
311     Integrator::setInitialSubStep(double n){
312     integrator_set_stepzero(blsys,n);
313     }
314    
315     void
316     Integrator::setMaxSubSteps(int n){
317     integrator_set_maxsubsteps(blsys,n);
318     }
319    
320     IntegratorSystem *
321     Integrator::getInternalType(){
322     return blsys;
323     }

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