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

Contents of /trunk/test.py

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1837 - (show annotations) (download) (as text)
Thu Aug 28 07:08:32 2008 UTC (15 years, 1 month ago) by jpye
File MIME type: text/x-python
File size: 45791 byte(s)
Changed Instance::as  to Instance::to in PyGTK interface, for compatibility
with Python 2.5 ('as' is a reserved word).
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 try:
32 import dl
33 _dlflags = dl.RTLD_GLOBAL|dl.RTLD_NOW
34 except:
35 # On platforms that unilaterally refuse to provide the 'dl' module
36 # we'll just set the value and see if it works.
37 print "Note: python 'dl' module not available on this system, guessing value of RTLD_* flags"
38 _dlflags = 258
39
40 sys.setdlopenflags(_dlflags)
41
42 class Ascend(unittest.TestCase):
43
44 def setUp(self):
45 import ascpy
46 self.L = ascpy.Library()
47
48 def tearDown(self):
49 self.L.clear()
50 del self.L
51
52 class AscendSelfTester(Ascend):
53
54 def _run(self,modelname,solvername="QRSlv",filename=None,parameters={}):
55 if filename==None:
56 filename = 'johnpye/%s.a4c' % modelname
57 self.L.load(filename)
58 T = self.L.findType(modelname)
59 M = T.getSimulation('sim')
60 M.setSolver(ascpy.Solver(solvername))
61 for k,v in parameters.iteritems():
62 M.setParameter(k,v)
63 M.solve(ascpy.Solver(solvername),ascpy.SolverReporter())
64 M.run(T.getMethod('self_test'))
65 return M
66
67 class TestCompiler(Ascend):
68
69 def _run(self,filen,modeln=""):
70 self.L.load('test/compiler/%s.a4c' % filen)
71 T = self.L.findType('%s%s' % (filen,modeln))
72 M = T.getSimulation('sim')
73 M.build()
74
75 def _runfail(self,filen,n,msg="failed"):
76 try:
77 self._run(filen,'fail%d' % n)
78 except Exception,e:
79 print "(EXPECTED) ERROR: %s" % e
80 return
81 self.fail(msg)
82
83
84 def testloading(self):
85 """library startup"""
86 pass
87
88 def testsystema4l(self):
89 """loading system.a4l?"""
90 self.L.load('system.a4l')
91
92 def testatomsa4l(self):
93 """loading atoms.a4l?"""
94 self.L.load('atoms.a4l')
95
96 def testmissingreq(self):
97 """flagging a missing REQUIRE"""
98 self._runfail('missingreq',1)
99
100 def testmissingsubreq(self):
101 """flagging a subsidiary missing REQUIRE"""
102 self._runfail('missingreq',1)
103
104 def defaultmethodstest(self,modelname):
105 self.L.load("test/defaultmethods.a4c")
106 T = self.L.findType(modelname)
107 M = T.getSimulation('sim')
108 M.run(T.getMethod('on_load'))
109 M.run(T.getMethod('self_test'))
110 return M
111
112 def testdefault1(self):
113 self.defaultmethodstest('testdefault1')
114
115 def testdefault2(self):
116 self.defaultmethodstest('testdefault2')
117
118 def testdefault3(self):
119 self.defaultmethodstest('testdefault3')
120
121 def testdefault4(self):
122 self.defaultmethodstest('testdefault4')
123
124 def testdefault5(self):
125 self.defaultmethodstest('testdefault5')
126
127 def testdefault6(self):
128 self.defaultmethodstest('testdefault6')
129
130 def testdefault7(self):
131 self.defaultmethodstest('testdefault7')
132
133 def testdefault8(self):
134 self.defaultmethodstest('testdefault8')
135
136 def testdefault9(self):
137 self.defaultmethodstest('testdefault9')
138
139 def testdefault10(self):
140 self.defaultmethodstest('testdefault10')
141
142 def testdefault11(self):
143 self.defaultmethodstest('testdefault11')
144
145 def testdefault12(self):
146 self.defaultmethodstest('testdefault12')
147
148 def testdefault13(self):
149 self.defaultmethodstest('testdefault13')
150
151 def testdefault14(self):
152 self.defaultmethodstest('testdefault14')
153
154 def testdefault15(self):
155 self.defaultmethodstest('testdefault15')
156
157 def testdefault16(self):
158 self.defaultmethodstest('testdefault16')
159
160 def testdefault17(self):
161 self.defaultmethodstest('testdefault17')
162 def testdefault18(self):
163 self.defaultmethodstest('testdefault18')
164
165 def testdefault19(self):
166 self.defaultmethodstest('testdefault19')
167
168 #def testdefault19fail(self):
169 # self.defaultmethodstest('testdefault19fail')
170
171 def testdefault20(self):
172 self.defaultmethodstest('testdefault20')
173
174 #def testdefault20fail(self):
175 # self.defaultmethodstest('testdefault20fail')
176
177 class TestSystem(AscendSelfTester):
178
179 def testwritegraph(self):
180 M = self._run('testlog10')
181 M.run(self.L.findType('testlog10').getMethod('on_load'))
182 if platform.system!="Windows":
183 f = file('temp.png','wb')
184 M.write(f,"dot")
185 f.close()
186 else:
187 self.fail("not implemented on windows")
188
189 class TestSolver(AscendSelfTester):
190
191 def testlog10(self):
192 M = self._run('testlog10')
193
194 def testrootsofpoly(self):
195 self._run('roots_of_poly',filename="roots_of_poly.a4c")
196
197 def testcollapsingcan(self):
198 self._run('collapsingcan',filename="collapsingcan.a4c")
199
200 def testdistancecalc(self):
201 self._run('distance_calc',filename="distance_calc.a4c")
202
203 def testconopt(self):
204 self._run('conopttest',"CONOPT",filename="conopttest.a4c")
205
206 def testcmslv2(self):
207 self._run('testcmslv2',"CMSlv")
208
209 def testsunpos1(self):
210 self._run('example_1_6_1',"QRSlv","johnpye/sunpos.a4c")
211
212 def testsunpos2(self):
213 self._run('example_1_6_2',"QRSlv","johnpye/sunpos.a4c")
214
215 def testsunpos3(self):
216 self._run('example_1_7_1',"QRSlv","johnpye/sunpos.a4c")
217
218 def testsunpos4(self):
219 self._run('example_1_7_2',"QRSlv","johnpye/sunpos.a4c")
220
221 def testsunpos5(self):
222 self._run('example_1_7_3',"QRSlv","johnpye/sunpos.a4c")
223
224 def testsunpos6(self):
225 self._run('example_1_8_1',"QRSlv","johnpye/sunpos.a4c")
226
227 def testinstanceas(self):
228 M = self._run('example_1_6_1',"QRSlv","johnpye/sunpos.a4c")
229 self.assertAlmostEqual( float(M.t_solar), M.t_solar.to("s"))
230 self.assertAlmostEqual( float(M.t_solar)/3600, M.t_solar.to("h"))
231
232 def testrelinclude(self):
233 self.L.load('test/relinclude.a4c')
234 T = self.L.findType('relinclude')
235 M = T.getSimulation('sim')
236 M.eq1.setIncluded(True)
237 M.eq2.setIncluded(False)
238 M.eq3.setIncluded(False)
239 M.solve(ascpy.Solver('QRSlv'),ascpy.SolverReporter())
240 self.assertAlmostEqual( float(M.z), 2.0)
241 M.eq1.setIncluded(False)
242 M.eq2.setIncluded(True)
243 M.eq3.setIncluded(False)
244 M.solve(ascpy.Solver('QRSlv'),ascpy.SolverReporter())
245 self.assertAlmostEqual( float(M.z), 4.0)
246 M.eq1.setIncluded(False)
247 M.eq2.setIncluded(False)
248 M.eq3.setIncluded(True)
249 M.solve(ascpy.Solver('QRSlv'),ascpy.SolverReporter())
250 self.assertAlmostEqual( float(M.z), 4.61043629206)
251
252
253 class TestLRSlv(AscendSelfTester):
254 def testonerel(self):
255 self._run('onerel',"LRSlv","test/lrslv/onerel.a4c")
256
257 # need to migrate to 'FIX boolvar', currently not supported...
258 # def testonerel(self):
259 # self._run('onerel',"LRSlv","test/lrslv/onerel.a4c")
260
261 def testsequencecrash(self):
262 try:
263 self._run('onerel',"LRSlv","test/lrslv/sequencecrash.a4c")
264 except:
265 pass
266 # it just has to not crash, that's all
267
268 def testsequence(self):
269 self._run('onerel',"LRSlv","test/lrslv/sequence.a4c")
270
271
272 class TestCMSlv(AscendSelfTester):
273 def testsonic(self):
274 M = self._run('sonic',"CMSlv","sonic.a4c")
275 assert(M.sonic_flow.getBoolValue())
276
277 # other side of boundary...
278 M.D.setRealValueWithUnits(4.,"cm")
279 T = self.L.findType('sonic')
280 M.solve(ascpy.Solver('CMSlv'),ascpy.SolverReporter())
281 M.run(T.getMethod('self_test'))
282 assert(not M.sonic_flow.getBoolValue())
283
284 def testheatex(self):
285 self._run('heatex',"CMSlv","heatex.a4c")
286 def testphaseeq(self):
287 self._run('phaseq',"CMSlv","phaseq.a4c")
288 def testpipeline(self):
289 self._run('pipeline',"CMSlv","pipeline.a4c"
290 ,{'infinity':3.2e9}
291 )
292 def testrachford(self):
293 self._run('rachford',"CMSlv","rachford.a4c")
294 def testlinmassbal(self):
295 self._run('linmassbal',"CMSlv","linmassbal.a4c")
296
297
298 class TestMatrix(AscendSelfTester):
299 def testlog10(self):
300 M = self._run('testlog10')
301 print "FETCHING MATRIX................."
302 X = M.getMatrix()
303 # this stuff crashes Windows because the FILE* structure used by Python is not the same
304 # as used by MinGW...
305 #print "GOT MATRIX"
306 #sys.stderr.flush()
307 #sys.stdout.flush()
308 #F = os.tmpfile()
309 #X.write(F.fileno,"mmio")
310 #F.seek(0)
311 #print F.read()
312
313 class TestIntegrator(Ascend):
314
315 def testListIntegrators(self):
316 I = ascpy.Integrator.getEngines()
317 s1 = sorted([str(i) for i in I])
318 assert 'IDA' in s1
319 assert 'LSODE' in s1
320
321 # this routine is reused by both testIDA and testLSODE
322 def _testIntegrator(self,integratorname):
323 self.L.load('johnpye/shm.a4c')
324 M = self.L.findType('shm').getSimulation('sim')
325 M.setSolver(ascpy.Solver('QRSlv'))
326 P = M.getParameters()
327 M.setParameter('feastol',1e-12)
328 print M.getChildren()
329 assert float(M.x) == 10.0
330 assert float(M.v) == 0.0
331 t_end = math.pi
332
333 I = ascpy.Integrator(M)
334 I.setReporter(ascpy.IntegratorReporterNull(I))
335 I.setEngine(integratorname);
336 I.setLinearTimesteps(ascpy.Units("s"), 0.0, t_end, 100);
337 I.setMinSubStep(0.0001); # these limits are required by IDA at present (numeric diff)
338 I.setMaxSubStep(0.1);
339 I.setInitialSubStep(0.001);
340 I.setMaxSubSteps(200);
341 if(integratorname=='IDA'):
342 I.setParameter('autodiff',False)
343 for p in M.getParameters():
344 print p.getName(),"=",p.getValue()
345 I.analyse();
346 I.solve();
347 print "At end of simulation,"
348 print "x = %f" % M.x
349 print "v = %f" % M.v
350 assert abs(float(M.x) + 10) < 1e-2
351 assert abs(float(M.v)) < 1e-2
352 assert I.getNumObservedVars() == 3
353
354 def testInvalidIntegrator(self):
355 self.L.load('johnpye/shm.a4c')
356 M = self.L.findType('shm').getSimulation('sim')
357 M.setSolver(ascpy.Solver('QRSlv'))
358 I = ascpy.Integrator(M)
359 try:
360 I.setEngine('___NONEXISTENT____')
361 except IndexError:
362 return
363 self.fail("setEngine did not raise error!")
364
365 def testLSODE(self):
366 self._testIntegrator('LSODE')
367
368 def testIDA(self):
369 self._testIntegrator('IDA')
370
371 def testparameters(self):
372 self.L.load('johnpye/shm.a4c')
373 M = self.L.findType('shm').getSimulation('sim')
374 M.build()
375 I = ascpy.Integrator(M)
376 I.setEngine('IDA')
377 P = I.getParameters()
378 for p in P:
379 print p.getName(),"=",p.getValue()
380 assert len(P)==12
381 assert P[0].isStr()
382 assert P[0].getName()=="linsolver"
383 assert P[0].getValue()=='DENSE'
384 assert P[2].getName()=="maxord"
385 assert P[3].getName()=="autodiff"
386 assert P[3].getValue()==True
387 assert P[8].getName()=="atolvect"
388 assert P[8].getBoolValue() == True
389 P[3].setBoolValue(False)
390 assert P[3].getBoolValue()==False
391 I.setParameters(P)
392 assert I.getParameterValue('autodiff')==False
393 I.setParameter('autodiff',True)
394 try:
395 v = I.getParameterValue('nonexist')
396 except KeyError:
397 pass
398 else:
399 self.fail('Failed to trip invalid Integrator parameter')
400
401 class TestLSODE(Ascend):
402
403 def testzill(self):
404 self.L.load('johnpye/zill.a4c')
405 T = self.L.findType('zill')
406 M = T.getSimulation('sim')
407 M.setSolver(ascpy.Solver('QRSlv'))
408 I = ascpy.Integrator(M)
409 I.setEngine('LSODE')
410 I.setMinSubStep(1e-7)
411 I.setMaxSubStep(0.001)
412 I.setMaxSubSteps(10000)
413 I.setReporter(ascpy.IntegratorReporterConsole(I))
414 I.setLinearTimesteps(ascpy.Units(), 1.0, 1.5, 5)
415 I.analyse()
416 I.solve()
417 M.run(T.getMethod('self_test'))
418
419 def testnewton(self):
420 sys.stderr.write("STARTING TESTNEWTON\n")
421 self.L.load('johnpye/newton.a4c')
422 T = self.L.findType('newton')
423 M = T.getSimulation('sim')
424 M.solve(ascpy.Solver("QRSlv"),ascpy.SolverReporter())
425 I = ascpy.Integrator(M)
426 I.setEngine('LSODE')
427 I.setParameter('rtolvect',False)
428 I.setParameter('rtol',1e-7)
429 I.setParameter('atolvect',False)
430 I.setParameter('atol',1e-7)
431 I.setMinSubStep(1e-7)
432 I.setMaxSubStep(0.001)
433 I.setMaxSubSteps(10000)
434
435 I.setReporter(ascpy.IntegratorReporterConsole(I))
436 I.setLinearTimesteps(ascpy.Units("s"), 0, 2*float(M.v)/float(M.g), 2)
437 I.analyse()
438 I.solve()
439 print "At end of simulation,"
440 print "x = %f" % M.x
441 print "v = %f" % M.v
442 M.run(T.getMethod('self_test'))
443
444 def testlotka(self):
445 self.L.load('johnpye/lotka.a4c')
446 M = self.L.findType('lotka').getSimulation('sim')
447 M.setSolver(ascpy.Solver("QRSlv"))
448 I = ascpy.Integrator(M)
449 I.setEngine('LSODE')
450 I.setReporter(ascpy.IntegratorReporterConsole(I))
451 I.setLinearTimesteps(ascpy.Units("s"), 0, 200, 5)
452 I.analyse()
453 print "Number of vars = %d" % I.getNumVars()
454 assert I.getNumVars()==2
455 I.solve()
456 assert I.getNumObservedVars() == 3;
457 assert abs(M.R - 832) < 1.0
458 assert abs(M.F - 21.36) < 0.1
459
460 def testwritegraph(self):
461 self.L.load('johnpye/lotka.a4c')
462 M = self.L.findType('lotka').getSimulation('sim')
463 F = file('lotka.png','w')
464 M.build()
465 M.write(F,"dot")
466
467
468 #-------------------------------------------------------------------------------
469 # Testing of a external blackbox functions
470
471 class TestBlackBox(AscendSelfTester):
472 def testparsefail0(self):
473 try:
474 self.L.load('test/blackbox/parsefail0.a4c')
475 self.fail("parsefail0 should not have loaded without errors")
476 except:
477 pass
478
479 def testparsefail1(self):
480 try:
481 self.L.load('test/blackbox/parsefail1.a4c')
482 self.fail("parsefail1 should not have loaded without errors")
483 except:
484 pass
485
486 def testparsefail2(self):
487 try:
488 self.L.load('test/blackbox/parsefail2.a4c')
489 self.fail("parsefail2 should not have loaded without errors")
490 except:
491 pass
492
493 def testparsefail3(self):
494 try:
495 self.L.load('test/blackbox/parsefail3.a4c')
496 self.fail("parsefail3 should not have loaded without errors")
497 except:
498 pass
499
500 def testparsefail4(self):
501 try:
502 self.L.load('test/blackbox/parsefail4.a4c')
503 self.fail("parsefail4 should not have loaded")
504 except:
505 pass
506
507 def testfail1(self):
508 """Mismatched arg counts check-- tests bbox, not ascend."""
509 self.L.load('test/blackbox/fail1.a4c')
510 try:
511 M = self.L.findType('fail1').getSimulation('sim')
512 self.fail("expected exception was not raised")
513 except RuntimeError,e:
514 print "Caught exception '%s', assumed ok" % e
515
516 def testfail2(self):
517 """Incorrect data arg check -- tests bbox, not ascend"""
518 self.L.load('test/blackbox/fail2.a4c')
519 try:
520 M = self.L.findType('fail2').getSimulation('sim')
521 self.fail("expected exception was not raised")
522 except RuntimeError,e:
523 print "Caught exception '%s', assumed ok (should mention errors during instantiation)" % e
524
525 def testpass1(self):
526 """simple single bbox forward solve"""
527 M = self._run('pass1',filename='test/blackbox/pass.a4c')
528
529 def testpass2(self):
530 """simple single bbox reverse solve"""
531 M = self._run('pass2',filename='test/blackbox/pass.a4c')
532
533 def testpass3(self):
534 """simple double bbox solve"""
535 M = self._run('pass3',filename='test/blackbox/pass3.a4c')
536
537 def testpass4(self):
538 """simple double bbox reverse solve"""
539 M = self._run('pass4',filename='test/blackbox/pass3.a4c')
540
541 def testpass5(self):
542 M = self._run('pass5',filename='test/blackbox/pass5.a4c')
543
544 def testpass6(self):
545 M = self._run('pass6',filename='test/blackbox/pass5.a4c')
546
547 def testpass7(self):
548 M = self._run('pass7',filename='test/blackbox/passmerge.a4c')
549
550 def testpass8(self):
551 M = self._run('pass8',filename='test/blackbox/passmerge.a4c')
552
553 def testpass9(self):
554 M = self._run('pass9',filename='test/blackbox/passmerge.a4c')
555
556 def testpass10(self):
557 M = self._run('pass10',filename='test/blackbox/passmerge.a4c')
558
559 def testpass11(self):
560 M = self._run('pass11',filename='test/blackbox/passmerge.a4c')
561
562 def testpass12(self):
563 M = self._run('pass12',filename='test/blackbox/passmerge.a4c')
564
565 # this test doesn't work: 'system is inconsistent' -- and structurally singular
566 # def testpass13(self):
567 # """cross-merged input/output solve"""
568 # M = self._run('pass13',filename='test/blackbox/passmerge.a4c')
569
570 def testpass14(self):
571 """cross-merged input/output reverse solve"""
572 M = self._run('pass14',filename='test/blackbox/passmerge.a4c')
573
574 def testpass20(self):
575 M = self._run('pass20',filename='test/blackbox/passarray.a4c')
576
577 def testparsefail21(self):
578 """dense array of black boxes wrong syntax"""
579 try:
580 self.L.load('test/blackbox/parsefail21.a4c')
581 self.fail("parsefail21 should not have loaded without errors")
582 except:
583 pass
584
585 def testpass22(self):
586 M = self._run('pass22',filename='test/blackbox/passarray.a4c')
587
588 def testpass23(self):
589 M = self._run('pass23',filename='test/blackbox/passarray.a4c')
590
591 def testpass61(self):
592 M = self._run('pass61',filename='test/blackbox/reinstantiate.a4c')
593
594 def testpass62(self):
595 M = self._run('pass62',filename='test/blackbox/reinstantiate.a4c')
596
597 def testpass64(self):
598 M = self._run('pass64',filename='test/blackbox/reinstantiate.a4c')
599
600 def testpass65(self):
601 M = self._run('pass65',filename='test/blackbox/reinstantiate.a4c')
602
603 def testpass66(self):
604 M = self._run('pass66',filename='test/blackbox/reinstantiate.a4c')
605
606 def testpass67(self):
607 M = self._run('pass67',filename='test/blackbox/reinstantiate.a4c')
608
609 class TestExtFn(AscendSelfTester):
610 def testextfntest(self):
611 M = self._run('extfntest',filename='johnpye/extfn/extfntest.a4c')
612 self.assertAlmostEqual(M.y, 2);
613 self.assertAlmostEqual(M.x, 1);
614 self.assertAlmostEqual(M.y, M.x + 1);
615
616 def testextrelfor(self):
617 M = self._run('extrelfor',filename='johnpye/extfn/extrelfor.a4c')
618
619 ## @TODO fix bug with badly-named bbox rel in a loop (Ben, maybe)
620 # def testextrelforbadnaming(self):
621 # self.L.load('johnpye/extfn/extrelforbadnaming.a4c')
622 # T = self.L.findType('extrelfor')
623 # M = T.getSimulation('sim')
624 # M.solve(ascpy.Solver('QRSlv'),ascpy.SolverReporter())
625 # print "x[1] = %f" % M.x[1]
626 # print "x[2] = %f" % M.x[2]
627 # print "x[3] = %f" % M.x[3]
628 # print "x[4] = %f" % M.x[4]
629 # print "x[5] = %f" % M.x[5]
630 # M.run(T.getMethod('self_test'))
631
632 def testextrelrepeat(self):
633 M = self._run('extrelrepeat',filename='johnpye/extfn/extrelrepeat.a4c')
634
635 #-------------------------------------------------------------------------------
636 # Testing of Sensitivity module
637
638 class TestSensitivity(AscendSelfTester):
639 def test1(self):
640 self.L.load('sensitivity_test.a4c')
641 T = self.L.findType('sensitivity_test')
642 M = T.getSimulation('sim',False)
643 M.run(T.getMethod('on_load'))
644 M.solve(ascpy.Solver('QRSlv'),ascpy.SolverReporter())
645 M.run(T.getMethod('analyse'))
646 M.run(T.getMethod('self_test'))
647
648 # def testall(self):
649 # self.L.load('sensitivity_test.a4c')
650 # T = self.L.findType('sensitivity_test_all')
651 # M = T.getSimulation('sim',False)
652 # M.run(T.getMethod('on_load'))
653 # M.solve(ascpy.Solver('QRSlv'),ascpy.SolverReporter())
654 # M.run(T.getMethod('analyse'))
655 # M.run(T.getMethod('self_test'))
656 # CAUSES CRASH
657
658 #-------------------------------------------------------------------------------
659 # Testing of a ExtPy - external python methods
660
661 class TestExtPy(AscendSelfTester):
662 def test1(self):
663 self.L.load('johnpye/extpy/extpytest.a4c')
664 T = self.L.findType('extpytest')
665 M = T.getSimulation('sim')
666 M.run(T.getMethod('self_test'))
667
668 def test2(self):
669 self.L.load('johnpye/extpy/extpytest.a4c')
670 T = self.L.findType('extpytest')
671 M = T.getSimulation('sim')
672 M.run(T.getMethod('pythonthing'))
673 M.run(T.getMethod('pythonthing'))
674 M.run(T.getMethod('pythonthing'))
675 M.run(T.getMethod('pythonthing'))
676 M.run(T.getMethod('pythonthing'))
677 M.run(T.getMethod('pythonthing'))
678 M.run(T.getMethod('pythonthing'))
679 # causes crash!
680
681 #-------------------------------------------------------------------------------
682 # Testing of saturated steam properties library (iapwssatprops.a4c)
683
684 class TestSteam(AscendSelfTester):
685
686 def testiapwssatprops1(self):
687 M = self._run('testiapwssatprops1',filename='steam/iapwssatprops.a4c')
688 def testiapwssatprops2(self):
689 M = self._run('testiapwssatprops2',filename='steam/iapwssatprops.a4c')
690 def testiapwssatprops3(self):
691 M = self._run('testiapwssatprops3',filename='steam/iapwssatprops.a4c')
692
693 # test the stream model basically works
694 def testsatsteamstream(self):
695 M = self._run('satsteamstream',filename='steam/satsteamstream.a4c')
696
697 # test that we can solve in terms of various (rho,u)
698 def testsatuv(self):
699 self.L.load('steam/iapwssat.a4c')
700 T = self.L.findType('testiapwssatuv')
701 M = T.getSimulation('sim',False)
702 M.run(T.getMethod('on_load'))
703 M.solve(ascpy.Solver('QRSlv'),ascpy.SolverReporter())
704 print "p = %f bar" % M.p.to('bar');
705 print "T = %f C" % (M.T.to('K') - 273.15);
706 print "x = %f" % M.x;
707 M.run(T.getMethod('self_test'))
708 M.run(T.getMethod('values2'))
709 # M.v.setRealValueWithUnits(1.0/450,"m^3/kg");
710 # M.solve(ascpy.Solver('QRSlv'),ascpy.SolverReporter())
711 M.solve(ascpy.Solver('QRSlv'),ascpy.SolverReporter())
712 print "p = %f bar" % M.p.to('bar');
713 print "T = %f C" % (M.T.to('K') - 273.15);
714 print "x = %f" % M.x;
715 M.run(T.getMethod('self_test2'))
716
717
718 ## @TODO fix error capture from bounds checking during initialisation
719 # def testiapwssat1(self):
720 # M = self._run('testiapwssat1',filename='steam/iapwssat.a4c')
721
722 def testdsgsat(self):
723 self.L.load('steam/dsgsat3.a4c')
724 T = self.L.findType('dsgsat3')
725 M = T.getSimulation('sim',False)
726 M.run(T.getMethod('on_load'))
727 M.solve(ascpy.Solver('QRSlv'),ascpy.SolverReporter())
728 self.assertAlmostEqual(M.dTw_dt[2],0.0);
729 M.run(T.getMethod('configure_dynamic'))
730 M.solve(ascpy.Solver('QRSlv'),ascpy.SolverReporter())
731 return M
732
733 def testdsgsatrepeat(self):
734 self.L.load('steam/dsgsat3.a4c')
735 T = self.L.findType('dsgsat3')
736 M = T.getSimulation('sim',False)
737 M.run(T.getMethod('on_load'))
738 M.solve(ascpy.Solver('QRSlv'),ascpy.SolverReporter())
739 M.run(T.getMethod('on_load'))
740 M.solve(ascpy.Solver('QRSlv'),ascpy.SolverReporter())
741 M.run(T.getMethod('on_load'))
742 M.solve(ascpy.Solver('QRSlv'),ascpy.SolverReporter())
743
744 def testvary(self):
745 self.L.load('steam/dsgsat3.a4c')
746 T = self.L.findType('dsgsat3')
747 M = T.getSimulation('sim',False)
748 M.run(T.getMethod('on_load'))
749 M.solve(ascpy.Solver('QRSlv'),ascpy.SolverReporter())
750 print "----- setting qdot_s -----"
751 M.qdot_s.setRealValueWithUnits(1000,"W/m")
752 M.solve(ascpy.Solver('QRSlv'),ascpy.SolverReporter())
753 print "----- setting qdot_s -----"
754 M.qdot_s.setRealValueWithUnits(2000,"W/m")
755 M.solve(ascpy.Solver('QRSlv'),ascpy.SolverReporter())
756
757 def teststeadylsode(self):
758 "test that steady conditions are stable with LSODE"
759 M = self.testdsgsat()
760 #M.qdot_s.setRealValueWithUnits(1000,"W/m")
761 M.solve(ascpy.Solver('QRSlv'),ascpy.SolverReporter())
762 #M.setParameter('
763 I = ascpy.Integrator(M)
764 I.setEngine('LSODE')
765 I.setReporter(ascpy.IntegratorReporterConsole(I))
766 I.setLinearTimesteps(ascpy.Units("s"), 0, 3600, 10)
767 I.analyse()
768 I.solve()
769
770 # def testpeturblsode(self):
771 # "test that steady conditions are stable with LSODE"
772 # M = self.testdsgsat()
773 # # here is the peturbation...
774 # M.qdot_s.setRealValueWithUnits(1000,"W/m")
775 # M.solve(ascpy.Solver('QRSlv'),ascpy.SolverReporter())
776 # I = ascpy.Integrator(M)
777 # I.setEngine('LSODE')
778 # I.setReporter(ascpy.IntegratorReporterConsole(I))
779 # I.setLinearTimesteps(ascpy.Units("s"), 0, 5, 1)
780 # I.analyse()
781 # I.solve()
782
783 def teststeadyida(self):
784 M = self.testdsgsat()
785 self.assertAlmostEqual(M.dTw_dt[2],0.0)
786 Tw1 = float(M.T_w[2])
787 T = self.L.findType('dsgsat3')
788 M.run(T.getMethod('free_states'))
789 I = ascpy.Integrator(M)
790 I.setEngine('IDA')
791 I.setParameter('linsolver','DENSE')
792 I.setParameter('safeeval',True)
793 I.setParameter('rtol',1e-4)
794 I.setParameter('atolvect',False)
795 I.setParameter('atol',1e-4)
796 I.setParameter('maxord',3)
797 I.setInitialSubStep(0.001)
798 I.setReporter(ascpy.IntegratorReporterConsole(I))
799 I.setLinearTimesteps(ascpy.Units("s"), 0, 3600, 10)
800 I.analyse()
801 I.solve()
802 self.assertAlmostEqual(float(M.T_w[2]),Tw1)
803 M.qdot_s.setRealValueWithUnits(1000,"W/m")
804 self.assertAlmostEqual(M.qdot_s.to("W/m"),1000)
805 M.solve(ascpy.Solver('QRSlv'),ascpy.SolverReporter())
806 print "dTw/dt = %f" % M.dTw_dt[2]
807 self.assertNotAlmostEqual(M.dTw_dt[2],0.0)
808 F=file('dsgsat.dot','w')
809 M.write(F,'dot')
810
811 def testpeturbida(self):
812 M = self.testdsgsat()
813 self.assertAlmostEqual(M.dTw_dt[2],0.0)
814 T = self.L.findType('dsgsat3')
815 M.run(T.getMethod('free_states'))
816 # here is the peturbation...
817 qdot_s = float(M.qdot_s)
818 print "OLD QDOT_S = %f" % qdot_s
819 M.qdot_s.setRealValueWithUnits(6000,"W/m")
820 # IDA has its own initial conditions solver, so no need to call QRSlv here
821 I = ascpy.Integrator(M)
822 I.setEngine('IDA')
823 I.setParameter('linsolver','DENSE')
824 I.setReporter(ascpy.IntegratorReporterConsole(I))
825 #I.setLinearTimesteps(ascpy.Units("s"), 0,300,300)
826 I.setLogTimesteps(ascpy.Units("s"), 0.009, 1200, 150)
827 I.analyse()
828 F = file('ga.mm','w')
829 I.writeMatrix(F,'dg/dz')
830 F = file('gd.mm','w')
831 I.writeMatrix(F,'dg/dx')
832 F = file('fa.mm','w')
833 I.writeMatrix(F,'df/dz')
834 F = file('fd.mm','w')
835 I.writeMatrix(F,'df/dx')
836 F = file('fdp.mm','w')
837 I.writeMatrix(F,"df/dx'")
838 I.solve()
839
840 #-------------------------------------------------------------------------------
841 # Testing of freesteam external steam properties functions
842
843 with_freesteam = True
844 try:
845 # we assume that if the freesteam python module is installed, the ASCEND
846 # external library will also be.
847 import freesteam
848 have_freesteam = True
849 except ImportError,e:
850 have_freesteam = False
851
852 if with_freesteam and have_freesteam:
853 class TestFreesteam(AscendSelfTester):
854 # def testfreesteamtest(self):
855 # """run the self-test cases bundled with freesteam"""
856 # self._run('testfreesteam',filename='testfreesteam.a4c')
857
858 def testload(self):
859 """check that we can load 'thermalequilibrium2' (IMPORT "freesteam", etc)"""
860 self.L.load('johnpye/thermalequilibrium2.a4c')
861
862 def testinstantiate(self):
863 """load an instantiate 'thermalequilibrium2'"""
864 self.testload()
865 M = self.L.findType('thermalequilibrium2').getSimulation('sim')
866 return M
867
868 def testintegrate(self):
869 """integrate transfer of heat from one mass of water/steam to another
870 according to Newton's law of cooling"""
871 M = self.testinstantiate()
872 M.setSolver(ascpy.Solver("QRSlv"))
873 I = ascpy.Integrator(M)
874 I.setEngine('LSODE')
875 I.setReporter(ascpy.IntegratorReporterConsole(I))
876 I.setLinearTimesteps(ascpy.Units("s"), 0, 3000, 30)
877 I.setMinSubStep(0.01)
878 I.setInitialSubStep(1)
879 I.analyse()
880 print "Number of vars = %d" % I.getNumVars()
881 assert I.getNumVars()==2
882 I.solve()
883 assert I.getNumObservedVars() == 3;
884 print "S[1].T = %f K" % M.S[1].T
885 print "S[2].T = %f K" % M.S[2].T
886 print "Q = %f W" % M.Q
887 self.assertAlmostEqual(float(M.S[1].T),506.77225109,4);
888 self.assertAlmostEqual(float(M.S[2].T),511.605173967,5);
889 self.assertAlmostEqual(float(M.Q),-48.32922877329,3);
890 self.assertAlmostEqual(float(M.t),3000);
891 print "Note that the above values have not been verified analytically"
892
893 def testcollapsingcan2(self):
894 """ solve the collapsing can model using IAPWS-IF97 steam props """
895 M = self._run("collapsingcan2",filename="collapsingcan2.a4c");
896
897 #-------------------------------------------------------------------------------
898 # Testing of the brent-solver EXTERNAL method
899
900 class TestBrent(AscendSelfTester):
901 def testbrent(self):
902 M = self._run('brent1',filename='test/brent.a4c')
903
904 #-------------------------------------------------------------------------------
905 # Testing of IDA's analysis module
906
907 class TestIDA(Ascend):
908 def _run(self,filen,modeln=""):
909 self.L.load('test/ida/%s.a4c' % filen)
910 T = self.L.findType('%s%s' % (filen,modeln))
911 M = T.getSimulation('sim')
912 M.build()
913 I = ascpy.Integrator(M)
914 I.setEngine('IDA')
915 I.analyse()
916 return M;
917
918 def _runfail(self,filen,n,msg="failed"):
919 try:
920 self._run(filen,'fail%d' % n)
921 except Exception,e:
922 print "(EXPECTED) ERROR: %s" % e
923 return
924 self.fail(msg)
925
926 def testsinglederiv(self):
927 self._run('singlederiv')
928
929 def testsinglederivfail1(self):
930 self._runfail('singlederiv',1
931 ,"t.ode_id=-1 did not trigger error")
932
933 def testsinglederivfail2(self):
934 self._runfail('singlederiv',2
935 ,"dy_dt.ode_id=2 did not trigger error")
936
937 def testsinglederivfail3(self):
938 self._runfail('singlederiv',3
939 ,"dy_dt.ode_type=3 did not trigger error")
940
941 def testsinglederivfail4(self):
942 self._runfail('singlederiv',4
943 ,"duplicate ode_type=1 did not trigger error")
944
945 def testsinglederivfail5(self):
946 self._runfail('singlederiv',5
947 ,"duplicate ode_type=1 did not trigger error")
948
949 def testsinglederivfail6(self):
950 self._runfail('singlederiv',6
951 ,"duplicate ode_type=1 did not trigger error")
952
953 def testtwoderiv(self):
954 self._run('twoderiv')
955
956 def testtwoderivfail1(self):
957 self._runfail('twoderiv',1)
958
959 def testtwoderivfail2(self):
960 self._runfail('twoderiv',2)
961
962 def testtwoderivfail3(self):
963 self._runfail('twoderiv',3)
964 def testtwoderivfail4(self):
965 self._runfail('twoderiv',4)
966 def testtwoderivfail5(self):
967 self._runfail('twoderiv',5)
968
969 def testnoderivs(self):
970 self._runfail('noderivs',1)
971
972 def testnoindeps(self):
973 self._runfail('indeps',1)
974
975 def testtwoindeps(self):
976 self._runfail('indeps',2)
977
978 def testfixedvars(self):
979 self._run('fixedvars')
980
981 def testfixedvars1(self):
982 self._run('fixedvars',1)
983
984 # fails the index check
985 # def testfixedvars2(self):
986 # self._run('fixedvars',2)
987
988 # fails the index check
989 # def testfixedvars3(self):
990 # self._run('fixedvars',3)
991
992 def testincidence(self):
993 self._run('incidence')
994
995 def testincidence1(self):
996 self._run('incidence',1)
997 def testincidence2(self):
998 self._run('incidence',2)
999 def testincidence3(self):
1000 M = self._run('incidence',3)
1001
1002 def testincidence4(self):
1003 self._run('incidence',4)
1004 def testincidencefail5(self):
1005 self._runfail('incidence',5)
1006
1007 def testwritematrix(self):
1008 self.L.load('test/ida/writematrix.a4c')
1009 T = self.L.findType('writematrix')
1010 M = T.getSimulation('sim')
1011 M.build()
1012 I = ascpy.Integrator(M)
1013 I.setEngine('IDA')
1014 I.analyse()
1015 # this stuff fails on Windows because FILE structure is different python vs mingw
1016 # F = os.tmpfile()
1017 # I.writeMatrix(F,"dF/dy")
1018 # F.seek(0)
1019 # print F.read()
1020 # F1 = os.tmpfile()
1021 # I.writeMatrix(F1,"dF/dy'")
1022 # F1.seek(0)
1023 # print F1.read()
1024 # F1 = os.tmpfile()
1025 # I.writeMatrix(F1,"dg/dx")
1026 # F1.seek(0)
1027 # print F1.read()
1028 # # for the moment you'll have to check these results manually.
1029
1030 def testwritematrix2(self):
1031 self.L.load('test/ida/writematrix.a4c')
1032 T = self.L.findType('writematrix2')
1033 M = T.getSimulation('sim')
1034 M.build()
1035 I = ascpy.Integrator(M)
1036 I.setEngine('IDA')
1037 I.analyse()
1038 # this stuff fails on Windows because FILE structure is different python vs mingw
1039 # F = os.tmpfile()
1040 # I.writeMatrix(F,"dF/dy")
1041 # F.seek(0)
1042 # print F.read()
1043 # F1 = os.tmpfile()
1044 # I.writeMatrix(F1,"dF/dy'")
1045 # F1.seek(0)
1046 # print F1.read()
1047 # F1 = os.tmpfile()
1048 # I.writeMatrix(F1,"dg/dx")
1049 # F1.seek(0)
1050 # print F1.read()
1051 #F1 = os.tmpfile()
1052 #I.writeMatrix(F1,"dydp/dyd")
1053 #F1.seek(0)
1054 #print F1.read()
1055 # for the moment you'll have to check these results manually.
1056
1057 def testindexproblem(self):
1058 self.L.load('test/ida/indexproblem.a4c')
1059 T = self.L.findType('indexproblem')
1060 M = T.getSimulation('sim')
1061 M.build()
1062 I = ascpy.Integrator(M)
1063 I.setEngine('IDA')
1064 I.analyse()
1065 pass
1066
1067 def testindexproblem2(self):
1068 self.L.load('test/ida/indexproblem.a4c')
1069 T = self.L.findType('indexproblem2')
1070 M = T.getSimulation('sim')
1071 M.build()
1072 I = ascpy.Integrator(M)
1073 I.setEngine('IDA')
1074 try:
1075 I.analyse()
1076 except Exception,e:
1077 return
1078 self.fail("Index problem not detected")
1079
1080 def testboundaries(self):
1081 self.L.load('test/ida/boundaries.a4c')
1082 T = self.L.findType('boundaries')
1083 M = T.getSimulation('sim')
1084 M.build()
1085 I = ascpy.Integrator(M)
1086 I.setEngine('IDA')
1087 I.analyse()
1088 I.setLogTimesteps(ascpy.Units("s"), 0.1, 20, 20)
1089 I.setParameter('linsolver','DENSE')
1090 I.setParameter('calcic','Y')
1091 I.setParameter('linsolver','DENSE')
1092 I.setParameter('safeeval',False)
1093 I.setReporter(ascpy.IntegratorReporterConsole(I))
1094 I.solve()
1095
1096 # doesn't work yet:
1097 # def testincidence5(self):
1098 # self._run('incidence',5)
1099
1100
1101 #-------------------------------------------------------------------------------
1102 # Testing of IDA models using DENSE linear solver
1103
1104 class TestIDADENSE(Ascend):
1105 """IDA DAE integrator, DENSE linear solver"""
1106
1107 def testlotka(self):
1108 self.L.load('johnpye/lotka.a4c')
1109 M = self.L.findType('lotka').getSimulation('sim')
1110 M.setSolver(ascpy.Solver("QRSlv"))
1111 I = ascpy.Integrator(M)
1112 I.setEngine('IDA')
1113 I.setReporter(ascpy.IntegratorReporterConsole(I))
1114 I.setLinearTimesteps(ascpy.Units("s"), 0, 200, 5);
1115 I.setParameter('linsolver','DENSE')
1116 I.setParameter('rtol',1e-8);
1117 I.analyse()
1118 assert I.getNumVars()==2
1119 assert abs(M.R - 1000) < 1e-300
1120 I.solve()
1121 assert I.getNumObservedVars() == 3
1122 assert abs(M.R - 832) < 1.0
1123 assert abs(M.F - 21.36) < 0.1
1124
1125 def testdenx(self):
1126 print "-----------------------------====="
1127 self.L.load('johnpye/idadenx.a4c')
1128 M = self.L.findType('idadenx').getSimulation('sim')
1129 M.setSolver(ascpy.Solver("QRSlv"))
1130 I = ascpy.Integrator(M)
1131 I.setEngine('IDA')
1132 I.setParameter('calcic','YA_YDP')
1133 I.setParameter('linsolver','DENSE')
1134 I.setParameter('safeeval',True)
1135 I.setReporter(ascpy.IntegratorReporterConsole(I))
1136 I.setLogTimesteps(ascpy.Units("s"), 0.4, 4e10, 11)
1137 I.setMaxSubStep(0);
1138 I.setInitialSubStep(0)
1139 I.setMaxSubSteps(0)
1140 I.setParameter('autodiff',True)
1141 I.analyse()
1142 I.solve()
1143 assert abs(float(M.y1) - 5.1091e-08) < 2e-9
1144 assert abs(float(M.y2) - 2.0437e-13) < 2e-14
1145 assert abs(float(M.y3) - 1.0) < 1e-5
1146
1147 def testhires(self):
1148 self.L.load('test/hires.a4c')
1149 T = self.L.findType('hires')
1150 M = T.getSimulation('sim')
1151 M.setSolver(ascpy.Solver('QRSlv'))
1152 I = ascpy.Integrator(M)
1153 I.setEngine('IDA')
1154 I.setParameter('linsolver','DENSE')
1155 I.setParameter('rtol',1.1e-15)
1156 I.setParameter('atolvect',0)
1157 I.setParameter('atol',1.1e-15)
1158 I.setReporter(ascpy.IntegratorReporterConsole(I))
1159 I.setLogTimesteps(ascpy.Units(""), 1, 321.8122, 5)
1160 I.setInitialSubStep(1e-5)
1161 I.setMaxSubSteps(10000)
1162 I.analyse()
1163 I.solve()
1164 for i in range(8):
1165 print "y[%d] = %.20g" % (i+1, M.y[i+1])
1166 M.run(T.getMethod('self_test'))
1167
1168 def testchemakzo(self):
1169 self.L.load('test/chemakzo.a4c')
1170 T = self.L.findType('chemakzo')
1171 M = T.getSimulation('sim')
1172 M.setSolver(ascpy.Solver('QRSlv'))
1173 I = ascpy.Integrator(M)
1174 I.setEngine('IDA')
1175 I.setParameter('linsolver','DENSE')
1176 I.setParameter('rtol',1e-15)
1177 I.setParameter('atolvect',0)
1178 I.setParameter('atol',1e-15)
1179 I.setReporter(ascpy.IntegratorReporterConsole(I))
1180 I.setLinearTimesteps(ascpy.Units("s"), 1, 180, 5)
1181 I.setInitialSubStep(1e-13)
1182 I.setMaxSubSteps(10000)
1183 I.analyse()
1184 I.solve()
1185 for i in range(6):
1186 print "y[%d] = %.20g" % (i+1, M.y[i+1])
1187 M.run(T.getMethod('self_test'))
1188
1189 def testtransamp(self):
1190 self.L.load('test/transamp.a4c')
1191 T = self.L.findType('transamp')
1192 M = T.getSimulation('sim')
1193 M.setSolver(ascpy.Solver('QRSlv'))
1194 I = ascpy.Integrator(M)
1195 I.setEngine('IDA')
1196 I.setParameter('linsolver','DENSE')
1197 I.setParameter('rtol',1e-7)
1198 I.setParameter('atolvect',0)
1199 I.setParameter('atol',1e-7)
1200 I.setReporter(ascpy.IntegratorReporterConsole(I))
1201 I.setLinearTimesteps(ascpy.Units("s"), 0.05, 0.2, 20)
1202 I.setInitialSubStep(0.00001)
1203 I.setMaxSubSteps(10000)
1204 I.analyse()
1205 I.solve()
1206 for i in range(6):
1207 print "y[%d] = %.20g" % (i+1, M.y[i+1])
1208 M.run(T.getMethod('self_test'))
1209
1210 # MODEL FAILS ANALYSIS: we need to add support for non-incident differential vars
1211 # def testpollution(self):
1212 # self.L.load('test/pollution.a4c')
1213 # T = self.L.findType('pollution')
1214 # M = T.getSimulation('sim')
1215 # M.setSolver(ascpy.Solver('QRSlv'))
1216 # I = ascpy.Integrator(M)
1217 # I.setEngine('IDA')
1218 # I.setParameter('linsolver','DENSE')
1219 # I.setParameter('rtol',1.1e-15)
1220 # I.setParameter('atolvect',0)
1221 # I.setParameter('atol',1.1e-15)
1222 # I.setReporter(ascpy.IntegratorReporterConsole(I))
1223 # I.setLogTimesteps(ascpy.Units("s"), 1, 60, 5)
1224 # I.setInitialSubStep(1e-5)
1225 # I.setMaxSubSteps(10000)
1226 # I.analyse()
1227 # I.solve()
1228 # for i in range(20):
1229 # print "y[%d] = %.20g" % (i+1, M.y[i+1])
1230 # M.run(T.getMethod('self_test'))
1231
1232 ## @TODO fails during IDACalcIC (model too big?)
1233 # def testkryx(self):
1234 # self.L.load('johnpye/idakryx.a4c')
1235 # ascpy.getCompiler().setUseRelationSharing(False)
1236 # M = self.L.findType('idakryx').getSimulation('sim')
1237 # M.setSolver(ascpy.Solver('QRSlv'))
1238 # M.build()
1239 # I = ascpy.Integrator(M)
1240 # I.setEngine('IDA')
1241 # I.setReporter(ascpy.IntegratorReporterConsole(I))
1242 # I.setParameter('linsolver','DENSE')
1243 # I.setParameter('maxl',8)
1244 # I.setParameter('gsmodified',False)
1245 # I.setParameter('autodiff',True)
1246 # I.setParameter('rtol',0)
1247 # I.setParameter('atol',1e-3);
1248 # I.setParameter('atolvect',False)
1249 # I.setParameter('calcic','YA_YDP')
1250 # I.analyse()
1251 # I.setLogTimesteps(ascpy.Units("s"), 0.01, 10.24, 11)
1252 # I.solve()
1253 # assert abs(M.u[2][2].getValue()) < 1e-5
1254
1255 #-------------------------------------------------------------------------------
1256 # Testing of IDA models using SPGMR linear solver (Krylov)
1257
1258 # these tests are disabled until SPGMR preconditioning has been implemented
1259 class TestIDASPGMR:#(Ascend):
1260 def testlotka(self):
1261 self.L.load('johnpye/lotka.a4c')
1262 M = self.L.findType('lotka').getSimulation('sim')
1263 M.setSolver(ascpy.Solver("QRSlv"))
1264 I = ascpy.Integrator(M)
1265 I.setEngine('IDA')
1266 I.setReporter(ascpy.IntegratorReporterConsole(I))
1267 I.setLinearTimesteps(ascpy.Units("s"), 0, 200, 5)
1268 I.setParameter('rtol',1e-8)
1269 I.analyse()
1270 assert I.getNumVars()==2
1271 assert abs(M.R - 1000) < 1e-300
1272 I.solve()
1273 assert I.getNumObservedVars() == 3
1274 assert abs(M.R - 832) < 1.0
1275 assert abs(M.F - 21.36) < 0.1
1276
1277
1278 def testkryx(self):
1279 self.L.load('johnpye/idakryx.a4c')
1280 M = self.L.findType('idakryx').getSimulation('sim')
1281 M.build()
1282 I = ascpy.Integrator(M)
1283 I.setEngine('IDA')
1284 I.setReporter(ascpy.IntegratorReporterConsole(I))
1285 I.setParameter('linsolver','SPGMR')
1286 I.setParameter('prec','JACOBI')
1287 I.setParameter('maxl',8)
1288 I.setParameter('gsmodified',False)
1289 I.setParameter('autodiff',True)
1290 I.setParameter('gsmodified',True)
1291 I.setParameter('rtol',0)
1292 I.setParameter('atol',1e-3);
1293 I.setParameter('atolvect',False)
1294 I.setParameter('calcic','Y')
1295 I.analyse()
1296 I.setLogTimesteps(ascpy.Units("s"), 0.01, 10.24, 10);
1297 print M.udot[1][3]
1298 I.solve()
1299 assert 0
1300
1301 def testzill(self):
1302 self.L.load('johnpye/zill.a4c')
1303 T = self.L.findType('zill')
1304 M = T.getSimulation('sim')
1305 M.setSolver(ascpy.Solver('QRSlv'))
1306 I = ascpy.Integrator(M)
1307 I.setEngine('IDA')
1308 I.setParameter('safeeval',False)
1309 I.setMinSubStep(1e-7)
1310 I.setMaxSubStep(0.001)
1311 I.setMaxSubSteps(10000)
1312 I.setReporter(ascpy.IntegratorReporterConsole(I))
1313 I.setLinearTimesteps(ascpy.Units(), 1.0, 1.5, 5)
1314 I.analyse()
1315 I.solve()
1316 M.run(T.getMethod('self_test'))
1317
1318 def testdenxSPGMR(self):
1319 self.L.load('johnpye/idadenx.a4c')
1320 M = self.L.findType('idadenx').getSimulation('sim')
1321 M.setSolver(ascpy.Solver('QRSlv'))
1322 I = ascpy.Integrator(M)
1323 I.setEngine('IDA')
1324 I.setReporter(ascpy.IntegratorReporterConsole(I))
1325 I.setLogTimesteps(ascpy.Units("s"), 0.4, 4e10, 11)
1326 I.setMaxSubStep(0);
1327 I.setInitialSubStep(0);
1328 I.setMaxSubSteps(0);
1329 I.setParameter('autodiff',True)
1330 I.setParameter('linsolver','SPGMR')
1331 I.setParameter('gsmodified',False)
1332 I.setParameter('maxncf',10)
1333 I.analyse()
1334 I.solve()
1335 assert abs(float(M.y1) - 5.1091e-08) < 1e-10
1336 assert abs(float(M.y2) - 2.0437e-13) < 1e-15
1337 assert abs(float(M.y3) - 1.0) < 1e-5
1338
1339 class TestDOPRI5(Ascend):
1340 def testlotka(self):
1341 self.L.load('test/dopri5/dopri5test.a4c')
1342 M = self.L.findType('dopri5test').getSimulation('sim')
1343 M.setSolver(ascpy.Solver("QRSlv"))
1344 M.solve(ascpy.Solver("QRSlv"),ascpy.SolverReporter())
1345 I = ascpy.Integrator(M)
1346 I.setEngine('DOPRI5')
1347 I.setReporter(ascpy.IntegratorReporterConsole(I))
1348 I.setLinearTimesteps(ascpy.Units("s"), 0, 200, 20)
1349 I.setParameter('rtol',1e-8)
1350 I.analyse()
1351 assert I.getNumVars()==1
1352 I.solve()
1353 def testaren(self):
1354 self.L.load('test/dopri5/aren.a4c')
1355 M = self.L.findType('aren').getSimulation('sim')
1356 M.setSolver(ascpy.Solver("QRSlv"))
1357 M.solve(ascpy.Solver("QRSlv"),ascpy.SolverReporter())
1358 I = ascpy.Integrator(M)
1359 I.setEngine('DOPRI5')
1360 I.setReporter(ascpy.IntegratorReporterConsole(I))
1361 #xend = 17.0652165601579625588917206249
1362 I.setLinearTimesteps(ascpy.Units("s"), 0, 17.0652165601579625588917206249, 10)
1363 I.setParameter('rtol',1e-7)
1364 I.setParameter('atol',1e-7)
1365 I.setParameter('tolvect',False)
1366 I.setMinSubStep(0);
1367 I.setMaxSubStep(0);
1368 I.setInitialSubStep(0);
1369 I.analyse()
1370 I.solve()
1371 print "y[0] = %f" % float(M.y[0])
1372 assert abs(float(M.y[0]) - 0.994) < 1e-5
1373 assert abs(float(M.y[1]) - 0.0) < 1e-5
1374
1375 class TestIPOPT(Ascend):
1376 def test1(self):
1377 self.L.load('test/ipopt/test1.a4c')
1378 M = self.L.findType('test1').getSimulation('sim')
1379 M.setSolver(ascpy.Solver("IPOPT"))
1380 M.setParameter('linear_solver','mumps')
1381 M.setParameter('hessian_approximation','limited-memory')
1382 M.solve(ascpy.Solver("IPOPT"),ascpy.SolverReporter())
1383
1384 # test some stuff for beam calculations
1385 class TestSection(Ascend):
1386 def test_compound3(self):
1387 self.L.load('johnpye/section.a4c')
1388 T = self.L.findType('compound_section_test3')
1389 M = T.getSimulation('sim')
1390 M.solve(ascpy.Solver("QRSlv"),ascpy.SolverReporter())
1391 M.run(T.getMethod('self_test'))
1392 def test_compound4(self):
1393 self.L.load('johnpye/section.a4c')
1394 T = self.L.findType('compound_section_test4')
1395 M = T.getSimulation('sim')
1396 M.solve(ascpy.Solver("QRSlv"),ascpy.SolverReporter())
1397 M.run(T.getMethod('self_test'))
1398 def test_compound2(self):
1399 self.L.load('johnpye/section.a4c')
1400 T = self.L.findType('compound_section_test2')
1401 M = T.getSimulation('sim')
1402 M.solve(ascpy.Solver("QRSlv"),ascpy.SolverReporter())
1403 M.run(T.getMethod('self_test'))
1404
1405 # move code above down here if you want to temporarily avoid testing it
1406 class NotToBeTested:
1407 def nothing(self):
1408 pass
1409
1410 def testnewton(self):
1411 sys.stderr.write("STARTING TESTNEWTON\n")
1412 self.L.load('johnpye/newton.a4c')
1413 T = self.L.findType('newton')
1414 M = T.getSimulation('sim')
1415 M.solve(ascpy.Solver("QRSlv"),ascpy.SolverReporter())
1416 I = ascpy.Integrator(M)
1417 I.setEngine('IDA')
1418 I.setParameter('linsolver','DENSE')
1419 I.setParameter('safeeval',True)
1420 I.setParameter('rtol',1e-8)
1421 I.setMaxSubStep(0.001)
1422 I.setMaxSubSteps(10000)
1423
1424 I.setReporter(ascpy.IntegratorReporterConsole(I))
1425 I.setLinearTimesteps(ascpy.Units("s"), 0, 2*float(M.v)/float(M.g), 2)
1426 I.analyse()
1427 I.solve()
1428 print "At end of simulation,"
1429 print "x = %f" % M.x
1430 print "v = %f" % M.v
1431 M.run(T.getMethod('self_test'))
1432
1433 def patchpath(VAR,SEP,addvals):
1434 restart = 0
1435 envpath = [os.path.abspath(i) for i in os.environ[VAR].split(SEP)]
1436 for l in addvals:
1437 if l in envpath[len(addvals):]:
1438 envpath.remove(l)
1439 restart = 1
1440 for l in addvals:
1441 if l not in envpath:
1442 envpath.insert(0,l)
1443 restart = 1
1444 os.environ[VAR] = SEP.join(envpath)
1445 return restart
1446
1447 if __name__=='__main__':
1448 # a whole bag of tricks to make sure we get the necessary dirs in our ascend, python and ld path vars
1449 restart = 0
1450
1451 if platform.system()=="Windows":
1452 LD_LIBRARY_PATH="PATH"
1453 SEP = ";"
1454 else:
1455 LD_LIBRARY_PATH="LD_LIBRARY_PATH"
1456 SEP = ":"
1457
1458 solverdir = os.path.abspath(os.path.join(sys.path[0],"solvers"))
1459 solverdirs = [os.path.join(solverdir,s) for s in "qrslv","cmslv","lrslv","conopt","ida","lsode","ipopt","dopri5"]
1460
1461 if not os.environ.get('ASCENDSOLVERS'):
1462 os.environ['ASCENDSOLVERS'] = SEP.join(solverdirs)
1463 restart = 1
1464 else:
1465 if patchpath('ASCENDSOLVERS',SEP,solverdirs):
1466 restart = 1
1467
1468 freesteamdir = os.path.expanduser("~/freesteam/ascend")
1469 modeldirs = [os.path.abspath(os.path.join(sys.path[0],"models")),os.path.abspath(freesteamdir)]
1470
1471 if not os.environ.get('ASCENDLIBRARY'):
1472 os.environ['ASCENDLIBRARY'] = SEP.join(modeldirs)
1473 restart = 1
1474 else:
1475 if patchpath('ASCENDLIBRARY',SEP,modeldirs):
1476 restart = 1
1477
1478 libdirs = ["pygtk","."]
1479 libdirs = [os.path.normpath(os.path.join(sys.path[0],l)) for l in libdirs]
1480 if not os.environ.get(LD_LIBRARY_PATH):
1481 os.environ[LD_LIBRARY_PATH]=SEP.join(libdirs)
1482 restart = 1
1483 else:
1484 envlibdirs = [os.path.normpath(i) for i in os.environ[LD_LIBRARY_PATH].split(SEP)]
1485 for l in libdirs:
1486 if l in envlibdirs[len(libdirs):]:
1487 envlibdirs.remove(l)
1488 restart = 1
1489 for l in libdirs:
1490 if l not in envlibdirs:
1491 envlibdirs.insert(0,l)
1492 restart = 1
1493 os.environ[LD_LIBRARY_PATH] = SEP.join(envlibdirs)
1494
1495 pypath = os.path.normpath(os.path.join(sys.path[0],"pygtk"))
1496 if not os.environ.get('PYTHONPATH'):
1497 os.environ['PYTHONPATH']=pypath
1498 else:
1499 envpypath = os.environ['PYTHONPATH'].split(SEP)
1500 if pypath not in envpypath:
1501 envpypath.insert(0,pypath)
1502 os.environ['PYTHONPATH']=SEP.join(envpypath)
1503 restart = 1
1504
1505 if restart and platform.system()!="Windows":
1506 script = os.path.join(sys.path[0],"test.py")
1507 sys.stderr.write("Restarting with...\n")
1508 sys.stderr.write(" export LD_LIBRARY_PATH=%s\n" % os.environ.get(LD_LIBRARY_PATH))
1509 sys.stderr.write(" export PYTHONPATH=%s\n" % os.environ.get('PYTHONPATH'))
1510 sys.stderr.write(" export ASCENDLIBRARY=%s\n" % os.environ.get('ASCENDLIBRARY'))
1511 sys.stderr.write(" export ASCENDSOLVERS=%s\n" % os.environ.get('ASCENDSOLVERS'))
1512 sys.stderr.flush()
1513 os.execvp("python",[script] + sys.argv)
1514 exit(1)
1515 else:
1516 sys.stderr.write("Got...\n")
1517 sys.stderr.write(" LD_LIBRARY_PATH=%s\n" % os.environ.get(LD_LIBRARY_PATH))
1518 sys.stderr.write(" PYTHONPATH=%s\n" % os.environ.get('PYTHONPATH'))
1519 sys.stderr.write(" ASCENDLIBRARY=%s\n" % os.environ.get('ASCENDLIBRARY'))
1520 sys.stderr.write(" ASCENDSOLVERS=%s\n" % os.environ.get('ASCENDSOLVERS'))
1521 sys.stderr.flush()
1522
1523 import ascpy
1524
1525 try:
1526 import cunit
1527 except:
1528 pass
1529
1530 atexit.register(ascpy.shutdown)
1531 #suite = unittest.TestSuite()
1532 #suite = unittest.defaultTestLoader.loadTestsFromName('__main__')
1533 #unittest.TextTestRunner(verbosity=2).run(suite)
1534 unittest.main()

Properties

Name Value
svn:executable *

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