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

Contents of /trunk/test.py

Parent Directory Parent Directory | Revision Log Revision Log


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

Properties

Name Value
svn:executable *

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