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

Contents of /trunk/test.py

Parent Directory Parent Directory | Revision Log Revision Log


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

Properties

Name Value
svn:executable *

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