/[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 2548 by jfenwick, Mon Jul 20 06:20:06 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-2009 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
 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, PowerLaw, IncompressibleIsotropicFlowCartesian
40  from esys.finley import Rectangle, Brick  from esys.finley import Rectangle, Brick
41  class Test_Simple2DModels(unittest.TestCase):  
42    from esys.escript.models import Mountains
43    from math import pi
44    
45    #====================================================================================================================
46    class Test_StokesProblemCartesian2D(unittest.TestCase):
47     def setUp(self):     def setUp(self):
48           NE=6
49           self.TOL=1e-3
50         self.domain=Rectangle(NE,NE,order=2,useFullElementOrder=True)         self.domain=Rectangle(NE,NE,order=2,useFullElementOrder=True)
51     def tearDown(self):     def tearDown(self):
52         del self.domain         del self.domain
53     def test_StokesProblemCartesian_PCG_P_0(self):     def test_PCG_P_0(self):
54         ETA=1.         ETA=1.
55         P1=0.         P1=0.
56    
# Line 56  class Test_Simple2DModels(unittest.TestC Line 65  class Test_Simple2DModels(unittest.TestC
65                
66         sp.initialize(f=F,fixed_u_mask=mask,eta=ETA)         sp.initialize(f=F,fixed_u_mask=mask,eta=ETA)
67         u0=(1-x[0])*x[0]*[0.,1.]         u0=(1-x[0])*x[0]*[0.,1.]
68         p0=Scalar(P1,ReducedSolution(self.domain))         p0=Scalar(-P1,ReducedSolution(self.domain))
69         u,p=sp.solve(u0,p0,show_details=VERBOSE, verbose=VERBOSE,max_iter=100,solver="PCG")         sp.setTolerance(self.TOL)
70           u,p=sp.solve(u0,p0,verbose=VERBOSE,max_iter=100,usePCG=True)
71                
72         error_v0=Lsup(u[0]-u0[0])         error_v0=Lsup(u[0]-u0[0])
73         error_v1=Lsup(u[1]-u0[1])/0.25         error_v1=Lsup(u[1]-u0[1])/0.25
74         zz=P1*x[0]*x[1]-p         error_p=Lsup(p+P1*x[0]*x[1])
75         error_p=Lsup(zz-integrate(zz))         self.failUnless(error_p<10*self.TOL, "pressure error too large.")
76         # print error_p, error_v0,error_v1         self.failUnless(error_v0<10*self.TOL, "0-velocity error too large.")
77         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.")  
78    
79     def test_StokesProblemCartesian_PCG_P_small(self):     def test_PCG_P_small(self):
80         ETA=1.         ETA=1.
81         P1=1.         P1=1.
82    
# Line 83  class Test_Simple2DModels(unittest.TestC Line 91  class Test_Simple2DModels(unittest.TestC
91                
92         sp.initialize(f=F,fixed_u_mask=mask,eta=ETA)         sp.initialize(f=F,fixed_u_mask=mask,eta=ETA)
93         u0=(1-x[0])*x[0]*[0.,1.]         u0=(1-x[0])*x[0]*[0.,1.]
94         p0=Scalar(P1,ReducedSolution(self.domain))         p0=Scalar(-P1,ReducedSolution(self.domain))
95         u,p=sp.solve(u0,p0,show_details=VERBOSE, verbose=VERBOSE,max_iter=100,solver="PCG")         sp.setTolerance(self.TOL*0.2)
96                 u,p=sp.solve(u0,p0, verbose=VERBOSE,max_iter=100,usePCG=True)
97         error_v0=Lsup(u[0]-u0[0])         error_v0=Lsup(u[0]-u0[0])
98         error_v1=Lsup(u[1]-u0[1])/0.25         error_v1=Lsup(u[1]-u0[1])/0.25
99         zz=P1*x[0]*x[1]-p         error_p=Lsup(P1*x[0]*x[1]+p)
100         error_p=Lsup(zz-integrate(zz))         self.failUnless(error_v0<10*self.TOL, "0-velocity error too large.")
101         # print error_p, error_v0,error_v1         self.failUnless(error_v1<10*self.TOL, "1-velocity error too large.")
102         self.failUnless(error_p<TOL,"pressure error too large.")         self.failUnless(error_p<10*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.")  
103    
104     def test_StokesProblemCartesian_PCG_P_large(self):     def test_PCG_P_large(self):
105         ETA=1.         ETA=1.
106         P1=1000.         P1=1000.
107    
# Line 110  class Test_Simple2DModels(unittest.TestC Line 116  class Test_Simple2DModels(unittest.TestC
116                
117         sp.initialize(f=F,fixed_u_mask=mask,eta=ETA)         sp.initialize(f=F,fixed_u_mask=mask,eta=ETA)
118         u0=(1-x[0])*x[0]*[0.,1.]         u0=(1-x[0])*x[0]*[0.,1.]
119         p0=Scalar(P1,ReducedSolution(self.domain))         p0=Scalar(-P1,ReducedSolution(self.domain))
120         u,p=sp.solve(u0,p0,show_details=VERBOSE, verbose=VERBOSE,max_iter=100,solver="PCG")         sp.setTolerance(self.TOL)
121           u,p=sp.solve(u0,p0, verbose=VERBOSE,max_iter=100,usePCG=True)
122                
123         error_v0=Lsup(u[0]-u0[0])         error_v0=Lsup(u[0]-u0[0])
124         error_v1=Lsup(u[1]-u0[1])/0.25         error_v1=Lsup(u[1]-u0[1])/0.25
125         zz=P1*x[0]*x[1]-p         error_p=Lsup(P1*x[0]*x[1]+p)/P1
126         error_p=Lsup(zz-integrate(zz))/P1         self.failUnless(error_p<10*self.TOL, "pressure error too large.")
127         # print error_p, error_v0,error_v1         self.failUnless(error_v0<10*self.TOL, "0-velocity error too large.")
128         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.")  
129    
130     def test_StokesProblemCartesian_GMRES_P_0(self):     def test_GMRES_P_0(self):
131         ETA=1.         ETA=1.
132         P1=0.         P1=0.
133    
# Line 137  class Test_Simple2DModels(unittest.TestC Line 142  class Test_Simple2DModels(unittest.TestC
142                
143         sp.initialize(f=F,fixed_u_mask=mask,eta=ETA)         sp.initialize(f=F,fixed_u_mask=mask,eta=ETA)
144         u0=(1-x[0])*x[0]*[0.,1.]         u0=(1-x[0])*x[0]*[0.,1.]
145         p0=Scalar(P1,ReducedSolution(self.domain))         p0=Scalar(-P1,ReducedSolution(self.domain))
146         u,p=sp.solve(u0,p0,show_details=VERBOSE, verbose=VERBOSE,max_iter=100,solver="GMRES",iter_restart=20)         sp.setTolerance(self.TOL)
147           u,p=sp.solve(u0,p0, verbose=VERBOSE,max_iter=50,usePCG=False,iter_restart=18)
148                
149         error_v0=Lsup(u[0]-u0[0])         error_v0=Lsup(u[0]-u0[0])
150         error_v1=Lsup(u[1]-u0[1])/0.25         error_v1=Lsup(u[1]-u0[1])/0.25
151         zz=P1*x[0]*x[1]-p         error_p=Lsup(P1*x[0]*x[1]+p)
152         error_p=Lsup(zz-integrate(zz))         self.failUnless(error_p<10*self.TOL, "pressure error too large.")
153           self.failUnless(error_v0<10*self.TOL, "0-velocity error too large.")
154           self.failUnless(error_v1<10*self.TOL, "1-velocity error too large.")
155    
156        # 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):  
157         ETA=1.         ETA=1.
158         P1=1.         P1=1.
159    
# Line 165  class Test_Simple2DModels(unittest.TestC Line 168  class Test_Simple2DModels(unittest.TestC
168                
169         sp.initialize(f=F,fixed_u_mask=mask,eta=ETA)         sp.initialize(f=F,fixed_u_mask=mask,eta=ETA)
170         u0=(1-x[0])*x[0]*[0.,1.]         u0=(1-x[0])*x[0]*[0.,1.]
171         p0=Scalar(P1,ReducedSolution(self.domain))         p0=Scalar(-P1,ReducedSolution(self.domain))
172         u,p=sp.solve(u0,p0,show_details=VERBOSE, verbose=VERBOSE,max_iter=100,solver="GMRES")         sp.setTolerance(self.TOL*0.1)
173           u,p=sp.solve(u0,p0, verbose=VERBOSE,max_iter=20,usePCG=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         error_p=Lsup(P1*x[0]*x[1]+p)
178         error_p=Lsup(zz-integrate(zz))         self.failUnless(error_p<10*self.TOL, "pressure error too large.")
179         # print error_p, error_v0,error_v1         self.failUnless(error_v0<10*self.TOL, "0-velocity error too large.")
180         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.")  
181    
182     def test_StokesProblemCartesian_GMRES_P_large(self):     def test_GMRES_P_large(self):
183         ETA=1.         ETA=1.
184         P1=1000.         P1=1000.
185    
# Line 192  class Test_Simple2DModels(unittest.TestC Line 194  class Test_Simple2DModels(unittest.TestC
194                
195         sp.initialize(f=F,fixed_u_mask=mask,eta=ETA)         sp.initialize(f=F,fixed_u_mask=mask,eta=ETA)
196         u0=(1-x[0])*x[0]*[0.,1.]         u0=(1-x[0])*x[0]*[0.,1.]
197         p0=Scalar(P1,ReducedSolution(self.domain))         p0=Scalar(-P1,ReducedSolution(self.domain))
198         u,p=sp.solve(u0,p0,show_details=VERBOSE, verbose=VERBOSE,max_iter=100,solver="GMRES")         sp.setTolerance(self.TOL)
199           u,p=sp.solve(u0,p0, verbose=VERBOSE,max_iter=100,usePCG=False)
200                
201         error_v0=Lsup(u[0]-u0[0])         error_v0=Lsup(u[0]-u0[0])
202         error_v1=Lsup(u[1]-u0[1])/0.25         error_v1=Lsup(u[1]-u0[1])/0.25
203         zz=P1*x[0]*x[1]-p         error_p=Lsup(P1*x[0]*x[1]+p)/P1
204         error_p=Lsup(zz-integrate(zz))/P1         self.failUnless(error_p<10*self.TOL, "pressure error too large.")
205         # print error_p, error_v0,error_v1         self.failUnless(error_v0<10*self.TOL, "0-velocity error too large.")
206         self.failUnless(error_p<TOL,"pressure error too large.")         self.failUnless(error_v1<10*self.TOL, "1-velocity error too large.")
207         self.failUnless(error_v0<TOL,"0-velocity error too large.")  #====================================================================================================================
208         self.failUnless(error_v1<TOL,"1-velocity error too large.")  class Test_StokesProblemCartesian3D(unittest.TestCase):
   
    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")  
         
        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_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):  
209     def setUp(self):     def setUp(self):
210           NE=6
211           self.TOL=1e-4
212         self.domain=Brick(NE,NE,NE,order=2,useFullElementOrder=True)         self.domain=Brick(NE,NE,NE,order=2,useFullElementOrder=True)
213     def tearDown(self):     def tearDown(self):
214         del self.domain         del self.domain
215     def test_StokesProblemCartesian_PCG_P_0(self):     def test_PCG_P_0(self):
216         ETA=1.         ETA=1.
217         P1=0.         P1=0.
218    
# Line 383  class Test_Simple3DModels(unittest.TestC Line 221  class Test_Simple3DModels(unittest.TestC
221         x=self.domain.getX()         x=self.domain.getX()
222         mask=whereZero(x[0])    * [1.,1.,1.] \         mask=whereZero(x[0])    * [1.,1.,1.] \
223                +whereZero(x[0]-1)  * [1.,1.,1.] \                +whereZero(x[0]-1)  * [1.,1.,1.] \
224                +whereZero(x[1])    * [1.,1.,1.] \                +whereZero(x[1])    * [1.,0.,1.] \
225                +whereZero(x[1]-1)  * [1.,1.,1.] \                +whereZero(x[1]-1)  * [1.,1.,1.] \
226                +whereZero(x[2])    * [1.,1.,0.] \                +whereZero(x[2])    * [1.,1.,0.] \
227                +whereZero(x[2]-1)  * [1.,1.,1.]                +whereZero(x[2]-1)  * [1.,1.,1.]
# Line 393  class Test_Simple3DModels(unittest.TestC Line 231  class Test_Simple3DModels(unittest.TestC
231                
232         sp.initialize(f=F,fixed_u_mask=mask,eta=ETA)         sp.initialize(f=F,fixed_u_mask=mask,eta=ETA)
233         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.]
234         p0=Scalar(P1,ReducedSolution(self.domain))         p0=Scalar(-P1,ReducedSolution(self.domain))
235         u,p=sp.solve(u0,p0,show_details=VERBOSE, verbose=VERBOSE,max_iter=100,solver="PCG")         sp.setTolerance(self.TOL)
236           u,p=sp.solve(u0,p0, verbose=VERBOSE ,max_iter=100,usePCG=True)
237                
238         error_v0=Lsup(u[0]-u0[0])         error_v0=Lsup(u[0]-u0[0])
239         error_v1=Lsup(u[1]-u0[1])         error_v1=Lsup(u[1]-u0[1])
240         error_v2=Lsup(u[2]-u0[2])/0.25**2         error_v2=Lsup(u[2]-u0[2])/0.25**2
241         zz=P1*x[0]*x[1]*x[2]-p         error_p=Lsup(P1*x[0]*x[1]*x[2]+p)
242         error_p=Lsup(zz-integrate(zz))         self.failUnless(error_p<10*self.TOL, "pressure error too large.")
243         # print error_p, error_v0,error_v1,error_v2         self.failUnless(error_v0<10*self.TOL, "0-velocity error too large.")
244         self.failUnless(error_p<TOL,"pressure error too large.")         self.failUnless(error_v1<10*self.TOL, "1-velocity error too large.")
245         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.")  
246    
247     def test_StokesProblemCartesian_PCG_P_small(self):     def test_PCG_P_small(self):
248         ETA=1.         ETA=1.
249         P1=1.         P1=1.
250    
# Line 415  class Test_Simple3DModels(unittest.TestC Line 252  class Test_Simple3DModels(unittest.TestC
252         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.]
253         mask=whereZero(x[0])    * [1.,1.,1.] \         mask=whereZero(x[0])    * [1.,1.,1.] \
254                +whereZero(x[0]-1)  * [1.,1.,1.] \                +whereZero(x[0]-1)  * [1.,1.,1.] \
255                +whereZero(x[1])    * [1.,1.,1.] \                +whereZero(x[1])    * [1.,0.,1.] \
256                +whereZero(x[1]-1)  * [1.,1.,1.] \                +whereZero(x[1]-1)  * [1.,1.,1.] \
257                +whereZero(x[2])    * [1.,1.,0.] \                +whereZero(x[2])    * [1.,1.,0.] \
258                +whereZero(x[2]-1)  * [1.,1.,1.]                +whereZero(x[2]-1)  * [1.,1.,1.]
# Line 425  class Test_Simple3DModels(unittest.TestC Line 262  class Test_Simple3DModels(unittest.TestC
262                
263         sp.initialize(f=F,fixed_u_mask=mask,eta=ETA)         sp.initialize(f=F,fixed_u_mask=mask,eta=ETA)
264         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.]
265         p0=Scalar(P1,ReducedSolution(self.domain))         p0=Scalar(-P1,ReducedSolution(self.domain))
266         u,p=sp.solve(u0,p0,show_details=VERBOSE, verbose=VERBOSE,max_iter=100,solver="PCG")         sp.setTolerance(self.TOL*0.1)
267                 u,p=sp.solve(u0,p0, verbose=VERBOSE ,max_iter=100,usePCG=True)
268         error_v0=Lsup(u[0]-u0[0])         error_v0=Lsup(u[0]-u0[0])
269         error_v1=Lsup(u[1]-u0[1])         error_v1=Lsup(u[1]-u0[1])
270         error_v2=Lsup(u[2]-u0[2])/0.25**2         error_v2=Lsup(u[2]-u0[2])/0.25**2
271         zz=P1*x[0]*x[1]*x[2]-p         error_p=Lsup(P1*x[0]*x[1]*x[2]+p)
272         error_p=Lsup(zz-integrate(zz))/P1         self.failUnless(error_p<10*self.TOL, "pressure error too large.")
273         # print error_p, error_v0,error_v1,error_v2         self.failUnless(error_v0<10*self.TOL, "0-velocity error too large.")
274         self.failUnless(error_p<TOL,"pressure error too large.")         self.failUnless(error_v1<10*self.TOL, "1-velocity error too large.")
275         self.failUnless(error_v0<TOL,"0-velocity error too large.")         self.failUnless(error_v2<10*self.TOL, "2-velocity error too large.")
276         self.failUnless(error_v1<TOL,"1-velocity error too large.")  
277         self.failUnless(error_v2<TOL,"2-velocity error too large.")     def test_PCG_P_large(self):
    def test_StokesProblemCartesian_PCG_P_large(self):  
