1 |
#!/usr/bin/env python |
2 |
import unittest |
3 |
import ascpy |
4 |
import math |
5 |
import os, subprocess, sys |
6 |
import atexit |
7 |
import cunit |
8 |
|
9 |
class Ascend(unittest.TestCase): |
10 |
|
11 |
def setUp(self): |
12 |
import ascpy |
13 |
self.L = ascpy.Library() |
14 |
|
15 |
def tearDown(self): |
16 |
self.L.clear() |
17 |
del self.L |
18 |
|
19 |
class AscendSelfTester(Ascend): |
20 |
|
21 |
def _run(self,modelname,solvername="QRSlv",filename=None): |
22 |
if filename==None: |
23 |
filename = 'johnpye/%s.a4c' % modelname |
24 |
self.L.load(filename) |
25 |
T = self.L.findType(modelname) |
26 |
M = T.getSimulation('sim') |
27 |
M.build() |
28 |
M.solve(ascpy.Solver(solvername),ascpy.SolverReporter()) |
29 |
M.run(T.getMethod('self_test')) |
30 |
return M |
31 |
|
32 |
class TestCompiler(Ascend): |
33 |
|
34 |
def testloading(self): |
35 |
pass |
36 |
|
37 |
def testsystema4l(self): |
38 |
self.L.load('system.a4l') |
39 |
|
40 |
def testatomsa4l(self): |
41 |
self.L.load('atoms.a4l') |
42 |
|
43 |
class TestSolver(AscendSelfTester): |
44 |
|
45 |
def testlog10(self): |
46 |
self._run('testlog10') |
47 |
|
48 |
def testconopt(self): |
49 |
self._run('testconopt',"CONOPT") |
50 |
|
51 |
def testcmslv2(self): |
52 |
self._run('testcmslv2',"CMSlv") |
53 |
|
54 |
def testsunpos1(self): |
55 |
self._run('example_1_6_1',"QRSlv","johnpye/sunpos.a4c") |
56 |
|
57 |
def testsunpos2(self): |
58 |
self._run('example_1_6_2',"QRSlv","johnpye/sunpos.a4c") |
59 |
|
60 |
def testsunpos3(self): |
61 |
self._run('example_1_7_1',"QRSlv","johnpye/sunpos.a4c") |
62 |
|
63 |
def testsunpos4(self): |
64 |
self._run('example_1_7_2',"QRSlv","johnpye/sunpos.a4c") |
65 |
|
66 |
def testsunpos5(self): |
67 |
self._run('example_1_7_3',"QRSlv","johnpye/sunpos.a4c") |
68 |
|
69 |
def testsunpos6(self): |
70 |
self._run('example_1_8_1',"QRSlv","johnpye/sunpos.a4c") |
71 |
|
72 |
class TestIntegrator(Ascend): |
73 |
|
74 |
def testListIntegrators(self): |
75 |
I = ascpy.Integrator.getEngines() |
76 |
s1 = sorted([str(i) for i in I.values()]) |
77 |
s2 = sorted(['IDA','LSODE','AWW']) |
78 |
assert s1==s2 |
79 |
|
80 |
# this routine is reused by both testIDA and testLSODE |
81 |
def _testIntegrator(self,integratorname): |
82 |
self.L.load('johnpye/shm.a4c') |
83 |
M = self.L.findType('shm').getSimulation('sim') |
84 |
M.setSolver(ascpy.Solver('QRSlv')) |
85 |
print M.getChildren() |
86 |
assert float(M.x) == 10.0 |
87 |
assert float(M.v) == 0.0 |
88 |
t_end = math.pi |
89 |
|
90 |
I = ascpy.Integrator(M) |
91 |
I.setReporter(ascpy.IntegratorReporterNull(I)) |
92 |
I.setEngine(integratorname); |
93 |
I.setLinearTimesteps(ascpy.Units("s"), 0.0, t_end, 100); |
94 |
I.setMinSubStep(0.0005); # these limits are required by IDA at present (numeric diff) |
95 |
I.setMaxSubStep(0.02); |
96 |
I.setInitialSubStep(0.001); |
97 |
I.setMaxSubSteps(200); |
98 |
if(integratorname=='IDA'): |
99 |
I.setParameter('autodiff',False) |
100 |
I.analyse(); |
101 |
I.solve(); |
102 |
print "At end of simulation," |
103 |
print "x = %f" % M.x |
104 |
print "v = %f" % M.v |
105 |
assert abs(float(M.x) + 10) < 1e-2 |
106 |
assert abs(float(M.v)) < 1e-2 |
107 |
assert I.getNumObservedVars() == 3 |
108 |
|
109 |
def testInvalidIntegrator(self): |
110 |
self.L.load('johnpye/shm.a4c') |
111 |
M = self.L.findType('shm').getSimulation('sim') |
112 |
M.setSolver(ascpy.Solver('QRSlv')) |
113 |
I = ascpy.Integrator(M) |
114 |
try: |
115 |
I.setEngine('___NONEXISTENT____') |
116 |
except RuntimeError: |
117 |
return |
118 |
self.fail("setEngine did not raise error!") |
119 |
|
120 |
def testLSODE(self): |
121 |
self._testIntegrator('LSODE') |
122 |
|
123 |
def testIDA(self): |
124 |
self._testIntegrator('IDA') |
125 |
|
126 |
def testparameters(self): |
127 |
self.L.load('johnpye/shm.a4c') |
128 |
M = self.L.findType('shm').getSimulation('sim') |
129 |
M.build() |
130 |
I = ascpy.Integrator(M) |
131 |
I.setEngine('IDA') |
132 |
P = I.getParameters() |
133 |
for p in P: |
134 |
print p.getName(),"=",p.getValue() |
135 |
assert len(P)==11 |
136 |
assert P[0].isStr() |
137 |
assert P[0].getName()=="linsolver" |
138 |
assert P[0].getValue()=='SPGMR' |
139 |
assert P[2].getName()=="autodiff" |
140 |
assert P[2].getValue()==True |
141 |
assert P[7].getName()=="atolvect" |
142 |
assert P[7].getBoolValue() == True |
143 |
P[2].setBoolValue(False) |
144 |
assert P[2].getBoolValue()==False |
145 |
I.setParameters(P) |
146 |
assert I.getParameterValue('autodiff')==False |
147 |
I.setParameter('autodiff',True) |
148 |
try: |
149 |
v = I.getParameterValue('nonexist') |
150 |
except KeyError: |
151 |
pass |
152 |
else: |
153 |
self.fail('Failed to trip invalid Integrator parameter') |
154 |
|
155 |
class TestLSODE(Ascend): |
156 |
|
157 |
def testzill(self): |
158 |
self.L.load('johnpye/zill.a4c') |
159 |
T = self.L.findType('zill') |
160 |
M = T.getSimulation('sim') |
161 |
M.setSolver(ascpy.Solver('QRSlv')) |
162 |
I = ascpy.Integrator(M) |
163 |
I.setEngine('LSODE') |
164 |
I.setMinSubStep(1e-7) |
165 |
I.setMaxSubStep(0.001) |
166 |
I.setMaxSubSteps(10000) |
167 |
I.setReporter(ascpy.IntegratorReporterConsole(I)) |
168 |
I.setLinearTimesteps(ascpy.Units(), 1.0, 1.5, 5) |
169 |
I.analyse() |
170 |
I.solve() |
171 |
M.run(T.getMethod('self_test')) |
172 |
|
173 |
def testnewton(self): |
174 |
sys.stderr.write("STARTING TESTNEWTON\n") |
175 |
self.L.load('johnpye/newton.a4c') |
176 |
T = self.L.findType('newton') |
177 |
M = T.getSimulation('sim') |
178 |
M.solve(ascpy.Solver("QRSlv"),ascpy.SolverReporter()) |
179 |
I = ascpy.Integrator(M) |
180 |
I.setEngine('LSODE') |
181 |
I.setParameter('rtolvect',False) |
182 |
I.setParameter('rtol',1e-7) |
183 |
I.setParameter('atolvect',False) |
184 |
I.setParameter('atol',1e-7) |
185 |
I.setMinSubStep(1e-7) |
186 |
I.setMaxSubStep(0.001) |
187 |
I.setMaxSubSteps(10000) |
188 |
|
189 |
I.setReporter(ascpy.IntegratorReporterConsole(I)) |
190 |
I.setLinearTimesteps(ascpy.Units("s"), 0, 2*float(M.v)/float(M.g), 2) |
191 |
I.analyse() |
192 |
I.solve() |
193 |
print "At end of simulation," |
194 |
print "x = %f" % M.x |
195 |
print "v = %f" % M.v |
196 |
M.run(T.getMethod('self_test')) |
197 |
|
198 |
def testlotka(self): |
199 |
self.L.load('johnpye/lotka.a4c') |
200 |
M = self.L.findType('lotka').getSimulation('sim') |
201 |
M.setSolver(ascpy.Solver("QRSlv")) |
202 |
I = ascpy.Integrator(M) |
203 |
I.setEngine('LSODE') |
204 |
I.setReporter(ascpy.IntegratorReporterConsole(I)) |
205 |
I.setLinearTimesteps(ascpy.Units("s"), 0, 200, 5) |
206 |
I.analyse() |
207 |
print "Number of vars = %d" % I.getNumVars() |
208 |
assert I.getNumVars()==2 |
209 |
I.solve() |
210 |
assert I.getNumObservedVars() == 3; |
211 |
assert abs(M.R - 832) < 1.0 |
212 |
assert abs(M.F - 21.36) < 0.1 |
213 |
|
214 |
#------------------------------------------------------------------------------- |
215 |
# Testing of a simple external function |
216 |
|
217 |
class TestExtFn(AscendSelfTester): |
218 |
def testextfntest(self): |
219 |
M = self._run('extfntest',filename='johnpye/extfn/extfntest.a4c') |
220 |
self.assertAlmostEqual(M.y, 2); |
221 |
self.assertAlmostEqual(M.x, 1); |
222 |
self.assertAlmostEqual(M.y, M.x + 1); |
223 |
|
224 |
# THIS TEST FAILS |
225 |
# def testextrelfor(self): |
226 |
# self.L.load('johnpye/extfn/extrelfor.a4c') |
227 |
# T = self.L.findType('extrelfor') |
228 |
# M = T.getSimulation('sim') |
229 |
# M.solve(ascpy.Solver('QRSlv'),ascpy.SolverReporter()) |
230 |
# print "x[1] = %f" % M.x[1] |
231 |
# print "x[2] = %f" % M.x[2] |
232 |
# print "x[3] = %f" % M.x[3] |
233 |
# M.run(T.getMethod('self_test')) |
234 |
|
235 |
def testextrelrepeat(self): |
236 |
M = self._run('extrelrepeat',filename='johnpye/extfn/extrelrepeat.a4c') |
237 |
|
238 |
#------------------------------------------------------------------------------- |
239 |
# Testing of a ExtPy - external python methods |
240 |
|
241 |
class TestExtPy(AscendSelfTester): |
242 |
def testextpytest(self): |
243 |
print "-------------------=--=-=-=-" |
244 |
M = self._run('extpytest',filename='johnpye/extpy/extpytest.a4c') |
245 |
|
246 |
#------------------------------------------------------------------------------- |
247 |
# Testing of freesteam external steam properties functions |
248 |
|
249 |
with_freesteam = False |
250 |
try: |
251 |
import freesteam |
252 |
have_freesteam = True |
253 |
except ImportError,e: |
254 |
have_freesteam = False |
255 |
|
256 |
if with_freesteam and have_freesteam: |
257 |
class TestFreesteam(Ascend): |
258 |
def testload(self): |
259 |
self.L.load('johnpye/thermalequilibrium2.a4c') |
260 |
|
261 |
def testinstantiate(self): |
262 |
self.testload() |
263 |
M = self.L.findType('thermalequilibrium2').getSimulation('sim') |
264 |
|
265 |
def testsolve(self): |
266 |
self.testinstantiate() |
267 |
M.setSolver(ascpy.Solver("QRSlv")) |
268 |
#I = ascpy.Integrator(M) |
269 |
#I.setEngine('LSODE') |
270 |
#I.setReporter(ascpy.IntegratorReporterConsole(I)) |
271 |
#I.setLinearTimesteps(ascpy.Units("s"), 0, 3000, 30) |
272 |
#I.setMinSubStep(0.01) |
273 |
#I.setInitialSubStep(0.1) |
274 |
#I.analyse() |
275 |
#print "Number of vars = %d" % I.getNumVars() |
276 |
#assert I.getNumVars()==2 |
277 |
#I.solve() |
278 |
#assert I.getNumObservedVars() == 3; |
279 |
#assert abs(M.R - 832) < 1.0 |
280 |
#assert abs(M.F - 21.36) < 0.1 |
281 |
|
282 |
|
283 |
#------------------------------------------------------------------------------- |
284 |
# Testing of IDA models using DENSE linear solver |
285 |
|
286 |
class TestIDADENSE(Ascend): |
287 |
"""IDA DAE integrator, DENSE linear solver""" |
288 |
|
289 |
def testnewton(self): |
290 |
sys.stderr.write("STARTING TESTNEWTON\n") |
291 |
self.L.load('johnpye/newton.a4c') |
292 |
T = self.L.findType('newton') |
293 |
M = T.getSimulation('sim') |
294 |
M.solve(ascpy.Solver("QRSlv"),ascpy.SolverReporter()) |
295 |
I = ascpy.Integrator(M) |
296 |
I.setEngine('IDA') |
297 |
I.setParameter('linsolver','DENSE') |
298 |
I.setParameter('safeeval',True) |
299 |
I.setParameter('rtol',1e-8) |
300 |
I.setMaxSubStep(0.001) |
301 |
I.setMaxSubSteps(10000) |
302 |
|
303 |
I.setReporter(ascpy.IntegratorReporterConsole(I)) |
304 |
I.setLinearTimesteps(ascpy.Units("s"), 0, 2*float(M.v)/float(M.g), 2) |
305 |
I.analyse() |
306 |
I.solve() |
307 |
print "At end of simulation," |
308 |
print "x = %f" % M.x |
309 |
print "v = %f" % M.v |
310 |
M.run(T.getMethod('self_test')) |
311 |
|
312 |
def testlotkaDENSE(self): |
313 |
self.L.load('johnpye/lotka.a4c') |
314 |
M = self.L.findType('lotka').getSimulation('sim') |
315 |
M.setSolver(ascpy.Solver("QRSlv")) |
316 |
I = ascpy.Integrator(M) |
317 |
I.setEngine('IDA') |
318 |
I.setReporter(ascpy.IntegratorReporterConsole(I)) |
319 |
I.setLinearTimesteps(ascpy.Units("s"), 0, 200, 5); |
320 |
I.setParameter('linsolver','DENSE') |
321 |
I.setParameter('rtol',1e-8); |
322 |
I.analyse() |
323 |
assert I.getNumVars()==2 |
324 |
assert abs(M.R - 1000) < 1e-300 |
325 |
I.solve() |
326 |
assert I.getNumObservedVars() == 3 |
327 |
assert abs(M.R - 832) < 1.0 |
328 |
assert abs(M.F - 21.36) < 0.1 |
329 |
|
330 |
def testdenx(self): |
331 |
print "-----------------------------=====" |
332 |
self.L.load('johnpye/idadenx.a4c') |
333 |
M = self.L.findType('idadenx').getSimulation('sim') |
334 |
M.setSolver(ascpy.Solver("QRSlv")) |
335 |
I = ascpy.Integrator(M) |
336 |
I.setEngine('IDA') |
337 |
I.setParameter('calcic','YA_YPD') |
338 |
I.setParameter('linsolver','DENSE') |
339 |
I.setReporter(ascpy.IntegratorReporterConsole(I)) |
340 |
I.setLogTimesteps(ascpy.Units("s"), 0.4, 4e10, 11) |
341 |
I.setMaxSubStep(0); |
342 |
I.setInitialSubStep(0) |
343 |
I.setMaxSubSteps(0); |
344 |
I.setParameter('autodiff',True) |
345 |
I.analyse() |
346 |
I.solve() |
347 |
assert abs(float(M.y1) - 5.1091e-08) < 2e-9 |
348 |
assert abs(float(M.y2) - 2.0437e-13) < 2e-14 |
349 |
assert abs(float(M.y3) - 1.0) < 1e-5 |
350 |
|
351 |
def testkryxDENSE(self): |
352 |
self.L.load('johnpye/idakryx.a4c') |
353 |
M = self.L.findType('idakryx').getSimulation('sim') |
354 |
M.setSolver(ascpy.Solver('QRSlv')) |
355 |
M.build() |
356 |
I = ascpy.Integrator(M) |
357 |
I.setEngine('IDA') |
358 |
I.setReporter(ascpy.IntegratorReporterConsole(I)) |
359 |
I.setParameter('linsolver','DENSE') |
360 |
I.setParameter('maxl',8) |
361 |
I.setParameter('gsmodified',False) |
362 |
I.setParameter('autodiff',True) |
363 |
I.setParameter('rtol',0) |
364 |
I.setParameter('atol',1e-3); |
365 |
I.setParameter('atolvect',False) |
366 |
I.setParameter('calcic','YA_YDP') |
367 |
I.analyse() |
368 |
I.setLogTimesteps(ascpy.Units("s"), 0.01, 10.24, 11) |
369 |
I.solve() |
370 |
assert abs(M.u[2][2].getValue()) < 1e-5 |
371 |
|
372 |
#------------------------------------------------------------------------------- |
373 |
# Testing of IDA models using SPGMR linear solver (Krylov) |
374 |
|
375 |
# these tests are disabled until SPGMR preconditioning has been implemented |
376 |
class TestIDASPGMR:#(Ascend): |
377 |
def testlotka(self): |
378 |
self.L.load('johnpye/lotka.a4c') |
379 |
M = self.L.findType('lotka').getSimulation('sim') |
380 |
M.setSolver(ascpy.Solver("QRSlv")) |
381 |
I = ascpy.Integrator(M) |
382 |
I.setEngine('IDA') |
383 |
I.setReporter(ascpy.IntegratorReporterConsole(I)) |
384 |
I.setLinearTimesteps(ascpy.Units("s"), 0, 200, 5) |
385 |
I.setParameter('rtol',1e-8) |
386 |
I.analyse() |
387 |
assert I.getNumVars()==2 |
388 |
assert abs(M.R - 1000) < 1e-300 |
389 |
I.solve() |
390 |
assert I.getNumObservedVars() == 3 |
391 |
assert abs(M.R - 832) < 1.0 |
392 |
assert abs(M.F - 21.36) < 0.1 |
393 |
|
394 |
|
395 |
def testkryx(self): |
396 |
self.L.load('johnpye/idakryx.a4c') |
397 |
M = self.L.findType('idakryx').getSimulation('sim') |
398 |
M.build() |
399 |
I = ascpy.Integrator(M) |
400 |
I.setEngine('IDA') |
401 |
I.setReporter(ascpy.IntegratorReporterConsole(I)) |
402 |
I.setParameter('linsolver','SPGMR') |
403 |
I.setParameter('prec','JACOBI') |
404 |
I.setParameter('maxl',8) |
405 |
I.setParameter('gsmodified',False) |
406 |
I.setParameter('autodiff',True) |
407 |
I.setParameter('gsmodified',True) |
408 |
I.setParameter('rtol',0) |
409 |
I.setParameter('atol',1e-3); |
410 |
I.setParameter('atolvect',False) |
411 |
I.setParameter('calcic','Y') |
412 |
I.analyse() |
413 |
I.setLogTimesteps(ascpy.Units("s"), 0.01, 10.24, 10); |
414 |
print M.udot[1][3] |
415 |
I.solve() |
416 |
assert 0 |
417 |
|
418 |
def testzill(self): |
419 |
self.L.load('johnpye/zill.a4c') |
420 |
T = self.L.findType('zill') |
421 |
M = T.getSimulation('sim') |
422 |
M.setSolver(ascpy.Solver('QRSlv')) |
423 |
I = ascpy.Integrator(M) |
424 |
I.setEngine('IDA') |
425 |
I.setParameter('safeeval',False) |
426 |
I.setMinSubStep(1e-7) |
427 |
I.setMaxSubStep(0.001) |
428 |
I.setMaxSubSteps(10000) |
429 |
I.setReporter(ascpy.IntegratorReporterConsole(I)) |
430 |
I.setLinearTimesteps(ascpy.Units(), 1.0, 1.5, 5) |
431 |
I.analyse() |
432 |
I.solve() |
433 |
M.run(T.getMethod('self_test')) |
434 |
|
435 |
def testdenxSPGMR(self): |
436 |
self.L.load('johnpye/idadenx.a4c') |
437 |
M = self.L.findType('idadenx').getSimulation('sim') |
438 |
M.setSolver(ascpy.Solver('QRSlv')) |
439 |
I = ascpy.Integrator(M) |
440 |
I.setEngine('IDA') |
441 |
I.setReporter(ascpy.IntegratorReporterConsole(I)) |
442 |
I.setLogTimesteps(ascpy.Units("s"), 0.4, 4e10, 11) |
443 |
I.setMaxSubStep(0); |
444 |
I.setInitialSubStep(0); |
445 |
I.setMaxSubSteps(0); |
446 |
I.setParameter('autodiff',True) |
447 |
I.setParameter('linsolver','SPGMR') |
448 |
I.setParameter('gsmodified',False) |
449 |
I.setParameter('maxncf',10) |
450 |
I.analyse() |
451 |
I.solve() |
452 |
assert abs(float(M.y1) - 5.1091e-08) < 1e-10 |
453 |
assert abs(float(M.y2) - 2.0437e-13) < 1e-15 |
454 |
assert abs(float(M.y3) - 1.0) < 1e-5 |
455 |
|
456 |
# move code above down here if you want to temporarily avoid testing it |
457 |
class NotToBeTested: |
458 |
def nothing(self): |
459 |
pass |
460 |
|
461 |
if __name__=='__main__': |
462 |
atexit.register(ascpy.shutdown) |
463 |
#suite = unittest.TestSuite() |
464 |
#suite = unittest.defaultTestLoader.loadTestsFromName('__main__') |
465 |
#unittest.TextTestRunner(verbosity=2).run(suite) |
466 |
unittest.main() |