/[escript]/trunk/finley/test/python/run_models.py
ViewVC logotype

Diff of /trunk/finley/test/python/run_models.py

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

revision 2099 by ksteube, Thu Sep 25 06:43:44 2008 UTC revision 2100 by gross, Wed Nov 26 08:13:00 2008 UTC
# Line 24  import tempfile Line 24  import tempfile
24                
25  from esys.escript import *  from esys.escript import *
26  from esys.finley import Rectangle  from esys.finley import Rectangle
27    from esys.escript.models import DarcyFlow
28  import sys  import sys
29  import os  import os
30  try:  try:
# Line 32  except KeyError: Line 33  except KeyError:
33       FINLEY_WORKDIR='.'       FINLEY_WORKDIR='.'
34    
35    
36  NE=6  VERBOSE=False # or True
 TOL=1.e-5  
 VERBOSE=False or True  
37    
38  from esys.escript import *  from esys.escript import *
39  from esys.escript.models import StokesProblemCartesian  from esys.escript.models import StokesProblemCartesian
40  from esys.finley import Rectangle, Brick  from esys.finley import Rectangle, Brick
41  class Test_Simple2DModels(unittest.TestCase):  
42    #====================================================================================================================
43    class Test_StokesProblemCartesian2D(unittest.TestCase):
44     def setUp(self):     def setUp(self):
45           NE=6
46           self.TOL=1.e-5
47         self.domain=Rectangle(NE,NE,order=2,useFullElementOrder=True)         self.domain=Rectangle(NE,NE,order=2,useFullElementOrder=True)
48     def tearDown(self):     def tearDown(self):
49         del self.domain         del self.domain
50     def test_StokesProblemCartesian_PCG_P_0(self):     def test_PCG_P_0(self):
51         ETA=1.         ETA=1.
52         P1=0.         P1=0.
53    
# Line 60  class Test_Simple2DModels(unittest.TestC Line 63  class Test_Simple2DModels(unittest.TestC
63         sp.initialize(f=F,fixed_u_mask=mask,eta=ETA)         sp.initialize(f=F,fixed_u_mask=mask,eta=ETA)
64         u0=(1-x[0])*x[0]*[0.,1.]         u0=(1-x[0])*x[0]*[0.,1.]
65         p0=Scalar(P1,ReducedSolution(self.domain))         p0=Scalar(P1,ReducedSolution(self.domain))
66         u,p=sp.solve(u0,p0,show_details=VERBOSE, verbose=VERBOSE,max_iter=100,solver="PCG")         sp.setTolerance(self.TOL)
67           u,p=sp.solve(u0,p0,show_details=VERBOSE, verbose=VERBOSE,max_iter=100,useUzawa=True)
68                
69         error_v0=Lsup(u[0]-u0[0])         error_v0=Lsup(u[0]-u0[0])
70         error_v1=Lsup(u[1]-u0[1])/0.25         error_v1=Lsup(u[1]-u0[1])/0.25
71         zz=P1*x[0]*x[1]-p         error_p=Lsup(p+P1*x[0]*x[1])
72         error_p=Lsup(zz-integrate(zz))         self.failUnless(error_p<10*self.TOL, "pressure error too large.")
73         # print error_p, error_v0,error_v1         self.failUnless(error_v0<10*self.TOL, "0-velocity error too large.")
74         self.failUnless(error_p<TOL,"pressure error too large.")         self.failUnless(error_v1<10*self.TOL, "1-velocity error too large.")
        self.failUnless(error_v0<TOL,"0-velocity error too large.")  
        self.failUnless(error_v1<TOL,"1-velocity error too large.")  
75    
76     def test_StokesProblemCartesian_PCG_P_small(self):     def test_PCG_P_small(self):
77         ETA=1.         ETA=1.
78         P1=1.         P1=1.
79    
# Line 87  class Test_Simple2DModels(unittest.TestC Line 89  class Test_Simple2DModels(unittest.TestC
89         sp.initialize(f=F,fixed_u_mask=mask,eta=ETA)         sp.initialize(f=F,fixed_u_mask=mask,eta=ETA)
90         u0=(1-x[0])*x[0]*[0.,1.]         u0=(1-x[0])*x[0]*[0.,1.]
91         p0=Scalar(P1,ReducedSolution(self.domain))         p0=Scalar(P1,ReducedSolution(self.domain))
92         u,p=sp.solve(u0,p0,show_details=VERBOSE, verbose=VERBOSE,max_iter=100,solver="PCG")         sp.setTolerance(self.TOL)
93                 u,p=sp.solve(u0,p0,show_details=VERBOSE, verbose=VERBOSE or True,max_iter=100,useUzawa=True)
94         error_v0=Lsup(u[0]-u0[0])         error_v0=Lsup(u[0]-u0[0])
95         error_v1=Lsup(u[1]-u0[1])/0.25         error_v1=Lsup(u[1]-u0[1])/0.25
96         zz=P1*x[0]*x[1]-p         error_p=Lsup(P1*x[0]*x[1]+p)
97         error_p=Lsup(zz-integrate(zz))         self.failUnless(error_v0<10*self.TOL, "0-velocity error too large.")
98         # print error_p, error_v0,error_v1         self.failUnless(error_v1<10*self.TOL, "1-velocity error too large.")
99         self.failUnless(error_p<TOL,"pressure error too large.")         self.failUnless(error_p<11*self.TOL, "pressure error too large.")
        self.failUnless(error_v0<TOL,"0-velocity error too large.")  
        self.failUnless(error_v1<TOL,"1-velocity error too large.")  
100    
101     def test_StokesProblemCartesian_PCG_P_large(self):     def test_PCG_P_large(self):
102         ETA=1.         ETA=1.
103         P1=1000.         P1=1000.
104    
# Line 114  class Test_Simple2DModels(unittest.TestC Line 114  class Test_Simple2DModels(unittest.TestC
114         sp.initialize(f=F,fixed_u_mask=mask,eta=ETA)         sp.initialize(f=F,fixed_u_mask=mask,eta=ETA)
115         u0=(1-x[0])*x[0]*[0.,1.]         u0=(1-x[0])*x[0]*[0.,1.]
116         p0=Scalar(P1,ReducedSolution(self.domain))         p0=Scalar(P1,ReducedSolution(self.domain))
117         u,p=sp.solve(u0,p0,show_details=VERBOSE, verbose=VERBOSE,max_iter=100,solver="PCG")         sp.setTolerance(self.TOL)
118           u,p=sp.solve(u0,p0,show_details=VERBOSE, verbose=VERBOSE,max_iter=100,useUzawa=True)
119                
120         error_v0=Lsup(u[0]-u0[0])         error_v0=Lsup(u[0]-u0[0])
121         error_v1=Lsup(u[1]-u0[1])/0.25         error_v1=Lsup(u[1]-u0[1])/0.25
122         zz=P1*x[0]*x[1]-p         error_p=Lsup(P1*x[0]*x[1]+p)/P1
123         error_p=Lsup(zz-integrate(zz))/P1         self.failUnless(error_p<10*self.TOL, "pressure error too large.")
124         # print error_p, error_v0,error_v1         self.failUnless(error_v0<10*self.TOL, "0-velocity error too large.")
125         self.failUnless(error_p<TOL,"pressure error too large.")         self.failUnless(error_v1<10*self.TOL, "1-velocity error too large.")
        self.failUnless(error_v0<TOL,"0-velocity error too large.")  
        self.failUnless(error_v1<TOL,"1-velocity error too large.")  
