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

Contents of /trunk/test.py

Parent Directory Parent Directory | Revision Log Revision Log


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

Properties

Name Value
svn:executable *

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