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

Diff of /branches/fprops2/test.py

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 1247 by johnpye, Sat Jan 27 00:11:34 2007 UTC revision 1562 by jpye, Mon Jul 30 13:22:43 2007 UTC
# Line 1  Line 1 
1  #!/usr/bin/env python  #!/usr/bin/env python
2  #   ASCEND modelling environment  #   ASCEND modelling environment
3  #   Copyright (C) 2006 Carnegie Mellon University  #   Copyright (C) 2006, 2007 Carnegie Mellon University
4  #  #
5  #   This program is free software; you can redistribute it and/or modify  #   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  #   it under the terms of the GNU General Public License as published by
# Line 28  import atexit Line 28  import atexit
28    
29  import platform  import platform
30  if platform.system() != "Windows":  if platform.system() != "Windows":
31      import dl      try:
32      sys.setdlopenflags(dl.RTLD_GLOBAL|dl.RTLD_NOW)          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):  class Ascend(unittest.TestCase):
43    
# Line 43  class Ascend(unittest.TestCase): Line 51  class Ascend(unittest.TestCase):
51    
52  class AscendSelfTester(Ascend):  class AscendSelfTester(Ascend):
53    
54      def _run(self,modelname,solvername="QRSlv",filename=None):      def _run(self,modelname,solvername="QRSlv",filename=None,parameters={}):
55          if filename==None:          if filename==None:
56              filename = 'johnpye/%s.a4c' % modelname              filename = 'johnpye/%s.a4c' % modelname
57          self.L.load(filename)          self.L.load(filename)
58          T = self.L.findType(modelname)          T = self.L.findType(modelname)
59          M = T.getSimulation('sim')          M = T.getSimulation('sim')
60          M.build()          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())              M.solve(ascpy.Solver(solvername),ascpy.SolverReporter())    
64          M.run(T.getMethod('self_test'))          M.run(T.getMethod('self_test'))
65          return M          return M
# Line 91  class TestCompiler(Ascend): Line 101  class TestCompiler(Ascend):
101          """flagging a subsidiary missing REQUIRE"""          """flagging a subsidiary missing REQUIRE"""
102          self._runfail('missingreq',1)          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):  class TestSolver(AscendSelfTester):
178            
179      def testlog10(self):      def testlog10(self):
180          self._run('testlog10')          M = self._run('testlog10')
181    
182      def testrootsofpoly(self):      def testrootsofpoly(self):
183          self._run('roots_of_poly',filename="roots_of_poly.a4c")          self._run('roots_of_poly',filename="roots_of_poly.a4c")
# Line 106  class TestSolver(AscendSelfTester): Line 189  class TestSolver(AscendSelfTester):
189          self._run('distance_calc',filename="distance_calc.a4c")          self._run('distance_calc',filename="distance_calc.a4c")
190    
191      def testconopt(self):      def testconopt(self):
192          self._run('testconopt',"CONOPT")                          self._run('conopttest',"CONOPT",filename="conopttest.a4c")              
193    
194      def testcmslv2(self):      def testcmslv2(self):
195          self._run('testcmslv2',"CMSlv")          self._run('testcmslv2',"CMSlv")
# Line 134  class TestSolver(AscendSelfTester): Line 217  class TestSolver(AscendSelfTester):
217          self.assertAlmostEqual( float(M.t_solar), M.t_solar.as("s"))          self.assertAlmostEqual( float(M.t_solar), M.t_solar.as("s"))
218          self.assertAlmostEqual( float(M.t_solar)/3600, M.t_solar.as("h"))          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):  class TestMatrix(AscendSelfTester):
294      def testlog10(self):      def testlog10(self):
295          M = self._run('testlog10')          M = self._run('testlog10')
296          print M.getMatrix().write(sys.stderr,"mmio")          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):  class TestIntegrator(Ascend):
309    
310      def testListIntegrators(self):      def testListIntegrators(self):
311          I = ascpy.Integrator.getEngines()          I = ascpy.Integrator.getEngines()
312          s1 = sorted([str(i) for i in I.values()])          s1 = sorted([str(i) for i in I])
313          s2 = sorted(['IDA','LSODE','AWW'])          assert 'IDA' in s1
314          assert s1==s2          assert 'LSODE' in s1
315    
316      # this routine is reused by both testIDA and testLSODE      # this routine is reused by both testIDA and testLSODE
317      def _testIntegrator(self,integratorname):      def _testIntegrator(self,integratorname):
# Line 188  class TestIntegrator(Ascend): Line 353  class TestIntegrator(Ascend):
353          I = ascpy.Integrator(M)          I = ascpy.Integrator(M)
354          try:          try:
355              I.setEngine('___NONEXISTENT____')              I.setEngine('___NONEXISTENT____')
356          except RuntimeError:          except IndexError:
357              return              return
358          self.fail("setEngine did not raise error!")          self.fail("setEngine did not raise error!")
359    
# Line 207  class TestIntegrator(Ascend): Line 372  class TestIntegrator(Ascend):
372          P = I.getParameters()          P = I.getParameters()
373          for p in P:          for p in P:
374              print p.getName(),"=",p.getValue()              print p.getName(),"=",p.getValue()
375          assert len(P)==11          assert len(P)==12
376          assert P[0].isStr()          assert P[0].isStr()
377          assert P[0].getName()=="linsolver"          assert P[0].getName()=="linsolver"
378          assert P[0].getValue()=='SPGMR'          assert P[0].getValue()=='DENSE'
379          assert P[2].getName()=="autodiff"          assert P[2].getName()=="maxord"
380          assert P[2].getValue()==True          assert P[3].getName()=="autodiff"
381          assert P[7].getName()=="atolvect"          assert P[3].getValue()==True
382          assert P[7].getBoolValue() == True          assert P[8].getName()=="atolvect"
383          P[2].setBoolValue(False)          assert P[8].getBoolValue() == True
384          assert P[2].getBoolValue()==False          P[3].setBoolValue(False)
385            assert P[3].getBoolValue()==False
386          I.setParameters(P)          I.setParameters(P)
387          assert I.getParameterValue('autodiff')==False          assert I.getParameterValue('autodiff')==False
388          I.setParameter('autodiff',True)          I.setParameter('autodiff',True)
# Line 286  class TestLSODE(Ascend): Line 452  class TestLSODE(Ascend):
452          assert abs(M.R - 832) < 1.0          assert abs(M.R - 832) < 1.0
453          assert abs(M.F - 21.36) < 0.1          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  # Testing of a external blackbox functions
465    
# Line 551  class TestSteam(AscendSelfTester): Line 725  class TestSteam(AscendSelfTester):
725          M.solve(ascpy.Solver('QRSlv'),ascpy.SolverReporter())          M.solve(ascpy.Solver('QRSlv'),ascpy.SolverReporter())
726          return M          return M
727    
728      def testintegLSODE(self):      def testdsgsatrepeat(self):
729          M = self.testdsgsat()          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")          M.qdot_s.setRealValueWithUnits(1000,"W/m")
747          M.solve(ascpy.Solver('QRSlv'),ascpy.SolverReporter())          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('          #M.setParameter('
758          I = ascpy.Integrator(M)          I = ascpy.Integrator(M)
759          I.setEngine('LSODE')          I.setEngine('LSODE')
         I.setParameter('meth','AM')  
         I.setParameter('maxord',12)  
760          I.setReporter(ascpy.IntegratorReporterConsole(I))          I.setReporter(ascpy.IntegratorReporterConsole(I))
761          I.setLinearTimesteps(ascpy.Units("s"), 0, 5, 1)          I.setLinearTimesteps(ascpy.Units("s"), 0, 3600, 10)
762          I.analyse()          I.analyse()
763          I.solve()          I.solve()
764    
765      def testintegIDA(self):  #   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()          M = self.testdsgsat()
780          self.assertAlmostEqual(M.dTw_dt[2],0.0)          self.assertAlmostEqual(M.dTw_dt[2],0.0)
781          Tw1 = float(M.T_w[2])          Tw1 = float(M.T_w[2])
# Line 575  class TestSteam(AscendSelfTester): Line 785  class TestSteam(AscendSelfTester):
785          I.setEngine('IDA')          I.setEngine('IDA')
786          I.setParameter('linsolver','DENSE')          I.setParameter('linsolver','DENSE')
787          I.setParameter('safeeval',True)          I.setParameter('safeeval',True)
788          I.setParameter('rtol',1e-8)          I.setParameter('rtol',1e-4)
789          I.setInitialSubStep(0.01)          I.setParameter('atolvect',False)
790          I.setMaxSubSteps(100)                I.setParameter('atol',1e-4)
791            I.setParameter('maxord',3)      
792            I.setInitialSubStep(0.001)
793          I.setReporter(ascpy.IntegratorReporterConsole(I))          I.setReporter(ascpy.IntegratorReporterConsole(I))
794          I.setLinearTimesteps(ascpy.Units("s"), 0, 3600, 100)          I.setLinearTimesteps(ascpy.Units("s"), 0, 3600, 10)
795          try:          I.analyse()
             I.analyse()  
         except Exception,e:  
             print "ERROR: %s" % e  
             I.writeDebug(sys.stdout)  
           
796          I.solve()          I.solve()
797          self.assertAlmostEqual(float(M.T_w[2]),Tw1)          self.assertAlmostEqual(float(M.T_w[2]),Tw1)
798          M.qdot_s.setRealValueWithUnits(1000,"W/m")          M.qdot_s.setRealValueWithUnits(1000,"W/m")
799          self.assertAlmostEqual(M.qdot_s.as("W/m"),1000)          self.assertAlmostEqual(M.qdot_s.as("W/m"),1000)
800          M.solve(ascpy.Solver('QRSlv'),ascpy.SolverReporter())          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)          self.assertNotAlmostEqual(M.dTw_dt[2],0.0)
803  #       I = ascpy.Integrator(M)          F=file('dsgsat.dot','w')
804  #       I.setEngine('LSODE')          M.write(F,'dot')
805  #       I.setReporter(ascpy.IntegratorReporterConsole(I))  
806  #       I.setReporter(ascpy.IntegratorReporterConsole(I))      def testpeturbida(self):    
807  #       I.setLinearTimesteps(ascpy.Units("s"), 0, 5, 100)          M = self.testdsgsat()
808  #       I.setMinSubStep(0.0001)          self.assertAlmostEqual(M.dTw_dt[2],0.0)
809  #       I.setMaxSubStep(100)          T = self.L.findType('dsgsat3')
810  #       I.setInitialSubStep(0.1)          M.run(T.getMethod('free_states'))
811  #       I.analyse()          # here is the peturbation...
812  #       I.solve()          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  # Testing of freesteam external steam properties functions
# Line 661  if with_freesteam and have_freesteam: Line 890  if with_freesteam and have_freesteam:
890              M = self._run("collapsingcan2",filename="collapsingcan2.a4c");              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  # Testing of IDA's analysis module
901    
902  class TestIDA(Ascend):  class TestIDA(Ascend):
# Line 740  class TestIDA(Ascend): Line 976  class TestIDA(Ascend):
976      def testfixedvars1(self):      def testfixedvars1(self):
977          self._run('fixedvars',1)          self._run('fixedvars',1)
978    
979      def testfixedvars2(self):  # fails the index check
980          self._run('fixedvars',2)  #   def testfixedvars2(self):
981    #       self._run('fixedvars',2)
982      def testfixedvars3(self):  
983          self._run('fixedvars',3)  # fails the index check
984    #   def testfixedvars3(self):
985    #       self._run('fixedvars',3)
986    
987      def testincidence(self):      def testincidence(self):
988          self._run('incidence')          self._run('incidence')
# Line 761  class TestIDA(Ascend): Line 999  class TestIDA(Ascend):
999      def testincidencefail5(self):      def testincidencefail5(self):
1000          self._runfail('incidence',5)          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:  # doesn't work yet:
1092  #   def testincidence5(self):  #   def testincidence5(self):
1093  #       self._run('incidence',5)  #       self._run('incidence',5)
# Line 772  class TestIDA(Ascend): Line 1099  class TestIDA(Ascend):
1099  class TestIDADENSE(Ascend):  class TestIDADENSE(Ascend):
1100      """IDA DAE integrator, DENSE linear solver"""      """IDA DAE integrator, DENSE linear solver"""
1101    
     def testnewton(self):  
         sys.stderr.write("STARTING TESTNEWTON\n")  
         self.L.load('johnpye/newton.a4c')  
         T = self.L.findType('newton')  
         M = T.getSimulation('sim')  
         M.solve(ascpy.Solver("QRSlv"),ascpy.SolverReporter())    
         I = ascpy.Integrator(M)  
         I.setEngine('IDA')  
         I.setParameter('linsolver','DENSE')  
         I.setParameter('safeeval',True)  
         I.setParameter('rtol',1e-8)  
         I.setMaxSubStep(0.001)  
         I.setMaxSubSteps(10000)  
           
         I.setReporter(ascpy.IntegratorReporterConsole(I))  
         I.setLinearTimesteps(ascpy.Units("s"), 0, 2*float(M.v)/float(M.g), 2)  
         I.analyse()  
         I.solve()  
         print "At end of simulation,"  
         print "x = %f" % M.x  
         print "v = %f" % M.v  
         M.run(T.getMethod('self_test'))  
   