126    
127     def test_StokesProblemCartesian_GMRES_P_0(self):     def test_GMRES_P_0(self):
128         ETA=1.         ETA=1.
129         P1=0.         P1=0.
130    
# Line 141  class Test_Simple2DModels(unittest.TestC Line 140  class Test_Simple2DModels(unittest.TestC
140         sp.initialize(f=F,fixed_u_mask=mask,eta=ETA)         sp.initialize(f=F,fixed_u_mask=mask,eta=ETA)
141         u0=(1-x[0])*x[0]*[0.,1.]         u0=(1-x[0])*x[0]*[0.,1.]
142         p0=Scalar(P1,ReducedSolution(self.domain))         p0=Scalar(P1,ReducedSolution(self.domain))
143         u,p=sp.solve(u0,p0,show_details=VERBOSE, verbose=VERBOSE,max_iter=100,solver="GMRES",iter_restart=20)         sp.setTolerance(self.TOL)
144           u,p=sp.solve(u0,p0,show_details=VERBOSE, verbose=VERBOSE,max_iter=50,useUzawa=False,iter_restart=18)
145           # u,p=sp.solve(u0,p0,show_details=VERBOSE, verbose=VERBOSE,max_iter=100,useUzawa=False,iter_restart=20)
146                
147         error_v0=Lsup(u[0]-u0[0])         error_v0=Lsup(u[0]-u0[0])
148         error_v1=Lsup(u[1]-u0[1])/0.25         error_v1=Lsup(u[1]-u0[1])/0.25
149         zz=P1*x[0]*x[1]-p         zz=P1*x[0]*x[1]+p
150         error_p=Lsup(zz-integrate(zz))         error_p=Lsup(zz-integrate(zz))
151           self.failUnless(error_p<10*self.TOL, "pressure error too large.")
152           self.failUnless(error_v0<10*self.TOL, "0-velocity error too large.")
153           self.failUnless(error_v1<10*self.TOL, "1-velocity error too large.")
154    
155        # print error_p, error_v0,error_v1     def test_GMRES_P_small(self):
        self.failUnless(error_p<TOL,"pressure error too large.")  
        self.failUnless(error_v0<TOL,"0-velocity error too large.")  
        self.failUnless(error_v1<TOL,"1-velocity error too large.")  
   
    def test_StokesProblemCartesian_GMRES_P_small(self):  
156         ETA=1.         ETA=1.
157         P1=1.         P1=1.
158    
# Line 169  class Test_Simple2DModels(unittest.TestC Line 168  class Test_Simple2DModels(unittest.TestC
168         sp.initialize(f=F,fixed_u_mask=mask,eta=ETA)         sp.initialize(f=F,fixed_u_mask=mask,eta=ETA)
169         u0=(1-x[0])*x[0]*[0.,1.]         u0=(1-x[0])*x[0]*[0.,1.]
170         p0=Scalar(P1,ReducedSolution(self.domain))         p0=Scalar(P1,ReducedSolution(self.domain))
171         u,p=sp.solve(u0,p0,show_details=VERBOSE, verbose=VERBOSE,max_iter=100,solver="GMRES")         sp.setTolerance(self.TOL*0+1.e-6)
172           # sp.setSubToleranceReductionFactor(0.1)
173           u,p=sp.solve(u0,p0,show_details=VERBOSE, verbose=VERBOSE,max_iter=20,useUzawa=False)
174                
175         error_v0=Lsup(u[0]-u0[0])         error_v0=Lsup(u[0]-u0[0])
176         error_v1=Lsup(u[1]-u0[1])/0.25         error_v1=Lsup(u[1]-u0[1])/0.25
177         zz=P1*x[0]*x[1]-p         zz=P1*x[0]*x[1]+p
178         error_p=Lsup(zz-integrate(zz))         error_p=Lsup(zz-integrate(zz))
179         # print error_p, error_v0,error_v1         # print error_p, error_v0,error_v1
180         self.failUnless(error_p<TOL,"pressure error too large.")         self.failUnless(error_p<10*self.TOL, "pressure error too large.")
181         self.failUnless(error_v0<TOL,"0-velocity error too large.")         self.failUnless(error_v0<10*self.TOL, "0-velocity error too large.")
182         self.failUnless(error_v1<TOL,"1-velocity error too large.")         self.failUnless(error_v1<10*self.TOL, "1-velocity error too large.")
183    
184     def test_StokesProblemCartesian_GMRES_P_large(self):     def test_GMRES_P_large(self):
185         ETA=1.         ETA=1.
186         P1=1000.         P1=1000.
187    
# Line 196  class Test_Simple2DModels(unittest.TestC Line 197  class Test_Simple2DModels(unittest.TestC
197         sp.initialize(f=F,fixed_u_mask=mask,eta=ETA)         sp.initialize(f=F,fixed_u_mask=mask,eta=ETA)
198         u0=(1-x[0])*x[0]*[0.,1.]         u0=(1-x[0])*x[0]*[0.,1.]
199         p0=Scalar(P1,ReducedSolution(self.domain))         p0=Scalar(P1,ReducedSolution(self.domain))
200         u,p=sp.solve(u0,p0,show_details=VERBOSE, verbose=VERBOSE,max_iter=100,solver="GMRES")         sp.setTolerance(self.TOL)
201                 u,p=sp.solve(u0,p0,show_details=VERBOSE, verbose=VERBOSE,max_iter=100,useUzawa=False)
        error_v0=Lsup(u[0]-u0[0])  
        error_v1=Lsup(u[1]-u0[1])/0.25  
        zz=P1*x[0]*x[1]-p  
        error_p=Lsup(zz-integrate(zz))/P1  
        # print error_p, error_v0,error_v1  
        self.failUnless(error_p<TOL,"pressure error too large.")  
        self.failUnless(error_v0<TOL,"0-velocity error too large.")  
        self.failUnless(error_v1<TOL,"1-velocity error too large.")  
   
    def test_StokesProblemCartesian_MINRES_P_0(self):  
        ETA=1.  
        P1=0.  
   
        x=self.domain.getX()  
        F=-P1*x[1]*[1.,0]+(2*ETA-P1*x[0])*[0.,1.]  
        mask=whereZero(x[0])    * [1.,1.] \  
               +whereZero(x[0]-1)  * [1.,1.] \  
               +whereZero(x[1])    * [1.,0.] \  
               +whereZero(x[1]-1)  * [1.,1.]  
         
        sp=StokesProblemCartesian(self.domain)  
         
        sp.initialize(f=F,fixed_u_mask=mask,eta=ETA)  
        u0=(1-x[0])*x[0]*[0.,1.]  
        p0=Scalar(P1,ReducedSolution(self.domain))  
        u,p=sp.solve(u0,p0,show_details=VERBOSE, verbose=VERBOSE,max_iter=100,solver="MINRES")  
         
        error_v0=Lsup(u[0]-u0[0])  
        error_v1=Lsup(u[1]-u0[1])/0.25  
        zz=P1*x[0]*x[1]-p  
        error_p=Lsup(zz-integrate(zz))  
        # print error_p, error_v0,error_v1  
        self.failUnless(error_p<TOL,"pressure error too large.")  
        self.failUnless(error_v0<TOL,"0-velocity error too large.")  
        self.failUnless(error_v1<TOL,"1-velocity error too large.")  
   
    def test_StokesProblemCartesian_MINRES_P_small(self):  
        ETA=1.  
        P1=1.  
   
        x=self.domain.getX()  
        F=-P1*x[1]*[1.,0]+(2*ETA-P1*x[0])*[0.,1.]  
        mask=whereZero(x[0])    * [1.,1.] \  
               +whereZero(x[0]-1)  * [1.,1.] \  
               +whereZero(x[1])    * [1.,0.] \  
               +whereZero(x[1]-1)  * [1.,1.]  
         
        sp=StokesProblemCartesian(self.domain)  
         
        sp.initialize(f=F,fixed_u_mask=mask,eta=ETA)  
        u0=(1-x[0])*x[0]*[0.,1.]  
        p0=Scalar(P1,ReducedSolution(self.domain))  
        u,p=sp.solve(u0,p0,show_details=VERBOSE, verbose=VERBOSE,max_iter=100,solver="MINRES")  
