/[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 2675 by gross, Mon Sep 21 02:32:17 2009 UTC revision 3527 by gross, Tue May 24 06:59:20 2011 UTC
# Line 1  Line 1 
1    # -*- coding: utf-8 -*-
2    
3  ########################################################  ########################################################
4  #  #
5  # Copyright (c) 2003-2009 by University of Queensland  # Copyright (c) 2003-2010 by University of Queensland
6  # Earth Systems Science Computational Center (ESSCC)  # Earth Systems Science Computational Center (ESSCC)
7  # http://www.uq.edu.au/esscc  # http://www.uq.edu.au/esscc
8  #  #
# Line 11  Line 12 
12  #  #
13  ########################################################  ########################################################
14    
15  __copyright__="""Copyright (c) 2003-2009 by University of Queensland  __copyright__="""Copyright (c) 2003-2010 by University of Queensland
16  Earth Systems Science Computational Center (ESSCC)  Earth Systems Science Computational Center (ESSCC)
17  http://www.uq.edu.au/esscc  http://www.uq.edu.au/esscc
18  Primary Business: Queensland, Australia"""  Primary Business: Queensland, Australia"""
# Line 24  import tempfile Line 25  import tempfile
25                
26    
27    
28  VERBOSE=False # or True  VERBOSE=False  and True
29    
30  from esys.escript import *  from esys.escript import *
31  from esys.escript.models import StokesProblemCartesian, PowerLaw, IncompressibleIsotropicFlowCartesian, FaultSystem, DarcyFlow  from esys.escript.models import StokesProblemCartesian, PowerLaw, IncompressibleIsotropicFlowCartesian, FaultSystem, DarcyFlow
# Line 46  class Test_StokesProblemCartesian2D(unit Line 47  class Test_StokesProblemCartesian2D(unit
47     def setUp(self):     def setUp(self):
48         NE=6         NE=6
49         self.TOL=1e-3         self.TOL=1e-3
50         self.domain=Rectangle(NE,NE,order=2,useFullElementOrder=True)         self.domain=Rectangle(NE,NE,order=-1)
51     def tearDown(self):     def tearDown(self):
52         del self.domain         del self.domain
53     def test_PCG_P_0(self):     def test_PCG_P_0(self):
# Line 66  class Test_StokesProblemCartesian2D(unit Line 67  class Test_StokesProblemCartesian2D(unit
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         sp.setTolerance(self.TOL)         sp.setTolerance(self.TOL)
70         u,p=sp.solve(u0,p0,verbose=VERBOSE,max_iter=100,usePCG=True)         u,p=sp.solve(u0*mask,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
# Line 91  class Test_StokesProblemCartesian2D(unit Line 92  class Test_StokesProblemCartesian2D(unit
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         sp.setTolerance(self.TOL*0.2)         sp.setTolerance(self.TOL)
96         u,p=sp.solve(u0,p0, verbose=VERBOSE,max_iter=100,usePCG=True)         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
# Line 118  class Test_StokesProblemCartesian2D(unit Line 119  class Test_StokesProblemCartesian2D(unit
119         p0=Scalar(-P1,ReducedSolution(self.domain))         p0=Scalar(-P1,ReducedSolution(self.domain))
120         sp.setTolerance(self.TOL)         sp.setTolerance(self.TOL)
121         u,p=sp.solve(u0,p0, verbose=VERBOSE,max_iter=100,usePCG=True)         u,p=sp.solve(u0,p0, verbose=VERBOSE,max_iter=100,usePCG=True)
122           # u,p=sp.solve(u0,p0, verbose=VERBOSE,max_iter=100,usePCG=True)
123                
124         error_v0=Lsup(u[0]-u0[0])         error_v0=Lsup(u[0]-u0[0])
125         error_v1=Lsup(u[1]-u0[1])/0.25         error_v1=Lsup(u[1]-u0[1])/0.25
# Line 168  class Test_StokesProblemCartesian2D(unit Line 170  class Test_StokesProblemCartesian2D(unit
170         sp.initialize(f=F,fixed_u_mask=mask,eta=ETA)         sp.initialize(f=F,fixed_u_mask=mask,eta=ETA)
171         u0=(1-x[0])*x[0]*[0.,1.]         u0=(1-x[0])*x[0]*[0.,1.]
172         p0=Scalar(-P1,ReducedSolution(self.domain))         p0=Scalar(-P1,ReducedSolution(self.domain))
173         sp.setTolerance(self.TOL*0.1)         sp.setTolerance(self.TOL)
174         u,p=sp.solve(u0,p0, verbose=VERBOSE,max_iter=20,usePCG=False)         u,p=sp.solve(u0,p0, verbose=VERBOSE,max_iter=20,usePCG=False)
175                
176         error_v0=Lsup(u[0]-u0[0])         error_v0=Lsup(u[0]-u0[0])
177         error_v1=Lsup(u[1]-u0[1])/0.25         error_v1=Lsup(u[1]-u0[1])/0.25
178         error_p=Lsup(P1*x[0]*x[1]+p)         error_p=Lsup(P1*x[0]*x[1]+p)
179    
180         self.failUnless(error_p<10*self.TOL, "pressure error too large.")         self.failUnless(error_p<10*self.TOL, "pressure error too large.")
181         self.failUnless(error_v0<10*self.TOL, "0-velocity error too large.")         self.failUnless(error_v0<10*self.TOL, "0-velocity error too large.")
182         self.failUnless(error_v1<10*self.TOL, "1-velocity error too large.")         self.failUnless(error_v1<10*self.TOL, "1-velocity error too large.")
# Line 208  class Test_StokesProblemCartesian3D(unit Line 211  class Test_StokesProblemCartesian3D(unit
211     def setUp(self):     def setUp(self):
212         NE=6         NE=6
213         self.TOL=1e-4         self.TOL=1e-4
214         self.domain=Brick(NE,NE,NE,order=2,useFullElementOrder=True)         self.domain=Brick(NE,NE,NE,order=-1)
215     def tearDown(self):     def tearDown(self):
216         del self.domain         del self.domain
217     def test_PCG_P_0(self):     def test_PCG_P_0(self):
# Line 262  class Test_StokesProblemCartesian3D(unit Line 265  class Test_StokesProblemCartesian3D(unit
265         sp.initialize(f=F,fixed_u_mask=mask,eta=ETA)         sp.initialize(f=F,fixed_u_mask=mask,eta=ETA)
266         u0=(1-x[0])*x[0]*(1-x[1])*x[1]*[0.,0.,1.]         u0=(1-x[0])*x[0]*(1-x[1])*x[1]*[0.,0.,1.]
267         p0=Scalar(-P1,ReducedSolution(self.domain))         p0=Scalar(-P1,ReducedSolution(self.domain))
268         sp.setTolerance(self.TOL*0.1)         sp.setTolerance(self.TOL)
269         u,p=sp.solve(u0,p0, verbose=VERBOSE ,max_iter=100,usePCG=True)         u,p=sp.solve(u0,p0, verbose=VERBOSE ,max_iter=100,usePCG=True)
270         error_v0=Lsup(u[0]-u0[0])         error_v0=Lsup(u[0]-u0[0])
271         error_v1=Lsup(u[1]-u0[1])         error_v1=Lsup(u[1]-u0[1])
# Line 293  class Test_StokesProblemCartesian3D(unit Line 296  class Test_StokesProblemCartesian3D(unit
296         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.]
297         p0=Scalar(-P1,ReducedSolution(self.domain))         p0=Scalar(-P1,ReducedSolution(self.domain))
298         sp.setTolerance(self.TOL)         sp.setTolerance(self.TOL)
299         u,p=sp.solve(u0,-p0, verbose=VERBOSE ,max_iter=100,usePCG=True)         u,p=sp.solve(u0,p0, verbose=VERBOSE ,max_iter=100,usePCG=True)
300                
301         error_v0=Lsup(u[0]-u0[0])         error_v0=Lsup(u[0]-u0[0])
302         error_v1=Lsup(u[1]-u0[1])         error_v1=Lsup(u[1]-u0[1])
# Line 354  class Test_StokesProblemCartesian3D(unit Line 357  class Test_StokesProblemCartesian3D(unit
357         sp.initialize(f=F,fixed_u_mask=mask,eta=ETA)         sp.initialize(f=F,fixed_u_mask=mask,eta=ETA)
358         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.]
359         p0=Scalar(-P1,ReducedSolution(self.domain))         p0=Scalar(-P1,ReducedSolution(self.domain))
360         sp.setTolerance(self.TOL*0.1)         sp.setTolerance(self.TOL/10)
361         u,p=sp.solve(u0,p0, verbose=VERBOSE,max_iter=100,usePCG=False)         u,p=sp.solve(u0,p0, verbose=VERBOSE,max_iter=100,usePCG=False)
362                
363         error_v0=Lsup(u[0]-u0[0])         error_v0=Lsup(u[0]-u0[0])
# Line 396  class Test_StokesProblemCartesian3D(unit Line 399  class Test_StokesProblemCartesian3D(unit
399         self.failUnless(error_v1<10*self.TOL, "1-velocity error too large.")         self.failUnless(error_v1<10*self.TOL, "1-velocity error too large.")
400         self.failUnless(error_v2<10*self.TOL, "2-velocity error too large.")         self.failUnless(error_v2<10*self.TOL, "2-velocity error too large.")
401  #====================================================================================================================  #====================================================================================================================
 class Test_Darcy(unittest.TestCase):  
     # this is a simple test for the darcy flux problem  
     #  
     #  
     #  p = 1/k * ( 1/2* (fb-f0)/xb* x **2 + f0 * x - ub*x ) +  p0  
     #  
     #  with f = (fb-f0)/xb* x + f0  
     #  
     #    u = f - k * p,x = ub  
     #  
     #  we prescribe pressure at x=x0=0 to p0  
     #  
     #  if we prescribe the pressure on the bottom x=xb we set  
     #  
     #  pb= 1/k * ( 1/2* (fb-f0)* xb + f0 * xb - ub*xb ) +  p0 = 1/k * ((fb+f0)/2 - ub ) * xb  + p0  
     #  
     #  which leads to ub = (fb+f0)/2-k*(pb-p0)/xb  
     #  
     def rescaleDomain(self):  
         x=self.dom.getX().copy()  
         for i in xrange(self.dom.getDim()):  
              x_inf=inf(x[i])  
              x_sup=sup(x[i])  
              if i == self.dom.getDim()-1:  
                 x[i]=-self.WIDTH*(x[i]-x_sup)/(x_inf-x_sup)  
              else:  
                 x[i]=self.WIDTH*(x[i]-x_inf)/(x_sup-x_inf)  
         self.dom.setX(x)  
     def getScalarMask(self,include_bottom=True):  
         x=self.dom.getX().copy()  
         x_inf=inf(x[self.dom.getDim()-1])  
         x_sup=sup(x[self.dom.getDim()-1])  
         out=whereZero(x[self.dom.getDim()-1]-x_sup)  
         if include_bottom: out+=whereZero(x[self.dom.getDim()-1]-x_inf)  
         return wherePositive(out)  
     def getVectorMask(self,include_bottom=True):  
         x=self.dom.getX().copy()  
         out=Vector(0.,Solution(self.dom))  
         for i in xrange(self.dom.getDim()):  
              x_inf=inf(x[i])  
              x_sup=sup(x[i])  
              if i != self.dom.getDim()-1: out[i]+=whereZero(x[i]-x_sup)  
              if i != self.dom.getDim()-1 or include_bottom: out[i]+=whereZero(x[i]-x_inf)  
         return wherePositive(out)  
   
     def setSolutionFixedBottom(self, p0, pb, f0, fb, k):  
          d=self.dom.getDim()  
          x=self.dom.getX()[d-1]  
          xb=inf(x)  
          u=Vector(0.,Solution(self.dom))+kronecker(d)[d-1]*((f0+fb)/2.-k*(pb-p0)/xb)  
          p=1./k*((fb-f0)/(xb*2.)* x**2 - (fb-f0)/2.*x)+(pb-p0)/xb*x +  p0  
          f= ((fb-f0)/xb* x + f0)*kronecker(Function(self.dom))[d-1]  
          return u,p,f  
           
     def testConstF_FixedBottom_smallK(self):  
         k=1.e-10  
         mp=self.getScalarMask(include_bottom=True)  
         mv=self.getVectorMask(include_bottom=False)  
         u_ref,p_ref,f=self.setSolutionFixedBottom(p0=2500,pb=4000.,f0=10.,fb=10.,k=k)  
         p=p_ref*mp  
         u=u_ref*mv  
         df=DarcyFlow(self.dom)  
         df.setValue(g=f,  
                       location_of_fixed_pressure=mp,  
                       location_of_fixed_flux=mv,  
                       permeability=Scalar(k,Function(self.dom)))  
         df.setTolerance(rtol=self.TOL)  
         v,p=df.solve(u_ref,p, max_iter=100, verbose=VERBOSE)  
         self.failUnless(Lsup(v-u_ref)<self.TOL*10.*Lsup(u_ref), "flux error too big.")  
         self.failUnless(Lsup(p-p_ref)<self.TOL*10.*Lsup(p_ref), "pressure error too big.")  
     def testConstF_FixedBottom_mediumK(self):  
         k=1.  
         mp=self.getScalarMask(include_bottom=True)  
         mv=self.getVectorMask(include_bottom=False)  
         u_ref,p_ref,f=self.setSolutionFixedBottom(p0=2500,pb=4000.,f0=10.,fb=10.,k=k)  
         p=p_ref*mp  
         u=u_ref*mv  
         df=DarcyFlow(self.dom)  
         df.setValue(g=f,  
                       location_of_fixed_pressure=mp,  
                       location_of_fixed_flux=mv,  
                       permeability=Scalar(k,Function(self.dom)))  
         df.setTolerance(rtol=self.TOL)  
         v,p=df.solve(u,p,max_iter=100, verbose=VERBOSE )  
         self.failUnless(Lsup(p-p_ref)<self.TOL*10.*Lsup(p_ref), "pressure error too big.")  
         self.failUnless(Lsup(v-u_ref)<self.TOL*10.*Lsup(u_ref), "flux error too big.")  
   
     def testConstF_FixedBottom_largeK(self):  
         k=1.e10  
         mp=self.getScalarMask(include_bottom=True)  
         mv=self.getVectorMask(include_bottom=False)  
         u_ref,p_ref,f=self.setSolutionFixedBottom(p0=2500,pb=4000.,f0=10.,fb=10.,k=k)  
         p=p_ref*mp  
         u=u_ref*mv  
         df=DarcyFlow(self.dom)  
         df.setValue(g=f,  
                       location_of_fixed_pressure=mp,  
                       location_of_fixed_flux=mv,  
                       permeability=Scalar(k,Function(self.dom)))  
         df.setTolerance(rtol=self.TOL)  
         v,p=df.solve(u,p, max_iter=100, verbose=VERBOSE)  
         self.failUnless(Lsup(v-u_ref)<self.TOL*10.*Lsup(u_ref), "flux error too big.")  
         self.failUnless(Lsup(p-p_ref)<self.TOL*10.*Lsup(p_ref), "pressure error too big.")  
   
     def testVarioF_FixedBottom_smallK(self):  
         k=1.e-10  
         mp=self.getScalarMask(include_bottom=True)  
         mv=self.getVectorMask(include_bottom=False)  
         u_ref,p_ref,f=self.setSolutionFixedBottom(p0=2500,pb=4000.,f0=10.,fb=30.,k=k)  
         p=p_ref*mp  
         u=u_ref*mv  
         df=DarcyFlow(self.dom)  
         df.setValue(g=f,  
                       location_of_fixed_pressure=mp,  
                       location_of_fixed_flux=mv,  
                       permeability=Scalar(k,Function(self.dom)))  
         df.setTolerance(rtol=self.TOL)  
         v,p=df.solve(u,p, max_iter=100, verbose=VERBOSE)  
         self.failUnless(Lsup(v-u_ref)<self.TOL*10.*Lsup(u_ref), "flux error too big.")  
         self.failUnless(Lsup(p-p_ref)<self.TOL*10.*Lsup(p_ref), "pressure error too big.")  
   
     def testVarioF_FixedBottom_mediumK(self):  
         k=1.  
         mp=self.getScalarMask(include_bottom=True)  
         mv=self.getVectorMask(include_bottom=False)  
         u_ref,p_ref,f=self.setSolutionFixedBottom(p0=2500,pb=4000.,f0=10.,fb=30.,k=k)  
         p=p_ref*mp  
         u=u_ref*mv  
         df=DarcyFlow(self.dom)  
         df.setValue(g=f,  
                       location_of_fixed_pressure=mp,  
                       location_of_fixed_flux=mv,  
                       permeability=Scalar(k,Function(self.dom)))  
         df.setTolerance(rtol=self.TOL)  
         v,p=df.solve(u,p, max_iter=100, verbose=VERBOSE)  
         self.failUnless(Lsup(v-u_ref)<self.TOL*10.*Lsup(u_ref), "flux error too big.")  
         self.failUnless(Lsup(p-p_ref)<self.TOL*10.*Lsup(p_ref), "pressure error too big.")  
   
     def testVarioF_FixedBottom_largeK(self):  
         k=1.e10  
         mp=self.getScalarMask(include_bottom=True)  
         mv=self.getVectorMask(include_bottom=False)  
         u_ref,p_ref,f=self.setSolutionFixedBottom(p0=2500,pb=4000.,f0=10.,fb=30.,k=k)  
         p=p_ref*mp  
         u=u_ref*mv  
         df=DarcyFlow(self.dom)  
         df.setValue(g=f,  
                       location_of_fixed_pressure=mp,  
                       location_of_fixed_flux=mv,  
                       permeability=Scalar(k,Function(self.dom)))  
         df.setTolerance(rtol=self.TOL)  
         v,p=df.solve(u,p, max_iter=100, verbose=VERBOSE)  
         self.failUnless(Lsup(v-u_ref)<self.TOL*10.*Lsup(u_ref), "flux error too big.")  
         self.failUnless(Lsup(p-p_ref)<self.TOL*10.*Lsup(p_ref), "pressure error too big.")  
   
     def testConstF_FreeBottom_smallK(self):  
         k=1.e-10  
         mp=self.getScalarMask(include_bottom=False)  
         mv=self.getVectorMask(include_bottom=True)  
         u_ref,p_ref,f=self.setSolutionFixedBottom(p0=2500,pb=4000.,f0=10.,fb=10.,k=k)  
         p=p_ref*mp  
         u=u_ref*mv  
         df=DarcyFlow(self.dom)  
         df.setValue(g=f,  
                     location_of_fixed_pressure=mp,  
                       location_of_fixed_flux=mv,  
                       permeability=Scalar(k,Function(self.dom)))  
         df.setTolerance(rtol=self.TOL)  
         v,p=df.solve(u,p, max_iter=100, verbose=VERBOSE)  
         self.failUnless(Lsup(v-u_ref)<self.TOL*10.*Lsup(u_ref), "flux error too big.")  
         self.failUnless(Lsup(p-p_ref)<self.TOL*10.*Lsup(p_ref), "pressure error too big.")  
   
     def testConstF_FreeBottom_mediumK(self):  
         k=1.  
         mp=self.getScalarMask(include_bottom=False)  
         mv=self.getVectorMask(include_bottom=True)  
         u_ref,p_ref,f=self.setSolutionFixedBottom(p0=2500,pb=4000.,f0=10.,fb=10.,k=k)  
         p=p_ref*mp  
         u=u_ref*mv  
         df=DarcyFlow(self.dom)  
         df.setValue(g=f,  
                       location_of_fixed_pressure=mp,  
                       location_of_fixed_flux=mv,  
                       permeability=Scalar(k,Function(self.dom)))  
         df.setTolerance(rtol=self.TOL)  
         v,p=df.solve(u,p, max_iter=100, verbose=VERBOSE)  
         self.failUnless(Lsup(v-u_ref)<self.TOL*10.*Lsup(u_ref), "flux error too big.")  
         self.failUnless(Lsup(p-p_ref)<self.TOL*10.*Lsup(p_ref), "pressure error too big.")  
   
     def testConstF_FreeBottom_largeK(self):  
         k=1.e10  
         mp=self.getScalarMask(include_bottom=False)  
         mv=self.getVectorMask(include_bottom=True)  
         u_ref,p_ref,f=self.setSolutionFixedBottom(p0=2500,pb=4000.,f0=10.,fb=10.,k=k)  
         p=p_ref*mp  
         u=u_ref*mv  
         df=DarcyFlow(self.dom)  
         df.setValue(g=f,  
                       location_of_fixed_pressure=mp,  
                       location_of_fixed_flux=mv,  
                       permeability=Scalar(k,Function(self.dom)))  
         df.setTolerance(rtol=self.TOL)  
         v,p=df.solve(u,p, max_iter=100, verbose=VERBOSE)  
         self.failUnless(Lsup(v-u_ref)<self.TOL*10.*Lsup(u_ref), "flux error too big.")  
         self.failUnless(Lsup(p-p_ref)<self.TOL*10.*Lsup(p_ref), "pressure error too big.")  
   
     def testVarioF_FreeBottom_smallK(self):  
         k=1.e-10  
         mp=self.getScalarMask(include_bottom=False)  
         mv=self.getVectorMask(include_bottom=True)  
         u_ref,p_ref,f=self.setSolutionFixedBottom(p0=2500,pb=4000.,f0=10.,fb=30.,k=k)  
         p=p_ref*mp  
         u=u_ref*mv  
         df=DarcyFlow(self.dom)  
         df.setValue(g=f,  
                       location_of_fixed_pressure=mp,  
                       location_of_fixed_flux=mv,  
                       permeability=Scalar(k,Function(self.dom)))  
         df.setTolerance(rtol=self.TOL)  
         v,p=df.solve(u,p, max_iter=100, verbose=VERBOSE)  
         self.failUnless(Lsup(v-u_ref)<self.TOL*25.*Lsup(u_ref), "flux error too big.")  # 25 because of disc. error.  
         self.failUnless(Lsup(p-p_ref)<self.TOL*25.*Lsup(p_ref), "pressure error too big.")  
   
     def testVarioF_FreeBottom_mediumK(self):  
         k=1.  
         mp=self.getScalarMask(include_bottom=False)  
         mv=self.getVectorMask(include_bottom=True)  
         u_ref,p_ref,f=self.setSolutionFixedBottom(p0=2500,pb=4000.,f0=10.,fb=30.,k=k)  
         p=p_ref*mp  
         u=u_ref*mv  
         df=DarcyFlow(self.dom)  
         df.setValue(g=f,  
                       location_of_fixed_pressure=mp,  
                       location_of_fixed_flux=mv,  
                       permeability=Scalar(k,Function(self.dom)))  
         df.setTolerance(rtol=self.TOL)  
         v,p=df.solve(u,p, max_iter=100, verbose=VERBOSE)  
         self.failUnless(Lsup(v-u_ref)<self.TOL*10.*Lsup(u_ref), "flux error too big.")  
         self.failUnless(Lsup(p-p_ref)<self.TOL*10.*Lsup(p_ref), "pressure error too big.")  
   
     def testVarioF_FreeBottom_largeK(self):  
         k=1.e10  
         mp=self.getScalarMask(include_bottom=False)  
         mv=self.getVectorMask(include_bottom=True)  
         u_ref,p_ref,f=self.setSolutionFixedBottom(p0=2500,pb=4000.,f0=10.,fb=30.,k=k)  
         p=p_ref*mp  
         u=u_ref*mv  
         df=DarcyFlow(self.dom)  
         df.setValue(g=f,  
                       location_of_fixed_pressure=mp,  
                       location_of_fixed_flux=mv,  
                       permeability=Scalar(k,Function(self.dom)))  
         df.setTolerance(rtol=self.TOL)  
         v,p=df.solve(u,p, max_iter=100, verbose=VERBOSE)  
         self.failUnless(Lsup(v-u_ref)<self.TOL*10.*Lsup(u_ref), "flux error too big.")  
         self.failUnless(Lsup(p-p_ref)<self.TOL*10.*Lsup(p_ref), "pressure error too big.")  
   
 class Test_Darcy2D(Test_Darcy):  
     TOL=1e-4  
     WIDTH=1.  
     def setUp(self):  
         NE=40  # wrning smaller NE may case a failure for VarioF tests due to discretization errors.  
         self.dom = Rectangle(NE,NE)  
         self.rescaleDomain()  
     def tearDown(self):  
         del self.dom  
 class Test_Darcy3D(Test_Darcy):  
     TOL=1e-4  
     WIDTH=1.  
     def setUp(self):  
         NE=25  # wrning smaller NE may case a failure for VarioF tests due to discretization errors.  
         self.dom = Brick(NE,NE,NE)  
         self.rescaleDomain()  
     def tearDown(self):  
         del self.dom  
   
