/[ascend]/trunk/test.py
ViewVC logotype

Contents of /trunk/test.py

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1280 - (show annotations) (download) (as text)
Sat Feb 10 10:47:23 2007 UTC (12 years, 4 months ago) by johnpye
File MIME type: text/x-python
File size: 35418 byte(s)
Fully operational test suite again (fixed TestIntegrator.testparameters)
1 #!/usr/bin/env python
2 # ASCEND modelling environment
3 # Copyright (C) 2006, 2007 Carnegie Mellon University
4 #
5 # This program is free software; you can redistribute it and/or modify
6 # it under the terms of the GNU General Public License as published by
7 # the Free Software Foundation; either version 2, or (at your option)
8 # any later version.
9 #
10 # This program is distributed in the hope that it will be useful,
11 # but WITHOUT ANY WARRANTY; without even the implied warranty of
12 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13 # GNU General Public License for more details.
14 #
15 # You should have received a copy of the GNU General Public License
16 # along with this program; if not, write to the Free Software
17 # Foundation, Inc., 59 Temple Place - Suite 330,
18 # Boston, MA 02111-1307, USA.
19
20 # This script gives a test suite for the high-level interface of ASCEND via
21 # Python. It is also planned to be a wrapper for the CUnit test suite, although
22 # this is still experimental.
23
24 import unittest
25 import os, sys
26 import math
27 import atexit
28
29 import platform
30 if platform.system() != "Windows":
31 import dl
32 sys.setdlopenflags(dl.RTLD_GLOBAL|dl.RTLD_NOW)
33
34 class Ascend(unittest.TestCase):
35
36 def setUp(self):
37 import ascpy
38 self.L = ascpy.Library()
39
40 def tearDown(self):
41 self.L.clear()
42 del self.L
43
44 class AscendSelfTester(Ascend):
45
46 def _run(self,modelname,solvername="QRSlv",filename=None,parameters={}):
47 if filename==None:
48 filename = 'johnpye/%s.a4c' % modelname
49 self.L.load(filename)
50 T = self.L.findType(modelname)
51 M = T.getSimulation('sim')
52 M.setSolver(ascpy.Solver(solvername))
53 for k,v in parameters.iteritems():
54 M.setParameter(k,v)
55 M.solve(ascpy.Solver(solvername),ascpy.SolverReporter())
56 M.run(T.getMethod('self_test'))
57 return M
58
59 class TestCompiler(Ascend):
60
61 def _run(self,filen,modeln=""):
62 self.L.load('test/compiler/%s.a4c' % filen)
63 T = self.L.findType('%s%s' % (filen,modeln))
64 M = T.getSimulation('sim')
65 M.build()
66
67 def _runfail(self,filen,n,msg="failed"):
68 try:
69 self._run(filen,'fail%d' % n)
70 except Exception,e:
71 print "(EXPECTED) ERROR: %s" % e
72 return
73 self.fail(msg)
74
75
76 def testloading(self):
77 """library startup"""
78 pass
79
80 def testsystema4l(self):
81 """loading system.a4l?"""
82 self.L.load('system.a4l')
83
84 def testatomsa4l(self):
85 """loading atoms.a4l?"""
86 self.L.load('atoms.a4l')
87
88 def testmissingreq(self):
89 """flagging a missing REQUIRE"""
90 self._runfail('missingreq',1)
91
92 def testmissingsubreq(self):
93 """flagging a subsidiary missing REQUIRE"""
94 self._runfail('missingreq',1)
95
96 class TestSolver(AscendSelfTester):
97
98 def testlog10(self):
99 self._run('testlog10')
100
101 def testrootsofpoly(self):
102 self._run('roots_of_poly',filename="roots_of_poly.a4c")
103
104 def testcollapsingcan(self):
105 self._run('collapsingcan',filename="collapsingcan.a4c")
106
107 def testdistancecalc(self):
108 self._run('distance_calc',filename="distance_calc.a4c")
109
110 def testconopt(self):
111 self._run('testconopt',"CONOPT")
112
113 def testcmslv2(self):
114 self._run('testcmslv2',"CMSlv")
115
116 def testsunpos1(self):
117 self._run('example_1_6_1',"QRSlv","johnpye/sunpos.a4c")
118
119 def testsunpos2(self):
120 self._run('example_1_6_2',"QRSlv","johnpye/sunpos.a4c")
121
122 def testsunpos3(self):
123 self._run('example_1_7_1',"QRSlv","johnpye/sunpos.a4c")
124
125 def testsunpos4(self):
126 self._run('example_1_7_2',"QRSlv","johnpye/sunpos.a4c")
127
128 def testsunpos5(self):
129 self._run('example_1_7_3',"QRSlv","johnpye/sunpos.a4c")
130
131 def testsunpos6(self):
132 self._run('example_1_8_1',"QRSlv","johnpye/sunpos.a4c")
133
134 def testinstanceas(self):
135 M = self._run('example_1_6_1',"QRSlv","johnpye/sunpos.a4c")
136 self.assertAlmostEqual( float(M.t_solar), M.t_solar.as("s"))
137 self.assertAlmostEqual( float(M.t_solar)/3600, M.t_solar.as("h"))
138
139
140 class TestCMSlv(AscendSelfTester):
141 def testheatex(self):
142 self._run('heatex',"CMSlv","heatex.a4c")
143 def testphaseeq(self):
144 self._run('phaseq',"CMSlv","phaseq.a4c")
145 def testpipeline(self):
146 self._run('pipeline',"CMSlv","pipeline.a4c"
147 ,{'infinity':3.2e9}
148 )
149 def testrachford(self):
150 self._run('rachford',"CMSlv","rachford.a4c")
151 def testlinmassbal(self):
152 self._run('linmassbal',"CMSlv","linmassbal.a4c")
153
154
155 class TestMatrix(AscendSelfTester):
156 def testlog10(self):
157 M = self._run('testlog10')
158 print M.getMatrix().write(sys.stderr,"mmio")
159
160
161 class TestIntegrator(Ascend):
162
163 def testListIntegrators(self):
164 I = ascpy.Integrator.getEngines()
165 s1 = sorted([str(i) for i in I.values()])
166 s2 = sorted(['IDA','LSODE','AWW'])
167 assert s1==s2
168
169 # this routine is reused by both testIDA and testLSODE
170 def _testIntegrator(self,integratorname):
171 self.L.load('johnpye/shm.a4c')
172 M = self.L.findType('shm').getSimulation('sim')
173 M.setSolver(ascpy.Solver('QRSlv'))
174 P = M.getParameters()
175 M.setParameter('feastol',1e-12)
176 print M.getChildren()
177 assert float(M.x) == 10.0
178 assert float(M.v) == 0.0
179 t_end = math.pi
180
181 I = ascpy.Integrator(M)
182 I.setReporter(ascpy.IntegratorReporterNull(I))
183 I.setEngine(integratorname);
184 I.setLinearTimesteps(ascpy.Units("s"), 0.0, t_end, 100);
185 I.setMinSubStep(0.0001); # these limits are required by IDA at present (numeric diff)
186 I.setMaxSubStep(0.1);
187 I.setInitialSubStep(0.001);
188 I.setMaxSubSteps(200);
189 if(integratorname=='IDA'):
190 I.setParameter('autodiff',False)
191 for p in M.getParameters():
192 print p.getName(),"=",p.getValue()
193 I.analyse();
194 I.solve();
195 print "At end of simulation,"
196 print "x = %f" % M.x
197 print "v = %f" % M.v
198 assert abs(float(M.x) + 10) < 1e-2
199 assert abs(float(M.v)) < 1e-2
200 assert I.getNumObservedVars() == 3
201
202 def testInvalidIntegrator(self):
203 self.L.load('johnpye/shm.a4c')
204 M = self.L.findType('shm').getSimulation('sim')
205 M.setSolver(ascpy.Solver('QRSlv'))
206 I = ascpy.Integrator(M)
207 try:
208 I.setEngine('___NONEXISTENT____')
209 except RuntimeError:
210 return
211 self.fail("setEngine did not raise error!")
212
213 def testLSODE(self):
214 self._testIntegrator('LSODE')
215
216 def testIDA(self):
217 self._testIntegrator('IDA')
218
219 def testparameters(self):
220 self.L.load('johnpye/shm.a4c')
221 M = self.L.findType('shm').getSimulation('sim')
222 M.build()
223 I = ascpy.Integrator(M)
224 I.setEngine('IDA')
225 P = I.getParameters()
226 for p in P:
227 print p.getName(),"=",p.getValue()
228 assert len(P)==12
229 assert P[0].isStr()
230 assert P[0].getName()=="linsolver"
231 assert P[0].getValue()=='SPGMR'
232 assert P[2].getName()=="maxord"
233 assert P[3].getName()=="autodiff"
234 assert P[3].getValue()==True
235 assert P[8].getName()=="atolvect"
236 assert P[8].getBoolValue() == True
237 P[3].setBoolValue(False)
238 assert P[3].getBoolValue()==False
239 I.setParameters(P)
240 assert I.getParameterValue('autodiff')==False
241 I.setParameter('autodiff',True)
242 try:
243 v = I.getParameterValue('nonexist')
244 except KeyError:
245 pass
246 else:
247 self.fail('Failed to trip invalid Integrator parameter')
248
249 class TestLSODE(Ascend):
250
251 def testzill(self):
252 self.L.load('johnpye/zill.a4c')
253 T = self.L.findType('zill')
254 M = T.getSimulation('sim')
255 M.setSolver(ascpy.Solver('QRSlv'))
256 I = ascpy.Integrator(M)
257 I.setEngine('LSODE')
258 I.setMinSubStep(1e-7)
259 I.setMaxSubStep(0.001)
260 I.setMaxSubSteps(10000)
261 I.setReporter(ascpy.IntegratorReporterConsole(I))
262 I.setLinearTimesteps(ascpy.Units(), 1.0, 1.5, 5)
263 I.analyse()
264 I.solve()
265 M.run(T.getMethod('self_test'))
266
267 def testnewton(self):
268 sys.stderr.write("STARTING TESTNEWTON\n")
269 self.L.load('johnpye/newton.a4c')
270 T = self.L.findType('newton')
271 M = T.getSimulation('sim')
272 M.solve(ascpy.Solver("QRSlv"),ascpy.SolverReporter())
273 I = ascpy.Integrator(M)
274 I.setEngine('LSODE')
275 I.setParameter('rtolvect',False)
276 I.setParameter('rtol',1e-7)
277 I.setParameter('atolvect',False)
278 I.setParameter('atol',1e-7)
279 I.setMinSubStep(1e-7)
280 I.setMaxSubStep(0.001)
281 I.setMaxSubSteps(10000)
282
283 I.setReporter(ascpy.IntegratorReporterConsole(I))
284 I.setLinearTimesteps(ascpy.Units("s"), 0, 2*float(M.v)/float(M.g), 2)
285 I.analyse()
286 I.solve()
287 print "At end of simulation,"
288 print "x = %f" % M.x
289 print "v = %f" % M.v
290 M.run(T.getMethod('self_test'))
291
292 def testlotka(self):
293 self.L.load('johnpye/lotka.a4c')
294 M = self.L.findType('lotka').getSimulation('sim')
295 M.setSolver(ascpy.Solver("QRSlv"))
296 I = ascpy.Integrator(M)
297 I.setEngine('LSODE')
298 I.setReporter(ascpy.IntegratorReporterConsole(I))
299 I.setLinearTimesteps(ascpy.Units("s"), 0, 200, 5)
300 I.analyse()
301 print "Number of vars = %d" % I.getNumVars()
302 assert I.getNumVars()==2
303 I.solve()
304 assert I.getNumObservedVars() == 3;
305 assert abs(M.R - 832) < 1.0
306 assert abs(M.F - 21.36) < 0.1
307
308 #-------------------------------------------------------------------------------
309 # Testing of a external blackbox functions
310
311 class TestBlackBox(AscendSelfTester):
312 def testparsefail0(self):
313 try:
314 self.L.load('test/blackbox/parsefail0.a4c')
315 self.fail("parsefail0 should not have loaded without errors")
316 except:
317 pass
318
319 def testparsefail1(self):
320 try:
321 self.L.load('test/blackbox/parsefail1.a4c')
322 self.fail("parsefail1 should not have loaded without errors")
323 except:
324 pass
325
326 def testparsefail2(self):
327 try:
328 self.L.load('test/blackbox/parsefail2.a4c')
329 self.fail("parsefail2 should not have loaded without errors")
330 except:
331 pass
332
333 def testparsefail3(self):
334 try:
335 self.L.load('test/blackbox/parsefail3.a4c')
336 self.fail("parsefail3 should not have loaded without errors")
337 except:
338 pass
339
340 def testparsefail4(self):
341 try:
342 self.L.load('test/blackbox/parsefail4.a4c')
343 self.fail("parsefail4 should not have loaded")
344 except:
345 pass
346
347 def testfail1(self):
348 """Mismatched arg counts check-- tests bbox, not ascend."""
349 self.L.load('test/blackbox/fail1.a4c')
350 try:
351 M = self.L.findType('fail1').getSimulation('sim')
352 self.fail("expected exception was not raised")
353 except RuntimeError,e:
354 print "Caught exception '%s', assumed ok" % e
355
356 def testfail2(self):
357 """Incorrect data arg check -- tests bbox, not ascend"""
358 self.L.load('test/blackbox/fail2.a4c')
359 try:
360 M = self.L.findType('fail2').getSimulation('sim')
361 self.fail("expected exception was not raised")
362 except RuntimeError,e:
363 print "Caught exception '%s', assumed ok (should mention errors during instantiation)" % e
364
365 def testpass1(self):
366 """simple single bbox forward solve"""
367 M = self._run('pass1',filename='test/blackbox/pass.a4c')
368
369 def testpass2(self):
370 """simple single bbox reverse solve"""
371 M = self._run('pass2',filename='test/blackbox/pass.a4c')
372
373 def testpass3(self):
374 """simple double bbox solve"""
375 M = self._run('pass3',filename='test/blackbox/pass3.a4c')
376
377 def testpass4(self):
378 """simple double bbox reverse solve"""
379 M = self._run('pass4',filename='test/blackbox/pass3.a4c')
380
381 def testpass5(self):
382 M = self._run('pass5',filename='test/blackbox/pass5.a4c')
383
384 def testpass6(self):
385 M = self._run('pass6',filename='test/blackbox/pass5.a4c')
386
387 def testpass7(self):
388 M = self._run('pass7',filename='test/blackbox/passmerge.a4c')
389
390 def testpass8(self):
391 M = self._run('pass8',filename='test/blackbox/passmerge.a4c')
392
393 def testpass9(self):
394 M = self._run('pass9',filename='test/blackbox/passmerge.a4c')
395
396 def testpass10(self):
397 M = self._run('pass10',filename='test/blackbox/passmerge.a4c')
398
399 def testpass11(self):
400 M = self._run('pass11',filename='test/blackbox/passmerge.a4c')
401
402 def testpass12(self):
403 M = self._run('pass12',filename='test/blackbox/passmerge.a4c')
404
405 # this test doesn't work: 'system is inconsistent' -- and structurally singular
406 # def testpass13(self):
407 # """cross-merged input/output solve"""
408 # M = self._run('pass13',filename='test/blackbox/passmerge.a4c')
409
410 def testpass14(self):
411 """cross-merged input/output reverse solve"""
412 M = self._run('pass14',filename='test/blackbox/passmerge.a4c')
413
414 def testpass20(self):
415 M = self._run('pass20',filename='test/blackbox/passarray.a4c')
416
417 def testparsefail21(self):
418 """dense array of black boxes wrong syntax"""
419 try:
420 self.L.load('test/blackbox/parsefail21.a4c')
421 self.fail("parsefail21 should not have loaded without errors")
422 except:
423 pass
424
425 def testpass22(self):
426 M = self._run('pass22',filename='test/blackbox/passarray.a4c')
427
428 def testpass23(self):
429 M = self._run('pass23',filename='test/blackbox/passarray.a4c')
430
431 def testpass61(self):
432 M = self._run('pass61',filename='test/blackbox/reinstantiate.a4c')
433
434 def testpass62(self):
435 M = self._run('pass62',filename='test/blackbox/reinstantiate.a4c')
436
437 def testpass64(self):
438 M = self._run('pass64',filename='test/blackbox/reinstantiate.a4c')
439
440 def testpass65(self):
441 M = self._run('pass65',filename='test/blackbox/reinstantiate.a4c')
442
443 def testpass66(self):
444 M = self._run('pass66',filename='test/blackbox/reinstantiate.a4c')
445
446 def testpass67(self):
447 M = self._run('pass67',filename='test/blackbox/reinstantiate.a4c')
448
449 class TestExtFn(AscendSelfTester):
450 def testextfntest(self):
451 M = self._run('extfntest',filename='johnpye/extfn/extfntest.a4c')
452 self.assertAlmostEqual(M.y, 2);
453 self.assertAlmostEqual(M.x, 1);
454 self.assertAlmostEqual(M.y, M.x + 1);
455
456 def testextrelfor(self):
457 M = self._run('extrelfor',filename='johnpye/extfn/extrelfor.a4c')
458
459 ## @TODO fix bug with badly-named bbox rel in a loop (Ben, maybe)
460 # def testextrelforbadnaming(self):
461 # self.L.load('johnpye/extfn/extrelforbadnaming.a4c')
462 # T = self.L.findType('extrelfor')
463 # M = T.getSimulation('sim')
464 # M.solve(ascpy.Solver('QRSlv'),ascpy.SolverReporter())
465 # print "x[1] = %f" % M.x[1]
466 # print "x[2] = %f" % M.x[2]
467 # print "x[3] = %f" % M.x[3]
468 # print "x[4] = %f" % M.x[4]
469 # print "x[5] = %f" % M.x[5]
470 # M.run(T.getMethod('self_test'))
471
472 def testextrelrepeat(self):
473 M = self._run('extrelrepeat',filename='johnpye/extfn/extrelrepeat.a4c')
474
475 #-------------------------------------------------------------------------------
476 # Testing of Sensitivity module
477
478 class TestSensitivity(AscendSelfTester):
479 def test1(self):
480 self.L.load('sensitivity_test.a4c')
481 T = self.L.findType('sensitivity_test')
482 M = T.getSimulation('sim',False)
483 M.run(T.getMethod('on_load'))
484 M.solve(ascpy.Solver('QRSlv'),ascpy.SolverReporter())
485 M.run(T.getMethod('analyse'))
486 M.run(T.getMethod('self_test'))
487
488 # def testall(self):
489 # self.L.load('sensitivity_test.a4c')
490 # T = self.L.findType('sensitivity_test_all')
491 # M = T.getSimulation('sim',False)
492 # M.run(T.getMethod('on_load'))
493 # M.solve(ascpy.Solver('QRSlv'),ascpy.SolverReporter())
494 # M.run(T.getMethod('analyse'))
495 # M.run(T.getMethod('self_test'))
496 # CAUSES CRASH
497
498 #-------------------------------------------------------------------------------
499 # Testing of a ExtPy - external python methods
500
501 class TestExtPy(AscendSelfTester):
502 def test1(self):
503 self.L.load('johnpye/extpy/extpytest.a4c')
504 T = self.L.findType('extpytest')
505 M = T.getSimulation('sim')
506 M.run(T.getMethod('self_test'))
507
508 def test2(self):
509 self.L.load('johnpye/extpy/extpytest.a4c')
510 T = self.L.findType('extpytest')
511 M = T.getSimulation('sim')
512 M.run(T.getMethod('pythonthing'))
513 M.run(T.getMethod('pythonthing'))
514 M.run(T.getMethod('pythonthing'))
515 M.run(T.getMethod('pythonthing'))
516 M.run(T.getMethod('pythonthing'))
517 M.run(T.getMethod('pythonthing'))
518 M.run(T.getMethod('pythonthing'))
519 # causes crash!
520
521 #-------------------------------------------------------------------------------
522 # Testing of saturated steam properties library (iapwssatprops.a4c)
523
524 class TestSteam(AscendSelfTester):
525
526 def testiapwssatprops1(self):
527 M = self._run('testiapwssatprops1',filename='steam/iapwssatprops.a4c')
528 def testiapwssatprops2(self):
529 M = self._run('testiapwssatprops2',filename='steam/iapwssatprops.a4c')
530 def testiapwssatprops3(self):
531 M = self._run('testiapwssatprops3',filename='steam/iapwssatprops.a4c')
532
533 # test the stream model basically works
534 def testsatsteamstream(self):
535 M = self._run('satsteamstream',filename='steam/satsteamstream.a4c')
536
537 # test that we can solve in terms of various (rho,u)
538 def testsatuv(self):
539 self.L.load('steam/iapwssat.a4c')
540 T = self.L.findType('testiapwssatuv')
541 M = T.getSimulation('sim',False)
542 M.run(T.getMethod('on_load'))
543 M.solve(ascpy.Solver('QRSlv'),ascpy.SolverReporter())
544 print "p = %f bar" % M.p.as('bar');
545 print "T = %f C" % (M.T.as('K') - 273.15);
546 print "x = %f" % M.x;
547 M.run(T.getMethod('self_test'))
548 M.run(T.getMethod('values2'))
549 # M.v.setRealValueWithUnits(1.0/450,"m^3/kg");
550 # M.solve(ascpy.Solver('QRSlv'),ascpy.SolverReporter())
551 M.solve(ascpy.Solver('QRSlv'),ascpy.SolverReporter())
552 print "p = %f bar" % M.p.as('bar');
553 print "T = %f C" % (M.T.as('K') - 273.15);
554 print "x = %f" % M.x;
555 M.run(T.getMethod('self_test2'))
556
557
558 ## @TODO fix error capture from bounds checking during initialisation
559 # def testiapwssat1(self):
560 # M = self._run('testiapwssat1',filename='steam/iapwssat.a4c')
561
562 def testdsgsat(self):
563 self.L.load('steam/dsgsat3.a4c')
564 T = self.L.findType('dsgsat3')
565 M = T.getSimulation('sim',False)
566 M.run(T.getMethod('on_load'))
567 M.solve(ascpy.Solver('QRSlv'),ascpy.SolverReporter())
568 self.assertAlmostEqual(M.dTw_dt[2],0.0);
569 M.run(T.getMethod('configure_dynamic'))
570 M.solve(ascpy.Solver('QRSlv'),ascpy.SolverReporter())
571 return M
572
573 def testdsgsatrepeat(self):
574 self.L.load('steam/dsgsat3.a4c')
575 T = self.L.findType('dsgsat3')
576 M = T.getSimulation('sim',False)
577 M.run(T.getMethod('on_load'))
578 M.solve(ascpy.Solver('QRSlv'),ascpy.SolverReporter())
579 M.run(T.getMethod('on_load'))
580 M.solve(ascpy.Solver('QRSlv'),ascpy.SolverReporter())
581 M.run(T.getMethod('on_load'))
582 M.solve(ascpy.Solver('QRSlv'),ascpy.SolverReporter())
583
584 def testvary(self):
585 self.L.load('steam/dsgsat3.a4c')
586 T = self.L.findType('dsgsat3')
587 M = T.getSimulation('sim',False)
588 M.run(T.getMethod('on_load'))
589 M.solve(ascpy.Solver('QRSlv'),ascpy.SolverReporter())
590 print "----- setting qdot_s -----"
591 M.qdot_s.setRealValueWithUnits(1000,"W/m")
592 M.solve(ascpy.Solver('QRSlv'),ascpy.SolverReporter())
593 print "----- setting qdot_s -----"
594 M.qdot_s.setRealValueWithUnits(2000,"W/m")
595 M.solve(ascpy.Solver('QRSlv'),ascpy.SolverReporter())
596
597 def teststeadylsode(self):
598 "test that steady conditions are stable with LSODE"
599 M = self.testdsgsat()
600 #M.qdot_s.setRealValueWithUnits(1000,"W/m")
601 M.solve(ascpy.Solver('QRSlv'),ascpy.SolverReporter())
602 #M.setParameter('
603 I = ascpy.Integrator(M)
604 I.setEngine('LSODE')
605 I.setReporter(ascpy.IntegratorReporterConsole(I))
606 I.setLinearTimesteps(ascpy.Units("s"), 0, 3600, 10)
607 I.analyse()
608 I.solve()
609
610 # DOESN'T CONVERGE
611 # def testpeturblsode(self):
612 # "test that steady conditions are stable with LSODE"
613 # M = self.testdsgsat()
614 # # here is the peturbation...
615 # M.qdot_s.setRealValueWithUnits(1000,"W/m")
616 # M.solve(ascpy.Solver('QRSlv'),ascpy.SolverReporter())
617 # I = ascpy.Integrator(M)
618 # I.setEngine('LSODE')
619 # I.setReporter(ascpy.IntegratorReporterConsole(I))
620 # I.setLinearTimesteps(ascpy.Units("s"), 0, 5, 1)
621 # I.analyse()
622 # I.solve()
623
624 def teststeadyida(self):
625 M = self.testdsgsat()
626 self.assertAlmostEqual(M.dTw_dt[2],0.0)
627 Tw1 = float(M.T_w[2])
628 T = self.L.findType('dsgsat3')
629 M.run(T.getMethod('free_states'))
630 I = ascpy.Integrator(M)
631 I.setEngine('IDA')
632 I.setParameter('linsolver','DENSE')
633 I.setParameter('safeeval',True)
634 I.setParameter('rtol',1e-4)
635 I.setParameter('atolvect',False)
636 I.setParameter('atol',1e-4)
637 I.setParameter('maxord',3)
638 I.setInitialSubStep(0.001)
639 I.setReporter(ascpy.IntegratorReporterConsole(I))
640 I.setLinearTimesteps(ascpy.Units("s"), 0, 3600, 10)
641 I.analyse()
642 I.solve()
643 self.assertAlmostEqual(float(M.T_w[2]),Tw1)
644 M.qdot_s.setRealValueWithUnits(1000,"W/m")
645 self.assertAlmostEqual(M.qdot_s.as("W/m"),1000)
646 M.solve(ascpy.Solver('QRSlv'),ascpy.SolverReporter())
647 print "dTw/dt = %f" % M.dTw_dt[2]
648 self.assertNotAlmostEqual(M.dTw_dt[2],0.0)
649
650 def testpeturbida(self):
651 M = self.testdsgsat()
652 self.assertAlmostEqual(M.dTw_dt[2],0.0)
653 T = self.L.findType('dsgsat3')
654 M.run(T.getMethod('free_states'))
655 # here is the peturbation...
656 M.qdot_s.setRealValueWithUnits(6000,"W/m")
657 # IDA has its own initial conditions solver, so no need to call QRSlv here
658 I = ascpy.Integrator(M)
659 I.setEngine('IDA')
660 I.setParameter('linsolver','DENSE')
661 I.setParameter('safeeval',True)
662 I.setParameter('rtol',1e-4)
663 I.setParameter('atolvect',False)
664 I.setParameter('atol',1e-4)
665 I.setParameter('maxord',3)
666 I.setInitialSubStep(0.001)
667 I.setReporter(ascpy.IntegratorReporterConsole(I))
668 I.setLogTimesteps(ascpy.Units("s"), 0.001, 3600, 40)
669 I.analyse()
670 I.solve()
671
672 #-------------------------------------------------------------------------------
673 # Testing of freesteam external steam properties functions
674
675 with_freesteam = True
676 try:
677 # we assume that if the freesteam python module is installed, the ASCEND
678 # external library will also be.
679 import freesteam
680 have_freesteam = True
681 except ImportError,e:
682 have_freesteam = False
683
684 if with_freesteam and have_freesteam:
685 class TestFreesteam(AscendSelfTester):
686 # def testfreesteamtest(self):
687 # """run the self-test cases bundled with freesteam"""
688 # self._run('testfreesteam',filename='testfreesteam.a4c')
689
690 def testload(self):
691 """check that we can load 'thermalequilibrium2' (IMPORT "freesteam", etc)"""
692 self.L.load('johnpye/thermalequilibrium2.a4c')
693
694 def testinstantiate(self):
695 """load an instantiate 'thermalequilibrium2'"""
696 self.testload()
697 M = self.L.findType('thermalequilibrium2').getSimulation('sim')
698 return M
699
700 def testintegrate(self):
701 """integrate transfer of heat from one mass of water/steam to another
702 according to Newton's law of cooling"""
703 M = self.testinstantiate()
704 M.setSolver(ascpy.Solver("QRSlv"))
705 I = ascpy.Integrator(M)
706 I.setEngine('LSODE')
707 I.setReporter(ascpy.IntegratorReporterConsole(I))
708 I.setLinearTimesteps(ascpy.Units("s"), 0, 3000, 30)
709 I.setMinSubStep(0.01)
710 I.setInitialSubStep(1)
711 I.analyse()
712 print "Number of vars = %d" % I.getNumVars()
713 assert I.getNumVars()==2
714 I.solve()
715 assert I.getNumObservedVars() == 3;
716 print "S[1].T = %f K" % M.S[1].T
717 print "S[2].T = %f K" % M.S[2].T
718 print "Q = %f W" % M.Q
719 self.assertAlmostEqual(float(M.S[1].T),506.77225109,4);
720 self.assertAlmostEqual(float(M.S[2].T),511.605173967,5);
721 self.assertAlmostEqual(float(M.Q),-48.32922877329,3);
722 self.assertAlmostEqual(float(M.t),3000);
723 print "Note that the above values have not been verified analytically"
724
725 def testcollapsingcan2(self):
726 """ solve the collapsing can model using IAPWS-IF97 steam props """
727 M = self._run("collapsingcan2",filename="collapsingcan2.a4c");
728
729 #-------------------------------------------------------------------------------
730 # Testing of IDA's analysis module
731
732 class TestIDA(Ascend):
733 def _run(self,filen,modeln=""):
734 self.L.load('test/ida/%s.a4c' % filen)
735 T = self.L.findType('%s%s' % (filen,modeln))
736 M = T.getSimulation('sim')
737 M.build()
738 I = ascpy.Integrator(M)
739 I.setEngine('IDA')
740 I.analyse()
741 return M;
742
743 def _runfail(self,filen,n,msg="failed"):
744 try:
745 self._run(filen,'fail%d' % n)
746 except Exception,e:
747 print "(EXPECTED) ERROR: %s" % e
748 return
749 self.fail(msg)
750
751 def testsinglederiv(self):
752 self._run('singlederiv')
753
754 def testsinglederivfail1(self):
755 self._runfail('singlederiv',1
756 ,"t.ode_id=-1 did not trigger error")
757
758 def testsinglederivfail2(self):
759 self._runfail('singlederiv',2
760 ,"dy_dt.ode_id=2 did not trigger error")
761
762 def testsinglederivfail3(self):
763 self._runfail('singlederiv',3
764 ,"dy_dt.ode_type=3 did not trigger error")
765
766 def testsinglederivfail4(self):
767 self._runfail('singlederiv',4
768 ,"duplicate ode_type=1 did not trigger error")
769
770 def testsinglederivfail5(self):
771 self._runfail('singlederiv',5
772 ,"duplicate ode_type=1 did not trigger error")
773
774 def testsinglederivfail6(self):
775 self._runfail('singlederiv',6
776 ,"duplicate ode_type=1 did not trigger error")
777
778 def testtwoderiv(self):
779 self._run('twoderiv')
780
781 def testtwoderivfail1(self):
782 self._runfail('twoderiv',1)
783
784 def testtwoderivfail2(self):
785 self._runfail('twoderiv',2)
786
787 def testtwoderivfail3(self):
788 self._runfail('twoderiv',3)
789 def testtwoderivfail4(self):
790 self._runfail('twoderiv',4)
791 def testtwoderivfail5(self):
792 self._runfail('twoderiv',5)
793
794 def testnoderivs(self):
795 self._runfail('noderivs',1)
796
797 def testnoindeps(self):
798 self._runfail('indeps',1)
799
800 def testtwoindeps(self):
801 self._runfail('indeps',2)
802
803 def testfixedvars(self):
804 self._run('fixedvars')
805
806 def testfixedvars1(self):
807 self._run('fixedvars',1)
808
809 def testfixedvars2(self):
810 self._run('fixedvars',2)
811
812 def testfixedvars3(self):
813 self._run('fixedvars',3)
814
815 def testincidence(self):
816 self._run('incidence')
817
818 def testincidence1(self):
819 self._run('incidence',1)
820 def testincidence2(self):
821 self._run('incidence',2)
822 def testincidence3(self):
823 M = self._run('incidence',3)
824
825 def testincidence4(self):
826 self._run('incidence',4)
827 def testincidencefail5(self):
828 self._runfail('incidence',5)
829
830 # doesn't work yet:
831 # def testincidence5(self):
832 # self._run('incidence',5)
833
834
835 #-------------------------------------------------------------------------------
836 # Testing of IDA models using DENSE linear solver
837
838 class TestIDADENSE(Ascend):
839 """IDA DAE integrator, DENSE linear solver"""
840
841 def testlotka(self):
842 self.L.load('johnpye/lotka.a4c')
843 M = self.L.findType('lotka').getSimulation('sim')
844 M.setSolver(ascpy.Solver("QRSlv"))
845 I = ascpy.Integrator(M)
846 I.setEngine('IDA')
847 I.setReporter(ascpy.IntegratorReporterConsole(I))
848 I.setLinearTimesteps(ascpy.Units("s"), 0, 200, 5);
849 I.setParameter('linsolver','DENSE')
850 I.setParameter('rtol',1e-8);
851 I.analyse()
852 assert I.getNumVars()==2
853 assert abs(M.R - 1000) < 1e-300
854 I.solve()
855 assert I.getNumObservedVars() == 3
856 assert abs(M.R - 832) < 1.0
857 assert abs(M.F - 21.36) < 0.1
858
859 def testdenx(self):
860 print "-----------------------------====="
861 self.L.load('johnpye/idadenx.a4c')
862 M = self.L.findType('idadenx').getSimulation('sim')
863 M.setSolver(ascpy.Solver("QRSlv"))
864 I = ascpy.Integrator(M)
865 I.setEngine('IDA')
866 I.setParameter('calcic','YA_YDP')
867 I.setParameter('linsolver','DENSE')
868 I.setParameter('safeeval',True)
869 I.setReporter(ascpy.IntegratorReporterConsole(I))
870 I.setLogTimesteps(ascpy.Units("s"), 0.4, 4e10, 11)
871 I.setMaxSubStep(0);
872 I.setInitialSubStep(0)
873 I.setMaxSubSteps(0)
874 I.setParameter('autodiff',True)
875 I.analyse()
876 I.solve()
877 assert abs(float(M.y1) - 5.1091e-08) < 2e-9
878 assert abs(float(M.y2) - 2.0437e-13) < 2e-14
879 assert abs(float(M.y3) - 1.0) < 1e-5
880
881 def testhires(self):
882 self.L.load('test/hires.a4c')
883 T = self.L.findType('hires')
884 M = T.getSimulation('sim')
885 M.setSolver(ascpy.Solver('QRSlv'))
886 I = ascpy.Integrator(M)
887 I.setEngine('IDA')
888 I.setParameter('linsolver','DENSE')
889 I.setParameter('rtol',1.1e-15)
890 I.setParameter('atolvect',0)
891 I.setParameter('atol',1.1e-15)
892 I.setReporter(ascpy.IntegratorReporterConsole(I))
893 I.setLogTimesteps(ascpy.Units(""), 1, 321.8122, 5)
894 I.setInitialSubStep(1e-5)
895 I.setMaxSubSteps(10000)
896 I.analyse()
897 I.solve()
898 for i in range(8):
899 print "y[%d] = %.20g" % (i+1, M.y[i+1])
900 M.run(T.getMethod('self_test'))
901
902 def testchemakzo(self):
903 self.L.load('test/chemakzo.a4c')
904 T = self.L.findType('chemakzo')
905 M = T.getSimulation('sim')
906 M.setSolver(ascpy.Solver('QRSlv'))
907 I = ascpy.Integrator(M)
908 I.setEngine('IDA')
909 I.setParameter('linsolver','DENSE')
910 I.setParameter('rtol',1e-15)
911 I.setParameter('atolvect',0)
912 I.setParameter('atol',1e-15)
913 I.setReporter(ascpy.IntegratorReporterConsole(I))
914 I.setLinearTimesteps(ascpy.Units("s"), 1, 180, 5)
915 I.setInitialSubStep(1e-13)
916 I.setMaxSubSteps(10000)
917 I.analyse()
918 I.solve()
919 for i in range(6):
920 print "y[%d] = %.20g" % (i+1, M.y[i+1])
921 M.run(T.getMethod('self_test'))
922
923 def testtransamp(self):
924 self.L.load('test/transamp.a4c')
925 T = self.L.findType('transamp')
926 M = T.getSimulation('sim')
927 M.setSolver(ascpy.Solver('QRSlv'))
928 I = ascpy.Integrator(M)
929 I.setEngine('IDA')
930 I.setParameter('linsolver','DENSE')
931 I.setParameter('rtol',1e-7)
932 I.setParameter('atolvect',0)
933 I.setParameter('atol',1e-7)
934 I.setReporter(ascpy.IntegratorReporterConsole(I))
935 I.setLinearTimesteps(ascpy.Units("s"), 0.05, 0.2, 20)
936 I.setInitialSubStep(0.00001)
937 I.setMaxSubSteps(10000)
938 I.analyse()
939 I.solve()
940 for i in range(6):
941 print "y[%d] = %.20g" % (i+1, M.y[i+1])
942 M.run(T.getMethod('self_test'))
943
944 # MODEL FAILS ANALYSIS: we need to add support for non-incident differential vars
945 # def testpollution(self):
946 # self.L.load('test/pollution.a4c')
947 # T = self.L.findType('pollution')
948 # M = T.getSimulation('sim')
949 # M.setSolver(ascpy.Solver('QRSlv'))
950 # I = ascpy.Integrator(M)
951 # I.setEngine('IDA')
952 # I.setParameter('linsolver','DENSE')
953 # I.setParameter('rtol',1.1e-15)
954 # I.setParameter('atolvect',0)
955 # I.setParameter('atol',1.1e-15)
956 # I.setReporter(ascpy.IntegratorReporterConsole(I))
957 # I.setLogTimesteps(ascpy.Units("s"), 1, 60, 5)
958 # I.setInitialSubStep(1e-5)
959 # I.setMaxSubSteps(10000)
960 # I.analyse()
961 # I.solve()
962 # for i in range(20):
963 # print "y[%d] = %.20g" % (i+1, M.y[i+1])
964 # M.run(T.getMethod('self_test'))
965
966 ## @TODO fails during IDACalcIC (model too big?)
967 # def testkryx(self):
968 # self.L.load('johnpye/idakryx.a4c')
969 # ascpy.getCompiler().setUseRelationSharing(False)
970 # M = self.L.findType('idakryx').getSimulation('sim')
971 # M.setSolver(ascpy.Solver('QRSlv'))
972 # M.build()
973 # I = ascpy.Integrator(M)
974 # I.setEngine('IDA')
975 # I.setReporter(ascpy.IntegratorReporterConsole(I))
976 # I.setParameter('linsolver','DENSE')
977 # I.setParameter('maxl',8)
978 # I.setParameter('gsmodified',False)
979 # I.setParameter('autodiff',True)
980 # I.setParameter('rtol',0)
981 # I.setParameter('atol',1e-3);
982 # I.setParameter('atolvect',False)
983 # I.setParameter('calcic','YA_YDP')
984 # I.analyse()
985 # I.setLogTimesteps(ascpy.Units("s"), 0.01, 10.24, 11)
986 # I.solve()
987 # assert abs(M.u[2][2].getValue()) < 1e-5
988
989 #-------------------------------------------------------------------------------
990 # Testing of IDA models using SPGMR linear solver (Krylov)
991
992 # these tests are disabled until SPGMR preconditioning has been implemented
993 class TestIDASPGMR:#(Ascend):
994 def testlotka(self):
995 self.L.load('johnpye/lotka.a4c')
996 M = self.L.findType('lotka').getSimulation('sim')
997 M.setSolver(ascpy.Solver("QRSlv"))
998 I = ascpy.Integrator(M)
999 I.setEngine('IDA')
1000 I.setReporter(ascpy.IntegratorReporterConsole(I))
1001 I.setLinearTimesteps(ascpy.Units("s"), 0, 200, 5)
1002 I.setParameter('rtol',1e-8)
1003 I.analyse()
1004 assert I.getNumVars()==2
1005 assert abs(M.R - 1000) < 1e-300
1006 I.solve()
1007 assert I.getNumObservedVars() == 3
1008 assert abs(M.R - 832) < 1.0
1009 assert abs(M.F - 21.36) < 0.1
1010
1011
1012 def testkryx(self):
1013 self.L.load('johnpye/idakryx.a4c')
1014 M = self.L.findType('idakryx').getSimulation('sim')
1015 M.build()
1016 I = ascpy.Integrator(M)
1017 I.setEngine('IDA')
1018 I.setReporter(ascpy.IntegratorReporterConsole(I))
1019 I.setParameter('linsolver','SPGMR')
1020 I.setParameter('prec','JACOBI')
1021 I.setParameter('maxl',8)
1022 I.setParameter('gsmodified',False)
1023 I.setParameter('autodiff',True)
1024 I.setParameter('gsmodified',True)
1025 I.setParameter('rtol',0)
1026 I.setParameter('atol',1e-3);
1027 I.setParameter('atolvect',False)
1028 I.setParameter('calcic','Y')
1029 I.analyse()
1030 I.setLogTimesteps(ascpy.Units("s"), 0.01, 10.24, 10);
1031 print M.udot[1][3]
1032 I.solve()
1033 assert 0
1034
1035 def testzill(self):
1036 self.L.load('johnpye/zill.a4c')
1037 T = self.L.findType('zill')
1038 M = T.getSimulation('sim')
1039 M.setSolver(ascpy.Solver('QRSlv'))
1040 I = ascpy.Integrator(M)
1041 I.setEngine('IDA')
1042 I.setParameter('safeeval',False)
1043 I.setMinSubStep(1e-7)
1044 I.setMaxSubStep(0.001)
1045 I.setMaxSubSteps(10000)
1046 I.setReporter(ascpy.IntegratorReporterConsole(I))
1047 I.setLinearTimesteps(ascpy.Units(), 1.0, 1.5, 5)
1048 I.analyse()
1049 I.solve()
1050 M.run(T.getMethod('self_test'))
1051
1052 def testdenxSPGMR(self):
1053 self.L.load('johnpye/idadenx.a4c')
1054 M = self.L.findType('idadenx').getSimulation('sim')
1055 M.setSolver(ascpy.Solver('QRSlv'))
1056 I = ascpy.Integrator(M)
1057 I.setEngine('IDA')
1058 I.setReporter(ascpy.IntegratorReporterConsole(I))
1059 I.setLogTimesteps(ascpy.Units("s"), 0.4, 4e10, 11)
1060 I.setMaxSubStep(0);
1061 I.setInitialSubStep(0);
1062 I.setMaxSubSteps(0);
1063 I.setParameter('autodiff',True)
1064 I.setParameter('linsolver','SPGMR')
1065 I.setParameter('gsmodified',False)
1066 I.setParameter('maxncf',10)
1067 I.analyse()
1068 I.solve()
1069 assert abs(float(M.y1) - 5.1091e-08) < 1e-10
1070 assert abs(float(M.y2) - 2.0437e-13) < 1e-15
1071 assert abs(float(M.y3) - 1.0) < 1e-5
1072
1073 # move code above down here if you want to temporarily avoid testing it
1074 class NotToBeTested:
1075 def nothing(self):
1076 pass
1077
1078 def testnewton(self):
1079 sys.stderr.write("STARTING TESTNEWTON\n")
1080 self.L.load('johnpye/newton.a4c')
1081 T = self.L.findType('newton')
1082 M = T.getSimulation('sim')
1083 M.solve(ascpy.Solver("QRSlv"),ascpy.SolverReporter())
1084 I = ascpy.Integrator(M)
1085 I.setEngine('IDA')
1086 I.setParameter('linsolver','DENSE')
1087 I.setParameter('safeeval',True)
1088 I.setParameter('rtol',1e-8)
1089 I.setMaxSubStep(0.001)
1090 I.setMaxSubSteps(10000)
1091
1092 I.setReporter(ascpy.IntegratorReporterConsole(I))
1093 I.setLinearTimesteps(ascpy.Units("s"), 0, 2*float(M.v)/float(M.g), 2)
1094 I.analyse()
1095 I.solve()
1096 print "At end of simulation,"
1097 print "x = %f" % M.x
1098 print "v = %f" % M.v
1099 M.run(T.getMethod('self_test'))
1100
1101 if __name__=='__main__':
1102 # a whole bag of tricks to make sure we get the necessary dirs in our ascend, python and ld path vars
1103 restart = 0
1104
1105 if platform.system()=="Windows":
1106 LD_LIBRARY_PATH="PATH"
1107 SEP = ";"
1108 else:
1109 LD_LIBRARY_PATH="LD_LIBRARY_PATH"
1110 SEP = ":"
1111
1112 freesteamdir = os.path.expanduser("~/freesteam/ascend")
1113 modeldirs = [os.path.abspath(os.path.join(sys.path[0],"models")),os.path.abspath(freesteamdir)]
1114 if not os.environ.get('ASCENDLIBRARY'):
1115 os.environ['ASCENDLIBRARY'] = SEP.join(modeldirs)
1116 restart = 1
1117 else:
1118 envmodelsdir = [os.path.abspath(i) for i in os.environ['ASCENDLIBRARY'].split(SEP)]
1119 for l in modeldirs:
1120 if l in envmodelsdir[len(modeldirs):]:
1121 envmodelsdir.remove(l)
1122 restart = 1
1123 for l in modeldirs:
1124 if l not in envmodelsdir:
1125 envmodelsdir.insert(0,l)
1126 restart = 1
1127 os.environ['ASCENDLIBRARY'] = SEP.join(envmodelsdir)
1128
1129 libdirs = ["pygtk","."]
1130 libdirs = [os.path.normpath(os.path.join(sys.path[0],l)) for l in libdirs]
1131 if not os.environ.get(LD_LIBRARY_PATH):
1132 os.environ[LD_LIBRARY_PATH]=SEP.join(libdirs)
1133 restart = 1
1134 else:
1135 envlibdirs = [os.path.normpath(i) for i in os.environ[LD_LIBRARY_PATH].split(SEP)]
1136 for l in libdirs:
1137 if l in envlibdirs[len(libdirs):]:
1138 envlibdirs.remove(l)
1139 restart = 1
1140 for l in libdirs:
1141 if l not in envlibdirs:
1142 envlibdirs.insert(0,l)
1143 restart = 1
1144 os.environ[LD_LIBRARY_PATH] = SEP.join(envlibdirs)
1145
1146 pypath = os.path.normpath(os.path.join(sys.path[0],"pygtk"))
1147 if not os.environ.get('PYTHONPATH'):
1148 os.environ['PYTHONPATH']=pypath
1149 else:
1150 envpypath = os.environ['PYTHONPATH'].split(SEP)
1151 if pypath not in envpypath:
1152 envpypath.insert(0,pypath)
1153 os.environ['PYTHONPATH']=SEP.join(envpypath)
1154 restart = 1
1155
1156 if restart:
1157 script = os.path.join(sys.path[0],"test.py")
1158 print "Restarting with..."
1159 print " export LD_LIBRARY_PATH=%s" % os.environ.get(LD_LIBRARY_PATH)
1160 print " export PYTHONPATH=%s" % os.environ.get('PYTHONPATH')
1161 print " export ASCENDLIBRARY=%s" % os.environ.get('ASCENDLIBRARY')
1162
1163 os.execvp("python",[script] + sys.argv)
1164
1165 import ascpy
1166
1167 try:
1168 import cunit
1169 except:
1170 pass
1171
1172 atexit.register(ascpy.shutdown)
1173 #suite = unittest.TestSuite()
1174 #suite = unittest.defaultTestLoader.loadTestsFromName('__main__')
1175 #unittest.TextTestRunner(verbosity=2).run(suite)
1176 unittest.main()

Properties

Name Value
svn:executable *

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