202                
203         error_v0=Lsup(u[0]-u0[0])         error_v0=Lsup(u[0]-u0[0])
204         error_v1=Lsup(u[1]-u0[1])/0.25         error_v1=Lsup(u[1]-u0[1])/0.25
205         zz=P1*x[0]*x[1]-p         zz=P1*x[0]*x[1]+p
206         error_p=Lsup(zz-integrate(zz))/P1         error_p=Lsup(zz-integrate(zz))/P1
207         # print error_p, error_v0,error_v1         self.failUnless(error_p<10*self.TOL, "pressure error too large.")
208         self.failUnless(error_p<TOL,"pressure error too large.")         self.failUnless(error_v0<10*self.TOL, "0-velocity error too large.")
209         self.failUnless(error_v0<TOL,"0-velocity error too large.")         self.failUnless(error_v1<10*self.TOL, "1-velocity error too large.")
210         self.failUnless(error_v1<TOL,"1-velocity error too large.")  #====================================================================================================================
211    class Test_StokesProblemCartesian3D(unittest.TestCase):
 #   def test_StokesProblemCartesian_MINRES_P_large(self):  
 #       ETA=1.  
 #       P1=1000.  
 #  
 #       x=self.domain.getX()  
 #       F=-P1*x[1]*[1.,0]+(2*ETA-P1*x[0])*[0.,1.]  
 #       mask=whereZero(x[0])    * [1.,1.] \  
 #              +whereZero(x[0]-1)  * [1.,1.] \  
 #              +whereZero(x[1])    * [1.,0.] \  
 #              +whereZero(x[1]-1)  * [1.,1.]  
         
 #       sp=StokesProblemCartesian(self.domain)  
         
 #       sp.initialize(f=F,fixed_u_mask=mask,eta=ETA)  
 #       u0=(1-x[0])*x[0]*[0.,1.]  
 #       p0=Scalar(P1,ReducedSolution(self.domain))  
 #       u,p=sp.solve(u0,p0,show_details=VERBOSE, verbose=VERBOSE,max_iter=100,solver="MINRES")  
         
 #       error_v0=Lsup(u[0]-u0[0])  
 #       error_v1=Lsup(u[1]-u0[1])/0.25  
 #       zz=P1*x[0]*x[1]-p  
 #       error_p=Lsup(zz-integrate(zz))/P1  
        # print error_p, error_v0,error_v1  
 #       self.failUnless(error_p<TOL,"pressure error too large.")  
 #       self.failUnless(error_v0<TOL,"0-velocity error too large.")  
 #       self.failUnless(error_v1<TOL,"1-velocity error too large.")  
   
   
 #   def test_StokesProblemCartesian_TFQMR_P_0(self):  
 #       ETA=1.  
 #       P1=0.  
   
 #       x=self.domain.getX()  
 #       F=-P1*x[1]*[1.,0]+(2*ETA-P1*x[0])*[0.,1.]  
 #       mask=whereZero(x[0])    * [1.,1.] \  
 #              +whereZero(x[0]-1)  * [1.,1.] \  
 #              +whereZero(x[1])    * [1.,0.] \  
 #              +whereZero(x[1]-1)  * [1.,1.]  
         
 #       sp=StokesProblemCartesian(self.domain)  
         
 #       sp.initialize(f=F,fixed_u_mask=mask,eta=ETA)  
 #       u0=(1-x[0])*x[0]*[0.,1.]  
 #       p0=Scalar(P1,ReducedSolution(self.domain))  
 #       u,p=sp.solve(u0,p0,show_details=VERBOSE, verbose=VERBOSE,max_iter=100,solver="TFQMR")  
         
 #       error_v0=Lsup(u[0]-u0[0])  
 #       error_v1=Lsup(u[1]-u0[1])/0.25  
 #       zz=P1*x[0]*x[1]-p  
 #       error_p=Lsup(zz-integrate(zz))  
        # print error_p, error_v0,error_v1  
 #       self.failUnless(error_p<TOL,"pressure error too large.")  
 #       self.failUnless(error_v0<TOL,"0-velocity error too large.")  
 #       self.failUnless(error_v1<TOL,"1-velocity error too large.")  
   
    def test_StokesProblemCartesian_TFQMR_P_small(self):  
        ETA=1.  
        P1=1.  
   
        x=self.domain.getX()  
        F=-P1*x[1]*[1.,0]+(2*ETA-P1*x[0])*[0.,1.]  
        mask=whereZero(x[0])    * [1.,1.] \  
               +whereZero(x[0]-1)  * [1.,1.] \  
               +whereZero(x[1])    * [1.,0.] \  
               +whereZero(x[1]-1)  * [1.,1.]  
         
        sp=StokesProblemCartesian(self.domain)  
         
        sp.initialize(f=F,fixed_u_mask=mask,eta=ETA)  
        u0=(1-x[0])*x[0]*[0.,1.]  
        p0=Scalar(P1,ReducedSolution(self.domain))  
        u,p=sp.solve(u0,p0,show_details=VERBOSE, verbose=VERBOSE,max_iter=100,solver="TFQMR")  
         
        error_v0=Lsup(u[0]-u0[0])  
        error_v1=Lsup(u[1]-u0[1])/0.25  
        zz=P1*x[0]*x[1]-p  
        error_p=Lsup(zz-integrate(zz))/P1  
        # print error_p, error_v0,error_v1  
        self.failUnless(error_p<TOL,"pressure error too large.")  
        self.failUnless(error_v0<TOL,"0-velocity error too large.")  
        self.failUnless(error_v1<TOL,"1-velocity error too large.")  
   
    def test_StokesProblemCartesian_TFQMR_P_large(self):  
        ETA=1.  
        P1=1000.  
   
        x=self.domain.getX()  
        F=-P1*x[1]*[1.,0]+(2*ETA-P1*x[0])*[0.,1.]  
        mask=whereZero(x[0])    * [1.,1.] \  
               +whereZero(x[0]-1)  * [1.,1.] \  
               +whereZero(x[1])    * [1.,0.] \  
               +whereZero(x[1]-1)  * [1.,1.]  
         
        sp=StokesProblemCartesian(self.domain)  
         
        sp.initialize(f=F,fixed_u_mask=mask,eta=ETA)  
        u0=(1-x[0])*x[0]*[0.,1.]  
        p0=Scalar(P1,ReducedSolution(self.domain))  
        u,p=sp.solve(u0,p0,show_details=VERBOSE, verbose=VERBOSE,max_iter=100,solver="TFQMR")  
         
        error_v0=Lsup(u[0]-u0[0])  
        error_v1=Lsup(u[1]-u0[1])/0.25  
        zz=P1*x[0]*x[1]-p  
        error_p=Lsup(zz-integrate(zz))/P1  
        # print error_p, error_v0,error_v1  
        self.failUnless(error_p<TOL,"pressure error too large.")  
        self.failUnless(error_v0<TOL,"0-velocity error too large.")  
        self.failUnless(error_v1<TOL,"1-velocity error too large.")  
   
   
   
 class Test_Simple3DModels(unittest.TestCase):  
