/[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 1562 by gross, Wed May 21 13:04:40 2008 UTC revision 2389 by gross, Wed Apr 15 06:46:01 2009 UTC
# Line 1  Line 1 
1  #######################################################  
2  #  ########################################################
 #           Copyright 2003-2007 by ACceSS MNRF  
 #       Copyright 2007 by University of Queensland  
3  #  #
4  #                http://esscc.uq.edu.au  # Copyright (c) 2003-2008 by University of Queensland
5  #        Primary Business: Queensland, Australia  # Earth Systems Science Computational Center (ESSCC)
6  #  Licensed under the Open Software License version 3.0  # http://www.uq.edu.au/esscc
 #     http://www.opensource.org/licenses/osl-3.0.php  
7  #  #
8  #######################################################  # Primary Business: Queensland, Australia
9    # Licensed under the Open Software License version 3.0
10    # http://www.opensource.org/licenses/osl-3.0.php
11  #  #
12    ########################################################
13    
14  __copyright__="""  Copyright (c) 2006 by ACcESS MNRF  __copyright__="""Copyright (c) 2003-2008 by University of Queensland
15                      http://www.access.edu.au  Earth Systems Science Computational Center (ESSCC)
16                  Primary Business: Queensland, Australia"""  http://www.uq.edu.au/esscc
17    Primary Business: Queensland, Australia"""
18  __license__="""Licensed under the Open Software License version 3.0  __license__="""Licensed under the Open Software License version 3.0
19               http://www.opensource.org/licenses/osl-3.0.php"""  http://www.opensource.org/licenses/osl-3.0.php"""
20    __url__="https://launchpad.net/escript-finley"
21    
22  import unittest  import unittest
23  import tempfile  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 29  except KeyError: Line 33  except KeyError:
33       FINLEY_WORKDIR='.'       FINLEY_WORKDIR='.'
34    
35    
36  NE=6  VERBOSE=False # or True
37  TOL=1.e-5  DETAIL_VERBOSE=False
 VERBOSE=False or True  
38    
39  from esys.escript import *  from esys.escript import *
40  from esys.escript.models import StokesProblemCartesian  from esys.escript.models import StokesProblemCartesian, PowerLaw
41  from esys.finley import Rectangle, Brick  from esys.finley import Rectangle, Brick
42  class Test_Simple2DModels(unittest.TestCase):  
43    from esys.escript.models import Mountains
44    from math import pi
45    
46    #====================================================================================================================
47    class Test_StokesProblemCartesian2D(unittest.TestCase):
48     def setUp(self):     def setUp(self):
49           NE=6
50           self.TOL=1e-3
51         self.domain=Rectangle(NE,NE,order=2,useFullElementOrder=True)         self.domain=Rectangle(NE,NE,order=2,useFullElementOrder=True)
52     def tearDown(self):     def tearDown(self):
53         del self.domain         del self.domain
54     def test_StokesProblemCartesian_PCG_P_0(self):     def test_PCG_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="PCG")  
         
        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_PCG_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="PCG")  
         
        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_PCG_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="PCG")  
         
        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_GMRES_P_0(self):  
55         ETA=1.         ETA=1.
56         P1=0.         P1=0.
57    
# Line 138  class Test_Simple2DModels(unittest.TestC Line 67  class Test_Simple2DModels(unittest.TestC
67         sp.initialize(f=F,fixed_u_mask=mask,eta=ETA)         sp.initialize(f=F,fixed_u_mask=mask,eta=ETA)
68         u0=(1-x[0])*x[0]*[0.,1.]         u0=(1-x[0])*x[0]*[0.,1.]
69         p0=Scalar(P1,ReducedSolution(self.domain))         p0=Scalar(P1,ReducedSolution(self.domain))
70         u,p=sp.solve(u0,p0,show_details=VERBOSE, verbose=VERBOSE,max_iter=100,solver="GMRES",iter_restart=20)         sp.setTolerance(self.TOL)
71           u,p=sp.solve(u0,p0,show_details=DETAIL_VERBOSE, verbose=VERBOSE,max_iter=100,usePCG=True)
72                
73         error_v0=Lsup(u[0]-u0[0])         error_v0=Lsup(u[0]-u0[0])
74         error_v1=Lsup(u[1]-u0[1])/0.25         error_v1=Lsup(u[1]-u0[1])/0.25
75         zz=P1*x[0]*x[1]-p         error_p=Lsup(p+P1*x[0]*x[1])
76         error_p=Lsup(zz-integrate(zz))         self.failUnless(error_p<10*self.TOL, "pressure error too large.")
77           self.failUnless(error_v0<10*self.TOL, "0-velocity error too large.")
78        # print error_p, error_v0,error_v1         self.failUnless(error_v1<10*self.TOL, "1-velocity error too large.")
        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.")  
79    
80     def test_StokesProblemCartesian_GMRES_P_small(self):     def test_PCG_P_small(self):
81         ETA=1.         ETA=1.
82         P1=1.         P1=1.
83    
# Line 166  class Test_Simple2DModels(unittest.TestC Line 93  class Test_Simple2DModels(unittest.TestC
93         sp.initialize(f=F,fixed_u_mask=mask,eta=ETA)         sp.initialize(f=F,fixed_u_mask=mask,eta=ETA)
94         u0=(1-x[0])*x[0]*[0.,1.]         u0=(1-x[0])*x[0]*[0.,1.]
95         p0=Scalar(P1,ReducedSolution(self.domain))         p0=Scalar(P1,ReducedSolution(self.domain))
96         u,p=sp.solve(u0,p0,show_details=VERBOSE, verbose=VERBOSE,max_iter=100,solver="GMRES")         sp.setTolerance(self.TOL*0.2)
97                 u,p=sp.solve(u0,p0,show_details=DETAIL_VERBOSE, verbose=VERBOSE,max_iter=100,usePCG=True)
98         error_v0=Lsup(u[0]-u0[0])         error_v0=Lsup(u[0]-u0[0])
99         error_v1=Lsup(u[1]-u0[1])/0.25         error_v1=Lsup(u[1]-u0[1])/0.25
100         zz=P1*x[0]*x[1]-p         error_p=Lsup(P1*x[0]*x[1]+p)
101         error_p=Lsup(zz-integrate(zz))         saveVTK("d.vtu",p=p, e=P1*x[0]*x[1]+p, p_ref=P1*x[0]*x[1])
102         # print error_p, error_v0,error_v1         self.failUnless(error_v0<10*self.TOL, "0-velocity error too large.")
103         self.failUnless(error_p<TOL,"pressure error too large.")         self.failUnless(error_v1<10*self.TOL, "1-velocity error too large.")
104         self.failUnless(error_v0<TOL,"0-velocity error too large.")         self.failUnless(error_p<10*self.TOL, "pressure error too large.")
        self.failUnless(error_v1<TOL,"1-velocity error too large.")  
105    
106     def test_StokesProblemCartesian_GMRES_P_large(self):     def test_PCG_P_large(self):
107         ETA=1.         ETA=1.
108         P1=1000.         P1=1000.
109    
# Line 193  class Test_Simple2DModels(unittest.TestC Line 119  class Test_Simple2DModels(unittest.TestC
119         sp.initialize(f=F,fixed_u_mask=mask,eta=ETA)         sp.initialize(f=F,fixed_u_mask=mask,eta=ETA)
120         u0=(1-x[0])*x[0]*[0.,1.]         u0=(1-x[0])*x[0]*[0.,1.]
121         p0=Scalar(P1,ReducedSolution(self.domain))         p0=Scalar(P1,ReducedSolution(self.domain))
122         u,p=sp.solve(u0,p0,show_details=VERBOSE, verbose=VERBOSE,max_iter=100,solver="GMRES")         sp.setTolerance(self.TOL)
123           u,p=sp.solve(u0,p0,show_details=DETAIL_VERBOSE, verbose=VERBOSE,max_iter=100,usePCG=True)
124                
125         error_v0=Lsup(u[0]-u0[0])         error_v0=Lsup(u[0]-u0[0])
126         error_v1=Lsup(u[1]-u0[1])/0.25         error_v1=Lsup(u[1]-u0[1])/0.25
127         zz=P1*x[0]*x[1]-p         error_p=Lsup(P1*x[0]*x[1]+p)/P1
128         error_p=Lsup(zz-integrate(zz))/P1         self.failUnless(error_p<10*self.TOL, "pressure error too large.")
129         # print error_p, error_v0,error_v1         self.failUnless(error_v0<10*self.TOL, "0-velocity error too large.")
130         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.")  
131    
132     def test_StokesProblemCartesian_MINRES_P_0(self):     def test_GMRES_P_0(self):
133         ETA=1.         ETA=1.
134         P1=0.         P1=0.
135    
# Line 220  class Test_Simple2DModels(unittest.TestC Line 145  class Test_Simple2DModels(unittest.TestC
145         sp.initialize(f=F,fixed_u_mask=mask,eta=ETA)         sp.initialize(f=F,fixed_u_mask=mask,eta=ETA)
146         u0=(1-x[0])*x[0]*[0.,1.]         u0=(1-x[0])*x[0]*[0.,1.]
147         p0=Scalar(P1,ReducedSolution(self.domain))         p0=Scalar(P1,ReducedSolution(self.domain))
148         u,p=sp.solve(u0,p0,show_details=VERBOSE, verbose=VERBOSE,max_iter=100,solver="MINRES")         sp.setTolerance(self.TOL)
149                 u,p=sp.solve(u0,p0,show_details=DETAIL_VERBOSE, verbose=VERBOSE,max_iter=50,usePCG=False,iter_restart=18)
        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")  
