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

revision 2960 by gross, Tue Mar 2 07:54:11 2010 UTC revision 3528 by gross, Tue May 24 07:02:08 2011 UTC
# Line 1  Line 1
1    # -*- coding: utf-8 -*-
2
3  ########################################################  ########################################################
4  #  #
# 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 398  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)
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)
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
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.TEST_TOL*Lsup(u_ref), "flux error too big.")
self.failUnless(Lsup(p-p_ref)<self.TEST_TOL*Lsup(p_ref), "pressure error too big.")
def testConstF_FixedBottom_mediumK(self):
k=1.
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.TEST_TOL*Lsup(p_ref), "pressure error too big.")
self.failUnless(Lsup(v-u_ref)<self.TEST_TOL*Lsup(u_ref), "flux error too big.")

def testConstF_FixedBottom_largeK(self):
k=1.e10
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.TEST_TOL*Lsup(u_ref), "flux error too big.")
self.failUnless(Lsup(p-p_ref)<self.TEST_TOL*Lsup(p_ref), "pressure error too big.")

def testVarioF_FixedBottom_smallK(self):
k=1.e-10
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.TEST_TOL*Lsup(u_ref), "flux error too big.")
self.failUnless(Lsup(p-p_ref)<self.TEST_TOL*Lsup(p_ref), "pressure error too big.")

def testVarioF_FixedBottom_mediumK(self):
k=1.
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.TEST_TOL*Lsup(u_ref), "flux error too big.")
self.failUnless(Lsup(p-p_ref)<self.TEST_TOL*Lsup(p_ref), "pressure error too big.")

def testVarioF_FixedBottom_largeK(self):
k=1.e10
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.TEST_TOL*Lsup(u_ref), "flux error too big.")
self.failUnless(Lsup(p-p_ref)<self.TEST_TOL*Lsup(p_ref), "pressure error too big.")

def testConstF_FreeBottom_smallK(self):
k=1.e-10
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.TEST_TOL*Lsup(u_ref), "flux error too big.")
self.failUnless(Lsup(p-p_ref)<self.TEST_TOL*Lsup(p_ref), "pressure error too big.")

def testConstF_FreeBottom_mediumK(self):
k=1.
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.TEST_TOL*Lsup(u_ref), "flux error too big.")
self.failUnless(Lsup(p-p_ref)<self.TEST_TOL*Lsup(p_ref), "pressure error too big.")

def testConstF_FreeBottom_largeK(self):
k=1.e10
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.TEST_TOL*Lsup(u_ref), "flux error too big.")
self.failUnless(Lsup(p-p_ref)<self.TEST_TOL*Lsup(p_ref), "pressure error too big.")

def testVarioF_FreeBottom_smallK(self):
k=1.e-10
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.TEST_TOL*Lsup(u_ref), "flux error too big.")
self.failUnless(Lsup(p-p_ref)<self.TEST_TOL*Lsup(p_ref), "pressure error too big.")

def testVarioF_FreeBottom_mediumK(self):
k=1.
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.TEST_TOL*Lsup(u_ref), "flux error too big.")
self.failUnless(Lsup(p-p_ref)<self.TEST_TOL*Lsup(p_ref), "pressure error too big.")

def testVarioF_FreeBottom_largeK(self):
k=1.e10
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.TEST_TOL*Lsup(u_ref), "flux error too big.")
self.failUnless(Lsup(p-p_ref)<self.TEST_TOL*Lsup(p_ref), "pressure error too big.")

class Test_Darcy2D(Test_Darcy):
TOL=1e-6
TEST_TOL=2.e-3
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-6
WIDTH=1.
TEST_TOL=4.e-3
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 1734  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()