212     def setUp(self):     def setUp(self):
213           NE=6
214           self.TOL=1.e-4
215         self.domain=Brick(NE,NE,NE,order=2,useFullElementOrder=True)         self.domain=Brick(NE,NE,NE,order=2,useFullElementOrder=True)
216     def tearDown(self):     def tearDown(self):
217         del self.domain         del self.domain
218     def test_StokesProblemCartesian_PCG_P_0(self):     def test_PCG_P_0(self):
219         ETA=1.         ETA=1.
220         P1=0.         P1=0.
221    
# Line 397  class Test_Simple3DModels(unittest.TestC Line 235  class Test_Simple3DModels(unittest.TestC
235         sp.initialize(f=F,fixed_u_mask=mask,eta=ETA)         sp.initialize(f=F,fixed_u_mask=mask,eta=ETA)
236         u0=(1-x[0])*x[0]*(1-x[1])*x[1]*[0.,0.,1.]         u0=(1-x[0])*x[0]*(1-x[1])*x[1]*[0.,0.,1.]
237         p0=Scalar(P1,ReducedSolution(self.domain))         p0=Scalar(P1,ReducedSolution(self.domain))
238         u,p=sp.solve(u0,p0,show_details=VERBOSE, verbose=VERBOSE,max_iter=100,solver="PCG")         sp.setTolerance(1.e-8)
239           u,p=sp.solve(u0,p0,show_details=VERBOSE, verbose=VERBOSE,max_iter=100,useUzawa=True)
240                
241         error_v0=Lsup(u[0]-u0[0])         error_v0=Lsup(u[0]-u0[0])
242         error_v1=Lsup(u[1]-u0[1])         error_v1=Lsup(u[1]-u0[1])
243         error_v2=Lsup(u[2]-u0[2])/0.25**2         error_v2=Lsup(u[2]-u0[2])/0.25**2
244         zz=P1*x[0]*x[1]*x[2]-p         zz=P1*x[0]*x[1]*x[2]+p
245         error_p=Lsup(zz-integrate(zz))         error_p=Lsup(zz-integrate(zz))
246         # print error_p, error_v0,error_v1,error_v2         # print error_p, error_v0,error_v1,error_v2
247         self.failUnless(error_p<TOL,"pressure error too large.")         self.failUnless(error_p<10*self.TOL, "pressure error too large.")
248         self.failUnless(error_v0<TOL,"0-velocity error too large.")         self.failUnless(error_v0<10*self.TOL, "0-velocity error too large.")
249         self.failUnless(error_v1<TOL,"1-velocity error too large.")         self.failUnless(error_v1<10*self.TOL, "1-velocity error too large.")
250         self.failUnless(error_v2<TOL,"2-velocity error too large.")         self.failUnless(error_v2<10*self.TOL, "2-velocity error too large.")
251    
252     def test_StokesProblemCartesian_PCG_P_small(self):     def test_PCG_P_small(self):
253         ETA=1.         ETA=1.
254         P1=1.         P1=1.
255    
# Line 429  class Test_Simple3DModels(unittest.TestC Line 268  class Test_Simple3DModels(unittest.TestC
268         sp.initialize(f=F,fixed_u_mask=mask,eta=ETA)         sp.initialize(f=F,fixed_u_mask=mask,eta=ETA)
269         u0=(1-x[0])*x[0]*(1-x[1])*x[1]*[0.,0.,1.]         u0=(1-x[0])*x[0]*(1-x[1])*x[1]*[0.,0.,1.]
270         p0=Scalar(P1,ReducedSolution(self.domain))         p0=Scalar(P1,ReducedSolution(self.domain))
271         u,p=sp.solve(u0,p0,show_details=VERBOSE, verbose=VERBOSE,max_iter=100,solver="PCG")         sp.setTolerance(1.e-8)
272           u,p=sp.solve(u0,p0,show_details=VERBOSE, verbose=VERBOSE,max_iter=100,useUzawa=True)
273                
274         error_v0=Lsup(u[0]-u0[0])         error_v0=Lsup(u[0]-u0[0])
275         error_v1=Lsup(u[1]-u0[1])         error_v1=Lsup(u[1]-u0[1])
276         error_v2=Lsup(u[2]-u0[2])/0.25**2         error_v2=Lsup(u[2]-u0[2])/0.25**2
277         zz=P1*x[0]*x[1]*x[2]-p         zz=P1*x[0]*x[1]*x[2]+p
278         error_p=Lsup(zz-integrate(zz))/P1         error_p=Lsup(zz-integrate(zz))/P1
279         # print error_p, error_v0,error_v1,error_v2         self.failUnless(error_p<10*self.TOL, "pressure error too large.")
280         self.failUnless(error_p<TOL,"pressure error too large.")         self.failUnless(error_v0<10*self.TOL, "0-velocity error too large.")
281         self.failUnless(error_v0<TOL,"0-velocity error too large.")         self.failUnless(error_v1<10*self.TOL, "1-velocity error too large.")
282         self.failUnless(error_v1<TOL,"1-velocity error too large.")         self.failUnless(error_v2<10*self.TOL, "2-velocity error too large.")
283         self.failUnless(error_v2<TOL,"2-velocity error too large.")     def test_PCG_P_large(self):
    def test_StokesProblemCartesian_PCG_P_large(self):  
284         ETA=1.         ETA=1.
285         P1=1000.         P1=1000.
286    
# Line 460  class Test_Simple3DModels(unittest.TestC Line 299  class Test_Simple3DModels(unittest.TestC
299         sp.initialize(f=F,fixed_u_mask=mask,eta=ETA)         sp.initialize(f=F,fixed_u_mask=mask,eta=ETA)
300         u0=(1-x[0])*x[0]*(1-x[1])*x[1]*[0.,0.,1.]         u0=(1-x[0])*x[0]*(1-x[1])*x[1]*[0.,0.,1.]
301         p0=Scalar(P1,ReducedSolution(self.domain))         p0=Scalar(P1,ReducedSolution(self.domain))
302         u,p=sp.solve(u0,p0,show_details=VERBOSE, verbose=VERBOSE,max_iter=100,solver="PCG")         sp.setTolerance(1.e-8)
303           u,p=sp.solve(u0,p0,show_details=VERBOSE, verbose=VERBOSE,max_iter=100,useUzawa=True)
304                
305         error_v0=Lsup(u[0]-u0[0])         error_v0=Lsup(u[0]-u0[0])
306         error_v1=Lsup(u[1]-u0[1])         error_v1=Lsup(u[1]-u0[1])
307         error_v2=Lsup(u[2]-u0[2])/0.25**2         error_v2=Lsup(u[2]-u0[2])/0.25**2
308         zz=P1*x[0]*x[1]*x[2]-p         zz=P1*x[0]*x[1]*x[2]+p
309         error_p=Lsup(zz-integrate(zz))/P1         error_p=Lsup(zz-integrate(zz))/P1
310         # print error_p, error_v0,error_v1,error_v2         # print error_p, error_v0,error_v1,error_v2
311         self.failUnless(error_p<TOL,"pressure error too large.")         self.failUnless(error_p<10*self.TOL, "pressure error too large.")
312         self.failUnless(error_v0<TOL,"0-velocity error too large.")         self.failUnless(error_v0<10*self.TOL, "0-velocity error too large.")
313         self.failUnless(error_v1<TOL,"1-velocity error too large.")         self.failUnless(error_v1<10*self.TOL, "1-velocity error too large.")
314         self.failUnless(error_v2<TOL,"2-velocity error too large.")         self.failUnless(error_v2<10*self.TOL, "2-velocity error too large.")
315    
316     def test_StokesProblemCartesian_GMRES_P_0(self):     def test_GMRES_P_0(self):
317         ETA=1.         ETA=1.
318         P1=0.         P1=0.
319    
# Line 493  class Test_Simple3DModels(unittest.TestC Line 333  class Test_Simple3DModels(unittest.TestC
333         sp.initialize(f=F,fixed_u_mask=mask,eta=ETA)         sp.initialize(f=F,fixed_u_mask=mask,eta=ETA)
334         u0=(1-x[0])*x[0]*(1-x[1])*x[1]*[0.,0.,1.]         u0=(1-x[0])*x[0]*(1-x[1])*x[1]*[0.,0.,1.]
335         p0=Scalar(P1,ReducedSolution(self.domain))         p0=Scalar(P1,ReducedSolution(self.domain))
336         u,p=sp.solve(u0,p0,show_details=VERBOSE, verbose=VERBOSE,max_iter=100,solver="GMRES",iter_restart=20)         sp.setTolerance(1.e-8)
337           u,p=sp.solve(u0,p0,show_details=VERBOSE, verbose=VERBOSE,max_iter=100,useUzawa=False,iter_restart=20)
338                
339         error_v0=Lsup(u[0]-u0[0])         error_v0=Lsup(u[0]-u0[0])
340         error_v1=Lsup(u[1]-u0[1])         error_v1=Lsup(u[1]-u0[1])
341         error_v2=Lsup(u[2]-u0[2])/0.25**2         error_v2=Lsup(u[2]-u0[2])/0.25**2
342         zz=P1*x[0]*x[1]*x[2]-p         zz=P1*x[0]*x[1]*x[2]+p
343         error_p=Lsup(zz-integrate(zz))         error_p=Lsup(zz-integrate(zz))
344         # print error_p, error_v0,error_v1,error_v2         # print error_p, error_v0,error_v1,error_v2
345         self.failUnless(error_p<TOL,"pressure error too large.")         self.failUnless(error_p<10*self.TOL, "pressure error too large.")
346         self.failUnless(error_v0<TOL,"0-velocity error too large.")         self.failUnless(error_v0<10*self.TOL, "0-velocity error too large.")
347         self.failUnless(error_v1<TOL,"1-velocity error too large.")         self.failUnless(error_v1<10*self.TOL, "1-velocity error too large.")
348         self.failUnless(error_v2<TOL,"2-velocity error too large.")         self.failUnless(error_v2<10*self.TOL, "2-velocity error too large.")
349     def test_StokesProblemCartesian_GMRES_P_small(self):     def test_GMRES_P_small(self):
350         ETA=1.         ETA=1.
351         P1=1.         P1=1.
352    
# Line 524  class Test_Simple3DModels(unittest.TestC Line 365  class Test_Simple3DModels(unittest.TestC
365         sp.initialize(f=F,fixed_u_mask=mask,eta=ETA)         sp.initialize(f=F,fixed_u_mask=mask,eta=ETA)
366         u0=(1-x[0])*x[0]*(1-x[1])*x[1]*[0.,0.,1.]         u0=(1-x[0])*x[0]*(1-x[1])*x[1]*[0.,0.,1.]
367         p0=Scalar(P1,ReducedSolution(self.domain))         p0=Scalar(P1,ReducedSolution(self.domain))
368         u,p=sp.solve(u0,p0,show_details=VERBOSE, verbose=VERBOSE,max_iter=100,solver="GMRES")         sp.setTolerance(1.e-8)
369           u,p=sp.solve(u0,p0,show_details=VERBOSE, verbose=VERBOSE,max_iter=100,useUzawa=False)
370                
371         error_v0=Lsup(u[0]-u0[0])         error_v0=Lsup(u[0]-u0[0])
372         error_v1=Lsup(u[1]-u0[1])         error_v1=Lsup(u[1]-u0[1])
373         error_v2=Lsup(u[2]-u0[2])/0.25**2         error_v2=Lsup(u[2]-u0[2])/0.25**2
374         zz=P1*x[0]*x[1]*x[2]-p         zz=P1*x[0]*x[1]*x[2]+p
375         error_p=Lsup(zz-integrate(zz))/P1         error_p=Lsup(zz-integrate(zz))/P1
376         # print error_p, error_v0,error_v1,error_v2         # print error_p, error_v0,error_v1,error_v2
377         self.failUnless(error_p<TOL,"pressure error too large.")         self.failUnless(error_p<10*self.TOL, "pressure error too large.")
378         self.failUnless(error_v0<TOL,"0-velocity error too large.")         self.failUnless(error_v0<10*self.TOL, "0-velocity error too large.")
379         self.failUnless(error_v1<TOL,"1-velocity error too large.")         self.failUnless(error_v1<10*self.TOL, "1-velocity error too large.")
380         self.failUnless(error_v2<TOL,"2-velocity error too large.")         self.failUnless(error_v2<10*self.TOL, "2-velocity error too large.")
381     def test_StokesProblemCartesian_GMRES_P_large(self):     def test_GMRES_P_large(self):
382         ETA=1.         ETA=1.
383         P1=1000.         P1=1000.
384    
# Line 555  class Test_Simple3DModels(unittest.TestC Line 397  class Test_Simple3DModels(unittest.TestC
397         sp.initialize(f=F,fixed_u_mask=mask,eta=ETA)         sp.initialize(f=F,fixed_u_mask=mask,eta=ETA)
398         u0=(1-x[0])*x[0]*(1-x[1])*x[1]*[0.,0.,1.]         u0=(1-x[0])*x[0]*(1-x[1])*x[1]*[0.,0.,1.]
399         p0=Scalar(P1,ReducedSolution(self.domain))         p0=Scalar(P1,ReducedSolution(self.domain))
400         u,p=sp.solve(u0,p0,show_details=VERBOSE, verbose=VERBOSE,max_iter=100,solver="GMRES")         sp.setTolerance(self.TOL+0*1.e-8)
401                 # sp.setSubToleranceReductionFactor(1.e-8/self.TOL)
402         error_v0=Lsup(u[0]-u0[0])         # sp.setSubToleranceReductionFactor(None)
403         error_v1=Lsup(u[1]-u0[1])         u,p=sp.solve(u0,p0,show_details=VERBOSE, verbose=VERBOSE,max_iter=100,useUzawa=False)
        error_v2=Lsup(u[2]-u0[2])/0.25**2  
        zz=P1*x[0]*x[1]*x[2]-p  
        error_p=Lsup(zz-integrate(zz))/P1  
        # print error_p, error_v0,error_v1,error_v2  
        self.failUnless(error_p<TOL,"pressure error too large.")  
        self.failUnless(error_v0<TOL,"0-velocity error too large.")  
        self.failUnless(error_v1<TOL,"1-velocity error too large.")  
        self.failUnless(error_v2<TOL,"2-velocity error too large.")  
   
    def test_StokesProblemCartesian_MINRES_P_0(self):  
        ETA=1.  
        P1=0.  
   
        x=self.domain.getX()  
        F=-P1*x[1]*x[2]*[1.,0.,0.]-P1*x[0]*x[2]*[0.,1.,0.]+(2*ETA*((1-x[0])*x[0]+(1-x[1])*x[1])-P1*x[0]*x[1])*[0.,0.,1.]  
        x=self.domain.getX()  
        mask=whereZero(x[0])    * [1.,1.,1.] \  
               +whereZero(x[0]-1)  * [1.,1.,1.] \  
               +whereZero(x[1])    * [1.,1.,1.] \  
               +whereZero(x[1]-1)  * [1.,1.,1.] \  
               +whereZero(x[2])    * [1.,1.,0.] \  
               +whereZero(x[2]-1)  * [1.,1.,1.]  
         
         
        sp=StokesProblemCartesian(self.domain)  
         
        sp.initialize(f=F,fixed_u_mask=mask,eta=ETA)  
        u0=(1-x[0])*x[0]*(1-x[1])*x[1]*[0.,0.,1.]  
        p0=Scalar(P1,ReducedSolution(self.domain))  
        u,p=sp.solve(u0,p0,show_details=VERBOSE, verbose=VERBOSE,max_iter=100,solver="MINRES")  
