1 |
# $Id: pdetools.py 867 2006-10-09 06:50:09Z gross $ |
2 |
|
3 |
""" |
4 |
Provides some tools related to PDEs. |
5 |
|
6 |
Currently includes: |
7 |
- Projector - to project a discontinuous |
8 |
- Locator - to trace values in data objects at a certain location |
9 |
- TimeIntegrationManager - to handel extraplotion in time |
10 |
- SaddlePointProblem - solver for Saddle point problems using the inexact uszawa scheme |
11 |
|
12 |
@var __author__: name of author |
13 |
@var __copyright__: copyrights |
14 |
@var __license__: licence agreement |
15 |
@var __url__: url entry point on documentation |
16 |
@var __version__: version |
17 |
@var __date__: date of the version |
18 |
""" |
19 |
|
20 |
__author__="Lutz Gross, l.gross@uq.edu.au" |
21 |
__copyright__=""" Copyright (c) 2006 by ACcESS MNRF |
22 |
http://www.access.edu.au |
23 |
Primary Business: Queensland, Australia""" |
24 |
__license__="""Licensed under the Open Software License version 3.0 |
25 |
http://www.opensource.org/licenses/osl-3.0.php""" |
26 |
__url__="http://www.iservo.edu.au/esys" |
27 |
__version__="$Revision$" |
28 |
__date__="$Date:$" |
29 |
|
30 |
from esys.escript import * |
31 |
from esys.escript.pdetools import SaddlePointProblem |
32 |
from esys.escript.linearPDEs import LinearPDE |
33 |
from esys.finley import Rectangle |
34 |
|
35 |
class SimpleStokesProblem(SaddlePointProblem): |
36 |
""" |
37 |
simple example of saddle point problem |
38 |
""" |
39 |
def __init__(self,domain): |
40 |
super(SimpleStokesProblem, self).__init__(self) |
41 |
|
42 |
self.__pde_u=LinearPDE(domain) |
43 |
self.__pde_u.setSymmetryOn() |
44 |
self.__pde_u.setValue(A=identityTensor4(dom)) |
45 |
|
46 |
self.__pde_p=LinearPDE(domain) |
47 |
self.__pde_p.setReducedOrderOn() |
48 |
self.__pde_p.setSymmetryOn() |
49 |
self.__pde_p.setValue(D=1.) |
50 |
|
51 |
def initialize(self,f=Data(),fixed_u_mask=Data()): |
52 |
self.__pde_u.setValue(q=fixed_u_mask,Y=f) |
53 |
def inner(self,p0,p1): |
54 |
return integrate(p0*p1,Function(self.__pde_p.getDomain())) |
55 |
|
56 |
def solve_f(self,u,p,tol=1.e-8): |
57 |
self.__pde_u.setTolerance(tol) |
58 |
self.__pde_u.setValue(X=grad(u)+p*kronecker(self.__pde_u.getDomain())) |
59 |
return self.__pde_u.getSolution() |
60 |
def solve_g(self,u,tol=1.e-8): |
61 |
self.__pde_p.setTolerance(tol) |
62 |
self.__pde_p.setValue(X=-u) |
63 |
dp=self.__pde_p.getSolution() |
64 |
return dp |
65 |
|
66 |
class StokesProblem(SaddlePointProblem): |
67 |
""" |
68 |
simple example of saddle point problem |
69 |
""" |
70 |
def __init__(self,domain): |
71 |
super(StokesProblem, self).__init__(self) |
72 |
self.domain=domain |
73 |
self.__pde_u=LinearPDE(domain,numEquations=self.domain.getDim(),numSolutions=self.domain.getDim()) |
74 |
self.__pde_u.setSymmetryOn() |
75 |
|
76 |
self.__pde_p=LinearPDE(domain) |
77 |
self.__pde_p.setReducedOrderOn() |
78 |
self.__pde_p.setSymmetryOn() |
79 |
|
80 |
def initialize(self,f=Data(),fixed_u_mask=Data(),eta=1): |
81 |
self.eta=eta |
82 |
A =self.__pde_u.createCoefficientOfGeneralPDE("A") |
83 |
for i in range(self.domain.getDim()): |
84 |
for j in range(self.domain.getDim()): |
85 |
A[i,j,j,i] += self.eta |
86 |
A[i,j,i,j] += self.eta |
87 |
self.__pde_p.setValue(D=1./self.eta) |
88 |
self.__pde_u.setValue(A=A,q=fixed_u_mask,Y=f) |
89 |
|
90 |
def inner(self,p0,p1): |
91 |
return integrate(p0*p1,Function(self.__pde_p.getDomain())) |
92 |
|
93 |
def solve_f(self,u,p,tol=1.e-8): |
94 |
self.__pde_u.setTolerance(tol) |
95 |
g=grad(u) |
96 |
self.__pde_u.setValue(X=self.eta*symmetric(g)+p*kronecker(self.__pde_u.getDomain())) |
97 |
return self.__pde_u.getSolution() |
98 |
|
99 |
def solve_g(self,u,tol=1.e-8): |
100 |
self.__pde_p.setTolerance(tol) |
101 |
self.__pde_p.setValue(X=-u) |
102 |
dp=self.__pde_p.getSolution() |
103 |
return dp |
104 |
|
105 |
NE=1 |
106 |
dom=Rectangle(NE,NE,order=2) |
107 |
# prop=SimpleStokesProblem(dom) |
108 |
prop=StokesProblem(dom) |
109 |
x=dom.getX() |
110 |
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) |
111 |
u0=Vector(0.,Solution(dom)) |
112 |
u0[0]=x[1]*whereZero(x[1]-1.) |
113 |
p0=Scalar(0,ReducedSolution(dom)) |
114 |
# prop.initialize(fixed_u_mask=mask) |
115 |
prop.initialize(fixed_u_mask=mask,eta=10.) |
116 |
u,p=prop.solve(u0,p0,tolerance=0.01) |
117 |
# saveVTK("stokes.xml",u=u,p=p,m=mask,u0=u0) |
118 |
|
119 |
eta=whereNegative(x[1]-0.5)*1.e6+whereNonNegative(x[1]-0.5) |
120 |
prop.initialize(fixed_u_mask=mask,eta=eta) |
121 |
u,p=prop.solve(u0,p0,tolerance=0.01,tolerance_u=0.1,relaxation=1.) |
122 |
saveVTK("stokes.xml",u=u,p=p,m=mask,u0=u0) |
123 |
|
124 |
# vim: expandtab shiftwidth=4: |