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

Diff of /trunk/test.py

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

revision 1221 by johnpye, Wed Jan 24 13:33:06 2007 UTC revision 1278 by johnpye, Fri Feb 9 02:35: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 43  class Ascend(unittest.TestCase): Line 43  class Ascend(unittest.TestCase):
43    
44  class AscendSelfTester(Ascend):  class AscendSelfTester(Ascend):
45    
46      def _run(self,modelname,solvername="QRSlv",filename=None):      def _run(self,modelname,solvername="QRSlv",filename=None,parameters={}):
47          if filename==None:          if filename==None:
48              filename = 'johnpye/%s.a4c' % modelname              filename = 'johnpye/%s.a4c' % modelname
49          self.L.load(filename)          self.L.load(filename)
50          T = self.L.findType(modelname)          T = self.L.findType(modelname)
51          M = T.getSimulation('sim')          M = T.getSimulation('sim')
52          M.build()          M.setSolver(ascpy.Solver(solvername))
53            for k,v in parameters.iteritems():
54                M.setParameter(k,v)
55          M.solve(ascpy.Solver(solvername),ascpy.SolverReporter())              M.solve(ascpy.Solver(solvername),ascpy.SolverReporter())    
56          M.run(T.getMethod('self_test'))          M.run(T.getMethod('self_test'))
57          return M          return M
# Line 134  class TestSolver(AscendSelfTester): Line 136  class TestSolver(AscendSelfTester):
136          self.assertAlmostEqual( float(M.t_solar), M.t_solar.as("s"))          self.assertAlmostEqual( float(M.t_solar), M.t_solar.as("s"))
137          self.assertAlmostEqual( float(M.t_solar)/3600, M.t_solar.as("h"))          self.assertAlmostEqual( float(M.t_solar)/3600, M.t_solar.as("h"))
138    
139    
140    class TestCMSlv(AscendSelfTester):
141        def testheatex(self):
142            self._run('heatex',"CMSlv","heatex.a4c")
143        def testphaseeq(self):
144            self._run('phaseq',"CMSlv","phaseq.a4c")
145        def testpipeline(self):
146            self._run('pipeline',"CMSlv","pipeline.a4c"
147                ,{'infinity':3.2e9}
148            )
149        def testrachford(self):
150            self._run('rachford',"CMSlv","rachford.a4c")
151        def testlinmassbal(self):
152            self._run('linmassbal',"CMSlv","linmassbal.a4c")
153    
154    
155  class TestMatrix(AscendSelfTester):  class TestMatrix(AscendSelfTester):
156      def testlog10(self):      def testlog10(self):
157          M = self._run('testlog10')          M = self._run('testlog10')
# Line 551  class TestSteam(AscendSelfTester): Line 569  class TestSteam(AscendSelfTester):
569          M.solve(ascpy.Solver('QRSlv'),ascpy.SolverReporter())          M.solve(ascpy.Solver('QRSlv'),ascpy.SolverReporter())
570          return M          return M
571    
572      def testintegLSODE(self):      def testdsgsatrepeat(self):
573          M = self.testdsgsat()          self.L.load('steam/dsgsat3.a4c')
574            T = self.L.findType('dsgsat3')
575            M = T.getSimulation('sim',False)
576            M.run(T.getMethod('on_load'))
577            M.solve(ascpy.Solver('QRSlv'),ascpy.SolverReporter())
578            M.run(T.getMethod('on_load'))
579            M.solve(ascpy.Solver('QRSlv'),ascpy.SolverReporter())
580            M.run(T.getMethod('on_load'))
581            M.solve(ascpy.Solver('QRSlv'),ascpy.SolverReporter())
582    
583        def testvary(self):
584            self.L.load('steam/dsgsat3.a4c')
585            T = self.L.findType('dsgsat3')
586            M = T.getSimulation('sim',False)
587            M.run(T.getMethod('on_load'))
588            M.solve(ascpy.Solver('QRSlv'),ascpy.SolverReporter())
589            print "----- setting qdot_s -----"
590          M.qdot_s.setRealValueWithUnits(1000,"W/m")          M.qdot_s.setRealValueWithUnits(1000,"W/m")
591          M.solve(ascpy.Solver('QRSlv'),ascpy.SolverReporter())          M.solve(ascpy.Solver('QRSlv'),ascpy.SolverReporter())
592            print "----- setting qdot_s -----"
593            M.qdot_s.setRealValueWithUnits(2000,"W/m")
594            M.solve(ascpy.Solver('QRSlv'),ascpy.SolverReporter())
595    
596        def teststeadylsode(self):
597            "test that steady conditions are stable with LSODE"
598            M = self.testdsgsat()
599            #M.qdot_s.setRealValueWithUnits(1000,"W/m")
600            M.solve(ascpy.Solver('QRSlv'),ascpy.SolverReporter())
601          #M.setParameter('          #M.setParameter('
602          I = ascpy.Integrator(M)          I = ascpy.Integrator(M)
603          I.setEngine('LSODE')          I.setEngine('LSODE')
         I.setParameter('meth','AM')  
         I.setParameter('maxord',12)  