150                
151         error_v0=Lsup(u[0]-u0[0])         error_v0=Lsup(u[0]-u0[0])
152         error_v1=Lsup(u[1]-u0[1])/0.25         error_v1=Lsup(u[1]-u0[1])/0.25
153         zz=P1*x[0]*x[1]-p         error_p=Lsup(P1*x[0]*x[1]+p)
154         error_p=Lsup(zz-integrate(zz))/P1         self.failUnless(error_p<10*self.TOL, "pressure error too large.")
155         # print error_p, error_v0,error_v1         self.failUnless(error_v0<10*self.TOL, "0-velocity error too large.")
156         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.")  
   
 #   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.")  
157    
158     def test_StokesProblemCartesian_TFQMR_P_small(self):     def test_GMRES_P_small(self):
159         ETA=1.         ETA=1.
160         P1=1.         P1=1.
161    
# Line 329  class Test_Simple2DModels(unittest.TestC Line 171  class Test_Simple2DModels(unittest.TestC
171         sp.initialize(f=F,fixed_u_mask=mask,eta=ETA)         sp.initialize(f=F,fixed_u_mask=mask,eta=ETA)
172         u0=(1-x[0])*x[0]*[0.,1.]         u0=(1-x[0])*x[0]*[0.,1.]
173         p0=Scalar(P1,ReducedSolution(self.domain))         p0=Scalar(P1,ReducedSolution(self.domain))
174         u,p=sp.solve(u0,p0,show_details=VERBOSE, verbose=VERBOSE,max_iter=100,solver="TFQMR")         sp.setTolerance(self.TOL*0.1)
175           u,p=sp.solve(u0,p0,show_details=DETAIL_VERBOSE, verbose=VERBOSE,max_iter=20,usePCG=False)
176                
177         error_v0=Lsup(u[0]-u0[0])         error_v0=Lsup(u[0]-u0[0])
178         error_v1=Lsup(u[1]-u0[1])/0.25         error_v1=Lsup(u[1]-u0[1])/0.25
179         zz=P1*x[0]*x[1]-p         error_p=Lsup(P1*x[0]*x[1]+p)
180         error_p=Lsup(zz-integrate(zz))/P1         self.failUnless(error_p<10*self.TOL, "pressure error too large.")
181         # print error_p, error_v0,error_v1         self.failUnless(error_v0<10*self.TOL, "0-velocity error too large.")
182         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.")  
183    
184     def test_StokesProblemCartesian_TFQMR_P_large(self):     def test_GMRES_P_large(self):
185         ETA=1.         ETA=1.
186         P1=1000.         P1=1000.
187    
# Line 356  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="TFQMR")         sp.setTolerance(self.TOL)
201           u,p=sp.solve(u0,p0,show_details=DETAIL_VERBOSE, verbose=VERBOSE,max_iter=100,usePCG=False)
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         error_p=Lsup(P1*x[0]*x[1]+p)/P1
206         error_p=Lsup(zz-integrate(zz))/P1         self.failUnless(error_p<10*self.TOL, "pressure error too large.")
207         # print error_p, error_v0,error_v1         self.failUnless(error_v0<10*self.TOL, "0-velocity error too large.")
208         self.failUnless(error_p<TOL,"pressure error too large.")         self.failUnless(error_v1<10*self.TOL, "1-velocity error too large.")
209         self.failUnless(error_v0<TOL,"0-velocity error too large.")  #====================================================================================================================
210         self.failUnless(error_v1<TOL,"1-velocity error too large.")  class Test_StokesProblemCartesian3D(unittest.TestCase):
   
   
   
 class Test_Simple3DModels(unittest.TestCase):  
211     def setUp(self):     def setUp(self):
212           NE=6
213           self.TOL=1e-4
214         self.domain=Brick(NE,NE,NE,order=2,useFullElementOrder=True)         self.domain=Brick(NE,NE,NE,order=2,useFullElementOrder=True)
215     def tearDown(self):     def tearDown(self):
216         del self.domain         del self.domain
217     def test_StokesProblemCartesian_PCG_P_0(self):     def test_PCG_P_0(self):
218         ETA=1.         ETA=1.
219         P1=0.         P1=0.
220    
# Line 383  class Test_Simple3DModels(unittest.TestC Line 223  class Test_Simple3DModels(unittest.TestC
223         x=self.domain.getX()         x=self.domain.getX()
224         mask=whereZero(x[0])    * [1.,1.,1.] \         mask=whereZero(x[0])    * [1.,1.,1.] \
225                +whereZero(x[0]-1)  * [1.,1.,1.] \                +whereZero(x[0]-1)  * [1.,1.,1.] \
226                +whereZero(x[1])    * [1.,1.,1.] \                +whereZero(x[1])    * [1.,0.,1.] \
227                +whereZero(x[1]-1)  * [1.,1.,1.] \                +whereZero(x[1]-1)  * [1.,1.,1.] \
228                +whereZero(x[2])    * [1.,1.,0.] \                +whereZero(x[2])    * [1.,1.,0.] \
229                +whereZero(x[2]-1)  * [1.,1.,1.]                +whereZero(x[2]-1)  * [1.,1.,1.]
# Line 394  class Test_Simple3DModels(unittest.TestC Line 234  class Test_Simple3DModels(unittest.TestC
234         sp.initialize(f=F,fixed_u_mask=mask,eta=ETA)         sp.initialize(f=F,fixed_u_mask=mask,eta=ETA)
235         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.]
236         p0=Scalar(P1,ReducedSolution(self.domain))         p0=Scalar(P1,ReducedSolution(self.domain))
237         u,p=sp.solve(u0,p0,show_details=VERBOSE, verbose=VERBOSE,max_iter=100,solver="PCG")         sp.setTolerance(self.TOL)
238           u,p=sp.solve(u0,p0,show_details=DETAIL_VERBOSE, verbose=VERBOSE ,max_iter=100,usePCG=True)
239                
240         error_v0=Lsup(u[0]-u0[0])         error_v0=Lsup(u[0]-u0[0])
241         error_v1=Lsup(u[1]-u0[1])         error_v1=Lsup(u[1]-u0[1])
242         error_v2=Lsup(u[2]-u0[2])/0.25**2         error_v2=Lsup(u[2]-u0[2])/0.25**2
243         zz=P1*x[0]*x[1]*x[2]-p         error_p=Lsup(P1*x[0]*x[1]*x[2]+p)
244         error_p=Lsup(zz-integrate(zz))         self.failUnless(error_p<10*self.TOL, "pressure error too large.")
245         # print error_p, error_v0,error_v1,error_v2         self.failUnless(error_v0<10*self.TOL, "0-velocity error too large.")
246         self.failUnless(error_p<TOL,"pressure error too large.")         self.failUnless(error_v1<10*self.TOL, "1-velocity error too large.")
247         self.failUnless(error_v0<TOL,"0-velocity error too large.")         self.failUnless(error_v2<10*self.TOL, "2-velocity error too large.")
        self.failUnless(error_v1<TOL,"1-velocity error too large.")  
        self.failUnless(error_v2<TOL,"2-velocity error too large.")  