404                
405         error_v0=Lsup(u[0]-u0[0])         error_v0=Lsup(u[0]-u0[0])
406         error_v1=Lsup(u[1]-u0[1])         error_v1=Lsup(u[1]-u0[1])
407         error_v2=Lsup(u[2]-u0[2])/0.25**2         error_v2=Lsup(u[2]-u0[2])/0.25**2
408         zz=P1*x[0]*x[1]*x[2]-p         zz=P1*x[0]*x[1]*x[2]+p
        error_p=Lsup(zz-integrate(zz))  
        # print error_p, error_v0,error_v1,error_v2  
        self.failUnless(error_p<TOL,"pressure error too large.")  
        self.failUnless(error_v0<TOL,"0-velocity error too large.")  
        self.failUnless(error_v1<TOL,"1-velocity error too large.")  
        self.failUnless(error_v2<TOL,"2-velocity error too large.")  
   
    def test_StokesProblemCartesian_MINRES_P_small(self):  
        ETA=1.  
        P1=1.  
   
        x=self.domain.getX()  
        F=-P1*x[1]*x[2]*[1.,0.,0.]-P1*x[0]*x[2]*[0.,1.,0.]+(2*ETA*((1-x[0])*x[0]+(1-x[1])*x[1])-P1*x[0]*x[1])*[0.,0.,1.]  
        mask=whereZero(x[0])    * [1.,1.,1.] \  
               +whereZero(x[0]-1)  * [1.,1.,1.] \  
               +whereZero(x[1])    * [1.,1.,1.] \  
               +whereZero(x[1]-1)  * [1.,1.,1.] \  
               +whereZero(x[2])    * [1.,1.,0.] \  
               +whereZero(x[2]-1)  * [1.,1.,1.]  
         
         
        sp=StokesProblemCartesian(self.domain)  
         
        sp.initialize(f=F,fixed_u_mask=mask,eta=ETA)  
        u0=(1-x[0])*x[0]*(1-x[1])*x[1]*[0.,0.,1.]  
        p0=Scalar(P1,ReducedSolution(self.domain))  
        u,p=sp.solve(u0,p0,show_details=VERBOSE, verbose=VERBOSE,max_iter=100,solver="MINRES")  
         
        error_v0=Lsup(u[0]-u0[0])  
        error_v1=Lsup(u[1]-u0[1])  
        error_v2=Lsup(u[2]-u0[2])/0.25**2  
        zz=P1*x[0]*x[1]*x[2]-p  
