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

revision 2155 by gross, Wed Nov 26 08:13:00 2008 UTC revision 2156 by gross, Mon Dec 15 05:09:02 2008 UTC
# Line 33  except KeyError: Line 33  except KeyError:
33       FINLEY_WORKDIR='.'       FINLEY_WORKDIR='.'
34
35
36  VERBOSE=False # or True  VERBOSE=False #
37
38  from esys.escript import *  from esys.escript import *
39  from esys.escript.models import StokesProblemCartesian  from esys.escript.models import StokesProblemCartesian
# Line 89  class Test_StokesProblemCartesian2D(unit Line 89  class Test_StokesProblemCartesian2D(unit
90         u0=(1-x[0])*x[0]*[0.,1.]         u0=(1-x[0])*x[0]*[0.,1.]
91         p0=Scalar(P1,ReducedSolution(self.domain))         p0=Scalar(P1,ReducedSolution(self.domain))
92         sp.setTolerance(self.TOL)         sp.setTolerance(self.TOL*0.2)
93         u,p=sp.solve(u0,p0,show_details=VERBOSE, verbose=VERBOSE or True,max_iter=100,useUzawa=True)         u,p=sp.solve(u0,p0,show_details=VERBOSE, verbose=VERBOSE,max_iter=100,useUzawa=True)
94         error_v0=Lsup(u[0]-u0[0])         error_v0=Lsup(u[0]-u0[0])
95         error_v1=Lsup(u[1]-u0[1])/0.25         error_v1=Lsup(u[1]-u0[1])/0.25
96         error_p=Lsup(P1*x[0]*x[1]+p)         error_p=Lsup(P1*x[0]*x[1]+p)
97           saveVTK("d.vtu",p=p, e=P1*x[0]*x[1]+p, p_ref=P1*x[0]*x[1])
98         self.failUnless(error_v0<10*self.TOL, "0-velocity error too large.")         self.failUnless(error_v0<10*self.TOL, "0-velocity error too large.")
99         self.failUnless(error_v1<10*self.TOL, "1-velocity error too large.")         self.failUnless(error_v1<10*self.TOL, "1-velocity error too large.")
100         self.failUnless(error_p<11*self.TOL, "pressure error too large.")         self.failUnless(error_p<10*self.TOL, "pressure error too large.")
101
102     def test_PCG_P_large(self):     def test_PCG_P_large(self):
103         ETA=1.         ETA=1.
# Line 142  class Test_StokesProblemCartesian2D(unit Line 143  class Test_StokesProblemCartesian2D(unit
143         p0=Scalar(P1,ReducedSolution(self.domain))         p0=Scalar(P1,ReducedSolution(self.domain))
144         sp.setTolerance(self.TOL)         sp.setTolerance(self.TOL)
145         u,p=sp.solve(u0,p0,show_details=VERBOSE, verbose=VERBOSE,max_iter=50,useUzawa=False,iter_restart=18)         u,p=sp.solve(u0,p0,show_details=VERBOSE, verbose=VERBOSE,max_iter=50,useUzawa=False,iter_restart=18)
# u,p=sp.solve(u0,p0,show_details=VERBOSE, verbose=VERBOSE,max_iter=100,useUzawa=False,iter_restart=20)
146
147         error_v0=Lsup(u[0]-u0[0])         error_v0=Lsup(u[0]-u0[0])
148         error_v1=Lsup(u[1]-u0[1])/0.25         error_v1=Lsup(u[1]-u0[1])/0.25
149         zz=P1*x[0]*x[1]+p         error_p=Lsup(P1*x[0]*x[1]+p)
error_p=Lsup(zz-integrate(zz))
150         self.failUnless(error_p<10*self.TOL, "pressure error too large.")         self.failUnless(error_p<10*self.TOL, "pressure error too large.")
151         self.failUnless(error_v0<10*self.TOL, "0-velocity error too large.")         self.failUnless(error_v0<10*self.TOL, "0-velocity error too large.")
152         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 168  class Test_StokesProblemCartesian2D(unit Line 167  class Test_StokesProblemCartesian2D(unit
168         u0=(1-x[0])*x[0]*[0.,1.]         u0=(1-x[0])*x[0]*[0.,1.]
169         p0=Scalar(P1,ReducedSolution(self.domain))         p0=Scalar(P1,ReducedSolution(self.domain))
170         sp.setTolerance(self.TOL*0+1.e-6)         sp.setTolerance(self.TOL*0.1)
# sp.setSubToleranceReductionFactor(0.1)
171         u,p=sp.solve(u0,p0,show_details=VERBOSE, verbose=VERBOSE,max_iter=20,useUzawa=False)         u,p=sp.solve(u0,p0,show_details=VERBOSE, verbose=VERBOSE,max_iter=20,useUzawa=False)
172
173         error_v0=Lsup(u[0]-u0[0])         error_v0=Lsup(u[0]-u0[0])
174         error_v1=Lsup(u[1]-u0[1])/0.25         error_v1=Lsup(u[1]-u0[1])/0.25
175         zz=P1*x[0]*x[1]+p         error_p=Lsup(P1*x[0]*x[1]+p)
error_p=Lsup(zz-integrate(zz))
# print error_p, error_v0,error_v1
176         self.failUnless(error_p<10*self.TOL, "pressure error too large.")         self.failUnless(error_p<10*self.TOL, "pressure error too large.")
177         self.failUnless(error_v0<10*self.TOL, "0-velocity error too large.")         self.failUnless(error_v0<10*self.TOL, "0-velocity error too large.")
178         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 202  class Test_StokesProblemCartesian2D(unit Line 198  class Test_StokesProblemCartesian2D(unit
198
199         error_v0=Lsup(u[0]-u0[0])         error_v0=Lsup(u[0]-u0[0])
200         error_v1=Lsup(u[1]-u0[1])/0.25         error_v1=Lsup(u[1]-u0[1])/0.25
201         zz=P1*x[0]*x[1]+p         error_p=Lsup(P1*x[0]*x[1]+p)/P1
error_p=Lsup(zz-integrate(zz))/P1
202         self.failUnless(error_p<10*self.TOL, "pressure error too large.")         self.failUnless(error_p<10*self.TOL, "pressure error too large.")
203         self.failUnless(error_v0<10*self.TOL, "0-velocity error too large.")         self.failUnless(error_v0<10*self.TOL, "0-velocity error too large.")
204         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 224  class Test_StokesProblemCartesian3D(unit Line 219  class Test_StokesProblemCartesian3D(unit
219         x=self.domain.getX()         x=self.domain.getX()
221                +whereZero(x[0]-1)  * [1.,1.,1.] \                +whereZero(x[0]-1)  * [1.,1.,1.] \
222                +whereZero(x[1])    * [1.,1.,1.] \                +whereZero(x[1])    * [1.,0.,1.] \
223                +whereZero(x[1]-1)  * [1.,1.,1.] \                +whereZero(x[1]-1)  * [1.,1.,1.] \
224                +whereZero(x[2])    * [1.,1.,0.] \                +whereZero(x[2])    * [1.,1.,0.] \
225                +whereZero(x[2]-1)  * [1.,1.,1.]                +whereZero(x[2]-1)  * [1.,1.,1.]
# Line 235  class Test_StokesProblemCartesian3D(unit Line 230  class Test_StokesProblemCartesian3D(unit
231         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.]
232         p0=Scalar(P1,ReducedSolution(self.domain))         p0=Scalar(P1,ReducedSolution(self.domain))
233         sp.setTolerance(1.e-8)         sp.setTolerance(self.TOL)
234         u,p=sp.solve(u0,p0,show_details=VERBOSE, verbose=VERBOSE,max_iter=100,useUzawa=True)         u,p=sp.solve(u0,p0,show_details=VERBOSE, verbose=VERBOSE ,max_iter=100,useUzawa=True)
235
236         error_v0=Lsup(u[0]-u0[0])         error_v0=Lsup(u[0]-u0[0])
237         error_v1=Lsup(u[1]-u0[1])         error_v1=Lsup(u[1]-u0[1])
238         error_v2=Lsup(u[2]-u0[2])/0.25**2         error_v2=Lsup(u[2]-u0[2])/0.25**2
239         zz=P1*x[0]*x[1]*x[2]+p         error_p=Lsup(P1*x[0]*x[1]*x[2]+p)
error_p=Lsup(zz-integrate(zz))
# print error_p, error_v0,error_v1,error_v2
240         self.failUnless(error_p<10*self.TOL, "pressure error too large.")         self.failUnless(error_p<10*self.TOL, "pressure error too large.")
241         self.failUnless(error_v0<10*self.TOL, "0-velocity error too large.")         self.failUnless(error_v0<10*self.TOL, "0-velocity error too large.")
242         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 257  class Test_StokesProblemCartesian3D(unit Line 250  class Test_StokesProblemCartesian3D(unit
250         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.]
252                +whereZero(x[0]-1)  * [1.,1.,1.] \                +whereZero(x[0]-1)  * [1.,1.,1.] \
253                +whereZero(x[1])    * [1.,1.,1.] \                +whereZero(x[1])    * [1.,0.,1.] \
254                +whereZero(x[1]-1)  * [1.,1.,1.] \                +whereZero(x[1]-1)  * [1.,1.,1.] \
255                +whereZero(x[2])    * [1.,1.,0.] \                +whereZero(x[2])    * [1.,1.,0.] \
256                +whereZero(x[2]-1)  * [1.,1.,1.]                +whereZero(x[2]-1)  * [1.,1.,1.]
# Line 268  class Test_StokesProblemCartesian3D(unit Line 261  class Test_StokesProblemCartesian3D(unit
262         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.]
263         p0=Scalar(P1,ReducedSolution(self.domain))         p0=Scalar(P1,ReducedSolution(self.domain))
264         sp.setTolerance(1.e-8)         sp.setTolerance(self.TOL)
265         u,p=sp.solve(u0,p0,show_details=VERBOSE, verbose=VERBOSE,max_iter=100,useUzawa=True)         u,p=sp.solve(u0,p0,show_details=VERBOSE, verbose=VERBOSE ,max_iter=100,useUzawa=True)
266
267         error_v0=Lsup(u[0]-u0[0])         error_v0=Lsup(u[0]-u0[0])
268         error_v1=Lsup(u[1]-u0[1])         error_v1=Lsup(u[1]-u0[1])
269         error_v2=Lsup(u[2]-u0[2])/0.25**2         error_v2=Lsup(u[2]-u0[2])/0.25**2
270         zz=P1*x[0]*x[1]*x[2]+p         error_p=Lsup(P1*x[0]*x[1]*x[2]+p)
error_p=Lsup(zz-integrate(zz))/P1
271         self.failUnless(error_p<10*self.TOL, "pressure error too large.")         self.failUnless(error_p<10*self.TOL, "pressure error too large.")
272         self.failUnless(error_v0<10*self.TOL, "0-velocity error too large.")         self.failUnless(error_v0<10*self.TOL, "0-velocity error too large.")
273         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 288  class Test_StokesProblemCartesian3D(unit Line 280  class Test_StokesProblemCartesian3D(unit
280         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.]
282                +whereZero(x[0]-1)  * [1.,1.,1.] \                +whereZero(x[0]-1)  * [1.,1.,1.] \
283                +whereZero(x[1])    * [1.,1.,1.] \                +whereZero(x[1])    * [1.,0.,1.] \
284                +whereZero(x[1]-1)  * [1.,1.,1.] \                +whereZero(x[1]-1)  * [1.,1.,1.] \
285                +whereZero(x[2])    * [1.,1.,0.] \                +whereZero(x[2])    * [1.,1.,0.] \
286                +whereZero(x[2]-1)  * [1.,1.,1.]                +whereZero(x[2]-1)  * [1.,1.,1.]
# Line 299  class Test_StokesProblemCartesian3D(unit Line 291  class Test_StokesProblemCartesian3D(unit
292         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.]
293         p0=Scalar(P1,ReducedSolution(self.domain))         p0=Scalar(P1,ReducedSolution(self.domain))
294         sp.setTolerance(1.e-8)         sp.setTolerance(self.TOL)
295         u,p=sp.solve(u0,p0,show_details=VERBOSE, verbose=VERBOSE,max_iter=100,useUzawa=True)         u,p=sp.solve(u0,-p0,show_details=VERBOSE, verbose=VERBOSE ,max_iter=100,useUzawa=True)
296
297         error_v0=Lsup(u[0]-u0[0])         error_v0=Lsup(u[0]-u0[0])
298         error_v1=Lsup(u[1]-u0[1])         error_v1=Lsup(u[1]-u0[1])
299         error_v2=Lsup(u[2]-u0[2])/0.25**2         error_v2=Lsup(u[2]-u0[2])/0.25**2
300         zz=P1*x[0]*x[1]*x[2]+p         error_p=Lsup(P1*x[0]*x[1]*x[2]+p)/P1
error_p=Lsup(zz-integrate(zz))/P1
# print error_p, error_v0,error_v1,error_v2
301         self.failUnless(error_p<10*self.TOL, "pressure error too large.")         self.failUnless(error_p<10*self.TOL, "pressure error too large.")
302         self.failUnless(error_v0<10*self.TOL, "0-velocity error too large.")         self.failUnless(error_v0<10*self.TOL, "0-velocity error too large.")
303         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 333  class Test_StokesProblemCartesian3D(unit Line 323  class Test_StokesProblemCartesian3D(unit
324         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.]
325         p0=Scalar(P1,ReducedSolution(self.domain))         p0=Scalar(P1,ReducedSolution(self.domain))
326         sp.setTolerance(1.e-8)         sp.setTolerance(self.TOL)
327         u,p=sp.solve(u0,p0,show_details=VERBOSE, verbose=VERBOSE,max_iter=100,useUzawa=False,iter_restart=20)         u,p=sp.solve(u0,p0,show_details=VERBOSE, verbose=VERBOSE,max_iter=100,useUzawa=False,iter_restart=20)
328
329         error_v0=Lsup(u[0]-u0[0])         error_v0=Lsup(u[0]-u0[0])
330         error_v1=Lsup(u[1]-u0[1])         error_v1=Lsup(u[1]-u0[1])
331         error_v2=Lsup(u[2]-u0[2])/0.25**2         error_v2=Lsup(u[2]-u0[2])/0.25**2
332         zz=P1*x[0]*x[1]*x[2]+p         error_p=Lsup(P1*x[0]*x[1]*x[2]+p)
error_p=Lsup(zz-integrate(zz))
333         # print error_p, error_v0,error_v1,error_v2         # print error_p, error_v0,error_v1,error_v2
334         self.failUnless(error_p<10*self.TOL, "pressure error too large.")         self.failUnless(error_p<10*self.TOL, "pressure error too large.")
335         self.failUnless(error_v0<10*self.TOL, "0-velocity error too large.")         self.failUnless(error_v0<10*self.TOL, "0-velocity error too large.")
# Line 365  class Test_StokesProblemCartesian3D(unit Line 354  class Test_StokesProblemCartesian3D(unit
355         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.]
356         p0=Scalar(P1,ReducedSolution(self.domain))         p0=Scalar(P1,ReducedSolution(self.domain))
357         sp.setTolerance(1.e-8)         sp.setTolerance(self.TOL)
358         u,p=sp.solve(u0,p0,show_details=VERBOSE, verbose=VERBOSE,max_iter=100,useUzawa=False)         u,p=sp.solve(u0,p0,show_details=VERBOSE, verbose=VERBOSE ,max_iter=100,useUzawa=False)
359
360         error_v0=Lsup(u[0]-u0[0])         error_v0=Lsup(u[0]-u0[0])
361         error_v1=Lsup(u[1]-u0[1])         error_v1=Lsup(u[1]-u0[1])
362         error_v2=Lsup(u[2]-u0[2])/0.25**2         error_v2=Lsup(u[2]-u0[2])/0.25**2
363         zz=P1*x[0]*x[1]*x[2]+p         error_p=Lsup(P1*x[0]*x[1]*x[2]+p)/P1
error_p=Lsup(zz-integrate(zz))/P1
# print error_p, error_v0,error_v1,error_v2
364         self.failUnless(error_p<10*self.TOL, "pressure error too large.")         self.failUnless(error_p<10*self.TOL, "pressure error too large.")
365         self.failUnless(error_v0<10*self.TOL, "0-velocity error too large.")         self.failUnless(error_v0<10*self.TOL, "0-velocity error too large.")
366         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 386  class Test_StokesProblemCartesian3D(unit Line 373  class Test_StokesProblemCartesian3D(unit
373         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                +whereZero(x[0]-1)  * [1.,1.,1.] \                +whereZero(x[0]-1)  * [1.,1.,1.] \
376                +whereZero(x[1])    * [1.,1.,1.] \                +whereZero(x[1])    * [1.,0.,1.] \
377                +whereZero(x[1]-1)  * [1.,1.,1.] \                +whereZero(x[1]-1)  * [1.,1.,1.] \
378                +whereZero(x[2])    * [1.,1.,0.] \                +whereZero(x[2])    * [1.,1.,0.] \
379                +whereZero(x[2]-1)  * [1.,1.,1.]                +whereZero(x[2]-1)  * [1.,1.,1.]
# Line 397  class Test_StokesProblemCartesian3D(unit Line 384  class Test_StokesProblemCartesian3D(unit
385         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.]
386         p0=Scalar(P1,ReducedSolution(self.domain))         p0=Scalar(P1,ReducedSolution(self.domain))
387         sp.setTolerance(self.TOL+0*1.e-8)         sp.setTolerance(self.TOL)
388         # sp.setSubToleranceReductionFactor(1.e-8/self.TOL)         u,p=sp.solve(u0,p0,show_details=VERBOSE, verbose=VERBOSE ,max_iter=100,useUzawa=False)
# sp.setSubToleranceReductionFactor(None)
u,p=sp.solve(u0,p0,show_details=VERBOSE, verbose=VERBOSE,max_iter=100,useUzawa=False)
389
390         error_v0=Lsup(u[0]-u0[0])         error_v0=Lsup(u[0]-u0[0])
391         error_v1=Lsup(u[1]-u0[1])         error_v1=Lsup(u[1]-u0[1])
392         error_v2=Lsup(u[2]-u0[2])/0.25**2         error_v2=Lsup(u[2]-u0[2])/0.25**2
393         zz=P1*x[0]*x[1]*x[2]+p         error_p=Lsup(P1*x[0]*x[1]*x[2]+p)/P1
error_p=Lsup(zz-integrate(zz))/P1
# print error_p, error_v0,error_v1,error_v2
394         self.failUnless(error_p<10*self.TOL, "pressure error too large.")         self.failUnless(error_p<10*self.TOL, "pressure error too large.")
395         self.failUnless(error_v0<10*self.TOL, "0-velocity error too large.")         self.failUnless(error_v0<10*self.TOL, "0-velocity error too large.")
396         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 511  class Test_Darcy(unittest.TestCase): Line 494  class Test_Darcy(unittest.TestCase):
494                        location_of_fixed_pressure=mp,                        location_of_fixed_pressure=mp,
495                        location_of_fixed_flux=mv,                        location_of_fixed_flux=mv,
496                        permeability=Scalar(k,Function(self.dom)))                        permeability=Scalar(k,Function(self.dom)))
497          v,p=df.solve(u,p,atol=0,rtol=self.TOL, max_iter=100, verbose=VERBOSE,sub_rtol=self.TOL/200)          v,p=df.solve(u,p,atol=0,rtol=self.TOL, max_iter=100, verbose=VERBOSE,sub_rtol=self.TOL*self.TOL)
498          self.failUnless(Lsup(v-u_ref)<self.TOL*10.*Lsup(u_ref), "flux error too big.")          self.failUnless(Lsup(v-u_ref)<self.TOL*10.*Lsup(u_ref), "flux error too big.")
499          self.failUnless(Lsup(p-p_ref)<self.TOL*10.*Lsup(p_ref), "pressure error too big.")          self.failUnless(Lsup(p-p_ref)<self.TOL*10.*Lsup(p_ref), "pressure error too big.")
500

Legend:
 Removed from v.2155 changed lines Added in v.2156