/[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 2563 by gross, Tue Jul 28 03:50:45 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-2009 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           mts=Mountains(self.domain,eps=EPS)
700           mts.setVelocity(v)
701           Z=mts.update()
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,eps=EPS)
727                 mts.setVelocity(v)
728         sp.initialize(f=F,fixed_u_mask=mask,eta=ETA)         Z=mts.update()
729         u0=(1-x[0])*x[0]*(1-x[1])*x[1]*[0.,0.,1.]        
730         p0=Scalar(P1,ReducedSolution(self.domain))         error_int=integrate(Z)
731         u,p=sp.solve(u0,p0,show_details=VERBOSE, verbose=VERBOSE,max_iter=100,solver="TFQMR")         self.failUnless(error_int<self.TOL, "Boundary intergral is too large.")
732                
733         error_v0=Lsup(u[0]-u0[0])  
734         error_v1=Lsup(u[1]-u0[1])  
735         error_v2=Lsup(u[2]-u0[2])/0.25**2  class Test_Rheologies(unittest.TestCase):
736         zz=P1*x[0]*x[1]*x[2]-p       """
737         error_p=Lsup(zz-integrate(zz))/P1       this is the program used to generate the powerlaw tests:
738         # print error_p, error_v0,error_v1,error_v2  
739         self.failUnless(error_p<TOL,"pressure error too large.")       TAU_Y=100.
740         self.failUnless(error_v0<TOL,"0-velocity error too large.")       N=10
741         self.failUnless(error_v1<TOL,"1-velocity error too large.")       M=5
742         self.failUnless(error_v2<TOL,"2-velocity error too large.")  
743         def getE(tau):
744             if tau<=TAU_Y:
745               return 1./(0.5+20*sqrt(tau))
746             else:
747               raise ValueError,"out of range."
748         tau=[ i* (TAU_Y/N) for i in range(N+1)] + [TAU_Y for i in range(M)]
749         e=[ tau[i]/getE(tau[i]) for i in range(N+1)]
750         e+= [ (i+1+j)* (max(e)/N) for j in range(M) ]
751    
752         print tau
753         print e
754         """
755         TOL=1.e-8
756         def test_PowerLaw_Init(self):
757             pl=PowerLaw(numMaterials=3,verbose=VERBOSE)
758    
759             self.failUnlessEqual(pl.getNumMaterials(),3,"num materials is wrong")
760             self.failUnlessEqual(pl.validMaterialId(0),True,"material id 0 not found")
761             self.failUnlessEqual(pl.validMaterialId(1),True,"material id 1 not found")
762             self.failUnlessEqual(pl.validMaterialId(2),True,"material id 2 not found")
763             self.failUnlessEqual(pl.validMaterialId(3),False,"material id 3 not found")
764    
765             self.failUnlessEqual(pl.getFriction(),None,"initial friction wrong.")
766             self.failUnlessEqual(pl.getTauY(),None,"initial tau y wrong.")
767             pl.setDruckerPragerLaw(tau_Y=10,friction=3)
768             self.failUnlessEqual(pl.getFriction(),3,"friction wrong.")
769             self.failUnlessEqual(pl.getTauY(),10,"tau y wrong.")
770    
771             self.failUnlessEqual(pl.getElasticShearModulus(),None,"initial shear modulus is wrong")
772             pl.setElasticShearModulus(1000)
773             self.failUnlessEqual(pl.getElasticShearModulus(),1000.,"shear modulus is wrong")
774    
775             e=pl.getEtaN()
776             self.failUnlessEqual(len(e),3,"initial length of etaN is wrong.")
777             self.failUnlessEqual(e,[None, None, None],"initial etaN are wrong.")
778             self.failUnlessEqual(pl.getEtaN(0),None,"initial material id 0 is wrong")
779             self.failUnlessEqual(pl.getEtaN(1),None,"initial material id 1 is wrong")
780             self.failUnlessEqual(pl.getEtaN(2),None,"initial material id 2 is wrong")
781             self.failUnlessRaises(ValueError, pl.getEtaN, 3)
782    
783             self.failUnlessRaises(ValueError,pl.setPowerLaws,eta_N=[10,20,30,40],tau_t=[2], power=[5,6,7])
784             self.failUnlessRaises(ValueError,pl.setPowerLaws,eta_N=[10,20,30],tau_t=[2,4,5,5], power=[5,6,7])
785             self.failUnlessRaises(ValueError,pl.setPowerLaws,eta_N=[10,20,30],tau_t=[2,1,2], power=[5,6,7,7])
786    
787             pl.setPowerLaws(eta_N=[10,20,30],tau_t=[100,200,300], power=[1,2,3])
788             self.failUnlessEqual(pl.getPower(),[1,2,3],"powers are wrong.")
789             self.failUnlessEqual(pl.getTauT(),[100,200,300],"tau t are wrong.")
790             self.failUnlessEqual(pl.getEtaN(),[10,20,30],"etaN are wrong.")
791    
792             pl.setPowerLaw(id=0,eta_N=40,tau_t=400, power=4)
793             self.failUnlessEqual(pl.getPower(),[4,2,3],"powers are wrong.")
794             self.failUnlessEqual(pl.getTauT(),[400,200,300],"tau t are wrong.")
795             self.failUnlessEqual(pl.getEtaN(),[40,20,30],"etaN are wrong.")
796    
797             self.failUnlessRaises(ValueError,pl.getPower,-1)
798             self.failUnlessEqual(pl.getPower(0),4,"power 0 is wrong.")
799             self.failUnlessEqual(pl.getPower(1),2,"power 1 is wrong.")
800             self.failUnlessEqual(pl.getPower(2),3,"power 2 is wrong.")
801             self.failUnlessRaises(ValueError,pl.getPower,3)
802    
803             self.failUnlessRaises(ValueError,pl.getTauT,-1)
804             self.failUnlessEqual(pl.getTauT(0),400,"tau t 0 is wrong.")
805             self.failUnlessEqual(pl.getTauT(1),200,"tau t 1 is wrong.")
806             self.failUnlessEqual(pl.getTauT(2),300,"tau t 2 is wrong.")
807             self.failUnlessRaises(ValueError,pl.getTauT,3)
808    
809             self.failUnlessRaises(ValueError,pl.getEtaN,-1)
810             self.failUnlessEqual(pl.getEtaN(0),40,"eta n 0 is wrong.")
811             self.failUnlessEqual(pl.getEtaN(1),20,"eta n 1 is wrong.")
812             self.failUnlessEqual(pl.getEtaN(2),30,"eta n 2 is wrong.")
813             self.failUnlessRaises(ValueError,pl.getEtaN,3)
814    
815         def checkResult(self,id,gamma_dot_, eta, tau_ref):
816             self.failUnless(eta>=0,"eta needs to be positive (test %s)"%id)
817             error=abs(gamma_dot_*eta-tau_ref)
818             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))
819            
820         def test_PowerLaw_Linear(self):
821             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]
822             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]
823             pl=PowerLaw(numMaterials=1,verbose=VERBOSE)
824             pl.setDruckerPragerLaw(tau_Y=100.)
825             pl.setPowerLaw(eta_N=2.)
826             pl.setEtaTolerance(self.TOL)
827             for i in xrange(len(taus)): self.checkResult(i,gamma_dot_s[i], pl.getEtaEff(gamma_dot_s[i]),taus[i])
828            
829         def test_PowerLaw_QuadLarge(self):
830             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]
831             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]
832             pl=PowerLaw(numMaterials=2,verbose=VERBOSE)
833             pl.setDruckerPragerLaw(tau_Y=100.)
834             pl.setPowerLaws(eta_N=[2.,0.01],tau_t=[1, 25.], power=[1,2])
835             pl.setEtaTolerance(self.TOL)
836             for i in xrange(len(taus)): self.checkResult(i,gamma_dot_s[i], pl.getEtaEff(gamma_dot_s[i]),taus[i])
837    
838         def test_PowerLaw_QuadSmall(self):
839             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]
840             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]
841             pl=PowerLaw(numMaterials=2,verbose=VERBOSE)
842             pl.setDruckerPragerLaw(tau_Y=100.)
843             pl.setPowerLaws(eta_N=[2.,10.],tau_t=[1, 25.], power=[1,2])
844             pl.setEtaTolerance(self.TOL)
845             for i in xrange(len(taus)): self.checkResult(i,gamma_dot_s[i], pl.getEtaEff(gamma_dot_s[i]),taus[i])
846    
847         def test_PowerLaw_CubeLarge(self):
848             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]
849             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]
850             pl=PowerLaw(numMaterials=2,verbose=VERBOSE)
851             pl.setDruckerPragerLaw(tau_Y=100.)
852             pl.setPowerLaws(eta_N=[2.,1./16.],tau_t=[1, 64.], power=[1,3])
853             pl.setEtaTolerance(self.TOL)
854             for i in xrange(len(taus)): self.checkResult(i,gamma_dot_s[i], pl.getEtaEff(gamma_dot_s[i]),taus[i])
855    
856         def test_PowerLaw_CubeSmall(self):
857             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]
858             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]
859             pl=PowerLaw(numMaterials=2,verbose=VERBOSE)
860             pl.setDruckerPragerLaw(tau_Y=100.)
861             pl.setPowerLaws(eta_N=[2.,25./4.],tau_t=[1, 64.], power=[1,3])
862             pl.setEtaTolerance(self.TOL)
863             for i in xrange(len(taus)): self.checkResult(i,gamma_dot_s[i], pl.getEtaEff(gamma_dot_s[i]),taus[i])
864    
865         def test_PowerLaw_QuadLarge_CubeLarge(self):
866             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]
867             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]
868    
869             pl=PowerLaw(numMaterials=3,verbose=VERBOSE)
870             pl.setDruckerPragerLaw(tau_Y=100.)
871             pl.setPowerLaws(eta_N=[2.,0.01,1./16.],tau_t=[1, 25.,64.], power=[1,2,3])
872             pl.setEtaTolerance(self.TOL)
873             for i in xrange(len(taus)): self.checkResult(i,gamma_dot_s[i], pl.getEtaEff(gamma_dot_s[i]),taus[i])
874    
875         def test_PowerLaw_QuadLarge_CubeSmall(self):
876             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]
877             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]
878    
879             pl=PowerLaw(numMaterials=3,verbose=VERBOSE)
880             pl.setDruckerPragerLaw(tau_Y=100.)
881             pl.setPowerLaws(eta_N=[2.,0.01,25./4.],tau_t=[1, 25.,64.], power=[1,2,3])
882             pl.setEtaTolerance(self.TOL)
883             for i in xrange(len(taus)): self.checkResult(i,gamma_dot_s[i], pl.getEtaEff(gamma_dot_s[i]),taus[i])
884    
885         def test_PowerLaw_QuadSmall_CubeLarge(self):
886             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]
887             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]
888             pl=PowerLaw(numMaterials=3,verbose=VERBOSE)
889             pl.setDruckerPragerLaw(tau_Y=100.)
890             pl.setPowerLaws(eta_N=[2.,10.,1./16.],tau_t=[1, 25.,64.], power=[1,2,3])
891             pl.setEtaTolerance(self.TOL)
892             for i in xrange(len(taus)): self.checkResult(i,gamma_dot_s[i], pl.getEtaEff(gamma_dot_s[i]),taus[i])
893    
894         def test_PowerLaw_QuadSmall_CubeSmall(self):
895             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]
896             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]
897             pl=PowerLaw(numMaterials=3,verbose=VERBOSE)
898             pl.setDruckerPragerLaw(tau_Y=100.)
899             pl.setPowerLaws(eta_N=[2.,10.,25./4.],tau_t=[1, 25.,64.], power=[1,2,3])
900             pl.setEtaTolerance(self.TOL)
901             for i in xrange(len(taus)): self.checkResult(i,gamma_dot_s[i], pl.getEtaEff(gamma_dot_s[i]),taus[i])
902    
903         def test_PowerLaw_withShear(self):
904             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]
905             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]
906             pl=PowerLaw(numMaterials=1,verbose=VERBOSE)
907             pl.setDruckerPragerLaw(tau_Y=100.)
908             pl.setPowerLaw(eta_N=2.)
909             pl.setElasticShearModulus(3.)
910             dt=1./3.
911             pl.setEtaTolerance(self.TOL)
912             self.failUnlessRaises(ValueError, pl.getEtaEff,gamma_dot_s[0])
913             for i in xrange(len(taus)): self.checkResult(i,gamma_dot_s[i], pl.getEtaEff(gamma_dot_s[i],dt=dt),taus[i])
914    
915    class Test_IncompressibleIsotropicFlowCartesian(unittest.TestCase):
916       TOL=1.e-6
917       VERBOSE=False
918       A=1.
919       P_max=100
920       NE=2*getMPISizeWorld()
921       tau_Y=10.
922       N_dt=10
923    
924       # material parameter:
925       tau_1=5.
926       tau_2=5.
927       eta_0=100.
928       eta_1=50.
929       eta_2=400.
930       N_1=2.
931       N_2=3.
932       def getReference(self, t):
933    
934          B=self.tau_Y/sqrt((self.dom.getDim()-1)*self.dom.getDim()*0.5)
935          x=self.dom.getX()
936    
937          s_00=min(self.A*t,B)
938          tau=sqrt((self.dom.getDim()-1)*self.dom.getDim()*0.5)*abs(s_00)
939          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.)
940    
941          alpha=0.5*inv_eta*s_00
942          if s_00 <= B and self.mu !=None: alpha+=1./(2*self.mu)*self.A
943          u_ref=x*alpha
944          u_ref[self.dom.getDim()-1]=(1.-x[self.dom.getDim()-1])*alpha*(self.dom.getDim()-1)
945          sigma_ref=kronecker(self.dom)*s_00
946          sigma_ref[self.dom.getDim()-1,self.dom.getDim()-1]=-s_00*(self.dom.getDim()-1)
947    
948          p_ref=self.P_max
949          for d in range(self.dom.getDim()): p_ref=p_ref*x[d]
950          p_ref-=integrate(p_ref)/vol(self.dom)
951          return u_ref, sigma_ref, p_ref
952    
953       def runIt(self, free=None):
954          x=self.dom.getX()
955          B=self.tau_Y/sqrt((self.dom.getDim()-1)*self.dom.getDim()*0.5)
956          dt=B/int(self.N_dt/2)
957          if self.VERBOSE: print "dt =",dt
958          if self.latestart:
959              t=dt
960          else:
961              t=0
962          v,s,p=self.getReference(t)
963    
964          mod=IncompressibleIsotropicFlowCartesian(self.dom, stress=s, v=v, p=p, t=t, numMaterials=3, verbose=self.VERBOSE)
965          mod.setDruckerPragerLaw(tau_Y=self.tau_Y,friction=None)
966          mod.setElasticShearModulus(self.mu)
967          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])
968          mod.setTolerance(self.TOL)
969          mod.setFlowTolerance(self.TOL)
970    
971          BF=Vector(self.P_max,Function(self.dom))
972          for d in range(self.dom.getDim()):
973              for d2 in range(self.dom.getDim()):
974                  if d!=d2: BF[d]*=x[d2]
975          v_mask=Vector(0,Solution(self.dom))
976          if free==None:
977             for d in range(self.dom.getDim()):
978                v_mask[d]=whereZero(x[d])+whereZero(x[d]-1.)
979          else:
980             for d in range(self.dom.getDim()):
981                if d == self.dom.getDim()-1:
982                   v_mask[d]=whereZero(x[d]-1.)
983                else:
984                   v_mask[d]=whereZero(x[d])
985          mod.setExternals(F=BF,fixed_v_mask=v_mask)
986          
987          n=self.dom.getNormal()
988          N_t=0
989          errors=[]
990          while N_t < self.N_dt:
991             t_ref=t+dt
992             v_ref, s_ref,p_ref=self.getReference(t_ref)
993             mod.setExternals(f=matrixmult(s_ref,n)-p_ref*n, v_boundary=v_ref)
994             mod.update(dt, iter_max=100, inner_iter_max=20, verbose=self.VERBOSE, usePCG=True)
995             self.check(N_t,mod,t_ref,v_ref, s_ref,p_ref)
996             t+=dt
997             N_t+=1
998    
999       def check(self,N_t,mod,t_ref,v_ref, s_ref,p_ref):
1000             p=mod.getPressure()
1001             p-=integrate(p)/vol(self.dom)
1002             error_p=Lsup(mod.getPressure()-p_ref)/Lsup(p_ref)
1003             error_s=Lsup(mod.getDeviatoricStress()-s_ref)/Lsup(s_ref)
1004             error_v=Lsup(mod.getVelocity()-v_ref)/Lsup(v_ref)
1005             error_t=abs(mod.getTime()-t_ref)/abs(t_ref)
1006             if self.VERBOSE: print "time step ",N_t,"time = ",mod.getTime(),"errors s,p,v = ",error_s, error_p, error_v
1007             self.failUnless( error_p <= 10*self.TOL, "time step %s: pressure error %s too high."%(N_t,error_p) )
1008             self.failUnless( error_s <= 10*self.TOL, "time step %s: stress error %s too high."%(N_t,error_s) )
1009             self.failUnless( error_v <= 10*self.TOL, "time step %s: velocity error %s too high."%(N_t,error_v) )
1010             self.failUnless( error_t <= 10*self.TOL, "time step %s: time marker error %s too high."%(N_t,error_t) )
1011       def tearDown(self):
1012            del self.dom
1013    
1014       def test_D2_Fixed_MuNone_LateStart(self):
1015           self.dom = Rectangle(self.NE,self.NE,order=2)
1016           self.mu=None
1017           self.latestart=True
1018           self.runIt()
1019       def test_D2_Fixed_Mu_LateStart(self):
1020           self.dom = Rectangle(self.NE,self.NE,order=2)
1021           self.mu=555.
1022           self.latestart=True
1023           self.runIt()
1024       def test_D2_Fixed_MuNone(self):
1025           self.dom = Rectangle(self.NE,self.NE,order=2)
1026           self.mu=None
1027           self.latestart=False
1028           self.runIt()
1029       def test_D2_Fixed_Mu(self):
1030           self.dom = Rectangle(self.NE,self.NE,order=2)
1031           self.mu=555.
1032           self.latestart=False
1033           self.runIt()
1034       def test_D2_Free_MuNone_LateStart(self):
1035           self.dom = Rectangle(self.NE,self.NE,order=2)
1036           self.mu=None
1037           self.latestart=True
1038           self.runIt(free=0)
1039       def test_D2_Free_Mu_LateStart(self):
1040           self.dom = Rectangle(self.NE,self.NE,order=2)
1041           self.mu=555.
1042           self.latestart=True
1043           self.runIt(free=0)
1044       def test_D2_Free_MuNone(self):
1045           self.dom = Rectangle(self.NE,self.NE,order=2)
1046           self.mu=None
1047           self.latestart=False
1048           self.runIt(free=0)
1049       def test_D2_Free_Mu(self):
1050           self.dom = Rectangle(self.NE,self.NE,order=2)
1051           self.mu=555.
1052           self.latestart=False
1053           self.runIt(free=0)
1054    
1055       def test_D3_Fixed_MuNone_LateStart(self):
1056           self.dom = Brick(self.NE,self.NE,self.NE,order=2)
1057           self.mu=None
1058           self.latestart=True
1059           self.runIt()
1060       def test_D3_Fixed_Mu_LateStart(self):
1061           self.dom = Brick(self.NE,self.NE,self.NE,order=2)
1062           self.mu=555.
1063           self.latestart=True
1064           self.runIt()
1065       def test_D3_Fixed_MuNone(self):
1066           self.dom = Brick(self.NE,self.NE,self.NE,order=2)
1067           self.mu=None
1068           self.latestart=False
1069           self.runIt()
1070       def test_D3_Fixed_Mu(self):
1071           self.dom = Brick(self.NE,self.NE,self.NE,order=2)
1072           self.mu=555.
1073           self.latestart=False
1074           self.runIt()
1075       def test_D3_Free_MuNone_LateStart(self):
1076           self.dom = Brick(self.NE,self.NE,self.NE,order=2)
1077           self.mu=None
1078           self.latestart=True
1079           self.runIt(free=0)
1080       def test_D3_Free_Mu_LateStart(self):
1081           self.dom = Brick(self.NE,self.NE,self.NE,order=2)
1082           self.mu=555.
1083           self.latestart=True
1084           self.runIt(free=0)
1085       def test_D3_Free_MuNone(self):
1086           self.dom = Brick(self.NE,self.NE,self.NE,order=2)
1087           self.mu=None
1088           self.latestart=False
1089           self.runIt(free=0)
1090       def test_D3_Free_Mu(self):
1091           self.dom = Brick(self.NE,self.NE,self.NE,order=2)
1092           self.mu=555.
1093           self.latestart=False
1094           self.runIt(free=0)
1095    
1096  if __name__ == '__main__':  if __name__ == '__main__':
1097     suite = unittest.TestSuite()     suite = unittest.TestSuite()
1098     suite.addTest(unittest.makeSuite(Test_Simple2DModels))     suite.addTest(unittest.makeSuite(Test_StokesProblemCartesian2D))
1099     # suite.addTest(unittest.makeSuite(Test_Simple3DModels))     suite.addTest(unittest.makeSuite(Test_Darcy3D))
1100       suite.addTest(unittest.makeSuite(Test_Darcy2D))
1101       suite.addTest(unittest.makeSuite(Test_StokesProblemCartesian3D))
1102       suite.addTest(unittest.makeSuite(Test_Mountains3D))
1103       suite.addTest(unittest.makeSuite(Test_Mountains2D))
1104       suite.addTest(unittest.makeSuite(Test_Rheologies))
1105       suite.addTest(unittest.makeSuite(Test_IncompressibleIsotropicFlowCartesian))
1106     s=unittest.TextTestRunner(verbosity=2).run(suite)     s=unittest.TextTestRunner(verbosity=2).run(suite)
1107     if not s.wasSuccessful(): sys.exit(1)     if not s.wasSuccessful(): sys.exit(1)
1108    

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

  ViewVC Help
Powered by ViewVC 1.1.26