409         error_p=Lsup(zz-integrate(zz))/P1         error_p=Lsup(zz-integrate(zz))/P1
410         # print error_p, error_v0,error_v1,error_v2         # print error_p, error_v0,error_v1,error_v2
411         self.failUnless(error_p<TOL,"pressure error too large.")         self.failUnless(error_p<10*self.TOL, "pressure error too large.")
412         self.failUnless(error_v0<TOL,"0-velocity error too large.")         self.failUnless(error_v0<10*self.TOL, "0-velocity error too large.")
413         self.failUnless(error_v1<TOL,"1-velocity error too large.")         self.failUnless(error_v1<10*self.TOL, "1-velocity error too large.")
414         self.failUnless(error_v2<TOL,"2-velocity error too large.")         self.failUnless(error_v2<10*self.TOL, "2-velocity error too large.")
415    #====================================================================================================================
416  #   def test_StokesProblemCartesian_MINRES_P_large(self):  class Test_Darcy(unittest.TestCase):
417  #       ETA=1.      # this is a simple test for the darcy flux problem
418  #       P1=1000.      #
419        #
420  #       x=self.domain.getX()      #  p = 1/k * ( 1/2* (fb-f0)/xb* x **2 + f0 * x - ub*x ) +  p0
421  #       F=-P1*x[1]*x[2]*[1.,0.,0.]-P1*x[0]*x[2]*[0.,1.,0.]+(2*ETA*((1-x[0])*x[0]+(1-x[1])*x[1])-P1*x[0]*x[1])*[0.,0.,1.]      #
422  #       mask=whereZero(x[0])    * [1.,1.,1.] \      #  with f = (fb-f0)/xb* x + f0
423  #              +whereZero(x[0]-1)  * [1.,1.,1.] \      #
424  #              +whereZero(x[1])    * [1.,1.,1.] \      #    u = f - k * p,x = ub
425  #              +whereZero(x[1]-1)  * [1.,1.,1.] \      #
426  #              +whereZero(x[2])    * [1.,1.,0.] \      #  we prescribe pressure at x=x0=0 to p0
427  #              +whereZero(x[2]-1)  * [1.,1.,1.]      #
428              #  if we prescribe the pressure on the bottom x=xb we set
429              #
430  #       sp=StokesProblemCartesian(self.domain)      #  pb= 1/k * ( 1/2* (fb-f0)* xb + f0 * xb - ub*xb ) +  p0 = 1/k * ((fb+f0)/2 - ub ) * xb  + p0
431              #
432  #       sp.initialize(f=F,fixed_u_mask=mask,eta=ETA)      #  which leads to ub = (fb+f0)/2-k*(pb-p0)/xb
433  #       u0=(1-x[0])*x[0]*(1-x[1])*x[1]*[0.,0.,1.]      #
434  #       p0=Scalar(P1,ReducedSolution(self.domain))      def rescaleDomain(self):
435  #       u,p=sp.solve(u0,p0,show_details=VERBOSE, verbose=VERBOSE,max_iter=100,solver="MINRES")          x=self.dom.getX().copy()
436                  for i in xrange(self.dom.getDim()):
437  #       error_v0=Lsup(u[0]-u0[0])               x_inf=inf(x[i])
438  #       error_v1=Lsup(u[1]-u0[1])               x_sup=sup(x[i])
439  #       error_v2=Lsup(u[2]-u0[2])/0.25**2               if i == self.dom.getDim()-1:
440  #       zz=P1*x[0]*x[1]*x[2]-p                  x[i]=-self.WIDTH*(x[i]-x_sup)/(x_inf-x_sup)
441  #       error_p=Lsup(zz-integrate(zz))/P1               else:
442         # print error_p, error_v0,error_v1,error_v2                  x[i]=self.WIDTH*(x[i]-x_inf)/(x_sup-x_inf)
443  #       self.failUnless(error_p<TOL,"pressure error too large.")          self.dom.setX(x)
444  #       self.failUnless(error_v0<TOL,"0-velocity error too large.")      def getScalarMask(self,include_bottom=True):
445  #       self.failUnless(error_v1<TOL,"1-velocity error too large.")          x=self.dom.getX().copy()
446  #       self.failUnless(error_v2<TOL,"2-velocity error too large.")          x_inf=inf(x[self.dom.getDim()-1])
447            x_sup=sup(x[self.dom.getDim()-1])
448  #   def test_StokesProblemCartesian_TFQMR_P_0(self):          out=whereZero(x[self.dom.getDim()-1]-x_sup)
449  #       ETA=1.          if include_bottom: out+=whereZero(x[self.dom.getDim()-1]-x_inf)
450  #       P1=0.          return wherePositive(out)
451        def getVectorMask(self,include_bottom=True):
452  #       x=self.domain.getX()          x=self.dom.getX().copy()
453  #       F=-P1*x[1]*x[2]*[1.,0.,0.]-P1*x[0]*x[2]*[0.,1.,0.]+(2*ETA*((1-x[0])*x[0]+(1-x[1])*x[1])-P1*x[0]*x[1])*[0.,0.,1.]          out=Vector(0.,Solution(self.dom))
454  #       x=self.domain.getX()          for i in xrange(self.dom.getDim()):
455  #       mask=whereZero(x[0])    * [1.,1.,1.] \               x_inf=inf(x[i])
456  #              +whereZero(x[0]-1)  * [1.,1.,1.] \               x_sup=sup(x[i])
457  #              +whereZero(x[1])    * [1.,1.,1.] \               if i != self.dom.getDim()-1: out[i]+=whereZero(x[i]-x_sup)
458  #              +whereZero(x[1]-1)  * [1.,1.,1.] \               if i != self.dom.getDim()-1 or include_bottom: out[i]+=whereZero(x[i]-x_inf)
459  #              +whereZero(x[2])    * [1.,1.,0.] \          return wherePositive(out)
460  #              +whereZero(x[2]-1)  * [1.,1.,1.]  
461              def setSolutionFixedBottom(self, p0, pb, f0, fb, k):
462                   d=self.dom.getDim()
463  #       sp=StokesProblemCartesian(self.domain)           x=self.dom.getX()[d-1]
464                   xb=inf(x)
465  #       sp.initialize(f=F,fixed_u_mask=mask,eta=ETA)           u=Vector(0.,Solution(self.dom))+kronecker(d)[d-1]*((f0+fb)/2.-k*(pb-p0)/xb)
466  #       u0=(1-x[0])*x[0]*(1-x[1])*x[1]*[0.,0.,1.]           p=1./k*((fb-f0)/(xb*2.)* x**2 - (fb-f0)/2.*x)+(pb-p0)/xb*x +  p0
467  #       p0=Scalar(P1,ReducedSolution(self.domain))           f= ((fb-f0)/xb* x + f0)*kronecker(Function(self.dom))[d-1]
468  #       u,p=sp.solve(u0,p0,show_details=VERBOSE, verbose=VERBOSE,max_iter=100,solver="TFQMR")           return u,p,f
469                  
470  #       error_v0=Lsup(u[0]-u0[0])      def testConstF_FixedBottom_smallK(self):
471  #       error_v1=Lsup(u[1]-u0[1])          k=1.e-10
472  #       error_v2=Lsup(u[2]-u0[2])/0.25**2          mp=self.getScalarMask(include_bottom=True)
473  #       zz=P1*x[0]*x[1]*x[2]-p          mv=self.getVectorMask(include_bottom=False)
474  #       error_p=Lsup(zz-integrate(zz))          u_ref,p_ref,f=self.setSolutionFixedBottom(p0=2500,pb=4000.,f0=10.,fb=10.,k=k)
475         # print error_p, error_v0,error_v1,error_v2          p=p_ref*mp
476  #       self.failUnless(error_p<TOL,"pressure error too large.")          u=u_ref*mv
477  #       self.failUnless(error_v0<TOL,"0-velocity error too large.")          df=DarcyFlow(self.dom)
478  #       self.failUnless(error_v1<TOL,"1-velocity error too large.")          df.setValue(g=f,
479  #       self.failUnless(error_v2<TOL,"2-velocity error too large.")                        location_of_fixed_pressure=mp,
480                          location_of_fixed_flux=mv,
481     def test_StokesProblemCartesian_TFQMR_P_small(self):                        permeability=Scalar(k,Function(self.dom)))
482         ETA=1.          v,p=df.solve(u,p,atol=0,rtol=self.TOL, max_iter=100, verbose=VERBOSE,sub_rtol=self.TOL/200)
483         P1=1.          self.failUnless(Lsup(v-u_ref)<self.TOL*10.*Lsup(u_ref), "flux error too big.")
484            self.failUnless(Lsup(p-p_ref)<self.TOL*10.*Lsup(p_ref), "pressure error too big.")
485         x=self.domain.getX()  
486         F=-P1*x[1]*x[2]*[1.,0.,0.]-P1*x[0]*x[2]*[0.,1.,0.]+(2*ETA*((1-x[0])*x[0]+(1-x[1])*x[1])-P1*x[0]*x[1])*[0.,0.,1.]      def testConstF_FixedBottom_mediumK(self):
487         mask=whereZero(x[0])    * [1.,1.,1.] \          k=1.
488                +whereZero(x[0]-1)  * [1.,1.,1.] \          mp=self.getScalarMask(include_bottom=True)
489                +whereZero(x[1])    * [1.,1.,1.] \          mv=self.getVectorMask(include_bottom=False)
490                +whereZero(x[1]-1)  * [1.,1.,1.] \          u_ref,p_ref,f=self.setSolutionFixedBottom(p0=2500,pb=4000.,f0=10.,fb=10.,k=k)
491                +whereZero(x[2])    * [1.,1.,0.] \          p=p_ref*mp
492                +whereZero(x[2]-1)  * [1.,1.,1.]          u=u_ref*mv
493                  df=DarcyFlow(self.dom)
494                  df.setValue(g=f,
495         sp=StokesProblemCartesian(self.domain)                        location_of_fixed_pressure=mp,
496                                location_of_fixed_flux=mv,
497         sp.initialize(f=F,fixed_u_mask=mask,eta=ETA)                        permeability=Scalar(k,Function(self.dom)))
498         u0=(1-x[0])*x[0]*(1-x[1])*x[1]*[0.,0.,1.]          v,p=df.solve(u,p,atol=0,rtol=self.TOL, max_iter=100, verbose=VERBOSE,sub_rtol=self.TOL/200)
499         p0=Scalar(P1,ReducedSolution(self.domain))          self.failUnless(Lsup(v-u_ref)<self.TOL*10.*Lsup(u_ref), "flux error too big.")
500         u,p=sp.solve(u0,p0,show_details=VERBOSE, verbose=VERBOSE,max_iter=100,solver="TFQMR")          self.failUnless(Lsup(p-p_ref)<self.TOL*10.*Lsup(p_ref), "pressure error too big.")
501          
502         error_v0=Lsup(u[0]-u0[0])      def testConstF_FixedBottom_largeK(self):
503         error_v1=Lsup(u[1]-u0[1])          k=1.e10
504         error_v2=Lsup(u[2]-u0[2])/0.25**2          mp=self.getScalarMask(include_bottom=True)
505         zz=P1*x[0]*x[1]*x[2]-p          mv=self.getVectorMask(include_bottom=False)
506         error_p=Lsup(zz-integrate(zz))/P1          u_ref,p_ref,f=self.setSolutionFixedBottom(p0=2500,pb=4000.,f0=10.,fb=10.,k=k)
507         # print error_p, error_v0,error_v1,error_v2          p=p_ref*mp
508         self.failUnless(error_p<TOL,"pressure error too large.")          u=u_ref*mv
509         self.failUnless(error_v0<TOL,"0-velocity error too large.")          df=DarcyFlow(self.dom)
510         self.failUnless(error_v1<TOL,"1-velocity error too large.")          df.setValue(g=f,
511         self.failUnless(error_v2<TOL,"2-velocity error too large.")                        location_of_fixed_pressure=mp,
512                          location_of_fixed_flux=mv,
513     def test_StokesProblemCartesian_TFQMR_P_large(self):                        permeability=Scalar(k,Function(self.dom)))
514         ETA=1.          v,p=df.solve(u,p,atol=0,rtol=self.TOL, max_iter=100, verbose=VERBOSE,sub_rtol=self.TOL/200)
515         P1=1000.          self.failUnless(Lsup(v-u_ref)<self.TOL*10.*Lsup(u_ref), "flux error too big.")
516            self.failUnless(Lsup(p-p_ref)<self.TOL*10.*Lsup(p_ref), "pressure error too big.")
517    
518        def testVarioF_FixedBottom_smallK(self):
519            k=1.e-10
520            mp=self.getScalarMask(include_bottom=True)
521            mv=self.getVectorMask(include_bottom=False)
522            u_ref,p_ref,f=self.setSolutionFixedBottom(p0=2500,pb=4000.,f0=10.,fb=30.,k=k)
523            p=p_ref*mp
524            u=u_ref*mv
525            df=DarcyFlow(self.dom)
526            df.setValue(g=f,
527                          location_of_fixed_pressure=mp,
528                          location_of_fixed_flux=mv,
529                          permeability=Scalar(k,Function(self.dom)))
530            v,p=df.solve(u,p,atol=0,rtol=self.TOL, max_iter=100, verbose=VERBOSE,sub_rtol=self.TOL/200)
531            self.failUnless(Lsup(v-u_ref)<self.TOL*10.*Lsup(u_ref), "flux error too big.")
532            self.failUnless(Lsup(p-p_ref)<self.TOL*10.*Lsup(p_ref), "pressure error too big.")
533    
534        def testVarioF_FixedBottom_mediumK(self):
535            k=1.
536            mp=self.getScalarMask(include_bottom=True)
537            mv=self.getVectorMask(include_bottom=False)
538            u_ref,p_ref,f=self.setSolutionFixedBottom(p0=2500,pb=4000.,f0=10.,fb=30.,k=k)
539            p=p_ref*mp
540            u=u_ref*mv
541            df=DarcyFlow(self.dom)
542            df.setValue(g=f,
543                          location_of_fixed_pressure=mp,
544                          location_of_fixed_flux=mv,
545                          permeability=Scalar(k,Function(self.dom)))
546            v,p=df.solve(u,p,atol=0,rtol=self.TOL, max_iter=100, verbose=VERBOSE,sub_rtol=self.TOL/200)
547            self.failUnless(Lsup(v-u_ref)<self.TOL*10.*Lsup(u_ref), "flux error too big.")
548            self.failUnless(Lsup(p-p_ref)<self.TOL*10.*Lsup(p_ref), "pressure error too big.")
549    
550        def testVarioF_FixedBottom_largeK(self):
551            k=1.e10
552            mp=self.getScalarMask(include_bottom=True)
553            mv=self.getVectorMask(include_bottom=False)
554            u_ref,p_ref,f=self.setSolutionFixedBottom(p0=2500,pb=4000.,f0=10.,fb=30.,k=k)
555            p=p_ref*mp
556            u=u_ref*mv
557            df=DarcyFlow(self.dom)
558            df.setValue(g=f,
559                          location_of_fixed_pressure=mp,
560                          location_of_fixed_flux=mv,
561                          permeability=Scalar(k,Function(self.dom)))
562            v,p=df.solve(u,p,atol=0,rtol=self.TOL, max_iter=100, verbose=VERBOSE,sub_rtol=self.TOL/200)
563            self.failUnless(Lsup(v-u_ref)<self.TOL*10.*Lsup(u_ref), "flux error too big.")
564            self.failUnless(Lsup(p-p_ref)<self.TOL*10.*Lsup(p_ref), "pressure error too big.")
565    
566        def testConstF_FreeBottom_smallK(self):
567            k=1.e-10
568            mp=self.getScalarMask(include_bottom=False)
569            mv=self.getVectorMask(include_bottom=True)
570            u_ref,p_ref,f=self.setSolutionFixedBottom(p0=2500,pb=4000.,f0=10.,fb=10.,k=k)
571            p=p_ref*mp
572            u=u_ref*mv
573            df=DarcyFlow(self.dom)
574            df.setValue(g=f,
575                          location_of_fixed_pressure=mp,
576                          location_of_fixed_flux=mv,
577                          permeability=Scalar(k,Function(self.dom)))
578            v,p=df.solve(u,p,atol=0,rtol=self.TOL, max_iter=100, verbose=VERBOSE,sub_rtol=self.TOL/200)
579            self.failUnless(Lsup(v-u_ref)<self.TOL*10.*Lsup(u_ref), "flux error too big.")
580            self.failUnless(Lsup(p-p_ref)<self.TOL*10.*Lsup(p_ref), "pressure error too big.")
581    
582        def testConstF_FreeBottom_mediumK(self):
583            k=1.
584            mp=self.getScalarMask(include_bottom=False)
585            mv=self.getVectorMask(include_bottom=True)
586            u_ref,p_ref,f=self.setSolutionFixedBottom(p0=2500,pb=4000.,f0=10.,fb=10.,k=k)
587            p=p_ref*mp
588            u=u_ref*mv
589            df=DarcyFlow(self.dom)
590            df.setValue(g=f,
591                          location_of_fixed_pressure=mp,
592                          location_of_fixed_flux=mv,
593                          permeability=Scalar(k,Function(self.dom)))
594            v,p=df.solve(u,p,atol=0,rtol=self.TOL, max_iter=100, verbose=VERBOSE,sub_rtol=self.TOL/200)
595            self.failUnless(Lsup(v-u_ref)<self.TOL*10.*Lsup(u_ref), "flux error too big.")
596            self.failUnless(Lsup(p-p_ref)<self.TOL*10.*Lsup(p_ref), "pressure error too big.")
597    
598        def testConstF_FreeBottom_largeK(self):
599            k=1.e10
600            mp=self.getScalarMask(include_bottom=False)
601            mv=self.getVectorMask(include_bottom=True)
602            u_ref,p_ref,f=self.setSolutionFixedBottom(p0=2500,pb=4000.,f0=10.,fb=10.,k=k)
603            p=p_ref*mp
604            u=u_ref*mv
605            df=DarcyFlow(self.dom)
606            df.setValue(g=f,
607                          location_of_fixed_pressure=mp,
608                          location_of_fixed_flux=mv,
609                          permeability=Scalar(k,Function(self.dom)))
610            v,p=df.solve(u,p,atol=0,rtol=self.TOL, max_iter=100, verbose=VERBOSE,sub_rtol=self.TOL/200)
611            self.failUnless(Lsup(v-u_ref)<self.TOL*10.*Lsup(u_ref), "flux error too big.")
612            self.failUnless(Lsup(p-p_ref)<self.TOL*10.*Lsup(p_ref), "pressure error too big.")
613    
614        def testVarioF_FreeBottom_smallK(self):
615            k=1.e-10
616            mp=self.getScalarMask(include_bottom=False)
617            mv=self.getVectorMask(include_bottom=True)
618            u_ref,p_ref,f=self.setSolutionFixedBottom(p0=2500,pb=4000.,f0=10.,fb=30.,k=k)
619            p=p_ref*mp
620            u=u_ref*mv
621            df=DarcyFlow(self.dom)
622            df.setValue(g=f,
623                          location_of_fixed_pressure=mp,
624                          location_of_fixed_flux=mv,
625                          permeability=Scalar(k,Function(self.dom)))
626            v,p=df.solve(u,p,atol=0,rtol=self.TOL, max_iter=100, verbose=VERBOSE,sub_rtol=self.TOL/200)
627            self.failUnless(Lsup(v-u_ref)<self.TOL*10.*Lsup(u_ref), "flux error too big.")
628            self.failUnless(Lsup(p-p_ref)<self.TOL*10.*Lsup(p_ref), "pressure error too big.")
629    
630        def testVarioF_FreeBottom_mediumK(self):
631            k=1.
632            mp=self.getScalarMask(include_bottom=False)
633            mv=self.getVectorMask(include_bottom=True)
634            u_ref,p_ref,f=self.setSolutionFixedBottom(p0=2500,pb=4000.,f0=10.,fb=30.,k=k)
635            p=p_ref*mp
636            u=u_ref*mv
637            df=DarcyFlow(self.dom)
638            df.setValue(g=f,
639                          location_of_fixed_pressure=mp,
640                          location_of_fixed_flux=mv,
641                          permeability=Scalar(k,Function(self.dom)))
642            v,p=df.solve(u,p,atol=0,rtol=self.TOL, max_iter=100, verbose=VERBOSE,sub_rtol=self.TOL/200)
643            self.failUnless(Lsup(v-u_ref)<self.TOL*10.*Lsup(u_ref), "flux error too big.")
644            self.failUnless(Lsup(p-p_ref)<self.TOL*10.*Lsup(p_ref), "pressure error too big.")
645    
646        def testVarioF_FreeBottom_largeK(self):
647            k=1.e10
648            mp=self.getScalarMask(include_bottom=False)
649            mv=self.getVectorMask(include_bottom=True)
650            u_ref,p_ref,f=self.setSolutionFixedBottom(p0=2500,pb=4000.,f0=10.,fb=30.,k=k)
651            p=p_ref*mp
652            u=u_ref*mv
653            df=DarcyFlow(self.dom)
654            df.setValue(g=f,
655                          location_of_fixed_pressure=mp,
656                          location_of_fixed_flux=mv,
657                          permeability=Scalar(k,Function(self.dom)))
658            v,p=df.solve(u,p,atol=0,rtol=self.TOL, max_iter=100, verbose=VERBOSE,sub_rtol=self.TOL/200)
659            self.failUnless(Lsup(v-u_ref)<self.TOL*10.*Lsup(u_ref), "flux error too big.")
660            self.failUnless(Lsup(p-p_ref)<self.TOL*10.*Lsup(p_ref), "pressure error too big.")
661    
662    class Test_Darcy2D(Test_Darcy):
663        TOL=1e-5
664        WIDTH=1.
665        def setUp(self):
666            NE=60  # wrning smaller NE may case a failure for VarioF tests due to discretization errors.
667            self.dom = Rectangle(NE,NE)
668            self.rescaleDomain()
669        def tearDown(self):
670            del self.dom
671    class Test_Darcy3D(Test_Darcy):
672        TOL=1e-4
673        WIDTH=1.
674        def setUp(self):
675            NE=25  # wrning smaller NE may case a failure for VarioF tests due to discretization errors.
676            self.dom = Brick(NE,NE,NE)
677            self.rescaleDomain()
678        def tearDown(self):
679            del self.dom
680    
        x=self.domain.getX()  
        F=-P1*x[1]*x[2]*[1.,0.,0.]-P1*x[0]*x[2]*[0.,1.,0.]+(2*ETA*((1-x[0])*x[0]+(1-x[1])*x[1])-P1*x[0]*x[1])*[0.,0.,1.]  
        mask=whereZero(x[0])    * [1.,1.,1.] \  
               +whereZero(x[0]-1)  * [1.,1.,1.] \  
               +whereZero(x[1])    * [1.,1.,1.] \  
               +whereZero(x[1]-1)  * [1.,1.,1.] \  
               +whereZero(x[2])    * [1.,1.,0.] \  
               +whereZero(x[2]-1)  * [1.,1.,1.]  
         
         
        sp=StokesProblemCartesian(self.domain)  
         
        sp.initialize(f=F,fixed_u_mask=mask,eta=ETA)  
        u0=(1-x[0])*x[0]*(1-x[1])*x[1]*[0.,0.,1.]  
        p0=Scalar(P1,ReducedSolution(self.domain))  
        u,p=sp.solve(u0,p0,show_details=VERBOSE, verbose=VERBOSE,max_iter=100,solver="TFQMR")  
         
        error_v0=Lsup(u[0]-u0[0])  
        error_v1=Lsup(u[1]-u0[1])  
        error_v2=Lsup(u[2]-u0[2])/0.25**2  
        zz=P1*x[0]*x[1]*x[2]-p  
        error_p=Lsup(zz-integrate(zz))/P1  
        # print error_p, error_v0,error_v1,error_v2  
        self.failUnless(error_p<TOL,"pressure error too large.")  
        self.failUnless(error_v0<TOL,"0-velocity error too large.")  
        self.failUnless(error_v1<TOL,"1-velocity error too large.")  
        self.failUnless(error_v2<TOL,"2-velocity error too large.")  
681    
682    
683  if __name__ == '__main__':  if __name__ == '__main__':
684     suite = unittest.TestSuite()     suite = unittest.TestSuite()
685     suite.addTest(unittest.makeSuite(Test_Simple2DModels))     suite.addTest(unittest.makeSuite(Test_StokesProblemCartesian2D))
686     # suite.addTest(unittest.makeSuite(Test_Simple3DModels))     suite.addTest(unittest.makeSuite(Test_Darcy3D))
687       suite.addTest(unittest.makeSuite(Test_Darcy2D))
688       suite.addTest(unittest.makeSuite(Test_StokesProblemCartesian3D))
689     s=unittest.TextTestRunner(verbosity=2).run(suite)     s=unittest.TextTestRunner(verbosity=2).run(suite)
690     if not s.wasSuccessful(): sys.exit(1)     if not s.wasSuccessful(): sys.exit(1)
691    

Legend:
Removed from v.2099  
changed lines
  Added in v.2100

  ViewVC Help
Powered by ViewVC 1.1.26