248    
249     def test_StokesProblemCartesian_PCG_P_small(self):     def test_PCG_P_small(self):
250         ETA=1.         ETA=1.
251         P1=1.         P1=1.
252    
# Line 415  class Test_Simple3DModels(unittest.TestC Line 254  class Test_Simple3DModels(unittest.TestC
254         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.]         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.]
255         mask=whereZero(x[0])    * [1.,1.,1.] \         mask=whereZero(x[0])    * [1.,1.,1.] \
256                +whereZero(x[0]-1)  * [1.,1.,1.] \                +whereZero(x[0]-1)  * [1.,1.,1.] \
257                +whereZero(x[1])    * [1.,1.,1.] \                +whereZero(x[1])    * [1.,0.,1.] \
258                +whereZero(x[1]-1)  * [1.,1.,1.] \                +whereZero(x[1]-1)  * [1.,1.,1.] \
259                +whereZero(x[2])    * [1.,1.,0.] \                +whereZero(x[2])    * [1.,1.,0.] \
260                +whereZero(x[2]-1)  * [1.,1.,1.]                +whereZero(x[2]-1)  * [1.,1.,1.]
# Line 426  class Test_Simple3DModels(unittest.TestC Line 265  class Test_Simple3DModels(unittest.TestC
265         sp.initialize(f=F,fixed_u_mask=mask,eta=ETA)         sp.initialize(f=F,fixed_u_mask=mask,eta=ETA)
266         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.]
267         p0=Scalar(P1,ReducedSolution(self.domain))         p0=Scalar(P1,ReducedSolution(self.domain))
268         u,p=sp.solve(u0,p0,show_details=VERBOSE, verbose=VERBOSE,max_iter=100,solver="PCG")         sp.setTolerance(self.TOL*0.1)
269                 u,p=sp.solve(u0,p0,show_details=DETAIL_VERBOSE, verbose=VERBOSE ,max_iter=100,usePCG=True)
270         error_v0=Lsup(u[0]-u0[0])         error_v0=Lsup(u[0]-u0[0])
271         error_v1=Lsup(u[1]-u0[1])         error_v1=Lsup(u[1]-u0[1])
272         error_v2=Lsup(u[2]-u0[2])/0.25**2         error_v2=Lsup(u[2]-u0[2])/0.25**2
273         zz=P1*x[0]*x[1]*x[2]-p         error_p=Lsup(P1*x[0]*x[1]*x[2]+p)
274         error_p=Lsup(zz-integrate(zz))/P1         self.failUnless(error_p<10*self.TOL, "pressure error too large.")
275         # print error_p, error_v0,error_v1,error_v2         self.failUnless(error_v0<10*self.TOL, "0-velocity error too large.")
276         self.failUnless(error_p<TOL,"pressure error too large.")         self.failUnless(error_v1<10*self.TOL, "1-velocity error too large.")
277         self.failUnless(error_v0<TOL,"0-velocity error too large.")         self.failUnless(error_v2<10*self.TOL, "2-velocity error too large.")
278         self.failUnless(error_v1<TOL,"1-velocity error too large.")     def test_PCG_P_large(self):
        self.failUnless(error_v2<TOL,"2-velocity error too large.")  
    def test_StokesProblemCartesian_PCG_P_large(self):  
279         ETA=1.         ETA=1.
280         P1=1000.         P1=1000.
281    
# Line 446  class Test_Simple3DModels(unittest.TestC Line 283  class Test_Simple3DModels(unittest.TestC
283         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.]         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.]
284         mask=whereZero(x[0])    * [1.,1.,1.] \         mask=whereZero(x[0])    * [1.,1.,1.] \
285                +whereZero(x[0]-1)  * [1.,1.,1.] \                +whereZero(x[0]-1)  * [1.,1.,1.] \
286                +whereZero(x[1])    * [1.,1.,1.] \                +whereZero(x[1])    * [1.,0.,1.] \
287                +whereZero(x[1]-1)  * [1.,1.,1.] \                +whereZero(x[1]-1)  * [1.,1.,1.] \
288                +whereZero(x[2])    * [1.,1.,0.] \                +whereZero(x[2])    * [1.,1.,0.] \
289                +whereZero(x[2]-1)  * [1.,1.,1.]                +whereZero(x[2]-1)  * [1.,1.,1.]
# Line 457  class Test_Simple3DModels(unittest.TestC Line 294  class Test_Simple3DModels(unittest.TestC
294         sp.initialize(f=F,fixed_u_mask=mask,eta=ETA)         sp.initialize(f=F,fixed_u_mask=mask,eta=ETA)
295         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.]
296         p0=Scalar(P1,ReducedSolution(self.domain))         p0=Scalar(P1,ReducedSolution(self.domain))
297         u,p=sp.solve(u0,p0,show_details=VERBOSE, verbose=VERBOSE,max_iter=100,solver="PCG")         sp.setTolerance(self.TOL)
298           u,p=sp.solve(u0,-p0,show_details=DETAIL_VERBOSE, verbose=VERBOSE ,max_iter=100,usePCG=True)
299                
300         error_v0=Lsup(u[0]-u0[0])         error_v0=Lsup(u[0]-u0[0])
301         error_v1=Lsup(u[1]-u0[1])         error_v1=Lsup(u[1]-u0[1])
302         error_v2=Lsup(u[2]-u0[2])/0.25**2         error_v2=Lsup(u[2]-u0[2])/0.25**2
303         zz=P1*x[0]*x[1]*x[2]-p         error_p=Lsup(P1*x[0]*x[1]*x[2]+p)/P1
304         error_p=Lsup(zz-integrate(zz))/P1         self.failUnless(error_p<10*self.TOL, "pressure error too large.")
305         # print error_p, error_v0,error_v1,error_v2         self.failUnless(error_v0<10*self.TOL, "0-velocity error too large.")
306         self.failUnless(error_p<TOL,"pressure error too large.")         self.failUnless(error_v1<10*self.TOL, "1-velocity error too large.")
307         self.failUnless(error_v0<TOL,"0-velocity error too large.")         self.failUnless(error_v2<10*self.TOL, "2-velocity error too large.")
        self.failUnless(error_v1<TOL,"1-velocity error too large.")  
        self.failUnless(error_v2<TOL,"2-velocity error too large.")  
308    
309     def test_StokesProblemCartesian_GMRES_P_0(self):     def test_GMRES_P_0(self):
310         ETA=1.         ETA=1.
311         P1=0.         P1=0.
312    
# Line 490  class Test_Simple3DModels(unittest.TestC Line 326  class Test_Simple3DModels(unittest.TestC
326         sp.initialize(f=F,fixed_u_mask=mask,eta=ETA)         sp.initialize(f=F,fixed_u_mask=mask,eta=ETA)
327         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.]
328         p0=Scalar(P1,ReducedSolution(self.domain))         p0=Scalar(P1,ReducedSolution(self.domain))
329         u,p=sp.solve(u0,p0,show_details=VERBOSE, verbose=VERBOSE,max_iter=100,solver="GMRES",iter_restart=20)         sp.setTolerance(self.TOL)
330           u,p=sp.solve(u0,p0,show_details=DETAIL_VERBOSE, verbose=VERBOSE,max_iter=100,usePCG=False,iter_restart=20)
331                
332         error_v0=Lsup(u[0]-u0[0])         error_v0=Lsup(u[0]-u0[0])
333         error_v1=Lsup(u[1]-u0[1])         error_v1=Lsup(u[1]-u0[1])
334         error_v2=Lsup(u[2]-u0[2])/0.25**2         error_v2=Lsup(u[2]-u0[2])/0.25**2
335         zz=P1*x[0]*x[1]*x[2]-p         error_p=Lsup(P1*x[0]*x[1]*x[2]+p)
336         error_p=Lsup(zz-integrate(zz))         self.failUnless(error_p<10*self.TOL, "pressure error too large.")
337         # print error_p, error_v0,error_v1,error_v2         self.failUnless(error_v0<10*self.TOL, "0-velocity error too large.")
338         self.failUnless(error_p<TOL,"pressure error too large.")         self.failUnless(error_v1<10*self.TOL, "1-velocity error too large.")
339         self.failUnless(error_v0<TOL,"0-velocity error too large.")         self.failUnless(error_v2<10*self.TOL, "2-velocity error too large.")
340         self.failUnless(error_v1<TOL,"1-velocity error too large.")     def test_GMRES_P_small(self):
        self.failUnless(error_v2<TOL,"2-velocity error too large.")  
    def test_StokesProblemCartesian_GMRES_P_small(self):  