1102      def testlotka(self):      def testlotka(self):
1103          self.L.load('johnpye/lotka.a4c')          self.L.load('johnpye/lotka.a4c')
1104          M = self.L.findType('lotka').getSimulation('sim')          M = self.L.findType('lotka').getSimulation('sim')
# Line 843  class TestIDADENSE(Ascend): Line 1147  class TestIDADENSE(Ascend):
1147          I = ascpy.Integrator(M)          I = ascpy.Integrator(M)
1148          I.setEngine('IDA')          I.setEngine('IDA')
1149          I.setParameter('linsolver','DENSE')          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)          I.setParameter('rtol',1e-7)
1193          I.setParameter('atolvect',0)          I.setParameter('atolvect',0)
1194          I.setParameter('atol',1e-7)          I.setParameter('atol',1e-7)
1195          I.setReporter(ascpy.IntegratorReporterConsole(I))          I.setReporter(ascpy.IntegratorReporterConsole(I))
1196          I.setLogTimesteps(ascpy.Units(""), 1, 321.8122, 5)          I.setLinearTimesteps(ascpy.Units("s"), 0.05, 0.2, 20)
1197          I.setMaxSubStep(1);          I.setInitialSubStep(0.00001)
1198          I.setInitialSubStep(1e-9)          I.setMaxSubSteps(10000)
         I.setMaxSubSteps(5000)  
