# Contents of /trunk/finley/test/python/run_darcy.py

Revision 3569 - (show annotations)
Thu Sep 1 02:42:36 2011 UTC (8 years, 1 month ago) by gross
File MIME type: text/x-python
File size: 14181 byte(s)
```Darcy flow solver and docu reviewed. Coupled solvers have been removed.

```
 1 # -*- coding: utf-8 -*- 2 3 ######################################################## 4 # 5 # Copyright (c) 2003-2010 by University of Queensland 6 # Earth Systems Science Computational Center (ESSCC) 7 8 # 9 # Primary Business: Queensland, Australia 10 # Licensed under the Open Software License version 3.0 11 12 # 13 ######################################################## 14 15 __copyright__="""Copyright (c) 2003-2010 by University of Queensland 16 Earth Systems Science Computational Center (ESSCC) 17 http://www.uq.edu.au/esscc 18 Primary Business: Queensland, Australia""" 19 __license__="""Licensed under the Open Software License version 3.0 20 21 __url__= 22 23 import unittest 24 import tempfile 25 26 27 28 VERBOSE=False and True 29 30 from esys.escript import * 31 from esys.escript.models import DarcyFlow 32 from esys.finley import Rectangle, Brick 33 34 from math import pi 35 import numpy 36 import sys 37 import os 38 #==================================================================================================================== 39 try: 40 FINLEY_WORKDIR=os.environ['FINLEY_WORKDIR'] 41 except KeyError: 42 FINLEY_WORKDIR='.' 43 44 class Test_Darcy(unittest.TestCase): 45 # this is a simple test for the darcy flux problem 46 # 47 # 48 # p = 1/k * ( 1/2* (fb-f0)/xb* x **2 + f0 * x - ub*x ) + p0 49 # 50 # with f = (fb-f0)/xb* x + f0 51 # 52 # u = f - k * p,x = ub 53 # 54 # we prescribe pressure at x=x0=0 to p0 55 # 56 # if we prescribe the pressure on the bottom x=xb we set 57 # 58 # pb= 1/k * ( 1/2* (fb-f0)* xb + f0 * xb - ub*xb ) + p0 = 1/k * ((fb+f0)/2 - ub ) * xb + p0 59 # 60 # which leads to ub = (fb+f0)/2-k*(pb-p0)/xb 61 # 62 def rescaleDomain(self): 63 x=self.dom.getX().copy() 64 for i in xrange(self.dom.getDim()): 65 x_inf=inf(x[i]) 66 x_sup=sup(x[i]) 67 if i == self.dom.getDim()-1: 68 x[i]=-self.WIDTH*(x[i]-x_sup)/(x_inf-x_sup) 69 else: 70 x[i]=self.WIDTH*(x[i]-x_inf)/(x_sup-x_inf) 71 self.dom.setX(x) 72 def getScalarMask(self,include_bottom=True): 73 x=self.dom.getX().copy() 74 x_inf=inf(x[self.dom.getDim()-1]) 75 x_sup=sup(x[self.dom.getDim()-1]) 76 out=whereZero(x[self.dom.getDim()-1]-x_sup) 77 if include_bottom: out+=whereZero(x[self.dom.getDim()-1]-x_inf) 78 return wherePositive(out) 79 def getVectorMask(self,include_bottom=True): 80 x=self.dom.getX().copy() 81 out=Vector(0.,Solution(self.dom)) 82 for i in xrange(self.dom.getDim()): 83 x_inf=inf(x[i]) 84 x_sup=sup(x[i]) 85 if i != self.dom.getDim()-1: out[i]+=whereZero(x[i]-x_sup) 86 if i != self.dom.getDim()-1 or include_bottom: out[i]+=whereZero(x[i]-x_inf) 87 return wherePositive(out) 88 89 def setSolutionFixedBottom(self, p0, pb, f0, fb, k): 90 d=self.dom.getDim() 91 x=self.dom.getX()[d-1] 92 xb=inf(x) 93 u=Vector(0.,Solution(self.dom))+kronecker(d)[d-1]*((f0+fb)/2.-k*(pb-p0)/xb) 94 p=1./k*((fb-f0)/(xb*2.)* x**2 - (fb-f0)/2.*x)+(pb-p0)/xb*x + p0 95 f= ((fb-f0)/xb* x + f0)*kronecker(Function(self.dom))[d-1] 96 return u,p,f 97 98 def testConstF_FixedBottom_smallK(self): 99 k=1.e-8 100 mp=self.getScalarMask(include_bottom=True) 101 mv=self.getVectorMask(include_bottom=False) 102 u_ref,p_ref,f=self.setSolutionFixedBottom(p0=2500,pb=4000.,f0=10.,fb=10.,k=k) 103 p=p_ref*mp 104 u=u_ref*mv 105 df=DarcyFlow(self.dom, verbose=VERBOSE, solver=self.SOLVER) 106 df.setValue(g=f, 107 location_of_fixed_pressure=mp, 108 location_of_fixed_flux=mv, 109 permeability=Scalar(k,Function(self.dom))) 110 v,p=df.solve(u_ref,p) 111 112 113 self.assertTrue(Lsup(v-u_ref)