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

Contents of /trunk/test.py

Parent Directory Parent Directory | Revision Log Revision Log


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

Properties

Name Value
svn:executable *

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