278         ETA=1.         ETA=1.
279         P1=1000.         P1=1000.
280    
# Line 446  class Test_Simple3DModels(unittest.TestC Line 282  class Test_Simple3DModels(unittest.TestC
282         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.]
283         mask=whereZero(x[0])    * [1.,1.,1.] \         mask=whereZero(x[0])    * [1.,1.,1.] \
284                +whereZero(x[0]-1)  * [1.,1.,1.] \                +whereZero(x[0]-1)  * [1.,1.,1.] \
285                +whereZero(x[1])    * [1.,1.,1.] \                +whereZero(x[1])    * [1.,0.,1.] \
286                +whereZero(x[1]-1)  * [1.,1.,1.] \                +whereZero(x[1]-1)  * [1.,1.,1.] \
287                +whereZero(x[2])    * [1.,1.,0.] \                +whereZero(x[2])    * [1.,1.,0.] \
288                +whereZero(x[2]-1)  * [1.,1.,1.]                +whereZero(x[2]-1)  * [1.,1.,1.]
# Line 456  class Test_Simple3DModels(unittest.TestC Line 292  class Test_Simple3DModels(unittest.TestC
292                
293         sp.initialize(f=F,fixed_u_mask=mask,eta=ETA)         sp.initialize(f=F,fixed_u_mask=mask,eta=ETA)
294         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.]
295         p0=Scalar(P1,ReducedSolution(self.domain))         p0=Scalar(-P1,ReducedSolution(self.domain))
296         u,p=sp.solve(u0,p0,show_details=VERBOSE, verbose=VERBOSE,max_iter=100,solver="PCG")         sp.setTolerance(self.TOL)
297           u,p=sp.solve(u0,-p0, verbose=VERBOSE ,max_iter=100,usePCG=True)
298                
299         error_v0=Lsup(u[0]-u0[0])         error_v0=Lsup(u[0]-u0[0])
300         error_v1=Lsup(u[1]-u0[1])         error_v1=Lsup(u[1]-u0[1])
301         error_v2=Lsup(u[2]-u0[2])/0.25**2         error_v2=Lsup(u[2]-u0[2])/0.25**2
302         zz=P1*x[0]*x[1]*x[2]-p         error_p=Lsup(P1*x[0]*x[1]*x[2]+p)/P1
303         error_p=Lsup(zz-integrate(zz))/P1         self.failUnless(error_p<10*self.TOL, "pressure error too large.")
304         # print error_p, error_v0,error_v1,error_v2         self.failUnless(error_v0<10*self.TOL, "0-velocity error too large.")
305         self.failUnless(error_p<TOL,"pressure error too large.")         self.failUnless(error_v1<10*self.TOL, "1-velocity error too large.")
306         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.")  
307    
308     def test_StokesProblemCartesian_GMRES_P_0(self):     def test_GMRES_P_0(self):
309         ETA=1.         ETA=1.
310         P1=0.         P1=0.
311    
# Line 489  class Test_Simple3DModels(unittest.TestC Line 324  class Test_Simple3DModels(unittest.TestC
324                
325         sp.initialize(f=F,fixed_u_mask=mask,eta=ETA)         sp.initialize(f=F,fixed_u_mask=mask,eta=ETA)
326         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.]
327         p0=Scalar(P1,ReducedSolution(self.domain))         p0=Scalar(-P1,ReducedSolution(self.domain))
328         u,p=sp.solve(u0,p0,show_details=VERBOSE, verbose=VERBOSE,max_iter=100,solver="GMRES",iter_restart=20)         sp.setTolerance(self.TOL)
329           u,p=sp.solve(u0,p0, verbose=VERBOSE,max_iter=100,usePCG=False,iter_restart=20)
330                
331         error_v0=Lsup(u[0]-u0[0])         error_v0=Lsup(u[0]-u0[0])
332         error_v1=Lsup(u[1]-u0[1])         error_v1=Lsup(u[1]-u0[1])
333         error_v2=Lsup(u[2]-u0[2])/0.25**2         error_v2=Lsup(u[2]-u0[2])/0.25**2
334         zz=P1*x[0]*x[1]*x[2]-p         error_p=Lsup(P1*x[0]*x[1]*x[2]+p)
335         error_p=Lsup(zz-integrate(zz))         self.failUnless(error_p<10*self.TOL, "pressure error too large.")
336         # print error_p, error_v0,error_v1,error_v2         self.failUnless(error_v0<10*self.TOL, "0-velocity error too large.")
337         self.failUnless(error_p<TOL,"pressure error too large.")         self.failUnless(error_v1<10*self.TOL, "1-velocity error too large.")
338         self.failUnless(error_v0<TOL,"0-velocity error too large.")         self.failUnless(error_v2<10*self.TOL, "2-velocity error too large.")
339         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):  
340         ETA=1.         ETA=1.
341         P1=1.         P1=1.
342    
# Line 520  class Test_Simple3DModels(unittest.TestC Line 354  class Test_Simple3DModels(unittest.TestC
354                
355         sp.initialize(f=F,fixed_u_mask=mask,eta=ETA)         sp.initialize(f=F,fixed_u_mask=mask,eta=ETA)
356         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.]
357         p0=Scalar(P1,ReducedSolution(self.domain))         p0=Scalar(-P1,ReducedSolution(self.domain))
358         u,p=sp.solve(u0,p0,show_details=VERBOSE, verbose=VERBOSE,max_iter=100,solver="GMRES")         sp.setTolerance(self.TOL*0.1)
359           u,p=sp.solve(u0,p0, verbose=VERBOSE,max_iter=100,usePCG=False)
360                
361         error_v0=Lsup(u[0]-u0[0])         error_v0=Lsup(u[0]-u0[0])
362         error_v1=Lsup(u[1]-u0[1])         error_v1=Lsup(u[1]-u0[1])
363         error_v2=Lsup(u[2]-u0[2])/0.25**2         error_v2=Lsup(u[2]-u0[2])/0.25**2
364         zz=P1*x[0]*x[1]*x[2]-p         error_p=Lsup(P1*x[0]*x[1]*x[2]+p)/P1
365         error_p=Lsup(zz-integrate(zz))/P1         self.failUnless(error_v0<10*self.TOL, "0-velocity error too large.")
366         # print error_p, error_v0,error_v1,error_v2         self.failUnless(error_v1<10*self.TOL, "1-velocity error too large.")
367         self.failUnless(error_p<TOL,"pressure error too large.")         self.failUnless(error_v2<10*self.TOL, "2-velocity error too large.")
368         self.failUnless(error_v0<TOL,"0-velocity error too large.")         self.failUnless(error_p<10*self.TOL, "pressure error too large.")
369         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):  
370         ETA=1.         ETA=1.
371         P1=1000.         P1=1000.
372    
# Line 541  class Test_Simple3DModels(unittest.TestC Line 374  class Test_Simple3DModels(unittest.TestC
374         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.]
375         mask=whereZero(x[0])    * [1.,1.,1.] \         mask=whereZero(x[0])    * [1.,1.,1.] \
376                +whereZero(x[0]-1)  * [1.,1.,1.] \                +whereZero(x[0]-1)  * [1.,1.,1.] \
377                +whereZero(x[1])    * [1.,1.,1.] \                +whereZero(x[1])    * [1.,0.,1.] \
378                +whereZero(x[1]-1)  * [1.,1.,1.] \                +whereZero(x[1]-1)  * [1.,1.,1.] \
379                +whereZero(x[2])    * [1.,1.,0.] \                +whereZero(x[2])    * [1.,1.,0.] \
380                +whereZero(x[2]-1)  * [1.,1.,1.]                +whereZero(x[2]-1)  * [1.,1.,1.]
# Line 551  class Test_Simple3DModels(unittest.TestC Line 384  class Test_Simple3DModels(unittest.TestC
384                
385         sp.initialize(f=F,fixed_u_mask=mask,eta=ETA)         sp.initialize(f=F,fixed_u_mask=mask,eta=ETA)
386         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.]
387         p0=Scalar(P1,ReducedSolution(self.domain))         p0=Scalar(-P1,ReducedSolution(self.domain))
388         u,p=sp.solve(u0,p0,show_details=VERBOSE, verbose=VERBOSE,max_iter=100,solver="GMRES")         sp.setTolerance(self.TOL)
389           u,p=sp.solve(u0,p0, verbose=VERBOSE ,max_iter=100,usePCG=False)
390                
391         error_v0=Lsup(u[0]-u0[0])         error_v0=Lsup(u[0]-u0[0])
392         error_v1=Lsup(u[1]-u0[1])         error_v1=Lsup(u[1]-u0[1])
393         error_v2=Lsup(u[2]-u0[2])/0.25**2         error_v2=Lsup(u[2]-u0[2])/0.25**2
394         zz=P1*x[0]*x[1]*x[2]-p         error_p=Lsup(P1*x[0]*x[1]*x[2]+p)/P1
395         error_p=Lsup(zz-integrate(zz))/P1         self.failUnless(error_p<10*self.TOL, "pressure error too large.")
396         # print error_p, error_v0,error_v1,error_v2         self.failUnless(error_v0<10*self.TOL, "0-velocity error too large.")
397         self.failUnless(error_p<TOL,"pressure error too large.")         self.failUnless(error_v1<10*self.TOL, "1-velocity error too large.")
398         self.failUnless(error_v0<TOL,"0-velocity error too large.")         self.failUnless(error_v2<10*self.TOL, "2-velocity error too large.")
399         self.failUnless(error_v1<TOL,"1-velocity error too large.")  #====================================================================================================================
400         self.failUnless(error_v2<TOL,"2-velocity error too large.")  class Test_Darcy(unittest.TestCase):
401        # this is a simple test for the darcy flux problem
402        #
403        #
404        #  p = 1/k * ( 1/2* (fb-f0)/xb* x **2 + f0 * x - ub*x ) +  p0
405        #
406        #  with f = (fb-f0)/xb* x + f0
407        #
408        #    u = f - k * p,x = ub
409        #
410        #  we prescribe pressure at x=x0=0 to p0
411        #
412        #  if we prescribe the pressure on the bottom x=xb we set
413        #
414        #  pb= 1/k * ( 1/2* (fb-f0)* xb + f0 * xb - ub*xb ) +  p0 = 1/k * ((fb+f0)/2 - ub ) * xb  + p0
415        #
416        #  which leads to ub = (fb+f0)/2-k*(pb-p0)/xb
417        #
418        def rescaleDomain(self):
419            x=self.dom.getX().copy()
420            for i in xrange(self.dom.getDim()):
421                 x_inf=inf(x[i])
422                 x_sup=sup(x[i])
423                 if i == self.dom.getDim()-1:
424                    x[i]=-self.WIDTH*(x[i]-x_sup)/(x_inf-x_sup)
425                 else:
426                    x[i]=self.WIDTH*(x[i]-x_inf)/(x_sup-x_inf)
427            self.dom.setX(x)
428        def getScalarMask(self,include_bottom=True):
429            x=self.dom.getX().copy()
430            x_inf=inf(x[self.dom.getDim()-1])
431            x_sup=sup(x[self.dom.getDim()-1])
432            out=whereZero(x[self.dom.getDim()-1]-x_sup)
433            if include_bottom: out+=whereZero(x[self.dom.getDim()-1]-x_inf)
434            return wherePositive(out)
435        def getVectorMask(self,include_bottom=True):
436            x=self.dom.getX().copy()
437            out=Vector(0.,Solution(self.dom))
438            for i in xrange(self.dom.getDim()):
439                 x_inf=inf(x[i])
440                 x_sup=sup(x[i])
441                 if i != self.dom.getDim()-1: out[i]+=whereZero(x[i]-x_sup)
442                 if i != self.dom.getDim()-1 or include_bottom: out[i]+=whereZero(x[i]-x_inf)
443            return wherePositive(out)
444    
445        def setSolutionFixedBottom(self, p0, pb, f0, fb, k):
446             d=self.dom.getDim()
447             x=self.dom.getX()[d-1]
448             xb=inf(x)
449             u=Vector(0.,Solution(self.dom))+kronecker(d)[d-1]*((f0+fb)/2.-k*(pb-p0)/xb)
450             p=1./k*((fb-f0)/(xb*2.)* x**2 - (fb-f0)/2.*x)+(pb-p0)/xb*x +  p0
451             f= ((fb-f0)/xb* x + f0)*kronecker(Function(self.dom))[d-1]
452             return u,p,f
453            
454        def testConstF_FixedBottom_smallK(self):
455            k=1.e-10
456            mp=self.getScalarMask(include_bottom=True)
457            mv=self.getVectorMask(include_bottom=False)
458            u_ref,p_ref,f=self.setSolutionFixedBottom(p0=2500,pb=4000.,f0=10.,fb=10.,k=k)
459            p=p_ref*mp
460            u=u_ref*mv
461            df=DarcyFlow(self.dom)
462            df.setValue(g=f,
463                          location_of_fixed_pressure=mp,
464                          location_of_fixed_flux=mv,
465                          permeability=Scalar(k,Function(self.dom)))
466            df.setTolerance(rtol=self.TOL)
467            v,p=df.solve(u_ref,p, max_iter=100, verbose=VERBOSE)
468            self.failUnless(Lsup(v-u_ref)<self.TOL*10.*Lsup(u_ref), "flux error too big.")
469            self.failUnless(Lsup(p-p_ref)<self.TOL*10.*Lsup(p_ref), "pressure error too big.")
470        def testConstF_FixedBottom_mediumK(self):
471            k=1.
472            mp=self.getScalarMask(include_bottom=True)
473            mv=self.getVectorMask(include_bottom=False)
474            u_ref,p_ref,f=self.setSolutionFixedBottom(p0=2500,pb=4000.,f0=10.,fb=10.,k=k)
475            p=p_ref*mp
476            u=u_ref*mv
477            df=DarcyFlow(self.dom)
478            df.setValue(g=f,
479                          location_of_fixed_pressure=mp,
480                          location_of_fixed_flux=mv,
481                          permeability=Scalar(k,Function(self.dom)))
482            df.setTolerance(rtol=self.TOL)
483            v,p=df.solve(u,p,max_iter=100, verbose=VERBOSE )
484            self.failUnless(Lsup(p-p_ref)<self.TOL*10.*Lsup(p_ref), "pressure error too big.")
485            self.failUnless(Lsup(v-u_ref)<self.TOL*10.*Lsup(u_ref), "flux error too big.")
486    
487        def testConstF_FixedBottom_largeK(self):
488            k=1.e10
489            mp=self.getScalarMask(include_bottom=True)
490            mv=self.getVectorMask(include_bottom=False)
491            u_ref,p_ref,f=self.setSolutionFixedBottom(p0=2500,pb=4000.,f0=10.,fb=10.,k=k)
492            p=p_ref*mp
493            u=u_ref*mv
494            df=DarcyFlow(self.dom)
495            df.setValue(g=f,
496                          location_of_fixed_pressure=mp,
497                          location_of_fixed_flux=mv,
498                          permeability=Scalar(k,Function(self.dom)))
499            df.setTolerance(rtol=self.TOL)
500            v,p=df.solve(u,p, max_iter=100, verbose=VERBOSE)
501            self.failUnless(Lsup(v-u_ref)<self.TOL*10.*Lsup(u_ref), "flux error too big.")
502            self.failUnless(Lsup(p-p_ref)<self.TOL*10.*Lsup(p_ref), "pressure error too big.")
503    
504        def testVarioF_FixedBottom_smallK(self):
505            k=1.e-10
506            mp=self.getScalarMask(include_bottom=True)
507            mv=self.getVectorMask(include_bottom=False)
508            u_ref,p_ref,f=self.setSolutionFixedBottom(p0=2500,pb=4000.,f0=10.,fb=30.,k=k)
509            p=p_ref*mp
510            u=u_ref*mv
511            df=DarcyFlow(self.dom)
512            df.setValue(g=f,
513                          location_of_fixed_pressure=mp,
514                          location_of_fixed_flux=mv,
515                          permeability=Scalar(k,Function(self.dom)))
516            df.setTolerance(rtol=self.TOL)
517            v,p=df.solve(u,p, max_iter=100, verbose=VERBOSE)
518            self.failUnless(Lsup(v-u_ref)<self.TOL*10.*Lsup(u_ref), "flux error too big.")
519            self.failUnless(Lsup(p-p_ref)<self.TOL*10.*Lsup(p_ref), "pressure error too big.")
520    
521        def testVarioF_FixedBottom_mediumK(self):
522            k=1.
523            mp=self.getScalarMask(include_bottom=True)
524            mv=self.getVectorMask(include_bottom=False)
525            u_ref,p_ref,f=self.setSolutionFixedBottom(p0=2500,pb=4000.,f0=10.,fb=30.,k=k)
526            p=p_ref*mp
527            u=u_ref*mv
528            df=DarcyFlow(self.dom)
529            df.setValue(g=f,
530                          location_of_fixed_pressure=mp,
531                          location_of_fixed_flux=mv,
532                          permeability=Scalar(k,Function(self.dom)))
533            df.setTolerance(rtol=self.TOL)
534            v,p=df.solve(u,p, max_iter=100, verbose=VERBOSE)
535            self.failUnless(Lsup(v-u_ref)<self.TOL*10.*Lsup(u_ref), "flux error too big.")
536            self.failUnless(Lsup(p-p_ref)<self.TOL*10.*Lsup(p_ref), "pressure error too big.")
537    
538        def testVarioF_FixedBottom_largeK(self):
539            k=1.e10
540            mp=self.getScalarMask(include_bottom=True)
541            mv=self.getVectorMask(include_bottom=False)
542            u_ref,p_ref,f=self.setSolutionFixedBottom(p0=2500,pb=4000.,f0=10.,fb=30.,k=k)
543            p=p_ref*mp
544            u=u_ref*mv
545            df=DarcyFlow(self.dom)
546            df.setValue(g=f,
547                          location_of_fixed_pressure=mp,
548                          location_of_fixed_flux=mv,
549                          permeability=Scalar(k,Function(self.dom)))
550            df.setTolerance(rtol=self.TOL)
551            v,p=df.solve(u,p, max_iter=100, verbose=VERBOSE)
552            self.failUnless(Lsup(v-u_ref)<self.TOL*10.*Lsup(u_ref), "flux error too big.")
553            self.failUnless(Lsup(p-p_ref)<self.TOL*10.*Lsup(p_ref), "pressure error too big.")
554    
555        def testConstF_FreeBottom_smallK(self):
556            k=1.e-10
557            mp=self.getScalarMask(include_bottom=False)
558            mv=self.getVectorMask(include_bottom=True)
559            u_ref,p_ref,f=self.setSolutionFixedBottom(p0=2500,pb=4000.,f0=10.,fb=10.,k=k)
560            p=p_ref*mp
561            u=u_ref*mv
562            df=DarcyFlow(self.dom)
563            df.setValue(g=f,
564                        location_of_fixed_pressure=mp,
565                          location_of_fixed_flux=mv,
566                          permeability=Scalar(k,Function(self.dom)))
567            df.setTolerance(rtol=self.TOL)
568            v,p=df.solve(u,p, max_iter=100, verbose=VERBOSE)
569            self.failUnless(Lsup(v-u_ref)<self.TOL*10.*Lsup(u_ref), "flux error too big.")
570            self.failUnless(Lsup(p-p_ref)<self.TOL*10.*Lsup(p_ref), "pressure error too big.")
571    
572        def testConstF_FreeBottom_mediumK(self):
573            k=1.
574            mp=self.getScalarMask(include_bottom=False)
575            mv=self.getVectorMask(include_bottom=True)
576            u_ref,p_ref,f=self.setSolutionFixedBottom(p0=2500,pb=4000.,f0=10.,fb=10.,k=k)
577            p=p_ref*mp
578            u=u_ref*mv
579            df=DarcyFlow(self.dom)
580            df.setValue(g=f,
581                          location_of_fixed_pressure=mp,
582                          location_of_fixed_flux=mv,
583                          permeability=Scalar(k,Function(self.dom)))
584            df.setTolerance(rtol=self.TOL)
585            v,p=df.solve(u,p, max_iter=100, verbose=VERBOSE)
586            self.failUnless(Lsup(v-u_ref)<self.TOL*10.*Lsup(u_ref), "flux error too big.")
587            self.failUnless(Lsup(p-p_ref)<self.TOL*10.*Lsup(p_ref), "pressure error too big.")
588    
589        def testConstF_FreeBottom_largeK(self):
590            k=1.e10
591            mp=self.getScalarMask(include_bottom=False)
592            mv=self.getVectorMask(include_bottom=True)
593            u_ref,p_ref,f=self.setSolutionFixedBottom(p0=2500,pb=4000.,f0=10.,fb=10.,k=k)
594            p=p_ref*mp
595            u=u_ref*mv
596            df=DarcyFlow(self.dom)
597            df.setValue(g=f,
598                          location_of_fixed_pressure=mp,
599                          location_of_fixed_flux=mv,
600                          permeability=Scalar(k,Function(self.dom)))
601            df.setTolerance(rtol=self.TOL)
602            v,p=df.solve(u,p, max_iter=100, verbose=VERBOSE)
603            self.failUnless(Lsup(v-u_ref)<self.TOL*10.*Lsup(u_ref), "flux error too big.")
604            self.failUnless(Lsup(p-p_ref)<self.TOL*10.*Lsup(p_ref), "pressure error too big.")
605    
606        def testVarioF_FreeBottom_smallK(self):
607            k=1.e-10
608            mp=self.getScalarMask(include_bottom=False)
609            mv=self.getVectorMask(include_bottom=True)
610            u_ref,p_ref,f=self.setSolutionFixedBottom(p0=2500,pb=4000.,f0=10.,fb=30.,k=k)
611            p=p_ref*mp
612            u=u_ref*mv
613            df=DarcyFlow(self.dom)
614            df.setValue(g=f,
615                          location_of_fixed_pressure=mp,
616                          location_of_fixed_flux=mv,
617                          permeability=Scalar(k,Function(self.dom)))
618            df.setTolerance(rtol=self.TOL)
619            v,p=df.solve(u,p, max_iter=100, verbose=VERBOSE)
620            self.failUnless(Lsup(v-u_ref)<self.TOL*25.*Lsup(u_ref), "flux error too big.")  # 25 because of disc. error.
621            self.failUnless(Lsup(p-p_ref)<self.TOL*25.*Lsup(p_ref), "pressure error too big.")
622    
623        def testVarioF_FreeBottom_mediumK(self):
624            k=1.
625            mp=self.getScalarMask(include_bottom=False)
626            mv=self.getVectorMask(include_bottom=True)
627            u_ref,p_ref,f=self.setSolutionFixedBottom(p0=2500,pb=4000.,f0=10.,fb=30.,k=k)
628            p=p_ref*mp
629            u=u_ref*mv
630            df=DarcyFlow(self.dom)
631            df.setValue(g=f,
632                          location_of_fixed_pressure=mp,
633                          location_of_fixed_flux=mv,
634                          permeability=Scalar(k,Function(self.dom)))
635            df.setTolerance(rtol=self.TOL)
636            v,p=df.solve(u,p, max_iter=100, verbose=VERBOSE)
637            self.failUnless(Lsup(v-u_ref)<self.TOL*10.*Lsup(u_ref), "flux error too big.")
638            self.failUnless(Lsup(p-p_ref)<self.TOL*10.*Lsup(p_ref), "pressure error too big.")
639    
640        def testVarioF_FreeBottom_largeK(self):
641            k=1.e10
642            mp=self.getScalarMask(include_bottom=False)
643            mv=self.getVectorMask(include_bottom=True)
644            u_ref,p_ref,f=self.setSolutionFixedBottom(p0=2500,pb=4000.,f0=10.,fb=30.,k=k)
645            p=p_ref*mp
646            u=u_ref*mv
647            df=DarcyFlow(self.dom)
648            df.setValue(g=f,
649                          location_of_fixed_pressure=mp,
650                          location_of_fixed_flux=mv,
651                          permeability=Scalar(k,Function(self.dom)))
652            df.setTolerance(rtol=self.TOL)
653            v,p=df.solve(u,p, max_iter=100, verbose=VERBOSE)
654            self.failUnless(Lsup(v-u_ref)<self.TOL*10.*Lsup(u_ref), "flux error too big.")
655            self.failUnless(Lsup(p-p_ref)<self.TOL*10.*Lsup(p_ref), "pressure error too big.")
656    
657    class Test_Darcy2D(Test_Darcy):
658        TOL=1e-4
659        WIDTH=1.
660        def setUp(self):
661            NE=40  # wrning smaller NE may case a failure for VarioF tests due to discretization errors.
662            self.dom = Rectangle(NE,NE)
663            self.rescaleDomain()
664        def tearDown(self):
665            del self.dom
666    class Test_Darcy3D(Test_Darcy):
667        TOL=1e-4
668        WIDTH=1.
669        def setUp(self):
670            NE=25  # wrning smaller NE may case a failure for VarioF tests due to discretization errors.
671            self.dom = Brick(NE,NE,NE)
672            self.rescaleDomain()
673        def tearDown(self):
674            del self.dom
675    
    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")  
         
        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.")  
   
    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  
        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.")  