402    
403  class Test_Mountains3D(unittest.TestCase):  class Test_Mountains3D(unittest.TestCase):
404     def setUp(self):     def setUp(self):
# Line 912  class Test_Rheologies(unittest.TestCase) Line 639  class Test_Rheologies(unittest.TestCase)
639           for i in xrange(len(taus)): self.checkResult(i,gamma_dot_s[i], pl.getEtaEff(gamma_dot_s[i],dt=dt),taus[i])           for i in xrange(len(taus)): self.checkResult(i,gamma_dot_s[i], pl.getEtaEff(gamma_dot_s[i],dt=dt),taus[i])
640    
641  class Test_IncompressibleIsotropicFlowCartesian(unittest.TestCase):  class Test_IncompressibleIsotropicFlowCartesian(unittest.TestCase):
642     TOL=1.e-6     TOL=1.e-5
643     VERBOSE=False     VERBOSE=False # or True
644     A=1.     A=1.
645     P_max=100     P_max=100
646     NE=2*getMPISizeWorld()     NE=2*getMPISizeWorld()
# Line 965  class Test_IncompressibleIsotropicFlowCa Line 692  class Test_IncompressibleIsotropicFlowCa
692        mod.setElasticShearModulus(self.mu)        mod.setElasticShearModulus(self.mu)
693        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])        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])
694        mod.setTolerance(self.TOL)        mod.setTolerance(self.TOL)
695        mod.setFlowTolerance(self.TOL)        mod.setEtaTolerance(self.TOL*0.1)
696    
697        BF=Vector(self.P_max,Function(self.dom))        BF=Vector(self.P_max,Function(self.dom))
698        for d in range(self.dom.getDim()):        for d in range(self.dom.getDim()):
# Line 990  class Test_IncompressibleIsotropicFlowCa Line 717  class Test_IncompressibleIsotropicFlowCa
717           t_ref=t+dt           t_ref=t+dt
718           v_ref, s_ref,p_ref=self.getReference(t_ref)           v_ref, s_ref,p_ref=self.getReference(t_ref)
719           mod.setExternals(f=matrixmult(s_ref,n)-p_ref*n, v_boundary=v_ref)           mod.setExternals(f=matrixmult(s_ref,n)-p_ref*n, v_boundary=v_ref)
720           mod.update(dt, iter_max=100, inner_iter_max=20, verbose=self.VERBOSE, usePCG=True)           # mod.update(dt, eta_iter_max=10, iter_max=50, verbose=self.VERBOSE, usePCG=True, max_correction_steps=30) new version
721             mod.update(dt, iter_max=50, verbose=self.VERBOSE, usePCG=True)
722           self.check(N_t,mod,t_ref,v_ref, s_ref,p_ref)           self.check(N_t,mod,t_ref,v_ref, s_ref,p_ref)
723           t+=dt           t+=dt
724           N_t+=1           N_t+=1
# Line 1004  class Test_IncompressibleIsotropicFlowCa Line 732  class Test_IncompressibleIsotropicFlowCa
732           error_t=abs(mod.getTime()-t_ref)/abs(t_ref)           error_t=abs(mod.getTime()-t_ref)/abs(t_ref)
733           if self.VERBOSE: print "time step ",N_t,"time = ",mod.getTime(),"errors s,p,v = ",error_s, error_p, error_v           if self.VERBOSE: print "time step ",N_t,"time = ",mod.getTime(),"errors s,p,v = ",error_s, error_p, error_v
734           self.failUnless( error_p <= 10*self.TOL, "time step %s: pressure error %s too high."%(N_t,error_p) )           self.failUnless( error_p <= 10*self.TOL, "time step %s: pressure error %s too high."%(N_t,error_p) )
          self.failUnless( error_s <= 10*self.TOL, "time step %s: stress error %s too high."%(N_t,error_s) )  