1199          I.analyse()          I.analyse()
1200          I.solve()          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'))          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?)  ## @TODO fails during IDACalcIC (model too big?)
1228  #   def testkryx(self):  #   def testkryx(self):
1229  #       self.L.load('johnpye/idakryx.a4c')  #       self.L.load('johnpye/idakryx.a4c')
# Line 962  class TestIDASPGMR:#(Ascend): Line 1331  class TestIDASPGMR:#(Ascend):
1331          assert abs(float(M.y2) - 2.0437e-13) < 1e-15          assert abs(float(M.y2) - 2.0437e-13) < 1e-15
1332          assert abs(float(M.y3) - 1.0) < 1e-5          assert abs(float(M.y3) - 1.0) < 1e-5
1333    
1334    class TestDOPRI5(Ascend):
1335        def testlotka(self):
1336            self.L.load('johnpye/dopri5/dopri5test.a4c')
1337            M = self.L.findType('dopri5test').getSimulation('sim')
1338            M.setSolver(ascpy.Solver("QRSlv"))
1339            M.solve(ascpy.Solver("QRSlv"),ascpy.SolverReporter())  
1340            I = ascpy.Integrator(M)
1341            I.setEngine('DOPRI5')
1342            I.setReporter(ascpy.IntegratorReporterConsole(I))
1343            I.setLinearTimesteps(ascpy.Units("s"), 0, 200, 20)
1344            I.setParameter('rtol',1e-8)
1345            I.analyse()
1346            assert I.getNumVars()==1
1347            I.solve()
1348        def testaren(self):
1349            self.L.load('johnpye/dopri5/aren.a4c')
1350            M = self.L.findType('aren').getSimulation('sim')
1351            M.setSolver(ascpy.Solver("QRSlv"))
1352            M.solve(ascpy.Solver("QRSlv"),ascpy.SolverReporter())
1353            I = ascpy.Integrator(M)
1354            I.setEngine('DOPRI5')
1355            I.setReporter(ascpy.IntegratorReporterConsole(I))
1356            #xend = 17.0652165601579625588917206249
1357            I.setLinearTimesteps(ascpy.Units("s"), 0, 17.0652165601579625588917206249, 10)
1358            I.setParameter('rtol',1e-7)
1359            I.setParameter('atol',1e-7)
1360            I.setParameter('tolvect',False)
1361            I.setMinSubStep(0);
1362            I.setMaxSubStep(0);
1363            I.setInitialSubStep(0);
1364            I.analyse()
1365            I.solve()
1366            assert abs(float(M.y[0]) - 0.994) < 1e-5
1367            assert abs(float(M.y[1]) - 0.0) < 1e-5
1368    
1369    class TestIPOPT(Ascend):
1370        def test1(self):
1371            self.L.load('test/ipopt/test1.a4c')
1372            M = self.L.findType('test1').getSimulation('sim')
1373            M.solve(ascpy.Solver("IPOPT"),ascpy.SolverReporter())
1374    
1375    # test some stuff for beam calculations
1376    class TestSection(Ascend):
1377        def test_compound3(self):
1378            self.L.load('johnpye/section.a4c')
1379            T = self.L.findType('compound_section_test3')
1380            M = T.getSimulation('sim')
1381            M.solve(ascpy.Solver("QRSlv"),ascpy.SolverReporter())
1382            M.run(T.getMethod('self_test'))
1383        def test_compound2(self):
1384            self.L.load('johnpye/section.a4c')
1385            T = self.L.findType('compound_section_test2')
1386            M = T.getSimulation('sim')
1387            M.solve(ascpy.Solver("QRSlv"),ascpy.SolverReporter())
1388            M.run(T.getMethod('self_test'))
1389    
1390  # move code above down here if you want to temporarily avoid testing it  # move code above down here if you want to temporarily avoid testing it
1391  class NotToBeTested:  class NotToBeTested:
1392      def nothing(self):      def nothing(self):
1393          pass          pass
1394    
1395        def testnewton(self):
1396            sys.stderr.write("STARTING TESTNEWTON\n")
1397            self.L.load('johnpye/newton.a4c')
1398            T = self.L.findType('newton')
1399            M = T.getSimulation('sim')
1400            M.solve(ascpy.Solver("QRSlv"),ascpy.SolverReporter())  
1401            I = ascpy.Integrator(M)
1402            I.setEngine('IDA')
1403            I.setParameter('linsolver','DENSE')
1404            I.setParameter('safeeval',True)
1405            I.setParameter('rtol',1e-8)
1406            I.setMaxSubStep(0.001)
1407            I.setMaxSubSteps(10000)
1408            
1409            I.setReporter(ascpy.IntegratorReporterConsole(I))
1410            I.setLinearTimesteps(ascpy.Units("s"), 0, 2*float(M.v)/float(M.g), 2)
1411            I.analyse()
1412            I.solve()
1413            print "At end of simulation,"
1414            print "x = %f" % M.x
1415            print "v = %f" % M.v
1416            M.run(T.getMethod('self_test'))
1417    
1418    def patchpath(VAR,SEP,addvals):
1419        restart = 0
1420        envpath = [os.path.abspath(i) for i in os.environ[VAR].split(SEP)]
1421        for l in addvals:
1422            if l in envpath[len(addvals):]:
1423                envpath.remove(l)
1424                restart = 1
1425        for l in addvals:
1426            if l not in envpath:
1427                envpath.insert(0,l)
1428                restart = 1
1429        os.environ[VAR] = SEP.join(envpath)
1430        return restart  
1431        
1432  if __name__=='__main__':  if __name__=='__main__':
1433      # a whole bag of tricks to make sure we get the necessary dirs in our ascend, python and ld path vars      # a whole bag of tricks to make sure we get the necessary dirs in our ascend, python and ld path vars
1434      restart = 0      restart = 0
# Line 978  if __name__=='__main__': Line 1440  if __name__=='__main__':
1440          LD_LIBRARY_PATH="LD_LIBRARY_PATH"          LD_LIBRARY_PATH="LD_LIBRARY_PATH"
1441          SEP = ":"          SEP = ":"
1442    
1443        solverdir = os.path.abspath(os.path.join(sys.path[0],"solvers"))
1444        solverdirs = [os.path.join(solverdir,s) for s in "qrslv","cmslv","lrslv","conopt","ida","lsode","ipopt"]
1445    
1446        if not os.environ.get('ASCENDSOLVERS'):
1447            os.environ['ASCENDSOLVERS'] = SEP.join(solverdirs)
1448            restart = 1
1449        else:
1450            if patchpath('ASCENDSOLVERS',SEP,solverdirs):
1451                restart = 1
1452        
1453      freesteamdir = os.path.expanduser("~/freesteam/ascend")      freesteamdir = os.path.expanduser("~/freesteam/ascend")
1454      modeldirs = [os.path.abspath(os.path.join(sys.path[0],"models")),os.path.abspath(freesteamdir)]      modeldirs = [os.path.abspath(os.path.join(sys.path[0],"models")),os.path.abspath(freesteamdir)]
1455        
1456      if not os.environ.get('ASCENDLIBRARY'):      if not os.environ.get('ASCENDLIBRARY'):
1457          os.environ['ASCENDLIBRARY'] = SEP.join(modeldirs)          os.environ['ASCENDLIBRARY'] = SEP.join(modeldirs)
1458          restart = 1          restart = 1
1459      else:      else:
1460          envmodelsdir = [os.path.abspath(i) for i in os.environ['ASCENDLIBRARY'].split(SEP)]          if patchpath('ASCENDLIBRARY',SEP,modeldirs):
1461          for l in modeldirs:              restart = 1
             if l in envmodelsdir[len(modeldirs):]:  
                 envmodelsdir.remove(l)  
                 restart = 1  
         for l in modeldirs:  
             if l not in envmodelsdir:  
                 envmodelsdir.insert(0,l)  
                 restart = 1  
         os.environ['ASCENDLIBRARY'] = SEP.join(envmodelsdir)      