341         ETA=1.         ETA=1.
342         P1=1.         P1=1.
343    
# Line 521  class Test_Simple3DModels(unittest.TestC Line 356  class Test_Simple3DModels(unittest.TestC
356         sp.initialize(f=F,fixed_u_mask=mask,eta=ETA)         sp.initialize(f=F,fixed_u_mask=mask,eta=ETA)
357         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.]
358         p0=Scalar(P1,ReducedSolution(self.domain))         p0=Scalar(P1,ReducedSolution(self.domain))
359         u,p=sp.solve(u0,p0,show_details=VERBOSE, verbose=VERBOSE,max_iter=100,solver="GMRES")         sp.setTolerance(self.TOL*0.1)
360           u,p=sp.solve(u0,p0,show_details=DETAIL_VERBOSE and False, verbose=VERBOSE,max_iter=100,usePCG=False)
361                
362         error_v0=Lsup(u[0]-u0[0])         error_v0=Lsup(u[0]-u0[0])
363         error_v1=Lsup(u[1]-u0[1])         error_v1=Lsup(u[1]-u0[1])
364         error_v2=Lsup(u[2]-u0[2])/0.25**2         error_v2=Lsup(u[2]-u0[2])/0.25**2
365         zz=P1*x[0]*x[1]*x[2]-p         error_p=Lsup(P1*x[0]*x[1]*x[2]+p)/P1
366         error_p=Lsup(zz-integrate(zz))/P1         self.failUnless(error_v0<10*self.TOL, "0-velocity error too large.")
367         # print error_p, error_v0,error_v1,error_v2         self.failUnless(error_v1<10*self.TOL, "1-velocity error too large.")
368         self.failUnless(error_p<TOL,"pressure error too large.")         self.failUnless(error_v2<10*self.TOL, "2-velocity error too large.")
369         self.failUnless(error_v0<TOL,"0-velocity error too large.")         self.failUnless(error_p<10*self.TOL, "pressure error too large.")
370         self.failUnless(error_v1<TOL,"1-velocity error too large.")     def test_GMRES_P_large(self):
        self.failUnless(error_v2<TOL,"2-velocity error too large.")  
    def test_StokesProblemCartesian_GMRES_P_large(self):  
371         ETA=1.         ETA=1.
372         P1=1000.         P1=1000.
373    
# Line 541  class Test_Simple3DModels(unittest.TestC Line 375  class Test_Simple3DModels(unittest.TestC
375         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.]         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.]
376         mask=whereZero(x[0])    * [1.,1.,1.] \         mask=whereZero(x[0])    * [1.,1.,1.] \
377                +whereZero(x[0]-1)  * [1.,1.,1.] \                +whereZero(x[0]-1)  * [1.,1.,1.] \
378                +whereZero(x[1])    * [1.,1.,1.] \                +whereZero(x[1])    * [1.,0.,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="GMRES")  
         
        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.")  
   
    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.] \  
