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