/[ascend]/branches/fprops2/test.py
ViewVC logotype

Contents of /branches/fprops2/test.py

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1452 - (show annotations) (download) (as text)
Mon May 28 14:01:51 2007 UTC (12 years, 1 month ago) by jpye
Original Path: trunk/test.py
File MIME type: text/x-python
File size: 42483 byte(s)
Integrators can now be dynamically loaded. DOPRI5 has progressed but still doesn't work.
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])
313 s2 = sorted(['IDA','LSODE'])
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 IndexError:
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(6000,"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.setReporter(ascpy.IntegratorReporterConsole(I))
820 #I.setLinearTimesteps(ascpy.Units("s"), 0,300,300)
821 I.setLogTimesteps(ascpy.Units("s"), 0.009, 1200, 150)
822 I.analyse()
823 F = file('ga.mm','w')
824 I.writeMatrix(F,'dg/dz')
825 F = file('gd.mm','w')
826 I.writeMatrix(F,'dg/dx')
827 F = file('fa.mm','w')
828 I.writeMatrix(F,'df/dz')
829 F = file('fd.mm','w')
830 I.writeMatrix(F,'df/dx')
831 F = file('fdp.mm','w')
832 I.writeMatrix(F,"df/dx'")
833 I.solve()
834
835 #-------------------------------------------------------------------------------
836 # Testing of freesteam external steam properties functions
837
838 with_freesteam = True
839 try:
840 # we assume that if the freesteam python module is installed, the ASCEND
841 # external library will also be.
842 import freesteam
843 have_freesteam = True
844 except ImportError,e:
845 have_freesteam = False
846
847 if with_freesteam and have_freesteam:
848 class TestFreesteam(AscendSelfTester):
849 # def testfreesteamtest(self):
850 # """run the self-test cases bundled with freesteam"""
851 # self._run('testfreesteam',filename='testfreesteam.a4c')
852
853 def testload(self):
854 """check that we can load 'thermalequilibrium2' (IMPORT "freesteam", etc)"""
855 self.L.load('johnpye/thermalequilibrium2.a4c')
856
857 def testinstantiate(self):
858 """load an instantiate 'thermalequilibrium2'"""
859 self.testload()
860 M = self.L.findType('thermalequilibrium2').getSimulation('sim')
861 return M
862
863 def testintegrate(self):
864 """integrate transfer of heat from one mass of water/steam to another
865 according to Newton's law of cooling"""
866 M = self.testinstantiate()
867 M.setSolver(ascpy.Solver("QRSlv"))
868 I = ascpy.Integrator(M)
869 I.setEngine('LSODE')
870 I.setReporter(ascpy.IntegratorReporterConsole(I))
871 I.setLinearTimesteps(ascpy.Units("s"), 0, 3000, 30)
872 I.setMinSubStep(0.01)
873 I.setInitialSubStep(1)
874 I.analyse()
875 print "Number of vars = %d" % I.getNumVars()
876 assert I.getNumVars()==2
877 I.solve()
878 assert I.getNumObservedVars() == 3;
879 print "S[1].T = %f K" % M.S[1].T
880 print "S[2].T = %f K" % M.S[2].T
881 print "Q = %f W" % M.Q
882 self.assertAlmostEqual(float(M.S[1].T),506.77225109,4);
883 self.assertAlmostEqual(float(M.S[2].T),511.605173967,5);
884 self.assertAlmostEqual(float(M.Q),-48.32922877329,3);
885 self.assertAlmostEqual(float(M.t),3000);
886 print "Note that the above values have not been verified analytically"
887
888 def testcollapsingcan2(self):
889 """ solve the collapsing can model using IAPWS-IF97 steam props """
890 M = self._run("collapsingcan2",filename="collapsingcan2.a4c");
891
892 #-------------------------------------------------------------------------------
893 # Testing of the brent-solver EXTERNAL method
894
895 class TestBrent(AscendSelfTester):
896 def testbrent(self):
897 M = self._run('brent1',filename='test/brent.a4c')
898
899 #-------------------------------------------------------------------------------
900 # Testing of IDA's analysis module
901
902 class TestIDA(Ascend):
903 def _run(self,filen,modeln=""):
904 self.L.load('test/ida/%s.a4c' % filen)
905 T = self.L.findType('%s%s' % (filen,modeln))
906 M = T.getSimulation('sim')
907 M.build()
908 I = ascpy.Integrator(M)
909 I.setEngine('IDA')
910 I.analyse()
911 return M;
912
913 def _runfail(self,filen,n,msg="failed"):
914 try:
915 self._run(filen,'fail%d' % n)
916 except Exception,e:
917 print "(EXPECTED) ERROR: %s" % e
918 return
919 self.fail(msg)
920
921 def testsinglederiv(self):
922 self._run('singlederiv')
923
924 def testsinglederivfail1(self):
925 self._runfail('singlederiv',1
926 ,"t.ode_id=-1 did not trigger error")
927
928 def testsinglederivfail2(self):
929 self._runfail('singlederiv',2
930 ,"dy_dt.ode_id=2 did not trigger error")
931
932 def testsinglederivfail3(self):
933 self._runfail('singlederiv',3
934 ,"dy_dt.ode_type=3 did not trigger error")
935
936 def testsinglederivfail4(self):
937 self._runfail('singlederiv',4
938 ,"duplicate ode_type=1 did not trigger error")
939
940 def testsinglederivfail5(self):
941 self._runfail('singlederiv',5
942 ,"duplicate ode_type=1 did not trigger error")
943
944 def testsinglederivfail6(self):
945 self._runfail('singlederiv',6
946 ,"duplicate ode_type=1 did not trigger error")
947
948 def testtwoderiv(self):
949 self._run('twoderiv')
950
951 def testtwoderivfail1(self):
952 self._runfail('twoderiv',1)
953
954 def testtwoderivfail2(self):
955 self._runfail('twoderiv',2)
956
957 def testtwoderivfail3(self):
958 self._runfail('twoderiv',3)
959 def testtwoderivfail4(self):
960 self._runfail('twoderiv',4)
961 def testtwoderivfail5(self):
962 self._runfail('twoderiv',5)
963
964 def testnoderivs(self):
965 self._runfail('noderivs',1)
966
967 def testnoindeps(self):
968 self._runfail('indeps',1)
969
970 def testtwoindeps(self):
971 self._runfail('indeps',2)
972
973 def testfixedvars(self):
974 self._run('fixedvars')
975
976 def testfixedvars1(self):
977 self._run('fixedvars',1)
978
979 # fails the index check
980 # def testfixedvars2(self):
981 # self._run('fixedvars',2)
982
983 # fails the index check
984 # def testfixedvars3(self):
985 # self._run('fixedvars',3)
986
987 def testincidence(self):
988 self._run('incidence')
989
990 def testincidence1(self):
991 self._run('incidence',1)
992 def testincidence2(self):
993 self._run('incidence',2)
994 def testincidence3(self):
995 M = self._run('incidence',3)
996
997 def testincidence4(self):
998 self._run('incidence',4)
999 def testincidencefail5(self):
1000 self._runfail('incidence',5)
1001
1002 def testwritematrix(self):
1003 self.L.load('test/ida/writematrix.a4c')
1004 T = self.L.findType('writematrix')
1005 M = T.getSimulation('sim')
1006 M.build()
1007 I = ascpy.Integrator(M)
1008 I.setEngine('IDA')
1009 I.analyse()
1010 # this stuff fails on Windows because FILE structure is different python vs mingw
1011 # F = os.tmpfile()
1012 # I.writeMatrix(F,"dF/dy")
1013 # F.seek(0)
1014 # print F.read()
1015 # F1 = os.tmpfile()
1016 # I.writeMatrix(F1,"dF/dy'")
1017 # F1.seek(0)
1018 # print F1.read()
1019 # F1 = os.tmpfile()
1020 # I.writeMatrix(F1,"dg/dx")
1021 # F1.seek(0)
1022 # print F1.read()
1023 # # for the moment you'll have to check these results manually.
1024
1025 def testwritematrix2(self):
1026 self.L.load('test/ida/writematrix.a4c')
1027 T = self.L.findType('writematrix2')
1028 M = T.getSimulation('sim')
1029 M.build()
1030 I = ascpy.Integrator(M)
1031 I.setEngine('IDA')
1032 I.analyse()
1033 # this stuff fails on Windows because FILE structure is different python vs mingw
1034 # F = os.tmpfile()
1035 # I.writeMatrix(F,"dF/dy")
1036 # F.seek(0)
1037 # print F.read()
1038 # F1 = os.tmpfile()
1039 # I.writeMatrix(F1,"dF/dy'")
1040 # F1.seek(0)
1041 # print F1.read()
1042 # F1 = os.tmpfile()
1043 # I.writeMatrix(F1,"dg/dx")
1044 # F1.seek(0)
1045 # print F1.read()
1046 #F1 = os.tmpfile()
1047 #I.writeMatrix(F1,"dydp/dyd")
1048 #F1.seek(0)
1049 #print F1.read()
1050 # for the moment you'll have to check these results manually.
1051
1052 def testindexproblem(self):
1053 self.L.load('test/ida/indexproblem.a4c')
1054 T = self.L.findType('indexproblem')
1055 M = T.getSimulation('sim')
1056 M.build()
1057 I = ascpy.Integrator(M)
1058 I.setEngine('IDA')
1059 I.analyse()
1060 pass
1061
1062 def testindexproblem2(self):
1063 self.L.load('test/ida/indexproblem.a4c')
1064 T = self.L.findType('indexproblem2')
1065 M = T.getSimulation('sim')
1066 M.build()
1067 I = ascpy.Integrator(M)
1068 I.setEngine('IDA')
1069 try:
1070 I.analyse()
1071 except Exception,e:
1072 return
1073 self.fail("Index problem not detected")
1074
1075 def testboundaries(self):
1076 self.L.load('test/ida/boundaries.a4c')
1077 T = self.L.findType('boundaries')
1078 M = T.getSimulation('sim')
1079 M.build()
1080 I = ascpy.Integrator(M)
1081 I.setEngine('IDA')
1082 I.analyse()
1083 I.setLogTimesteps(ascpy.Units("s"), 0.1, 20, 20)
1084 I.setParameter('linsolver','DENSE')
1085 I.setParameter('calcic','Y')
1086 I.setParameter('linsolver','DENSE')
1087 I.setParameter('safeeval',False)
1088 I.setReporter(ascpy.IntegratorReporterConsole(I))
1089 I.solve()
1090
1091 # doesn't work yet:
1092 # def testincidence5(self):
1093 # self._run('incidence',5)
1094
1095
1096 #-------------------------------------------------------------------------------
1097 # Testing of IDA models using DENSE linear solver
1098
1099 class TestIDADENSE(Ascend):
1100 """IDA DAE integrator, DENSE linear solver"""
1101
1102 def testlotka(self):
1103 self.L.load('johnpye/lotka.a4c')
1104 M = self.L.findType('lotka').getSimulation('sim')
1105 M.setSolver(ascpy.Solver("QRSlv"))
1106 I = ascpy.Integrator(M)
1107 I.setEngine('IDA')
1108 I.setReporter(ascpy.IntegratorReporterConsole(I))
1109 I.setLinearTimesteps(ascpy.Units("s"), 0, 200, 5);
1110 I.setParameter('linsolver','DENSE')
1111 I.setParameter('rtol',1e-8);
1112 I.analyse()
1113 assert I.getNumVars()==2
1114 assert abs(M.R - 1000) < 1e-300
1115 I.solve()
1116 assert I.getNumObservedVars() == 3
1117 assert abs(M.R - 832) < 1.0
1118 assert abs(M.F - 21.36) < 0.1
1119
1120 def testdenx(self):
1121 print "-----------------------------====="
1122 self.L.load('johnpye/idadenx.a4c')
1123 M = self.L.findType('idadenx').getSimulation('sim')
1124 M.setSolver(ascpy.Solver("QRSlv"))
1125 I = ascpy.Integrator(M)
1126 I.setEngine('IDA')
1127 I.setParameter('calcic','YA_YDP')
1128 I.setParameter('linsolver','DENSE')
1129 I.setParameter('safeeval',True)
1130 I.setReporter(ascpy.IntegratorReporterConsole(I))
1131 I.setLogTimesteps(ascpy.Units("s"), 0.4, 4e10, 11)
1132 I.setMaxSubStep(0);
1133 I.setInitialSubStep(0)
1134 I.setMaxSubSteps(0)
1135 I.setParameter('autodiff',True)
1136 I.analyse()
1137 I.solve()
1138 assert abs(float(M.y1) - 5.1091e-08) < 2e-9
1139 assert abs(float(M.y2) - 2.0437e-13) < 2e-14
1140 assert abs(float(M.y3) - 1.0) < 1e-5
1141
1142 def testhires(self):
1143 self.L.load('test/hires.a4c')
1144 T = self.L.findType('hires')
1145 M = T.getSimulation('sim')
1146 M.setSolver(ascpy.Solver('QRSlv'))
1147 I = ascpy.Integrator(M)
1148 I.setEngine('IDA')
1149 I.setParameter('linsolver','DENSE')
1150 I.setParameter('rtol',1.1e-15)
1151 I.setParameter('atolvect',0)
1152 I.setParameter('atol',1.1e-15)
1153 I.setReporter(ascpy.IntegratorReporterConsole(I))
1154 I.setLogTimesteps(ascpy.Units(""), 1, 321.8122, 5)
1155 I.setInitialSubStep(1e-5)
1156 I.setMaxSubSteps(10000)
1157 I.analyse()
1158 I.solve()
1159 for i in range(8):
1160 print "y[%d] = %.20g" % (i+1, M.y[i+1])
1161 M.run(T.getMethod('self_test'))
1162
1163 def testchemakzo(self):
1164 self.L.load('test/chemakzo.a4c')
1165 T = self.L.findType('chemakzo')
1166 M = T.getSimulation('sim')
1167 M.setSolver(ascpy.Solver('QRSlv'))
1168 I = ascpy.Integrator(M)
1169 I.setEngine('IDA')
1170 I.setParameter('linsolver','DENSE')
1171 I.setParameter('rtol',1e-15)
1172 I.setParameter('atolvect',0)
1173 I.setParameter('atol',1e-15)
1174 I.setReporter(ascpy.IntegratorReporterConsole(I))
1175 I.setLinearTimesteps(ascpy.Units("s"), 1, 180, 5)
1176 I.setInitialSubStep(1e-13)
1177 I.setMaxSubSteps(10000)
1178 I.analyse()
1179 I.solve()
1180 for i in range(6):
1181 print "y[%d] = %.20g" % (i+1, M.y[i+1])
1182 M.run(T.getMethod('self_test'))
1183
1184 def testtransamp(self):
1185 self.L.load('test/transamp.a4c')
1186 T = self.L.findType('transamp')
1187 M = T.getSimulation('sim')
1188 M.setSolver(ascpy.Solver('QRSlv'))
1189 I = ascpy.Integrator(M)
1190 I.setEngine('IDA')
1191 I.setParameter('linsolver','DENSE')
1192 I.setParameter('rtol',1e-7)
1193 I.setParameter('atolvect',0)
1194 I.setParameter('atol',1e-7)
1195 I.setReporter(ascpy.IntegratorReporterConsole(I))
1196 I.setLinearTimesteps(ascpy.Units("s"), 0.05, 0.2, 20)
1197 I.setInitialSubStep(0.00001)
1198 I.setMaxSubSteps(10000)
1199 I.analyse()
1200 I.solve()
1201 for i in range(6):
1202 print "y[%d] = %.20g" % (i+1, M.y[i+1])
1203 M.run(T.getMethod('self_test'))
1204
1205 # MODEL FAILS ANALYSIS: we need to add support for non-incident differential vars
1206 # def testpollution(self):
1207 # self.L.load('test/pollution.a4c')
1208 # T = self.L.findType('pollution')
1209 # M = T.getSimulation('sim')
1210 # M.setSolver(ascpy.Solver('QRSlv'))
1211 # I = ascpy.Integrator(M)
1212 # I.setEngine('IDA')
1213 # I.setParameter('linsolver','DENSE')
1214 # I.setParameter('rtol',1.1e-15)
1215 # I.setParameter('atolvect',0)
1216 # I.setParameter('atol',1.1e-15)
1217 # I.setReporter(ascpy.IntegratorReporterConsole(I))
1218 # I.setLogTimesteps(ascpy.Units("s"), 1, 60, 5)
1219 # I.setInitialSubStep(1e-5)
1220 # I.setMaxSubSteps(10000)
1221 # I.analyse()
1222 # I.solve()
1223 # for i in range(20):
1224 # print "y[%d] = %.20g" % (i+1, M.y[i+1])
1225 # M.run(T.getMethod('self_test'))
1226
1227 ## @TODO fails during IDACalcIC (model too big?)
1228 # def testkryx(self):
1229 # self.L.load('johnpye/idakryx.a4c')
1230 # ascpy.getCompiler().setUseRelationSharing(False)
1231 # M = self.L.findType('idakryx').getSimulation('sim')
1232 # M.setSolver(ascpy.Solver('QRSlv'))
1233 # M.build()
1234 # I = ascpy.Integrator(M)
1235 # I.setEngine('IDA')
1236 # I.setReporter(ascpy.IntegratorReporterConsole(I))
1237 # I.setParameter('linsolver','DENSE')
1238 # I.setParameter('maxl',8)
1239 # I.setParameter('gsmodified',False)
1240 # I.setParameter('autodiff',True)
1241 # I.setParameter('rtol',0)
1242 # I.setParameter('atol',1e-3);
1243 # I.setParameter('atolvect',False)
1244 # I.setParameter('calcic','YA_YDP')
1245 # I.analyse()
1246 # I.setLogTimesteps(ascpy.Units("s"), 0.01, 10.24, 11)
1247 # I.solve()
1248 # assert abs(M.u[2][2].getValue()) < 1e-5
1249
1250 #-------------------------------------------------------------------------------
1251 # Testing of IDA models using SPGMR linear solver (Krylov)
1252
1253 # these tests are disabled until SPGMR preconditioning has been implemented
1254 class TestIDASPGMR:#(Ascend):
1255 def testlotka(self):
1256 self.L.load('johnpye/lotka.a4c')
1257 M = self.L.findType('lotka').getSimulation('sim')
1258 M.setSolver(ascpy.Solver("QRSlv"))
1259 I = ascpy.Integrator(M)
1260 I.setEngine('IDA')
1261 I.setReporter(ascpy.IntegratorReporterConsole(I))
1262 I.setLinearTimesteps(ascpy.Units("s"), 0, 200, 5)
1263 I.setParameter('rtol',1e-8)
1264 I.analyse()
1265 assert I.getNumVars()==2
1266 assert abs(M.R - 1000) < 1e-300
1267 I.solve()
1268 assert I.getNumObservedVars() == 3
1269 assert abs(M.R - 832) < 1.0
1270 assert abs(M.F - 21.36) < 0.1
1271
1272
1273 def testkryx(self):
1274 self.L.load('johnpye/idakryx.a4c')
1275 M = self.L.findType('idakryx').getSimulation('sim')
1276 M.build()
1277 I = ascpy.Integrator(M)
1278 I.setEngine('IDA')
1279 I.setReporter(ascpy.IntegratorReporterConsole(I))
1280 I.setParameter('linsolver','SPGMR')
1281 I.setParameter('prec','JACOBI')
1282 I.setParameter('maxl',8)
1283 I.setParameter('gsmodified',False)
1284 I.setParameter('autodiff',True)
1285 I.setParameter('gsmodified',True)
1286 I.setParameter('rtol',0)
1287 I.setParameter('atol',1e-3);
1288 I.setParameter('atolvect',False)
1289 I.setParameter('calcic','Y')
1290 I.analyse()
1291 I.setLogTimesteps(ascpy.Units("s"), 0.01, 10.24, 10);
1292 print M.udot[1][3]
1293 I.solve()
1294 assert 0
1295
1296 def testzill(self):
1297 self.L.load('johnpye/zill.a4c')
1298 T = self.L.findType('zill')
1299 M = T.getSimulation('sim')
1300 M.setSolver(ascpy.Solver('QRSlv'))
1301 I = ascpy.Integrator(M)
1302 I.setEngine('IDA')
1303 I.setParameter('safeeval',False)
1304 I.setMinSubStep(1e-7)
1305 I.setMaxSubStep(0.001)
1306 I.setMaxSubSteps(10000)
1307 I.setReporter(ascpy.IntegratorReporterConsole(I))
1308 I.setLinearTimesteps(ascpy.Units(), 1.0, 1.5, 5)
1309 I.analyse()
1310 I.solve()
1311 M.run(T.getMethod('self_test'))
1312
1313 def testdenxSPGMR(self):
1314 self.L.load('johnpye/idadenx.a4c')
1315 M = self.L.findType('idadenx').getSimulation('sim')
1316 M.setSolver(ascpy.Solver('QRSlv'))
1317 I = ascpy.Integrator(M)
1318 I.setEngine('IDA')
1319 I.setReporter(ascpy.IntegratorReporterConsole(I))
1320 I.setLogTimesteps(ascpy.Units("s"), 0.4, 4e10, 11)
1321 I.setMaxSubStep(0);
1322 I.setInitialSubStep(0);
1323 I.setMaxSubSteps(0);
1324 I.setParameter('autodiff',True)
1325 I.setParameter('linsolver','SPGMR')
1326 I.setParameter('gsmodified',False)
1327 I.setParameter('maxncf',10)
1328 I.analyse()
1329 I.solve()
1330 assert abs(float(M.y1) - 5.1091e-08) < 1e-10
1331 assert abs(float(M.y2) - 2.0437e-13) < 1e-15
1332 assert abs(float(M.y3) - 1.0) < 1e-5
1333
1334 # move code above down here if you want to temporarily avoid testing it
1335 class NotToBeTested:
1336 def nothing(self):
1337 pass
1338
1339 def testnewton(self):
1340 sys.stderr.write("STARTING TESTNEWTON\n")
1341 self.L.load('johnpye/newton.a4c')
1342 T = self.L.findType('newton')
1343 M = T.getSimulation('sim')
1344 M.solve(ascpy.Solver("QRSlv"),ascpy.SolverReporter())
1345 I = ascpy.Integrator(M)
1346 I.setEngine('IDA')
1347 I.setParameter('linsolver','DENSE')
1348 I.setParameter('safeeval',True)
1349 I.setParameter('rtol',1e-8)
1350 I.setMaxSubStep(0.001)
1351 I.setMaxSubSteps(10000)
1352
1353 I.setReporter(ascpy.IntegratorReporterConsole(I))
1354 I.setLinearTimesteps(ascpy.Units("s"), 0, 2*float(M.v)/float(M.g), 2)
1355 I.analyse()
1356 I.solve()
1357 print "At end of simulation,"
1358 print "x = %f" % M.x
1359 print "v = %f" % M.v
1360 M.run(T.getMethod('self_test'))
1361
1362 if __name__=='__main__':
1363 # a whole bag of tricks to make sure we get the necessary dirs in our ascend, python and ld path vars
1364 restart = 0
1365
1366 if platform.system()=="Windows":
1367 LD_LIBRARY_PATH="PATH"
1368 SEP = ";"
1369 else:
1370 LD_LIBRARY_PATH="LD_LIBRARY_PATH"
1371 SEP = ":"
1372
1373 freesteamdir = os.path.expanduser("~/freesteam/ascend")
1374 modeldirs = [os.path.abspath(os.path.join(sys.path[0],"models")),os.path.abspath(freesteamdir)]
1375 if not os.environ.get('ASCENDLIBRARY'):
1376 os.environ['ASCENDLIBRARY'] = SEP.join(modeldirs)
1377 restart = 1
1378 else:
1379 envmodelsdir = [os.path.abspath(i) for i in os.environ['ASCENDLIBRARY'].split(SEP)]
1380 for l in modeldirs:
1381 if l in envmodelsdir[len(modeldirs):]:
1382 envmodelsdir.remove(l)
1383 restart = 1
1384 for l in modeldirs:
1385 if l not in envmodelsdir:
1386 envmodelsdir.insert(0,l)
1387 restart = 1
1388 os.environ['ASCENDLIBRARY'] = SEP.join(envmodelsdir)
1389
1390 libdirs = ["pygtk","."]
1391 libdirs = [os.path.normpath(os.path.join(sys.path[0],l)) for l in libdirs]
1392 if not os.environ.get(LD_LIBRARY_PATH):
1393 os.environ[LD_LIBRARY_PATH]=SEP.join(libdirs)
1394 restart = 1
1395 else:
1396 envlibdirs = [os.path.normpath(i) for i in os.environ[LD_LIBRARY_PATH].split(SEP)]
1397 for l in libdirs:
1398 if l in envlibdirs[len(libdirs):]:
1399 envlibdirs.remove(l)
1400 restart = 1
1401 for l in libdirs:
1402 if l not in envlibdirs:
1403 envlibdirs.insert(0,l)
1404 restart = 1
1405 os.environ[LD_LIBRARY_PATH] = SEP.join(envlibdirs)
1406
1407 pypath = os.path.normpath(os.path.join(sys.path[0],"pygtk"))
1408 if not os.environ.get('PYTHONPATH'):
1409 os.environ['PYTHONPATH']=pypath
1410 else:
1411 envpypath = os.environ['PYTHONPATH'].split(SEP)
1412 if pypath not in envpypath:
1413 envpypath.insert(0,pypath)
1414 os.environ['PYTHONPATH']=SEP.join(envpypath)
1415 restart = 1
1416
1417 if restart:
1418 if platform.system()=="Windows":
1419 pass
1420 else:
1421 script = os.path.join(sys.path[0],"test.py")
1422 sys.stderr.write("Restarting with...\n")
1423 sys.stderr.write(" export LD_LIBRARY_PATH=%s\n" % os.environ.get(LD_LIBRARY_PATH))
1424 sys.stderr.write(" export PYTHONPATH=%s\n" % os.environ.get('PYTHONPATH'))
1425 sys.stderr.write(" export ASCENDLIBRARY=%s\n" % os.environ.get('ASCENDLIBRARY'))
1426 sys.stderr.flush()
1427 os.execvp("python",[script] + sys.argv)
1428
1429 import ascpy
1430
1431 try:
1432 import cunit
1433 except:
1434 pass
1435
1436 atexit.register(ascpy.shutdown)
1437 #suite = unittest.TestSuite()
1438 #suite = unittest.defaultTestLoader.loadTestsFromName('__main__')
1439 #unittest.TextTestRunner(verbosity=2).run(suite)
1440 unittest.main()

Properties

Name Value
svn:executable *

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