379                +whereZero(x[1]-1)  * [1.,1.,1.] \                +whereZero(x[1]-1)  * [1.,1.,1.] \
380                +whereZero(x[2])    * [1.,1.,0.] \                +whereZero(x[2])    * [1.,1.,0.] \
381                +whereZero(x[2]-1)  * [1.,1.,1.]                +whereZero(x[2]-1)  * [1.,1.,1.]
# Line 585  class Test_Simple3DModels(unittest.TestC Line 386  class Test_Simple3DModels(unittest.TestC
386         sp.initialize(f=F,fixed_u_mask=mask,eta=ETA)         sp.initialize(f=F,fixed_u_mask=mask,eta=ETA)
387         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.]
388         p0=Scalar(P1,ReducedSolution(self.domain))         p0=Scalar(P1,ReducedSolution(self.domain))
389         u,p=sp.solve(u0,p0,show_details=VERBOSE, verbose=VERBOSE,max_iter=100,solver="MINRES")         sp.setTolerance(self.TOL)
390           u,p=sp.solve(u0,p0,show_details=DETAIL_VERBOSE, verbose=VERBOSE ,max_iter=100,usePCG=False)
391                
392         error_v0=Lsup(u[0]-u0[0])         error_v0=Lsup(u[0]-u0[0])
393         error_v1=Lsup(u[1]-u0[1])         error_v1=Lsup(u[1]-u0[1])
394         error_v2=Lsup(u[2]-u0[2])/0.25**2         error_v2=Lsup(u[2]-u0[2])/0.25**2
395         zz=P1*x[0]*x[1]*x[2]-p         error_p=Lsup(P1*x[0]*x[1]*x[2]+p)/P1
396         error_p=Lsup(zz-integrate(zz))         self.failUnless(error_p<10*self.TOL, "pressure error too large.")
397         # print error_p, error_v0,error_v1,error_v2         self.failUnless(error_v0<10*self.TOL, "0-velocity error too large.")
398         self.failUnless(error_p<TOL,"pressure error too large.")         self.failUnless(error_v1<10*self.TOL, "1-velocity error too large.")
399         self.failUnless(error_v0<TOL,"0-velocity error too large.")         self.failUnless(error_v2<10*self.TOL, "2-velocity error too large.")
400         self.failUnless(error_v1<TOL,"1-velocity error too large.")  #====================================================================================================================
401         self.failUnless(error_v2<TOL,"2-velocity error too large.")  class Test_Darcy(unittest.TestCase):
402        # this is a simple test for the darcy flux problem
403     def test_StokesProblemCartesian_MINRES_P_small(self):      #
404         ETA=1.      #
405         P1=1.      #  p = 1/k * ( 1/2* (fb-f0)/xb* x **2 + f0 * x - ub*x ) +  p0
406        #
407        #  with f = (fb-f0)/xb* x + f0
408        #
409        #    u = f - k * p,x = ub
410        #
411        #  we prescribe pressure at x=x0=0 to p0
412        #
413        #  if we prescribe the pressure on the bottom x=xb we set
414        #
415        #  pb= 1/k * ( 1/2* (fb-f0)* xb + f0 * xb - ub*xb ) +  p0 = 1/k * ((fb+f0)/2 - ub ) * xb  + p0
416        #
417        #  which leads to ub = (fb+f0)/2-k*(pb-p0)/xb
418        #
419        def rescaleDomain(self):
420            x=self.dom.getX().copy()
421            for i in xrange(self.dom.getDim()):
422                 x_inf=inf(x[i])
423                 x_sup=sup(x[i])
424                 if i == self.dom.getDim()-1:
425                    x[i]=-self.WIDTH*(x[i]-x_sup)/(x_inf-x_sup)
426                 else:
427                    x[i]=self.WIDTH*(x[i]-x_inf)/(x_sup-x_inf)
428            self.dom.setX(x)
429        def getScalarMask(self,include_bottom=True):
430            x=self.dom.getX().copy()
431            x_inf=inf(x[self.dom.getDim()-1])
432            x_sup=sup(x[self.dom.getDim()-1])
433            out=whereZero(x[self.dom.getDim()-1]-x_sup)
434            if include_bottom: out+=whereZero(x[self.dom.getDim()-1]-x_inf)
435            return wherePositive(out)
436        def getVectorMask(self,include_bottom=True):
437            x=self.dom.getX().copy()
438            out=Vector(0.,Solution(self.dom))
439            for i in xrange(self.dom.getDim()):
440                 x_inf=inf(x[i])
441                 x_sup=sup(x[i])
442                 if i != self.dom.getDim()-1: out[i]+=whereZero(x[i]-x_sup)
443                 if i != self.dom.getDim()-1 or include_bottom: out[i]+=whereZero(x[i]-x_inf)
444            return wherePositive(out)
445    
446        def setSolutionFixedBottom(self, p0, pb, f0, fb, k):
447             d=self.dom.getDim()
448             x=self.dom.getX()[d-1]
449             xb=inf(x)
450             u=Vector(0.,Solution(self.dom))+kronecker(d)[d-1]*((f0+fb)/2.-k*(pb-p0)/xb)
451             p=1./k*((fb-f0)/(xb*2.)* x**2 - (fb-f0)/2.*x)+(pb-p0)/xb*x +  p0
452             f= ((fb-f0)/xb* x + f0)*kronecker(Function(self.dom))[d-1]
453             return u,p,f
454            
455        def testConstF_FixedBottom_smallK(self):
456            k=1.e-10
457            mp=self.getScalarMask(include_bottom=True)
458            mv=self.getVectorMask(include_bottom=False)
459            u_ref,p_ref,f=self.setSolutionFixedBottom(p0=2500,pb=4000.,f0=10.,fb=10.,k=k)
460            p=p_ref*mp
461            u=u_ref*mv
462            df=DarcyFlow(self.dom)
463            df.setValue(g=f,
464                          location_of_fixed_pressure=mp,
465                          location_of_fixed_flux=mv,
466                          permeability=Scalar(k,Function(self.dom)))
467            df.setTolerance(rtol=self.TOL)
468            df.setSubProblemTolerance()
469            v,p=df.solve(u_ref,p, max_iter=100, verbose=VERBOSE)
470            self.failUnless(Lsup(v-u_ref)<self.TOL*10.*Lsup(u_ref), "flux error too big.")
471            self.failUnless(Lsup(p-p_ref)<self.TOL*10.*Lsup(p_ref), "pressure error too big.")
472        def testConstF_FixedBottom_mediumK(self):
473            k=1.
474            mp=self.getScalarMask(include_bottom=True)
475            mv=self.getVectorMask(include_bottom=False)
476            u_ref,p_ref,f=self.setSolutionFixedBottom(p0=2500,pb=4000.,f0=10.,fb=10.,k=k)
477            p=p_ref*mp
478            u=u_ref*mv
479            df=DarcyFlow(self.dom)
480            df.setValue(g=f,
481                          location_of_fixed_pressure=mp,
482                          location_of_fixed_flux=mv,
483                          permeability=Scalar(k,Function(self.dom)))
484            df.setTolerance(rtol=self.TOL)
485            df.setSubProblemTolerance()
486            v,p=df.solve(u,p,max_iter=100, verbose=VERBOSE )
487            self.failUnless(Lsup(p-p_ref)<self.TOL*10.*Lsup(p_ref), "pressure error too big.")
488            self.failUnless(Lsup(v-u_ref)<self.TOL*10.*Lsup(u_ref), "flux error too big.")
489    
490        def testConstF_FixedBottom_largeK(self):
491            k=1.e10
492            mp=self.getScalarMask(include_bottom=True)
493            mv=self.getVectorMask(include_bottom=False)
494            u_ref,p_ref,f=self.setSolutionFixedBottom(p0=2500,pb=4000.,f0=10.,fb=10.,k=k)
495            p=p_ref*mp
496            u=u_ref*mv
497            df=DarcyFlow(self.dom)
498            df.setValue(g=f,
499                          location_of_fixed_pressure=mp,
500                          location_of_fixed_flux=mv,
501                          permeability=Scalar(k,Function(self.dom)))
502            df.setTolerance(rtol=self.TOL)
503            df.setSubProblemTolerance()
504            v,p=df.solve(u,p, max_iter=100, verbose=VERBOSE)
505            self.failUnless(Lsup(v-u_ref)<self.TOL*10.*Lsup(u_ref), "flux error too big.")
506            self.failUnless(Lsup(p-p_ref)<self.TOL*10.*Lsup(p_ref), "pressure error too big.")
507    
508        def testVarioF_FixedBottom_smallK(self):
509            k=1.e-10
510            mp=self.getScalarMask(include_bottom=True)
511            mv=self.getVectorMask(include_bottom=False)
512            u_ref,p_ref,f=self.setSolutionFixedBottom(p0=2500,pb=4000.,f0=10.,fb=30.,k=k)
513            p=p_ref*mp
514            u=u_ref*mv
515            df=DarcyFlow(self.dom)
516            df.setValue(g=f,
517                          location_of_fixed_pressure=mp,
518                          location_of_fixed_flux=mv,
519                          permeability=Scalar(k,Function(self.dom)))
520            df.setTolerance(rtol=self.TOL)
521            df.setSubProblemTolerance()
522            v,p=df.solve(u,p, max_iter=100, verbose=VERBOSE)
523            self.failUnless(Lsup(v-u_ref)<self.TOL*10.*Lsup(u_ref), "flux error too big.")
524            self.failUnless(Lsup(p-p_ref)<self.TOL*10.*Lsup(p_ref), "pressure error too big.")
525    
526        def testVarioF_FixedBottom_mediumK(self):
527            k=1.
528            mp=self.getScalarMask(include_bottom=True)
529            mv=self.getVectorMask(include_bottom=False)
530            u_ref,p_ref,f=self.setSolutionFixedBottom(p0=2500,pb=4000.,f0=10.,fb=30.,k=k)
531            p=p_ref*mp
532            u=u_ref*mv
533            df=DarcyFlow(self.dom)
534            df.setValue(g=f,
535                          location_of_fixed_pressure=mp,
536                          location_of_fixed_flux=mv,
537                          permeability=Scalar(k,Function(self.dom)))
538            df.setTolerance(rtol=self.TOL)
539            df.setSubProblemTolerance()
540            v,p=df.solve(u,p, max_iter=100, verbose=VERBOSE)
541            self.failUnless(Lsup(v-u_ref)<self.TOL*10.*Lsup(u_ref), "flux error too big.")
542            self.failUnless(Lsup(p-p_ref)<self.TOL*10.*Lsup(p_ref), "pressure error too big.")
543    
544        def testVarioF_FixedBottom_largeK(self):
545            k=1.e10
546            mp=self.getScalarMask(include_bottom=True)
547            mv=self.getVectorMask(include_bottom=False)
548            u_ref,p_ref,f=self.setSolutionFixedBottom(p0=2500,pb=4000.,f0=10.,fb=30.,k=k)
549            p=p_ref*mp
550            u=u_ref*mv
551            df=DarcyFlow(self.dom)
552            df.setValue(g=f,
553                          location_of_fixed_pressure=mp,
554                          location_of_fixed_flux=mv,
555                          permeability=Scalar(k,Function(self.dom)))
556            df.setTolerance(rtol=self.TOL)
557            df.setSubProblemTolerance()
558            v,p=df.solve(u,p, max_iter=100, verbose=VERBOSE)
559            self.failUnless(Lsup(v-u_ref)<self.TOL*10.*Lsup(u_ref), "flux error too big.")
560            self.failUnless(Lsup(p-p_ref)<self.TOL*10.*Lsup(p_ref), "pressure error too big.")
561    
562        def testConstF_FreeBottom_smallK(self):
563            k=1.e-10
564            mp=self.getScalarMask(include_bottom=False)
565            mv=self.getVectorMask(include_bottom=True)
566            u_ref,p_ref,f=self.setSolutionFixedBottom(p0=2500,pb=4000.,f0=10.,fb=10.,k=k)
567            p=p_ref*mp
568            u=u_ref*mv
569            df=DarcyFlow(self.dom)
570            df.setValue(g=f,
571                        location_of_fixed_pressure=mp,
572                          location_of_fixed_flux=mv,
573                          permeability=Scalar(k,Function(self.dom)))
574            df.setTolerance(rtol=self.TOL)
575            df.setSubProblemTolerance()
576            v,p=df.solve(u,p, max_iter=100, verbose=VERBOSE)
577            self.failUnless(Lsup(v-u_ref)<self.TOL*10.*Lsup(u_ref), "flux error too big.")
578            self.failUnless(Lsup(p-p_ref)<self.TOL*10.*Lsup(p_ref), "pressure error too big.")
579    
580        def testConstF_FreeBottom_mediumK(self):
581            k=1.
582            mp=self.getScalarMask(include_bottom=False)
583            mv=self.getVectorMask(include_bottom=True)
584            u_ref,p_ref,f=self.setSolutionFixedBottom(p0=2500,pb=4000.,f0=10.,fb=10.,k=k)
585            p=p_ref*mp
586            u=u_ref*mv
587            df=DarcyFlow(self.dom)
588            df.setValue(g=f,
589                          location_of_fixed_pressure=mp,
590                          location_of_fixed_flux=mv,
591                          permeability=Scalar(k,Function(self.dom)))
592            df.setTolerance(rtol=self.TOL)
593            df.setSubProblemTolerance()
594            v,p=df.solve(u,p, max_iter=100, verbose=VERBOSE)
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            df.setTolerance(rtol=self.TOL)
611            df.setSubProblemTolerance()
612            v,p=df.solve(u,p, max_iter=100, verbose=VERBOSE)
613            self.failUnless(Lsup(v-u_ref)<self.TOL*10.*Lsup(u_ref), "flux error too big.")
614            self.failUnless(Lsup(p-p_ref)<self.TOL*10.*Lsup(p_ref), "pressure error too big.")
615    
616        def testVarioF_FreeBottom_smallK(self):
617            k=1.e-10
618            mp=self.getScalarMask(include_bottom=False)
619            mv=self.getVectorMask(include_bottom=True)
620            u_ref,p_ref,f=self.setSolutionFixedBottom(p0=2500,pb=4000.,f0=10.,fb=30.,k=k)
621            p=p_ref*mp
622            u=u_ref*mv
623            df=DarcyFlow(self.dom)
624            df.setValue(g=f,
625                          location_of_fixed_pressure=mp,
626                          location_of_fixed_flux=mv,
627                          permeability=Scalar(k,Function(self.dom)))
628            df.setTolerance(rtol=self.TOL)
629            df.setSubProblemTolerance()
630            v,p=df.solve(u,p, max_iter=100, verbose=VERBOSE)
631            self.failUnless(Lsup(v-u_ref)<self.TOL*25.*Lsup(u_ref), "flux error too big.")  # 25 because of disc. error.
632            self.failUnless(Lsup(p-p_ref)<self.TOL*25.*Lsup(p_ref), "pressure error too big.")
633    
634        def testVarioF_FreeBottom_mediumK(self):
635            k=1.
636            mp=self.getScalarMask(include_bottom=False)
637            mv=self.getVectorMask(include_bottom=True)
638            u_ref,p_ref,f=self.setSolutionFixedBottom(p0=2500,pb=4000.,f0=10.,fb=30.,k=k)
639            p=p_ref*mp
640            u=u_ref*mv
641            df=DarcyFlow(self.dom)
642            df.setValue(g=f,
643                          location_of_fixed_pressure=mp,
644                          location_of_fixed_flux=mv,
645                          permeability=Scalar(k,Function(self.dom)))
646            df.setTolerance(rtol=self.TOL)
647            df.setSubProblemTolerance()
648            v,p=df.solve(u,p, max_iter=100, verbose=VERBOSE)
649            self.failUnless(Lsup(v-u_ref)<self.TOL*10.*Lsup(u_ref), "flux error too big.")
650            self.failUnless(Lsup(p-p_ref)<self.TOL*10.*Lsup(p_ref), "pressure error too big.")
651    
652        def testVarioF_FreeBottom_largeK(self):
653            k=1.e10
654            mp=self.getScalarMask(include_bottom=False)
655            mv=self.getVectorMask(include_bottom=True)
656            u_ref,p_ref,f=self.setSolutionFixedBottom(p0=2500,pb=4000.,f0=10.,fb=30.,k=k)
657            p=p_ref*mp
658            u=u_ref*mv
659            df=DarcyFlow(self.dom)
660            df.setValue(g=f,
661                          location_of_fixed_pressure=mp,
662                          location_of_fixed_flux=mv,
663                          permeability=Scalar(k,Function(self.dom)))
664            df.setTolerance(rtol=self.TOL)
665            df.setSubProblemTolerance()
666            v,p=df.solve(u,p, max_iter=100, verbose=VERBOSE)
667            self.failUnless(Lsup(v-u_ref)<self.TOL*10.*Lsup(u_ref), "flux error too big.")
668            self.failUnless(Lsup(p-p_ref)<self.TOL*10.*Lsup(p_ref), "pressure error too big.")
669    
670    class Test_Darcy2D(Test_Darcy):
671        TOL=1e-4
672        WIDTH=1.
673        def setUp(self):
674            NE=40  # wrning smaller NE may case a failure for VarioF tests due to discretization errors.
675            self.dom = Rectangle(NE,NE)
676            self.rescaleDomain()
677        def tearDown(self):
678            del self.dom
679    class Test_Darcy3D(Test_Darcy):
680        TOL=1e-4
681        WIDTH=1.
682        def setUp(self):
683            NE=25  # wrning smaller NE may case a failure for VarioF tests due to discretization errors.
684            self.dom = Brick(NE,NE,NE)
685            self.rescaleDomain()
686        def tearDown(self):
687            del self.dom
688    
        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  
        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_large(self):  
 #       ETA=1.  
 #       P1=1000.  
   
 #       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  
 #       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_TFQMR_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="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))  
        # 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.")  