735           self.failUnless( error_v <= 10*self.TOL, "time step %s: velocity error %s too high."%(N_t,error_v) )           self.failUnless( error_v <= 10*self.TOL, "time step %s: velocity error %s too high."%(N_t,error_v) )
736           self.failUnless( error_t <= 10*self.TOL, "time step %s: time marker error %s too high."%(N_t,error_t) )           self.failUnless( error_t <= 10*self.TOL, "time step %s: time marker error %s too high."%(N_t,error_t) )
737             self.failUnless( error_s <= 10*self.TOL, "time step %s: stress error %s too high."%(N_t,error_s) )
738     def tearDown(self):     def tearDown(self):
739          del self.dom          del self.dom
740    
# Line 1103  class Test_FaultSystem(unittest.TestCase Line 831  class Test_FaultSystem(unittest.TestCase
831        f.addFault(V0=[0.5,0.], strikes=[3.*pi/4], ls=[0.70710678118654757], tag=1)        f.addFault(V0=[0.5,0.], strikes=[3.*pi/4], ls=[0.70710678118654757], tag=1)
832        f.addFault(V0=[1.,0.5], strikes=[pi, pi/2], ls=[0.5,0.5], tag=2)        f.addFault(V0=[1.,0.5], strikes=[pi, pi/2], ls=[0.5,0.5], tag=2)
833    
834        m, t, l=f.getMaxValue(x[0]*(1.-x[0])*(1-x[1]))        u=x[0]*(1.-x[0])*(1-x[1])
835          t, loc=f.getMaxValue(u)
836          p=f.getParametrization(x,t)[0]
837          m, l=loc(u), loc(p)
838        self.failUnless(  m == 0.25, "wrong max value")        self.failUnless(  m == 0.25, "wrong max value")
839        self.failUnless(  t == 1, "wrong max tag")        self.failUnless(  t == 1, "wrong max tag")
840        self.failUnless(  l == 0., "wrong max location")        self.failUnless(  l == 0., "wrong max location")
841        m, t, l=f.getMaxValue(x[1]*(1.-x[1])*(1-x[0])*x[0])  
842          u=x[1]*(1.-x[1])*(1-x[0])*x[0]
843          t, loc=f.getMaxValue(u)
844          p=f.getParametrization(x,t)[0]
845          m, l=loc(u), loc(p)
846        self.failUnless(  m == 0.0625, "wrong max value")        self.failUnless(  m == 0.0625, "wrong max value")
847        self.failUnless(  t == 2, "wrong max tag")        self.failUnless(  t == 2, "wrong max tag")
848        self.failUnless(  l == 0.5, "wrong max location")        self.failUnless(  l == 0.5, "wrong max location")
849        m, t, l=f.getMaxValue(x[0]*(1.-x[0])*x[1])  
850          u=x[0]*(1.-x[0])*x[1]
851          t, loc=f.getMaxValue(u)
852          p=f.getParametrization(x,t)[0]
853          m, l=loc(u), loc(p)
854        self.failUnless(  m == 0.25, "wrong max value")        self.failUnless(  m == 0.25, "wrong max value")
855        self.failUnless(  t == 2, "wrong max tag")        self.failUnless(  t == 2, "wrong max tag")
856        self.failUnless(  l == 1.0, "wrong max location")        self.failUnless(  l == 1.0, "wrong max location")
857        m, t, l= f.getMaxValue(x[1]*(1.-x[1])*x[0])  
858          u=x[1]*(1.-x[1])*x[0]
859          t, loc=f.getMaxValue(u)
860          p=f.getParametrization(x,t)[0]
861          m, l=loc(u), loc(p)
862        self.failUnless(  m == 0.25, "wrong max value")        self.failUnless(  m == 0.25, "wrong max value")
863        self.failUnless(  t == 2, "wrong max tag")        self.failUnless(  t == 2, "wrong max tag")
864        self.failUnless(  l == 0., "wrong max location")        self.failUnless(  l == 0., "wrong max location")
865        m, t, l= f.getMaxValue(x[1]*(1.-x[1])*(1.-x[0]))  
866          u=x[1]*(1.-x[1])*(1.-x[0])
867          t, loc=f.getMaxValue(u)
868          p=f.getParametrization(x,t)[0]
869          m, l=loc(u), loc(p)
870        self.failUnless(  m == 0.25, "wrong max value")        self.failUnless(  m == 0.25, "wrong max value")
871        self.failUnless(  t == 1, "wrong max tag")        self.failUnless(  t == 1, "wrong max tag")
872        self.failUnless(  abs(l-0.70710678118654) <= self.EPS,  "wrong max location")        self.failUnless(  abs(l-0.70710678118654) <= self.EPS,  "wrong max location")
# Line 1130  class Test_FaultSystem(unittest.TestCase Line 877  class Test_FaultSystem(unittest.TestCase
877        f.addFault(V0=[0.5,0.], strikes=[3.*pi/4], ls=[0.70710678118654757], tag=1)        f.addFault(V0=[0.5,0.], strikes=[3.*pi/4], ls=[0.70710678118654757], tag=1)
878        f.addFault(V0=[1.,0.5], strikes=[pi, pi/2], ls=[0.5,0.5], tag=2)        f.addFault(V0=[1.,0.5], strikes=[pi, pi/2], ls=[0.5,0.5], tag=2)
879    
880        m, t, l=f.getMinValue(-x[0]*(1.-x[0])*(1-x[1]))        u=-x[0]*(1.-x[0])*(1-x[1])
881          t, loc=f.getMinValue(u)
882          p=f.getParametrization(x,t)[0]
883          m, l=loc(u), loc(p)
884        self.failUnless(  m == -0.25, "wrong min value")        self.failUnless(  m == -0.25, "wrong min value")
885        self.failUnless(  t == 1, "wrong min tag")        self.failUnless(  t == 1, "wrong min tag")
886        self.failUnless(  l == 0., "wrong min location")        self.failUnless(  l == 0., "wrong min location")
887        m, t, l=f.getMinValue(-x[1]*(1.-x[1])*(1-x[0])*x[0])        u=-x[1]*(1.-x[1])*(1-x[0])*x[0]
888          t, loc=f.getMinValue(u)
889          p=f.getParametrization(x,t)[0]
890          m, l=loc(u), loc(p)
891        self.failUnless(  m == -0.0625, "wrong min value")        self.failUnless(  m == -0.0625, "wrong min value")
892        self.failUnless(  t == 2, "wrong min tag")        self.failUnless(  t == 2, "wrong min tag")
893        self.failUnless(  l == 0.5, "wrong min location")        self.failUnless(  l == 0.5, "wrong min location")
894        m, t, l=f.getMinValue(-x[0]*(1.-x[0])*x[1])        u=-x[0]*(1.-x[0])*x[1]
895          t, loc=f.getMinValue(u)
896          p=f.getParametrization(x,t)[0]
897          m, l=loc(u), loc(p)
898        self.failUnless(  m == -0.25, "wrong min value")        self.failUnless(  m == -0.25, "wrong min value")
899        self.failUnless(  t == 2, "wrong min tag")        self.failUnless(  t == 2, "wrong min tag")
900        self.failUnless(  l == 1.0, "wrong min location")        self.failUnless(  l == 1.0, "wrong min location")
901        m, t, l= f.getMinValue(-x[1]*(1.-x[1])*x[0])        u=-x[1]*(1.-x[1])*x[0]
902          t, loc=f.getMinValue(u)
903          p=f.getParametrization(x,t)[0]
904          m, l=loc(u), loc(p)
905        self.failUnless(  m == -0.25, "wrong min value")        self.failUnless(  m == -0.25, "wrong min value")
906        self.failUnless(  t == 2, "wrong min tag")        self.failUnless(  t == 2, "wrong min tag")
907        self.failUnless(  l == 0., "wrong min location")        self.failUnless(  l == 0., "wrong min location")
908        m, t, l= f.getMinValue(-x[1]*(1.-x[1])*(1.-x[0]))        u=-x[1]*(1.-x[1])*(1.-x[0])
909          t, loc=f.getMinValue(u)
910          p=f.getParametrization(x,t)[0]
911          m, l=loc(u), loc(p)
912        self.failUnless(  m == -0.25, "wrong min value")        self.failUnless(  m == -0.25, "wrong min value")
913        self.failUnless(  t == 1, "wrong min tag")        self.failUnless(  t == 1, "wrong min tag")
914        self.failUnless(  abs(l-0.70710678118654) <= self.EPS,  "wrong min location")        self.failUnless(  abs(l-0.70710678118654) <= self.EPS,  "wrong min location")
# Line 1695  class Test_FaultSystem(unittest.TestCase Line 1457  class Test_FaultSystem(unittest.TestCase
1457        s,d=f.getSideAndDistance([5.,12.,-4], tag=2)        s,d=f.getSideAndDistance([5.,12.,-4], tag=2)
1458        self.failUnless( s<0, "wrong side.")        self.failUnless( s<0, "wrong side.")
1459        self.failUnless( abs(d-2.*0.70710678118654757)<self.EPS, "wrong distance.")        self.failUnless( abs(d-2.*0.70710678118654757)<self.EPS, "wrong distance.")
1460    
1461  if __name__ == '__main__':  if __name__ == '__main__':
1462     suite = unittest.TestSuite()     suite = unittest.TestSuite()
1463     suite.addTest(unittest.makeSuite(Test_FaultSystem))     # suite.addTest(unittest.makeSuite(Test_FaultSystem))
1464     suite.addTest(unittest.makeSuite(Test_StokesProblemCartesian2D))     # suite.addTest(unittest.makeSuite(Test_StokesProblemCartesian2D))
    suite.addTest(unittest.makeSuite(Test_Darcy3D))  
    suite.addTest(unittest.makeSuite(Test_Darcy2D))  
1465     suite.addTest(unittest.makeSuite(Test_StokesProblemCartesian3D))     suite.addTest(unittest.makeSuite(Test_StokesProblemCartesian3D))
1466     suite.addTest(unittest.makeSuite(Test_Mountains3D))     suite.addTest(unittest.makeSuite(Test_Mountains3D))
1467     suite.addTest(unittest.makeSuite(Test_Mountains2D))     suite.addTest(unittest.makeSuite(Test_Mountains2D))
# Line 1709  if __name__ == '__main__': Line 1469  if __name__ == '__main__':
1469     suite.addTest(unittest.makeSuite(Test_IncompressibleIsotropicFlowCartesian))     suite.addTest(unittest.makeSuite(Test_IncompressibleIsotropicFlowCartesian))
1470     s=unittest.TextTestRunner(verbosity=2).run(suite)     s=unittest.TextTestRunner(verbosity=2).run(suite)
1471     if not s.wasSuccessful(): sys.exit(1)     if not s.wasSuccessful(): sys.exit(1)
   

Legend:
Removed from v.2675  
changed lines
  Added in v.3527

  ViewVC Help
Powered by ViewVC 1.1.26