604          I.setReporter(ascpy.IntegratorReporterConsole(I))          I.setReporter(ascpy.IntegratorReporterConsole(I))
605          I.setLinearTimesteps(ascpy.Units("s"), 0, 5, 1)          I.setLinearTimesteps(ascpy.Units("s"), 0, 5, 1)
606          I.analyse()          I.analyse()
607          I.solve()          I.solve()
608    
609      def testintegIDA(self):  # DOESN'T CONVERGE
610    #   def testpeturblsode(self):
611    #       "test that steady conditions are stable with LSODE"
612    #       M = self.testdsgsat()
613    #       # here is the peturbation...
614    #       M.qdot_s.setRealValueWithUnits(1000,"W/m")
615    #       M.solve(ascpy.Solver('QRSlv'),ascpy.SolverReporter())
616    #       I = ascpy.Integrator(M)
617    #       I.setEngine('LSODE')
618    #       I.setReporter(ascpy.IntegratorReporterConsole(I))
619    #       I.setLinearTimesteps(ascpy.Units("s"), 0, 5, 1)
620    #       I.analyse()
621    #       I.solve()
622    
623        def teststeadyida(self):    
624          M = self.testdsgsat()          M = self.testdsgsat()
625          self.assertAlmostEqual(M.dTw_dt[2],0.0)          self.assertAlmostEqual(M.dTw_dt[2],0.0)
626          Tw1 = float(M.T_w[2])          Tw1 = float(M.T_w[2])
# Line 575  class TestSteam(AscendSelfTester): Line 630  class TestSteam(AscendSelfTester):
630          I.setEngine('IDA')          I.setEngine('IDA')
631          I.setParameter('linsolver','DENSE')          I.setParameter('linsolver','DENSE')
632          I.setParameter('safeeval',True)          I.setParameter('safeeval',True)
633          I.setParameter('rtol',1e-8)          I.setParameter('rtol',1e-5)
634          I.setInitialSubStep(0.01)          I.setInitialSubStep(0.01)
635          I.setMaxSubSteps(100)                I.setMaxSubSteps(100)      
636          I.setReporter(ascpy.IntegratorReporterConsole(I))          I.setReporter(ascpy.IntegratorReporterConsole(I))
637          I.setLinearTimesteps(ascpy.Units("s"), 0, 3600, 100)          I.setLinearTimesteps(ascpy.Units("s"), 0, 3600, 5)
638          try:          I.analyse()
             I.analyse()  
         except Exception,e:  
             print "ERROR: %s" % e  
             I.writeDebug(sys.stdout)  
           
