/[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 6651 - (show annotations)
Wed Feb 7 02:12:08 2018 UTC (20 months, 1 week ago) by jfenwick
File MIME type: text/x-python
File size: 4658 byte(s)
Make everyone sad by touching all the files

Copyright dates update

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

  ViewVC Help
Powered by ViewVC 1.1.26