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

Revision 873 - (hide annotations)
Mon Oct 16 04:07:33 2006 UTC (16 years, 1 month ago) by gross
File MIME type: text/x-python
File size: 4368 byte(s)
```uszawa scheme runs with variable viscosity
```
 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 23 Primary Business: Queensland, Australia""" 24 __license__="""Licensed under the Open Software License version 3.0 25 26 __url__= 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=60 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=0. 114 prop.initialize(fixed_u_mask=mask,eta=10.) 115 u,p=prop.solve(u0,p0,relaxation=1.,iter_max=50,tolerance=0.01) 116 saveVTK("stokes.xml",u=u,p=p,m=mask,u0=u0) 117 118 eta=whereNegative(x[1]-0.5)*1.e6+whereNonNegative(x[1]-0.5) 119 prop.initialize(fixed_u_mask=mask,eta=eta) 120 u,p=prop.solve(u0,p0,relaxation=1.,iter_max=50,tolerance=0.01) 121 saveVTK("stokes.xml",u=u,p=p,m=mask,u0=u0) 122 123 # vim: expandtab shiftwidth=4:

 ViewVC Help Powered by ViewVC 1.1.26