1462    
1463      libdirs = ["pygtk","."]      libdirs = ["pygtk","."]
1464      libdirs = [os.path.normpath(os.path.join(sys.path[0],l)) for l in libdirs]      libdirs = [os.path.normpath(os.path.join(sys.path[0],l)) for l in libdirs]
# Line 1024  if __name__=='__main__': Line 1489  if __name__=='__main__':
1489    
1490      if restart:      if restart:
1491          script = os.path.join(sys.path[0],"test.py")                              script = os.path.join(sys.path[0],"test.py")                    
1492          print "Restarting with..."          sys.stderr.write("Restarting with...\n")
1493          print "  export LD_LIBRARY_PATH=%s" % os.environ.get(LD_LIBRARY_PATH)          sys.stderr.write("  export LD_LIBRARY_PATH=%s\n" % os.environ.get(LD_LIBRARY_PATH))
1494          print "  export PYTHONPATH=%s" % os.environ.get('PYTHONPATH')          sys.stderr.write("  export PYTHONPATH=%s\n" % os.environ.get('PYTHONPATH'))
1495          print "  export ASCENDLIBRARY=%s" % os.environ.get('ASCENDLIBRARY')          sys.stderr.write("  export ASCENDLIBRARY=%s\n" % os.environ.get('ASCENDLIBRARY'))
1496            sys.stderr.write("  export ASCENDSOLVERS=%s\n" % os.environ.get('ASCENDSOLVERS'))
1497            sys.stderr.flush()
1498          os.execvp("python",[script] + sys.argv)          os.execvp("python",[script] + sys.argv)
1499            exit(1)
1500        else:
1501            sys.stderr.write("Got...\n")
1502            sys.stderr.write("  LD_LIBRARY_PATH=%s\n" % os.environ.get(LD_LIBRARY_PATH))
1503            sys.stderr.write("  PYTHONPATH=%s\n" % os.environ.get('PYTHONPATH'))
1504            sys.stderr.write("  ASCENDLIBRARY=%s\n" % os.environ.get('ASCENDLIBRARY'))
1505            sys.stderr.write("  ASCENDSOLVERS=%s\n" % os.environ.get('ASCENDSOLVERS'))
1506            sys.stderr.flush()
1507    
1508      import ascpy      import ascpy
1509    

Legend:
Removed from v.1247  
changed lines
  Added in v.1562

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