676    
677     def test_StokesProblemCartesian_TFQMR_P_small(self):  class Test_Mountains3D(unittest.TestCase):
678         ETA=1.     def setUp(self):
679         P1=1.         NE=16
680           self.TOL=1e-4
681           self.domain=Brick(NE,NE,NE,order=1,useFullElementOrder=True)
682       def tearDown(self):
683           del self.domain
684       def test_periodic(self):
685           EPS=0.01
686    
687         x=self.domain.getX()         x=self.domain.getX()
688         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))
689         mask=whereZero(x[0])    * [1.,1.,1.] \         a0=1
690                +whereZero(x[0]-1)  * [1.,1.,1.] \         a1=1
691                +whereZero(x[1])    * [1.,1.,1.] \         n0=2
692                +whereZero(x[1]-1)  * [1.,1.,1.] \         n1=2
693                +whereZero(x[2])    * [1.,1.,0.] \         n2=0.5
694                +whereZero(x[2]-1)  * [1.,1.,1.]         a2=-(a0*n0+a1*n1)/n2
695                 v[0]=a0*sin(pi*n0*x[0])* cos(pi*n1*x[1])* cos(pi*n2*x[2])
696                 v[1]=a1*cos(pi*n0*x[0])* sin(pi*n1*x[1])* cos(pi*n2*x[2])
697         sp=StokesProblemCartesian(self.domain)         v[2]=a2*cos(pi*n0*x[0])* cos(pi*n1*x[1])* sin(pi*n2*x[2])
698    
699           H_t=Scalar(0.0, Solution(self.domain))
700           mts=Mountains(self.domain,v,eps=EPS,z=1)
701           u,Z=mts.update(u=v,H_t=H_t)
702                
703         sp.initialize(f=F,fixed_u_mask=mask,eta=ETA)         error_int=integrate(Z)
704         u0=(1-x[0])*x[0]*(1-x[1])*x[1]*[0.,0.,1.]         self.failUnless(error_int<self.TOL, "Boundary intergral is too large.")
        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.")  