689    
690     def test_StokesProblemCartesian_TFQMR_P_small(self):  class Test_Mountains3D(unittest.TestCase):
691         ETA=1.     def setUp(self):
692         P1=1.         NE=16
693           self.TOL=1e-4
694           self.domain=Brick(NE,NE,NE,order=1,useFullElementOrder=True)
695       def tearDown(self):
696           del self.domain
697       def test_periodic(self):
698           EPS=0.01
699    
700         x=self.domain.getX()         x=self.domain.getX()
701         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.]         v = Vector(0.0, Solution(self.domain))
702         mask=whereZero(x[0])    * [1.,1.,1.] \         a0=1
703                +whereZero(x[0]-1)  * [1.,1.,1.] \         a1=1
704                +whereZero(x[1])    * [1.,1.,1.] \         n0=2
705                +whereZero(x[1]-1)  * [1.,1.,1.] \         n1=2
706                +whereZero(x[2])    * [1.,1.,0.] \         n2=0.5
707                +whereZero(x[2]-1)  * [1.,1.,1.]         a2=-(a0*n0+a1*n1)/n2
708           v[0]=a0*sin(pi*n0*x[0])* cos(pi*n1*x[1])* cos(pi*n2*x[2])
709           v[1]=a1*cos(pi*n0*x[0])* sin(pi*n1*x[1])* cos(pi*n2*x[2])
710           v[2]=a2*cos(pi*n0*x[0])* cos(pi*n1*x[1])* sin(pi*n2*x[2])
711    
712           H_t=Scalar(0.0, Solution(self.domain))
713           mts=Mountains(self.domain,v,eps=EPS,z=1)
714           u,Z=mts.update(u=v,H_t=H_t)
715                
716                 error_int=integrate(Z)
717         sp=StokesProblemCartesian(self.domain)         self.failUnless(error_int<self.TOL, "Boundary intergral is too large.")
         
        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.")  