639          I.solve()          I.solve()
640          self.assertAlmostEqual(float(M.T_w[2]),Tw1)          self.assertAlmostEqual(float(M.T_w[2]),Tw1)
641          M.qdot_s.setRealValueWithUnits(1000,"W/m")          M.qdot_s.setRealValueWithUnits(1000,"W/m")
642          self.assertAlmostEqual(M.qdot_s.as("W/m"),1000)          self.assertAlmostEqual(M.qdot_s.as("W/m"),1000)
643          M.solve(ascpy.Solver('QRSlv'),ascpy.SolverReporter())          M.solve(ascpy.Solver('QRSlv'),ascpy.SolverReporter())
644            print "dTw/dt = %f" % M.dTw_dt[2]
645          self.assertNotAlmostEqual(M.dTw_dt[2],0.0)          self.assertNotAlmostEqual(M.dTw_dt[2],0.0)
646  #       I = ascpy.Integrator(M)  
647  #       I.setEngine('LSODE')      def testpeturbida(self):    
648  #       I.setReporter(ascpy.IntegratorReporterConsole(I))          M = self.testdsgsat()
649  #       I.setReporter(ascpy.IntegratorReporterConsole(I))          self.assertAlmostEqual(M.dTw_dt[2],0.0)
650  #       I.setLinearTimesteps(ascpy.Units("s"), 0, 5, 100)          T = self.L.findType('dsgsat3')
651  #       I.setMinSubStep(0.0001)          M.run(T.getMethod('free_states'))
652  #       I.setMaxSubStep(100)          # here is the peturbation...
653  #       I.setInitialSubStep(0.1)          M.qdot_s.setRealValueWithUnits(6000,"W/m")
654  #       I.analyse()          # IDA has its own initial conditions solver, so no need to call QRSlv here
655  #       I.solve()          I = ascpy.Integrator(M)
656            I.setEngine('IDA')
657            I.setParameter('linsolver','DENSE')
658            I.setParameter('safeeval',True)
659            I.setParameter('rtol',1e-5)
660            I.setParameter('atolvect',False)
661            I.setParameter('atol',1e-5)
662            I.setInitialSubStep(0.1)
663            I.setReporter(ascpy.IntegratorReporterConsole(I))
664            I.setLogTimesteps(ascpy.Units("s"), 0.01, 3600, 20)
665            I.analyse()
666            I.solve()
667                    
668  #-------------------------------------------------------------------------------  #-------------------------------------------------------------------------------
669  # Testing of freesteam external steam properties functions  # Testing of freesteam external steam properties functions
# Line 672  class TestIDA(Ascend): Line 734  class TestIDA(Ascend):
734          I = ascpy.Integrator(M)          I = ascpy.Integrator(M)
735          I.setEngine('IDA')          I.setEngine('IDA')
736          I.analyse()          I.analyse()
737            return M;
738    
739      def _runfail(self,filen,n,msg="failed"):      def _runfail(self,filen,n,msg="failed"):
740          try:          try:
# Line 739  class TestIDA(Ascend): Line 802  class TestIDA(Ascend):
802      def testfixedvars1(self):      def testfixedvars1(self):
803          self._run('fixedvars',1)          self._run('fixedvars',1)
804    
805  # CAUSES A CRASH !?!?!!      def testfixedvars2(self):
806  #   def testfixedvars2(self):          self._run('fixedvars',2)
807  #       self._run('fixedvars',2)  
808        def testfixedvars3(self):
809  # CAUSES A CRASH !!!?!?!?          self._run('fixedvars',3)
 #   def testfixedvars3(self):  
 #       self._run('fixedvars',3)  
810    
811      def testincidence(self):      def testincidence(self):
812          self._run('incidence')          self._run('incidence')
# Line 755  class TestIDA(Ascend): Line 816  class TestIDA(Ascend):
816      def testincidence2(self):      def testincidence2(self):
817          self._run('incidence',2)          self._run('incidence',2)
818      def testincidence3(self):      def testincidence3(self):
819          self._run('incidence',3)          M = self._run('incidence',3)
820    
821      def testincidence4(self):      def testincidence4(self):
822          self._run('incidence',4)          self._run('incidence',4)
823      def testincidencefail5(self):      def testincidencefail5(self):
# Line 772  class TestIDA(Ascend): Line 834  class TestIDA(Ascend):
834  class TestIDADENSE(Ascend):  class TestIDADENSE(Ascend):
835      """IDA DAE integrator, DENSE linear solver"""      """IDA DAE integrator, DENSE linear solver"""
836    
     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'))  
   