705    
706     def test_StokesProblemCartesian_TFQMR_P_large(self):  class Test_Mountains2D(unittest.TestCase):
707         ETA=1.     def setUp(self):
708         P1=1000.         NE=16
709           self.TOL=1e-4
710           self.domain=Rectangle(NE,NE,order=1,useFullElementOrder=True)
711       def tearDown(self):
712           del self.domain
713       def test_periodic(self):
714           EPS=0.01
715    
716         x=self.domain.getX()         x=self.domain.getX()
717         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))
718         mask=whereZero(x[0])    * [1.,1.,1.] \         a0=1
719                +whereZero(x[0]-1)  * [1.,1.,1.] \         n0=1
720                +whereZero(x[1])    * [1.,1.,1.] \         n1=0.5
721                +whereZero(x[1]-1)  * [1.,1.,1.] \         a1=-(a0*n0)/n1
722                +whereZero(x[2])    * [1.,1.,0.] \         v[0]=a0*sin(pi*n0*x[0])* cos(pi*n1*x[1])
723                +whereZero(x[2]-1)  * [1.,1.,1.]         v[1]=a1*cos(pi*n0*x[0])* sin(pi*n1*x[1])
724                
725                 H_t=Scalar(0.0, Solution(self.domain))
726         sp=StokesProblemCartesian(self.domain)         mts=Mountains(self.domain,v,eps=EPS,z=1)
727                 u,Z=mts.update(u=v,H_t=H_t)
728         sp.initialize(f=F,fixed_u_mask=mask,eta=ETA)        
729         u0=(1-x[0])*x[0]*(1-x[1])*x[1]*[0.,0.,1.]         error_int=integrate(Z)
730         p0=Scalar(P1,ReducedSolution(self.domain))         self.failUnless(error_int<self.TOL, "Boundary intergral is too large.")
731         u,p=sp.solve(u0,p0,show_details=VERBOSE, verbose=VERBOSE,max_iter=100,solver="TFQMR")        
732          
733         error_v0=Lsup(u[0]-u0[0])  
734         error_v1=Lsup(u[1]-u0[1])  class Test_Rheologies(unittest.TestCase):
735         error_v2=Lsup(u[2]-u0[2])/0.25**2       """
736         zz=P1*x[0]*x[1]*x[2]-p       this is the program used to generate the powerlaw tests:
737         error_p=Lsup(zz-integrate(zz))/P1  
738         # print error_p, error_v0,error_v1,error_v2       TAU_Y=100.
739         self.failUnless(error_p<TOL,"pressure error too large.")       N=10
740         self.failUnless(error_v0<TOL,"0-velocity error too large.")       M=5
741         self.failUnless(error_v1<TOL,"1-velocity error too large.")  
742         self.failUnless(error_v2<TOL,"2-velocity error too large.")       def getE(tau):
743             if tau<=TAU_Y:
744               return 1./(0.5+20*sqrt(tau))
745             else:
746               raise ValueError,"out of range."
747         tau=[ i* (TAU_Y/N) for i in range(N+1)] + [TAU_Y for i in range(M)]
748         e=[ tau[i]/getE(tau[i]) for i in range(N+1)]
749         e+= [ (i+1+j)* (max(e)/N) for j in range(M) ]
750    
751         print tau
752         print e
753         """
754         TOL=1.e-8
755         def test_PowerLaw_Init(self):
756             pl=PowerLaw(numMaterials=3,verbose=VERBOSE)
757    
758             self.failUnlessEqual(pl.getNumMaterials(),3,"num materials is wrong")
759             self.failUnlessEqual(pl.validMaterialId(0),True,"material id 0 not found")
760             self.failUnlessEqual(pl.validMaterialId(1),True,"material id 1 not found")
761             self.failUnlessEqual(pl.validMaterialId(2),True,"material id 2 not found")
762             self.failUnlessEqual(pl.validMaterialId(3),False,"material id 3 not found")
763    
764             self.failUnlessEqual(pl.getFriction(),None,"initial friction wrong.")
765             self.failUnlessEqual(pl.getTauY(),None,"initial tau y wrong.")
766             pl.setDruckerPragerLaw(tau_Y=10,friction=3)
767             self.failUnlessEqual(pl.getFriction(),3,"friction wrong.")
768             self.failUnlessEqual(pl.getTauY(),10,"tau y wrong.")
769    
770             self.failUnlessEqual(pl.getElasticShearModulus(),None,"initial shear modulus is wrong")
771             pl.setElasticShearModulus(1000)
772             self.failUnlessEqual(pl.getElasticShearModulus(),1000.,"shear modulus is wrong")
773    
774             e=pl.getEtaN()
775             self.failUnlessEqual(len(e),3,"initial length of etaN is wrong.")
776             self.failUnlessEqual(e,[None, None, None],"initial etaN are wrong.")
777             self.failUnlessEqual(pl.getEtaN(0),None,"initial material id 0 is wrong")
778             self.failUnlessEqual(pl.getEtaN(1),None,"initial material id 1 is wrong")
779             self.failUnlessEqual(pl.getEtaN(2),None,"initial material id 2 is wrong")
780             self.failUnlessRaises(ValueError, pl.getEtaN, 3)
781    
782             self.failUnlessRaises(ValueError,pl.setPowerLaws,eta_N=[10,20,30,40],tau_t=[2], power=[5,6,7])
783             self.failUnlessRaises(ValueError,pl.setPowerLaws,eta_N=[10,20,30],tau_t=[2,4,5,5], power=[5,6,7])
784             self.failUnlessRaises(ValueError,pl.setPowerLaws,eta_N=[10,20,30],tau_t=[2,1,2], power=[5,6,7,7])
785    
786             pl.setPowerLaws(eta_N=[10,20,30],tau_t=[100,200,300], power=[1,2,3])
787             self.failUnlessEqual(pl.getPower(),[1,2,3],"powers are wrong.")
788             self.failUnlessEqual(pl.getTauT(),[100,200,300],"tau t are wrong.")
789             self.failUnlessEqual(pl.getEtaN(),[10,20,30],"etaN are wrong.")
790    
791             pl.setPowerLaw(id=0,eta_N=40,tau_t=400, power=4)
792             self.failUnlessEqual(pl.getPower(),[4,2,3],"powers are wrong.")
793             self.failUnlessEqual(pl.getTauT(),[400,200,300],"tau t are wrong.")
794             self.failUnlessEqual(pl.getEtaN(),[40,20,30],"etaN are wrong.")
795    
796             self.failUnlessRaises(ValueError,pl.getPower,-1)
797             self.failUnlessEqual(pl.getPower(0),4,"power 0 is wrong.")
798             self.failUnlessEqual(pl.getPower(1),2,"power 1 is wrong.")
799             self.failUnlessEqual(pl.getPower(2),3,"power 2 is wrong.")
800             self.failUnlessRaises(ValueError,pl.getPower,3)
801    
802             self.failUnlessRaises(ValueError,pl.getTauT,-1)
803             self.failUnlessEqual(pl.getTauT(0),400,"tau t 0 is wrong.")
804             self.failUnlessEqual(pl.getTauT(1),200,"tau t 1 is wrong.")
805             self.failUnlessEqual(pl.getTauT(2),300,"tau t 2 is wrong.")
806             self.failUnlessRaises(ValueError,pl.getTauT,3)
807    
808             self.failUnlessRaises(ValueError,pl.getEtaN,-1)
809             self.failUnlessEqual(pl.getEtaN(0),40,"eta n 0 is wrong.")
810             self.failUnlessEqual(pl.getEtaN(1),20,"eta n 1 is wrong.")
811             self.failUnlessEqual(pl.getEtaN(2),30,"eta n 2 is wrong.")
812             self.failUnlessRaises(ValueError,pl.getEtaN,3)
813    
814         def checkResult(self,id,gamma_dot_, eta, tau_ref):
815             self.failUnless(eta>=0,"eta needs to be positive (test %s)"%id)
816             error=abs(gamma_dot_*eta-tau_ref)
817             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))
818            
819         def test_PowerLaw_Linear(self):
820             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]
821             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]
822             pl=PowerLaw(numMaterials=1,verbose=VERBOSE)
823             pl.setDruckerPragerLaw(tau_Y=100.)
824             pl.setPowerLaw(eta_N=2.)
825             pl.setEtaTolerance(self.TOL)
826             for i in xrange(len(taus)): self.checkResult(i,gamma_dot_s[i], pl.getEtaEff(gamma_dot_s[i]),taus[i])
827            
828         def test_PowerLaw_QuadLarge(self):
829             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]
830             gamma_dot_s=[0.0, 405.0, 1610.0, 3615.0, 6420.0, 10025.0, 14430.0, 19635.0, 25640.0, 32445.0, 40050.0, 44055.0, 48060.0, 52065.0, 56070.0, 60075.0]
831             pl=PowerLaw(numMaterials=2,verbose=VERBOSE)
832             pl.setDruckerPragerLaw(tau_Y=100.)
833             pl.setPowerLaws(eta_N=[2.,0.01],tau_t=[1, 25.], power=[1,2])
834             pl.setEtaTolerance(self.TOL)
835             for i in xrange(len(taus)): self.checkResult(i,gamma_dot_s[i], pl.getEtaEff(gamma_dot_s[i]),taus[i])
836    
837         def test_PowerLaw_QuadSmall(self):
838             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]
839             gamma_dot_s=[0.0, 5.4, 11.6, 18.6, 26.4, 35.0, 44.4, 54.6, 65.6, 77.4, 90.0, 99.0, 108.0, 117.0, 126.0, 135.0]
840             pl=PowerLaw(numMaterials=2,verbose=VERBOSE)
841             pl.setDruckerPragerLaw(tau_Y=100.)
842             pl.setPowerLaws(eta_N=[2.,10.],tau_t=[1, 25.], power=[1,2])
843             pl.setEtaTolerance(self.TOL)
844             for i in xrange(len(taus)): self.checkResult(i,gamma_dot_s[i], pl.getEtaEff(gamma_dot_s[i]),taus[i])
845    
846         def test_PowerLaw_CubeLarge(self):
847             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]
848             gamma_dot_s=[0.0, 8.90625, 41.25, 120.46875, 270.0, 513.28125, 873.75, 1374.84375, 2040.0, 2892.65625, 3956.25, 4351.875, 4747.5, 5143.125, 5538.75, 5934.375]
849             pl=PowerLaw(numMaterials=2,verbose=VERBOSE)
850             pl.setDruckerPragerLaw(tau_Y=100.)
851             pl.setPowerLaws(eta_N=[2.,1./16.],tau_t=[1, 64.], power=[1,3])
852             pl.setEtaTolerance(self.TOL)
853             for i in xrange(len(taus)): self.checkResult(i,gamma_dot_s[i], pl.getEtaEff(gamma_dot_s[i]),taus[i])
854    
855         def test_PowerLaw_CubeSmall(self):
856             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]
857             gamma_dot_s=[0.0, 5.0390625, 10.3125, 16.0546875, 22.5, 29.8828125, 38.4375, 48.3984375, 60.0, 73.4765625, 89.0625, 97.96875, 106.875, 115.78125, 124.6875, 133.59375]
858             pl=PowerLaw(numMaterials=2,verbose=VERBOSE)
859             pl.setDruckerPragerLaw(tau_Y=100.)
860             pl.setPowerLaws(eta_N=[2.,25./4.],tau_t=[1, 64.], power=[1,3])
861             pl.setEtaTolerance(self.TOL)
862             for i in xrange(len(taus)): self.checkResult(i,gamma_dot_s[i], pl.getEtaEff(gamma_dot_s[i]),taus[i])
863    
864         def test_PowerLaw_QuadLarge_CubeLarge(self):
865             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]
866             gamma_dot_s=[0.0, 408.90625, 1641.25, 3720.46875, 6670.0, 10513.28125, 15273.75, 20974.84375, 27640.000000000004, 35292.65625, 43956.25, 48351.875, 52747.5, 57143.125, 61538.75, 65934.375]
867    
868             pl=PowerLaw(numMaterials=3,verbose=VERBOSE)
869             pl.setDruckerPragerLaw(tau_Y=100.)
870             pl.setPowerLaws(eta_N=[2.,0.01,1./16.],tau_t=[1, 25.,64.], power=[1,2,3])
871             pl.setEtaTolerance(self.TOL)
872             for i in xrange(len(taus)): self.checkResult(i,gamma_dot_s[i], pl.getEtaEff(gamma_dot_s[i]),taus[i])
873    
874         def test_PowerLaw_QuadLarge_CubeSmall(self):
875             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]
876             gamma_dot_s=[0.0, 405.0390625, 1610.3125, 3616.0546875, 6422.5, 10029.8828125, 14438.4375, 19648.3984375, 25660.0, 32473.4765625, 40089.0625, 44097.96875, 48106.875, 52115.78125, 56124.6875, 60133.59375]
877    
878             pl=PowerLaw(numMaterials=3,verbose=VERBOSE)
879             pl.setDruckerPragerLaw(tau_Y=100.)
880             pl.setPowerLaws(eta_N=[2.,0.01,25./4.],tau_t=[1, 25.,64.], power=[1,2,3])
881             pl.setEtaTolerance(self.TOL)
882             for i in xrange(len(taus)): self.checkResult(i,gamma_dot_s[i], pl.getEtaEff(gamma_dot_s[i]),taus[i])
883    
884         def test_PowerLaw_QuadSmall_CubeLarge(self):
885             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]
886             gamma_dot_s=[0.0, 9.30625, 42.85, 124.06875, 276.4, 523.28125, 888.15, 1394.44375, 2065.6, 2925.05625, 3996.25, 4395.875, 4795.5, 5195.125, 5594.75, 5994.375]
887             pl=PowerLaw(numMaterials=3,verbose=VERBOSE)
888             pl.setDruckerPragerLaw(tau_Y=100.)
889             pl.setPowerLaws(eta_N=[2.,10.,1./16.],tau_t=[1, 25.,64.], power=[1,2,3])
890             pl.setEtaTolerance(self.TOL)
891             for i in xrange(len(taus)): self.checkResult(i,gamma_dot_s[i], pl.getEtaEff(gamma_dot_s[i]),taus[i])
892    
893         def test_PowerLaw_QuadSmall_CubeSmall(self):
894             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]
895             gamma_dot_s=[0.0, 5.4390625, 11.9125, 19.6546875, 28.9, 39.8828125, 52.8375, 67.9984375, 85.6, 105.8765625, 129.0625, 141.96875, 154.875, 167.78125, 180.6875, 193.59375]
896             pl=PowerLaw(numMaterials=3,verbose=VERBOSE)
897             pl.setDruckerPragerLaw(tau_Y=100.)
898             pl.setPowerLaws(eta_N=[2.,10.,25./4.],tau_t=[1, 25.,64.], power=[1,2,3])
899             pl.setEtaTolerance(self.TOL)
900             for i in xrange(len(taus)): self.checkResult(i,gamma_dot_s[i], pl.getEtaEff(gamma_dot_s[i]),taus[i])
901    
902         def test_PowerLaw_withShear(self):
903             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]
904             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]
905             pl=PowerLaw(numMaterials=1,verbose=VERBOSE)
906             pl.setDruckerPragerLaw(tau_Y=100.)
907             pl.setPowerLaw(eta_N=2.)
908             pl.setElasticShearModulus(3.)
909             dt=1./3.
910             pl.setEtaTolerance(self.TOL)
911             self.failUnlessRaises(ValueError, pl.getEtaEff,gamma_dot_s[0])
912             for i in xrange(len(taus)): self.checkResult(i,gamma_dot_s[i], pl.getEtaEff(gamma_dot_s[i],dt=dt),taus[i])
913    
914    class Test_IncompressibleIsotropicFlowCartesian(unittest.TestCase):
915       TOL=1.e-6
916       VERBOSE=False
917       A=1.
918       P_max=100
919       NE=2*getMPISizeWorld()
920       tau_Y=10.
921       N_dt=10
922    
923       # material parameter:
924       tau_1=5.
925       tau_2=5.
926       eta_0=100.
927       eta_1=50.
928       eta_2=400.
929       N_1=2.
930       N_2=3.
931       def getReference(self, t):
932    
933          B=self.tau_Y/sqrt((self.dom.getDim()-1)*self.dom.getDim()*0.5)
934          x=self.dom.getX()
935    
936          s_00=min(self.A*t,B)
937          tau=sqrt((self.dom.getDim()-1)*self.dom.getDim()*0.5)*abs(s_00)
938          inv_eta= 1./self.eta_0 + 1./self.eta_1*(tau/self.tau_1)**(self.N_1-1.) + 1./self.eta_2*(tau/self.tau_2)**(self.N_2-1.)
939    
940          alpha=0.5*inv_eta*s_00
941          if s_00 <= B and self.mu !=None: alpha+=1./(2*self.mu)*self.A
942          u_ref=x*alpha
943          u_ref[self.dom.getDim()-1]=(1.-x[self.dom.getDim()-1])*alpha*(self.dom.getDim()-1)
944          sigma_ref=kronecker(self.dom)*s_00
945          sigma_ref[self.dom.getDim()-1,self.dom.getDim()-1]=-s_00*(self.dom.getDim()-1)
946    
947          p_ref=self.P_max
948          for d in range(self.dom.getDim()): p_ref=p_ref*x[d]
949          p_ref-=integrate(p_ref)/vol(self.dom)
950          return u_ref, sigma_ref, p_ref
951    
952       def runIt(self, free=None):
953          x=self.dom.getX()
954          B=self.tau_Y/sqrt((self.dom.getDim()-1)*self.dom.getDim()*0.5)
955          dt=B/int(self.N_dt/2)
956          if self.VERBOSE: print "dt =",dt
957          if self.latestart:
958              t=dt
959          else:
960              t=0
961          v,s,p=self.getReference(t)
962    
963          mod=IncompressibleIsotropicFlowCartesian(self.dom, stress=s, v=v, p=p, t=t, numMaterials=3, verbose=self.VERBOSE)
964          mod.setDruckerPragerLaw(tau_Y=self.tau_Y,friction=None)
965          mod.setElasticShearModulus(self.mu)
966          mod.setPowerLaws([self.eta_0, self.eta_1, self.eta_2], [ 1., self.tau_1, self.tau_2],  [1.,self.N_1,self.N_2])
967          mod.setTolerance(self.TOL)
968          mod.setFlowTolerance(self.TOL)
969    
970          BF=Vector(self.P_max,Function(self.dom))
971          for d in range(self.dom.getDim()):
972              for d2 in range(self.dom.getDim()):
973                  if d!=d2: BF[d]*=x[d2]
974          v_mask=Vector(0,Solution(self.dom))
975          if free==None:
976             for d in range(self.dom.getDim()):
977                v_mask[d]=whereZero(x[d])+whereZero(x[d]-1.)
978          else:
979             for d in range(self.dom.getDim()):
980                if d == self.dom.getDim()-1:
981                   v_mask[d]=whereZero(x[d]-1.)
982                else:
983                   v_mask[d]=whereZero(x[d])
984          mod.setExternals(F=BF,fixed_v_mask=v_mask)
985          
986          n=self.dom.getNormal()
987          N_t=0
988          errors=[]
989          while N_t < self.N_dt:
990             t_ref=t+dt
991             v_ref, s_ref,p_ref=self.getReference(t_ref)
992             mod.setExternals(f=matrixmult(s_ref,n)-p_ref*n, v_boundary=v_ref)
993             mod.update(dt, iter_max=100, inner_iter_max=20, verbose=self.VERBOSE, usePCG=True)
994             self.check(N_t,mod,t_ref,v_ref, s_ref,p_ref)
995             t+=dt
996             N_t+=1
997    
998       def check(self,N_t,mod,t_ref,v_ref, s_ref,p_ref):
999             p=mod.getPressure()
1000             p-=integrate(p)/vol(self.dom)
1001             error_p=Lsup(mod.getPressure()-p_ref)/Lsup(p_ref)
1002             error_s=Lsup(mod.getDeviatoricStress()-s_ref)/Lsup(s_ref)
1003             error_v=Lsup(mod.getVelocity()-v_ref)/Lsup(v_ref)
1004             error_t=abs(mod.getTime()-t_ref)/abs(t_ref)
1005             if self.VERBOSE: print "time step ",N_t,"time = ",mod.getTime(),"errors s,p,v = ",error_s, error_p, error_v
1006             self.failUnless( error_p <= 10*self.TOL, "time step %s: pressure error %s too high."%(N_t,error_p) )
1007             self.failUnless( error_s <= 10*self.TOL, "time step %s: stress error %s too high."%(N_t,error_s) )
1008             self.failUnless( error_v <= 10*self.TOL, "time step %s: velocity error %s too high."%(N_t,error_v) )
1009             self.failUnless( error_t <= 10*self.TOL, "time step %s: time marker error %s too high."%(N_t,error_t) )
1010       def tearDown(self):
1011            del self.dom
1012    
1013       def test_D2_Fixed_MuNone_LateStart(self):
1014           self.dom = Rectangle(self.NE,self.NE,order=2)
1015           self.mu=None
1016           self.latestart=True
1017           self.runIt()
1018       def test_D2_Fixed_Mu_LateStart(self):
1019           self.dom = Rectangle(self.NE,self.NE,order=2)
1020           self.mu=555.
1021           self.latestart=True
1022           self.runIt()
1023       def test_D2_Fixed_MuNone(self):
1024           self.dom = Rectangle(self.NE,self.NE,order=2)
1025           self.mu=None
1026           self.latestart=False
1027           self.runIt()
1028       def test_D2_Fixed_Mu(self):
1029           self.dom = Rectangle(self.NE,self.NE,order=2)
1030           self.mu=555.
1031           self.latestart=False
1032           self.runIt()
1033       def test_D2_Free_MuNone_LateStart(self):
1034           self.dom = Rectangle(self.NE,self.NE,order=2)
1035           self.mu=None
1036           self.latestart=True
1037           self.runIt(free=0)
1038       def test_D2_Free_Mu_LateStart(self):
1039           self.dom = Rectangle(self.NE,self.NE,order=2)
1040           self.mu=555.
1041           self.latestart=True
1042           self.runIt(free=0)
1043       def test_D2_Free_MuNone(self):
1044           self.dom = Rectangle(self.NE,self.NE,order=2)
1045           self.mu=None
1046           self.latestart=False
1047           self.runIt(free=0)
1048       def test_D2_Free_Mu(self):
1049           self.dom = Rectangle(self.NE,self.NE,order=2)
1050           self.mu=555.
1051           self.latestart=False
1052           self.runIt(free=0)
1053    
1054       def test_D3_Fixed_MuNone_LateStart(self):
1055           self.dom = Brick(self.NE,self.NE,self.NE,order=2)
1056           self.mu=None
1057           self.latestart=True
1058           self.runIt()
1059       def test_D3_Fixed_Mu_LateStart(self):
1060           self.dom = Brick(self.NE,self.NE,self.NE,order=2)
1061           self.mu=555.
1062           self.latestart=True
1063           self.runIt()
1064       def test_D3_Fixed_MuNone(self):
1065           self.dom = Brick(self.NE,self.NE,self.NE,order=2)
1066           self.mu=None
1067           self.latestart=False
1068           self.runIt()
1069       def test_D3_Fixed_Mu(self):
1070           self.dom = Brick(self.NE,self.NE,self.NE,order=2)
1071           self.mu=555.
1072           self.latestart=False
1073           self.runIt()
1074       def test_D3_Free_MuNone_LateStart(self):
1075           self.dom = Brick(self.NE,self.NE,self.NE,order=2)
1076           self.mu=None
1077           self.latestart=True
1078           self.runIt(free=0)
1079       def test_D3_Free_Mu_LateStart(self):
1080           self.dom = Brick(self.NE,self.NE,self.NE,order=2)
1081           self.mu=555.
1082           self.latestart=True
1083           self.runIt(free=0)
1084       def test_D3_Free_MuNone(self):
1085           self.dom = Brick(self.NE,self.NE,self.NE,order=2)
1086           self.mu=None
1087           self.latestart=False
1088           self.runIt(free=0)
1089       def test_D3_Free_Mu(self):
1090           self.dom = Brick(self.NE,self.NE,self.NE,order=2)
1091           self.mu=555.
1092           self.latestart=False
1093           self.runIt(free=0)
1094    
1095  if __name__ == '__main__':  if __name__ == '__main__':
1096     suite = unittest.TestSuite()     suite = unittest.TestSuite()
1097     suite.addTest(unittest.makeSuite(Test_Simple2DModels))     suite.addTest(unittest.makeSuite(Test_StokesProblemCartesian2D))
1098     # suite.addTest(unittest.makeSuite(Test_Simple3DModels))     suite.addTest(unittest.makeSuite(Test_Darcy3D))
1099       suite.addTest(unittest.makeSuite(Test_Darcy2D))
1100       suite.addTest(unittest.makeSuite(Test_StokesProblemCartesian3D))
1101       suite.addTest(unittest.makeSuite(Test_Mountains3D))
1102       suite.addTest(unittest.makeSuite(Test_Mountains2D))
1103       suite.addTest(unittest.makeSuite(Test_Rheologies))
1104       suite.addTest(unittest.makeSuite(Test_IncompressibleIsotropicFlowCartesian))
1105     s=unittest.TextTestRunner(verbosity=2).run(suite)     s=unittest.TextTestRunner(verbosity=2).run(suite)
1106     if not s.wasSuccessful(): sys.exit(1)     if not s.wasSuccessful(): sys.exit(1)
1107    

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

  ViewVC Help
Powered by ViewVC 1.1.26