718    
719     def test_StokesProblemCartesian_TFQMR_P_large(self):  class Test_Mountains2D(unittest.TestCase):
720         ETA=1.     def setUp(self):
721         P1=1000.         NE=16
722           self.TOL=1e-4
723           self.domain=Rectangle(NE,NE,order=1,useFullElementOrder=True)
724       def tearDown(self):
725           del self.domain
726       def test_periodic(self):
727           EPS=0.01
728    
729         x=self.domain.getX()         x=self.domain.getX()
730         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.]         v = Vector(0.0, Solution(self.domain))
731         mask=whereZero(x[0])    * [1.,1.,1.] \         a0=1
732                +whereZero(x[0]-1)  * [1.,1.,1.] \         n0=1
733                +whereZero(x[1])    * [1.,1.,1.] \         n1=0.5
734                +whereZero(x[1]-1)  * [1.,1.,1.] \         a1=-(a0*n0)/n1
735                +whereZero(x[2])    * [1.,1.,0.] \         v[0]=a0*sin(pi*n0*x[0])* cos(pi*n1*x[1])
736                +whereZero(x[2]-1)  * [1.,1.,1.]         v[1]=a1*cos(pi*n0*x[0])* sin(pi*n1*x[1])
737                
738                 H_t=Scalar(0.0, Solution(self.domain))
739         sp=StokesProblemCartesian(self.domain)         mts=Mountains(self.domain,v,eps=EPS,z=1)
740                 u,Z=mts.update(u=v,H_t=H_t)
741         sp.initialize(f=F,fixed_u_mask=mask,eta=ETA)        
742         u0=(1-x[0])*x[0]*(1-x[1])*x[1]*[0.,0.,1.]         error_int=integrate(Z)
743         p0=Scalar(P1,ReducedSolution(self.domain))         self.failUnless(error_int<self.TOL, "Boundary intergral is too large.")
744         u,p=sp.solve(u0,p0,show_details=VERBOSE, verbose=VERBOSE,max_iter=100,solver="TFQMR")        
745          
746         error_v0=Lsup(u[0]-u0[0])  
747         error_v1=Lsup(u[1]-u0[1])  class Test_Rheologies(unittest.TestCase):
748         error_v2=Lsup(u[2]-u0[2])/0.25**2  
749         zz=P1*x[0]*x[1]*x[2]-p       """
750         error_p=Lsup(zz-integrate(zz))/P1       this is the program used to generate the powerlaw tests:
751         # print error_p, error_v0,error_v1,error_v2  
752         self.failUnless(error_p<TOL,"pressure error too large.")       TAU_Y=100.
753         self.failUnless(error_v0<TOL,"0-velocity error too large.")       N=10
754         self.failUnless(error_v1<TOL,"1-velocity error too large.")       M=5
755         self.failUnless(error_v2<TOL,"2-velocity error too large.")  
756         def getE(tau):
757             if tau<=TAU_Y:
758               return 1./(0.5+20*sqrt(tau))
759             else:
760               raise ValueError,"out of range."
761         tau=[ i* (TAU_Y/N) for i in range(N+1)] + [TAU_Y for i in range(M)]
762         e=[ tau[i]/getE(tau[i]) for i in range(N+1)]
763         e+= [ (i+1+j)* (max(e)/N) for j in range(M) ]
764    
765         print tau
766         print e
767         """
768         TOL=1.e-8
769         def test_PowerLaw_Init(self):
770             pl=PowerLaw(numMaterials=3,verbose=VERBOSE)
771    
772             self.failUnlessEqual(pl.getNumMaterials(),3,"num materials is wrong")
773             self.failUnlessEqual(pl.validMaterialId(0),True,"material id 0 not found")
774             self.failUnlessEqual(pl.validMaterialId(1),True,"material id 1 not found")
775             self.failUnlessEqual(pl.validMaterialId(2),True,"material id 2 not found")
776             self.failUnlessEqual(pl.validMaterialId(3),False,"material id 3 not found")
777    
778             self.failUnlessEqual(pl.getFriction(),None,"initial friction wrong.")
779             self.failUnlessEqual(pl.getTauY(),None,"initial tau y wrong.")
780             pl.setDruckerPragerLaw(tau_Y=10,friction=3)
781             self.failUnlessEqual(pl.getFriction(),3,"friction wrong.")
782             self.failUnlessEqual(pl.getTauY(),10,"tau y wrong.")
783    
784             self.failUnlessEqual(pl.getElasticShearModulus(),None,"initial shear modulus is wrong")
785             pl.setElasticShearModulus(1000)
786             self.failUnlessEqual(pl.getElasticShearModulus(),1000.,"shear modulus is wrong")
787    
788             e=pl.getEtaN()
789             self.failUnlessEqual(len(e),3,"initial length of etaN is wrong.")
790             self.failUnlessEqual(e,[None, None, None],"initial etaN are wrong.")
791             self.failUnlessEqual(pl.getEtaN(0),None,"initial material id 0 is wrong")
792             self.failUnlessEqual(pl.getEtaN(1),None,"initial material id 1 is wrong")
793             self.failUnlessEqual(pl.getEtaN(2),None,"initial material id 2 is wrong")
794             self.failUnlessRaises(ValueError, pl.getEtaN, 3)
795    
796             self.failUnlessRaises(ValueError,pl.setPowerLaws,eta_N=[10,20,30,40],tau_t=[2], power=[5,6,7])
797             self.failUnlessRaises(ValueError,pl.setPowerLaws,eta_N=[10,20,30],tau_t=[2,4,5,5], power=[5,6,7])
798             self.failUnlessRaises(ValueError,pl.setPowerLaws,eta_N=[10,20,30],tau_t=[2,1,2], power=[5,6,7,7])
799    
800             pl.setPowerLaws(eta_N=[10,20,30],tau_t=[100,200,300], power=[1,2,3])
801             self.failUnlessEqual(pl.getPower(),[1,2,3],"powers are wrong.")
802             self.failUnlessEqual(pl.getTauT(),[100,200,300],"tau t are wrong.")
803             self.failUnlessEqual(pl.getEtaN(),[10,20,30],"etaN are wrong.")
804    
805             pl.setPowerLaw(id=0,eta_N=40,tau_t=400, power=4)
806             self.failUnlessEqual(pl.getPower(),[4,2,3],"powers are wrong.")
807             self.failUnlessEqual(pl.getTauT(),[400,200,300],"tau t are wrong.")
808             self.failUnlessEqual(pl.getEtaN(),[40,20,30],"etaN are wrong.")
809    
810             self.failUnlessRaises(ValueError,pl.getPower,-1)
811             self.failUnlessEqual(pl.getPower(0),4,"power 0 is wrong.")
812             self.failUnlessEqual(pl.getPower(1),2,"power 1 is wrong.")
813             self.failUnlessEqual(pl.getPower(2),3,"power 2 is wrong.")
814             self.failUnlessRaises(ValueError,pl.getPower,3)
815    
816             self.failUnlessRaises(ValueError,pl.getTauT,-1)
817             self.failUnlessEqual(pl.getTauT(0),400,"tau t 0 is wrong.")
818             self.failUnlessEqual(pl.getTauT(1),200,"tau t 1 is wrong.")
819             self.failUnlessEqual(pl.getTauT(2),300,"tau t 2 is wrong.")
820             self.failUnlessRaises(ValueError,pl.getTauT,3)
821    
822             self.failUnlessRaises(ValueError,pl.getEtaN,-1)
823             self.failUnlessEqual(pl.getEtaN(0),40,"eta n 0 is wrong.")
824             self.failUnlessEqual(pl.getEtaN(1),20,"eta n 1 is wrong.")
825             self.failUnlessEqual(pl.getEtaN(2),30,"eta n 2 is wrong.")
826             self.failUnlessRaises(ValueError,pl.getEtaN,3)
827    
828         def checkResult(self,id,gamma_dot_, eta, tau_ref):
829             self.failUnless(eta>=0,"eta needs to be positive (test %s)"%id)
830             error=abs(gamma_dot_*eta-tau_ref)
831             self.failUnless(error<=self.TOL*tau_ref,"eta is wrong: error = gamma_dot_*eta-tau_ref = %s * %s - %s = %s (test %s)"%(gamma_dot_,eta,tau_ref,error,id))
832            
833         def test_PowerLaw_Linear(self):
834             taus= [0.0, 10.0, 20.0, 30.0, 40.0, 50.0, 60.0, 70.0, 80.0, 90.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0]
835             gamma_dot_s=[0.0, 5.0, 10.0, 15.0, 20.0, 25.0, 30.0, 35.0, 40.0, 45.0, 50.0, 55.0, 60.0, 65.0, 70.0, 75.0]
836             pl=PowerLaw(numMaterials=1,verbose=VERBOSE)
837             pl.setDruckerPragerLaw(tau_Y=100.)
838             pl.setPowerLaw(eta_N=2.)
839             pl.setEtaTolerance(self.TOL)
840             for i in xrange(len(taus)): self.checkResult(i,gamma_dot_s[i], pl.getEtaEff(gamma_dot_s[i]),taus[i])
841            
842         def test_PowerLaw_QuadLarge(self):
843             taus=[0.0, 10.0, 20.0, 30.0, 40.0, 50.0, 60.0, 70.0, 80.0, 90.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0]
844             gamma_dot_s=[0.0, 637.45553203367592, 1798.8543819998317, 3301.3353450309969, 5079.6442562694074, 7096.067811865476, 9325.1600308977995, 11748.240371477059, 14350.835055998656, 17121.29936490925, 20050.0, 22055.0, 24060.0, 26065.0, 28070.0, 30075.0]
845             pl=PowerLaw(numMaterials=2,verbose=VERBOSE)
846             pl.setDruckerPragerLaw(tau_Y=100.)
847             pl.setPowerLaws(eta_N=[2.,0.01],tau_t=[1, 25.], power=[1,2])
848             pl.setEtaTolerance(self.TOL)
849             for i in xrange(len(taus)): self.checkResult(i,gamma_dot_s[i], pl.getEtaEff(gamma_dot_s[i]),taus[i])
850    
851         def test_PowerLaw_QuadSmall(self):
852             taus=[0.0, 10.0, 20.0, 30.0, 40.0, 50.0, 60.0, 70.0, 80.0, 90.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0]
853             gamma_dot_s=[0.0, 5.632455532033676, 11.788854381999831, 18.286335345030995, 25.059644256269408, 32.071067811865476, 39.295160030897804, 46.713240371477056, 54.310835055998652, 62.076299364909254, 70.0, 77.0, 84.0, 91.0, 98.0, 105.0]
854             pl=PowerLaw(numMaterials=2,verbose=VERBOSE)
855             pl.setDruckerPragerLaw(tau_Y=100.)
856             pl.setPowerLaws(eta_N=[2.,10.],tau_t=[1, 25.], power=[1,2])
857             pl.setEtaTolerance(self.TOL)
858             for i in xrange(len(taus)): self.checkResult(i,gamma_dot_s[i], pl.getEtaEff(gamma_dot_s[i]),taus[i])
859    
860         def test_PowerLaw_CubeLarge(self):
861             taus=[0.0, 10.0, 20.0, 30.0, 40.0, 50.0, 60.0, 70.0, 80.0, 90.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0]
862             gamma_dot_s=[0.0, 51.415888336127786, 157.36125994561547, 304.6468153816889, 487.84283811405851, 703.60440414872664, 949.57131887226353, 1223.9494765692673, 1525.3084267560891, 1852.4689652218574, 2204.4346900318833, 2424.8781590350718, 2645.3216280382603, 2865.7650970414484, 3086.2085660446369, 3306.6520350478249]
863             pl=PowerLaw(numMaterials=2,verbose=VERBOSE)
864             pl.setDruckerPragerLaw(tau_Y=100.)
865             pl.setPowerLaws(eta_N=[2.,1./16.],tau_t=[1, 64.], power=[1,3])
866             pl.setEtaTolerance(self.TOL)
867             for i in xrange(len(taus)): self.checkResult(i,gamma_dot_s[i], pl.getEtaEff(gamma_dot_s[i]),taus[i])
868    
869         def test_PowerLaw_CubeSmall(self):
870             taus=[0.0, 10.0, 20.0, 30.0, 40.0, 50.0, 60.0, 70.0, 80.0, 90.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0]
871             gamma_dot_s=[0.0, 5.4641588833612778, 11.473612599456157, 17.89646815381689, 24.678428381140588, 31.786044041487269, 39.195713188722635, 46.889494765692675, 54.853084267560895, 63.074689652218574, 71.544346900318828, 78.698781590350706, 85.853216280382583, 93.007650970414474, 100.16208566044635, 107.316520350478]
872             pl=PowerLaw(numMaterials=2,verbose=VERBOSE)
873             pl.setDruckerPragerLaw(tau_Y=100.)
874             pl.setPowerLaws(eta_N=[2.,25./4.],tau_t=[1, 64.], power=[1,3])
875             pl.setEtaTolerance(self.TOL)
876             for i in xrange(len(taus)): self.checkResult(i,gamma_dot_s[i], pl.getEtaEff(gamma_dot_s[i]),taus[i])
877    
878         def test_PowerLaw_QuadLarge_CubeLarge(self):
879             taus=[0.0, 10.0, 20.0, 30.0, 40.0, 50.0, 60.0, 70.0, 80.0, 90.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0]
880             gamma_dot_s=[0.0, 683.87142036980367, 1946.2156419454475, 3590.982160412686, 5547.4870943834667, 7774.6722160142008, 10244.731349770063, 12937.189848046326, 15836.143482754744, 18928.768330131104, 22204.434690031885, 24424.878159035074, 26645.321628038262, 28865.765097041451, 31086.208566044639, 33306.652035047824]
881             pl=PowerLaw(numMaterials=3,verbose=VERBOSE)
882             pl.setDruckerPragerLaw(tau_Y=100.)
883             pl.setPowerLaws(eta_N=[2.,0.01,1./16.],tau_t=[1, 25.,64.], power=[1,2,3])
884             pl.setEtaTolerance(self.TOL)
885             for i in xrange(len(taus)): self.checkResult(i,gamma_dot_s[i], pl.getEtaEff(gamma_dot_s[i]),taus[i])
886    
887         def test_PowerLaw_QuadLarge_CubeSmall(self):
888             taus=[0.0, 10.0, 20.0, 30.0, 40.0, 50.0, 60.0, 70.0, 80.0, 90.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0]
889             gamma_dot_s=[0.0, 637.9196909170372, 1800.3279945992881, 3304.2318131848137, 5084.3226846505486, 7102.853855906963, 9334.3557440865225, 11760.129866242751, 14365.688140266215, 17139.374054561471, 20071.544346900318, 22078.698781590349, 24085.853216280382, 26093.007650970416, 28100.162085660446, 30107.316520350476]
890             pl=PowerLaw(numMaterials=3,verbose=VERBOSE)
891             pl.setDruckerPragerLaw(tau_Y=100.)
892             pl.setPowerLaws(eta_N=[2.,0.01,25./4.],tau_t=[1, 25.,64.], power=[1,2,3])
893             pl.setEtaTolerance(self.TOL)
894             for i in xrange(len(taus)): self.checkResult(i,gamma_dot_s[i], pl.getEtaEff(gamma_dot_s[i]),taus[i])
895    
896         def test_PowerLaw_QuadSmall_CubeLarge(self):
897             taus=[0.0, 10.0, 20.0, 30.0, 40.0, 50.0, 60.0, 70.0, 80.0, 90.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0]
898             gamma_dot_s=[0.0, 52.04834386816146, 159.15011432761528, 307.93315072671987, 492.9024823703279, 710.67547196059206, 958.86647890316135, 1235.6627169407443, 1539.6192618120876, 1869.5452645867665, 2224.4346900318833, 2446.8781590350718, 2669.3216280382603, 2891.7650970414484, 3114.2085660446369, 3336.6520350478249]
899    
900             pl=PowerLaw(numMaterials=3,verbose=VERBOSE)
901             pl.setDruckerPragerLaw(tau_Y=100.)
902             pl.setPowerLaws(eta_N=[2.,10.,1./16.],tau_t=[1, 25.,64.], power=[1,2,3])
903             pl.setEtaTolerance(self.TOL)
904             for i in xrange(len(taus)): self.checkResult(i,gamma_dot_s[i], pl.getEtaEff(gamma_dot_s[i]),taus[i])
905    
906         def test_PowerLaw_QuadSmall_CubeSmall(self):
907             taus=[0.0, 10.0, 20.0, 30.0, 40.0, 50.0, 60.0, 70.0, 80.0, 90.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0]
908             gamma_dot_s=[0.0, 6.0966144153949529, 13.262466981455987, 21.182803498847885, 29.738072637409996, 38.857111853352741, 48.49087321962044, 58.602735137169738, 69.16391932355954, 80.150989017127827, 91.544346900318828, 100.69878159035071, 109.85321628038258, 119.00765097041447, 128.16208566044637, 137.31652035047824]
909             pl=PowerLaw(numMaterials=3,verbose=VERBOSE)
910             pl.setDruckerPragerLaw(tau_Y=100.)
911             pl.setPowerLaws(eta_N=[2.,10.,25./4.],tau_t=[1, 25.,64.], power=[1,2,3])
912             pl.setEtaTolerance(self.TOL)
913             for i in xrange(len(taus)): self.checkResult(i,gamma_dot_s[i], pl.getEtaEff(gamma_dot_s[i]),taus[i])
914    
915         def test_PowerLaw_withShear(self):
916             taus=[0.0, 10.0, 20.0, 30.0, 40.0, 50.0, 60.0, 70.0, 80.0, 90.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0]
917             gamma_dot_s=[0.0, 15.0, 30.0, 45.0, 60.0, 75.0, 90.0, 105.0, 120.0, 135.0, 150.0, 165.0, 180.0, 195.0, 210.0, 225.0]
918             pl=PowerLaw(numMaterials=1,verbose=VERBOSE)
919             pl.setDruckerPragerLaw(tau_Y=100.)
920             pl.setPowerLaw(eta_N=2.)
921             pl.setElasticShearModulus(3.)
922             dt=1./3.
923             pl.setEtaTolerance(self.TOL)
924             self.failUnlessRaises(ValueError, pl.getEtaEff,gamma_dot_s[0])
925             for i in xrange(len(taus)): self.checkResult(i,gamma_dot_s[i], pl.getEtaEff(gamma_dot_s[i],dt=dt),taus[i])
926    
927  if __name__ == '__main__':  if __name__ == '__main__':
928     suite = unittest.TestSuite()     suite = unittest.TestSuite()
929     suite.addTest(unittest.makeSuite(Test_Simple2DModels))    
930     # suite.addTest(unittest.makeSuite(Test_Simple3DModels))     suite.addTest(unittest.makeSuite(Test_StokesProblemCartesian2D))
931       suite.addTest(unittest.makeSuite(Test_Darcy3D))
932       suite.addTest(unittest.makeSuite(Test_Darcy2D))
933       suite.addTest(unittest.makeSuite(Test_StokesProblemCartesian3D))
934       suite.addTest(unittest.makeSuite(Test_Mountains3D))
935       suite.addTest(unittest.makeSuite(Test_Mountains2D))
936       suite.addTest(unittest.makeSuite(Test_Rheologies))
937     s=unittest.TextTestRunner(verbosity=2).run(suite)     s=unittest.TextTestRunner(verbosity=2).run(suite)
938     if not s.wasSuccessful(): sys.exit(1)     if not s.wasSuccessful(): sys.exit(1)
939    

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

  ViewVC Help
Powered by ViewVC 1.1.26