/[escript]/trunk/finley/test/python/stokes_problems.py
ViewVC logotype

Contents of /trunk/finley/test/python/stokes_problems.py

Parent Directory Parent Directory | Revision Log Revision Log


Revision 4657 - (show annotations)
Thu Feb 6 06:12:20 2014 UTC (5 years, 4 months ago) by jfenwick
File MIME type: text/x-python
File size: 4621 byte(s)
I changed some files.
Updated copyright notices, added GeoComp.



1
2 ##############################################################################
3 #
4 # Copyright (c) 2003-2014 by University of Queensland
5 # http://www.uq.edu.au
6 #
7 # Primary Business: Queensland, Australia
8 # Licensed under the Open Software License version 3.0
9 # http://www.opensource.org/licenses/osl-3.0.php
10 #
11 # Development until 2012 by Earth Systems Science Computational Center (ESSCC)
12 # Development 2012-2013 by School of Earth Sciences
13 # Development from 2014 by Centre for Geoscience Computing (GeoComp)
14 #
15 ##############################################################################
16
17 __copyright__="""Copyright (c) 2003-2014 by University of Queensland
18 http://www.uq.edu.au
19 Primary Business: Queensland, Australia"""
20 __license__="""Licensed under the Open Software License version 3.0
21 http://www.opensource.org/licenses/osl-3.0.php"""
22 __url__="https://launchpad.net/escript-finley"
23
24 """
25 solvers for the stokes problem
26
27 :var __author__: name of author
28 :var __copyright__: copyrights
29 :var __license__: licence agreement
30 :var __url__: url entry point on documentation
31 :var __version__: version
32 :var __date__: date of the version
33 """
34
35 __author__="Lutz Gross, l.gross@uq.edu.au"
36
37 from esys.escript import *
38 from esys.escript.pdetools import SaddlePointProblem
39 from esys.escript.linearPDEs import LinearPDE
40 from esys.finley import Rectangle
41 from esys.weipa import saveVTK
42
43 class SimpleStokesProblem(SaddlePointProblem):
44 """
45 simple example of saddle point problem
46 """
47 def __init__(self,domain):
48 super(SimpleStokesProblem, self).__init__(self)
49
50 self.__pde_u=LinearPDE(domain)
51 self.__pde_u.setSymmetryOn()
52 self.__pde_u.setValue(A=identityTensor4(dom))
53
54 self.__pde_p=LinearPDE(domain)
55 self.__pde_p.setReducedOrderOn()
56 self.__pde_p.setSymmetryOn()
57 self.__pde_p.setValue(D=1.)
58
59 def initialize(self,f=Data(),fixed_u_mask=Data()):
60 self.__pde_u.setValue(q=fixed_u_mask,Y=f)
61 def inner(self,p0,p1):
62 return integrate(p0*p1,Function(self.__pde_p.getDomain()))
63
64 def solve_f(self,u,p,tol=1.e-8):
65 self.__pde_u.setTolerance(tol)
66 self.__pde_u.setValue(X=grad(u)+p*kronecker(self.__pde_u.getDomain()))
67 return self.__pde_u.getSolution()
68 def solve_g(self,u,tol=1.e-8):
69 self.__pde_p.setTolerance(tol)
70 self.__pde_p.setValue(X=-u)
71 dp=self.__pde_p.getSolution()
72 return dp
73
74 class StokesProblem(SaddlePointProblem):
75 """
76 simple example of saddle point problem
77 """
78 def __init__(self,domain):
79 super(StokesProblem, self).__init__(self)
80 self.domain=domain
81 self.__pde_u=LinearPDE(domain,numEquations=self.domain.getDim(),numSolutions=self.domain.getDim())
82 self.__pde_u.setSymmetryOn()
83
84 self.__pde_p=LinearPDE(domain)
85 self.__pde_p.setReducedOrderOn()
86 self.__pde_p.setSymmetryOn()
87
88 def initialize(self,f=Data(),fixed_u_mask=Data(),eta=1):
89 self.eta=eta
90 A =self.__pde_u.createCoefficientOfGeneralPDE("A")
91 for i in range(self.domain.getDim()):
92 for j in range(self.domain.getDim()):
93 A[i,j,j,i] += self.eta
94 A[i,j,i,j] += self.eta
95 self.__pde_p.setValue(D=1./self.eta)
96 self.__pde_u.setValue(A=A,q=fixed_u_mask,Y=f)
97
98 def inner(self,p0,p1):
99 return integrate(p0*p1,Function(self.__pde_p.getDomain()))
100
101 def solve_f(self,u,p,tol=1.e-8):
102 self.__pde_u.setTolerance(tol)
103 g=grad(u)
104 self.__pde_u.setValue(X=self.eta*symmetric(g)+p*kronecker(self.__pde_u.getDomain()))
105 return self.__pde_u.getSolution()
106
107 def solve_g(self,u,tol=1.e-8):
108 self.__pde_p.setTolerance(tol)
109 self.__pde_p.setValue(X=-u)
110 dp=self.__pde_p.getSolution()
111 return dp
112
113 NE=50
114 dom=Rectangle(NE,NE,order=2)
115 # prop=SimpleStokesProblem(dom)
116 prop=StokesProblem(dom)
117 x=dom.getX()
118 mask=(whereZero(x[0])+whereZero(x[0]-1.)+whereZero(x[1]-1.))*unitVector(0,dom)+(whereZero(x[1]-1.)+whereZero(x[1]))*unitVector(1,dom)
119 u0=Vector(0.,Solution(dom))
120 u0[0]=x[1]*whereZero(x[1]-1.)
121 p0=Scalar(0,ReducedSolution(dom))
122 # prop.initialize(fixed_u_mask=mask)
123 prop.initialize(fixed_u_mask=mask,eta=10.)
124 u,p=prop.solve(u0,p0,tolerance=0.01)
125 # saveVTK("stokes.vtu",u=u,p=p,m=mask,u0=u0)
126
127 eta=whereNegative(x[1]-0.5)*1.e6+whereNonNegative(x[1]-0.5)
128 prop.initialize(fixed_u_mask=mask,eta=eta)
129 u,p=prop.solve(u0,p0,tolerance=0.01,tolerance_u=0.1,accepted_reduction=0.8)
130 saveVTK("stokes.vtu",u=u,p=p,m=mask,u0=u0)
131
132 # vim: expandtab shiftwidth=4:

  ViewVC Help
Powered by ViewVC 1.1.26