837      def testlotka(self):      def testlotka(self):
838          self.L.load('johnpye/lotka.a4c')          self.L.load('johnpye/lotka.a4c')
839          M = self.L.findType('lotka').getSimulation('sim')          M = self.L.findType('lotka').getSimulation('sim')
# Line 820  class TestIDADENSE(Ascend): Line 859  class TestIDADENSE(Ascend):
859          M.setSolver(ascpy.Solver("QRSlv"))          M.setSolver(ascpy.Solver("QRSlv"))
860          I = ascpy.Integrator(M)          I = ascpy.Integrator(M)
861          I.setEngine('IDA')          I.setEngine('IDA')
862          I.setParameter('calcic','YA_YPD')          I.setParameter('calcic','YA_YDP')
863          I.setParameter('linsolver','DENSE')          I.setParameter('linsolver','DENSE')
864            I.setParameter('safeeval',True)
865          I.setReporter(ascpy.IntegratorReporterConsole(I))          I.setReporter(ascpy.IntegratorReporterConsole(I))
866          I.setLogTimesteps(ascpy.Units("s"), 0.4, 4e10, 11)          I.setLogTimesteps(ascpy.Units("s"), 0.4, 4e10, 11)
867          I.setMaxSubStep(0);          I.setMaxSubStep(0);
868          I.setInitialSubStep(0)          I.setInitialSubStep(0)
869          I.setMaxSubSteps(0);          I.setMaxSubSteps(0)
870          I.setParameter('autodiff',True)          I.setParameter('autodiff',True)
871          I.analyse()          I.analyse()
872          I.solve()          I.solve()
# Line 834  class TestIDADENSE(Ascend): Line 874  class TestIDADENSE(Ascend):
874          assert abs(float(M.y2) - 2.0437e-13) < 2e-14          assert abs(float(M.y2) - 2.0437e-13) < 2e-14
875          assert abs(float(M.y3) - 1.0) < 1e-5          assert abs(float(M.y3) - 1.0) < 1e-5
876    
877        def testhires(self):
878            self.L.load('test/hires.a4c')
879            T = self.L.findType('hires')
880            M = T.getSimulation('sim')
881            M.setSolver(ascpy.Solver('QRSlv'))
882            I = ascpy.Integrator(M)
883            I.setEngine('IDA')
884            I.setParameter('linsolver','DENSE')
885            I.setParameter('rtol',1.1e-15)
886            I.setParameter('atolvect',0)
887            I.setParameter('atol',1.1e-15)
888            I.setReporter(ascpy.IntegratorReporterConsole(I))
889            I.setLogTimesteps(ascpy.Units(""), 1, 321.8122, 5)
890            I.setInitialSubStep(1e-5)
891            I.setMaxSubSteps(10000)
892            I.analyse()
893            I.solve()
894            for i in range(8):
895                print "y[%d] = %.20g" % (i+1, M.y[i+1])
896            M.run(T.getMethod('self_test'))
897    
898        def testchemakzo(self):
899            self.L.load('test/chemakzo.a4c')
900            T = self.L.findType('chemakzo')
901            M = T.getSimulation('sim')
902            M.setSolver(ascpy.Solver('QRSlv'))
903            I = ascpy.Integrator(M)
904            I.setEngine('IDA')
905            I.setParameter('linsolver','DENSE')
906            I.setParameter('rtol',1e-15)
907            I.setParameter('atolvect',0)
908            I.setParameter('atol',1e-15)
909            I.setReporter(ascpy.IntegratorReporterConsole(I))
910            I.setLinearTimesteps(ascpy.Units("s"), 1, 180, 5)
911            I.setInitialSubStep(1e-13)
912            I.setMaxSubSteps(10000)
913            I.analyse()
914            I.solve()
915            for i in range(6):
916                print "y[%d] = %.20g" % (i+1, M.y[i+1])
917            M.run(T.getMethod('self_test'))
918    
919        def testtransamp(self):
920            self.L.load('test/transamp.a4c')
921            T = self.L.findType('transamp')
922            M = T.getSimulation('sim')
923            M.setSolver(ascpy.Solver('QRSlv'))
924            I = ascpy.Integrator(M)
925            I.setEngine('IDA')
926            I.setParameter('linsolver','DENSE')
927            I.setParameter('rtol',1e-7)
928            I.setParameter('atolvect',0)
929            I.setParameter('atol',1e-7)
930            I.setReporter(ascpy.IntegratorReporterConsole(I))
931            I.setLinearTimesteps(ascpy.Units("s"), 0.05, 0.2, 20)
932            I.setInitialSubStep(0.00001)
933            I.setMaxSubSteps(10000)
934            I.analyse()
935            I.solve()
936            for i in range(6):
937                print "y[%d] = %.20g" % (i+1, M.y[i+1])
938            M.run(T.getMethod('self_test'))
939    
940    # MODEL FAILS ANALYSIS: we need to add support for non-incident differential vars
941    #   def testpollution(self):
942    #       self.L.load('test/pollution.a4c')
943    #       T = self.L.findType('pollution')
944    #       M = T.getSimulation('sim')
945    #       M.setSolver(ascpy.Solver('QRSlv'))
946    #       I = ascpy.Integrator(M)
947    #       I.setEngine('IDA')
948    #       I.setParameter('linsolver','DENSE')
949    #       I.setParameter('rtol',1.1e-15)
950    #       I.setParameter('atolvect',0)
951    #       I.setParameter('atol',1.1e-15)
952    #       I.setReporter(ascpy.IntegratorReporterConsole(I))
953    #       I.setLogTimesteps(ascpy.Units("s"), 1, 60, 5)
954    #       I.setInitialSubStep(1e-5)
955    #       I.setMaxSubSteps(10000)
956    #       I.analyse()
957    #       I.solve()
958    #       for i in range(20):
959    #           print "y[%d] = %.20g" % (i+1, M.y[i+1])
960    #       M.run(T.getMethod('self_test'))
961    
962  ## @TODO fails during IDACalcIC (model too big?)  ## @TODO fails during IDACalcIC (model too big?)
963  #   def testkryx(self):  #   def testkryx(self):
964  #       self.L.load('johnpye/idakryx.a4c')  #       self.L.load('johnpye/idakryx.a4c')
# Line 946  class NotToBeTested: Line 1071  class NotToBeTested:
1071      def nothing(self):      def nothing(self):
1072          pass          pass
1073    
1074        def testnewton(self):
1075            sys.stderr.write("STARTING TESTNEWTON\n")
1076            self.L.load('johnpye/newton.a4c')
1077            T = self.L.findType('newton')
1078            M = T.getSimulation('sim')
1079            M.solve(ascpy.Solver("QRSlv"),ascpy.SolverReporter())  
1080            I = ascpy.Integrator(M)
1081            I.setEngine('IDA')
1082            I.setParameter('linsolver','DENSE')
1083            I.setParameter('safeeval',True)
1084            I.setParameter('rtol',1e-8)
1085            I.setMaxSubStep(0.001)
1086            I.setMaxSubSteps(10000)
1087            
1088            I.setReporter(ascpy.IntegratorReporterConsole(I))
1089            I.setLinearTimesteps(ascpy.Units("s"), 0, 2*float(M.v)/float(M.g), 2)
1090            I.analyse()
1091            I.solve()
1092            print "At end of simulation,"
1093            print "x = %f" % M.x
1094            print "v = %f" % M.v
1095            M.run(T.getMethod('self_test'))
1096    
1097  if __name__=='__main__':  if __name__=='__main__':
1098      # 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
1099      restart = 0      restart = 0

Legend:
Removed from v.1221  
changed lines
  Added in v.1278

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