/[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 877 - (hide annotations)
Wed Oct 25 03:06:58 2006 UTC (14 years, 9 months ago) by gross
File MIME type: text/x-python
File size: 4412 byte(s)
backtraking in the saddlepoint problem (not perfect yet)
1 gross 873 # $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 gross 877 NE=1
106 gross 873 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 gross 877 p0=Scalar(0,ReducedSolution(dom))
114     # prop.initialize(fixed_u_mask=mask)
115 gross 873 prop.initialize(fixed_u_mask=mask,eta=10.)
116 gross 877 u,p=prop.solve(u0,p0,tolerance=0.01)
117     # saveVTK("stokes.xml",u=u,p=p,m=mask,u0=u0)
118 gross 873
119     eta=whereNegative(x[1]-0.5)*1.e6+whereNonNegative(x[1]-0.5)
120     prop.initialize(fixed_u_mask=mask,eta=eta)
121 gross 877 u,p=prop.solve(u0,p0,tolerance=0.01,tolerance_u=0.1,relaxation=1.)
122 gross 873 saveVTK("stokes.xml",u=u,p=p,m=mask,u0=u0)
123    
124     # vim: expandtab shiftwidth=4:

  ViewVC Help
Powered by ViewVC 1.1.26