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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 6651 - (hide annotations)
Wed Feb 7 02:12:08 2018 UTC (3 years, 4 months 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 ksteube 1809
2 jfenwick 3981 ##############################################################################
3 ksteube 1312 #
4 jfenwick 6651 # Copyright (c) 2003-2018 by The University of Queensland
5 jfenwick 3981 # http://www.uq.edu.au
6 ksteube 1312 #
7 ksteube 1809 # Primary Business: Queensland, Australia
8 jfenwick 6112 # Licensed under the Apache License, version 2.0
9     # http://www.apache.org/licenses/LICENSE-2.0
10 ksteube 1312 #
11 jfenwick 3981 # Development until 2012 by Earth Systems Science Computational Center (ESSCC)
12 jfenwick 4657 # Development 2012-2013 by School of Earth Sciences
13     # Development from 2014 by Centre for Geoscience Computing (GeoComp)
14 jfenwick 3981 #
15     ##############################################################################
16 gross 873
17 sshaw 5706 from __future__ import print_function, division
18    
19 jfenwick 6651 __copyright__="""Copyright (c) 2003-2018 by The University of Queensland
20 jfenwick 3981 http://www.uq.edu.au
21 ksteube 1809 Primary Business: Queensland, Australia"""
22 jfenwick 6112 __license__="""Licensed under the Apache License, version 2.0
23     http://www.apache.org/licenses/LICENSE-2.0"""
24 jfenwick 2344 __url__="https://launchpad.net/escript-finley"
25 ksteube 1809
26 gross 873 """
27 gross 883 solvers for the stokes problem
28 gross 873
29 jfenwick 2625 :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 gross 873 """
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 caltinay 3346 from esys.weipa import saveVTK
44 gross 873
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 gross 878 NE=50
116 gross 873 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 gross 877 p0=Scalar(0,ReducedSolution(dom))
124     # prop.initialize(fixed_u_mask=mask)
125 gross 873 prop.initialize(fixed_u_mask=mask,eta=10.)
126 gross 877 u,p=prop.solve(u0,p0,tolerance=0.01)
127 caltinay 2534 # saveVTK("stokes.vtu",u=u,p=p,m=mask,u0=u0)
128 gross 873
129     eta=whereNegative(x[1]-0.5)*1.e6+whereNonNegative(x[1]-0.5)
130     prop.initialize(fixed_u_mask=mask,eta=eta)
131 gross 878 u,p=prop.solve(u0,p0,tolerance=0.01,tolerance_u=0.1,accepted_reduction=0.8)
132 caltinay 2534 saveVTK("stokes.vtu",u=u,p=p,m=mask,u0=u0)
133 gross 873
134     # vim: expandtab shiftwidth=4:

  ViewVC Help
Powered by ViewVC 1.1.26