36 |
from escript import * |
from escript import * |
37 |
import util |
import util |
38 |
from linearPDEs import LinearPDE |
from linearPDEs import LinearPDE |
39 |
from pdetools import HomogeneousSaddlePointProblem |
from pdetools import HomogeneousSaddlePointProblem,Projector |
40 |
|
|
41 |
class StokesProblemCartesian(HomogeneousSaddlePointProblem): |
class StokesProblemCartesian(HomogeneousSaddlePointProblem): |
42 |
""" |
""" |
63 |
self.vol=util.integrate(1.,Function(self.domain)) |
self.vol=util.integrate(1.,Function(self.domain)) |
64 |
self.__pde_u=LinearPDE(domain,numEquations=self.domain.getDim(),numSolutions=self.domain.getDim()) |
self.__pde_u=LinearPDE(domain,numEquations=self.domain.getDim(),numSolutions=self.domain.getDim()) |
65 |
self.__pde_u.setSymmetryOn() |
self.__pde_u.setSymmetryOn() |
66 |
self.__pde_u.setSolverMethod(preconditioner=LinearPDE.ILU0) |
self.__pde_u.setSolverMethod(preconditioner=LinearPDE.RILU) |
67 |
|
|
68 |
self.__pde_prec=LinearPDE(domain) |
self.__pde_prec=LinearPDE(domain) |
69 |
self.__pde_prec.setReducedOrderOn() |
self.__pde_prec.setReducedOrderOn() |
82 |
for j in range(self.domain.getDim()): |
for j in range(self.domain.getDim()): |
83 |
A[i,j,j,i] += 1. |
A[i,j,j,i] += 1. |
84 |
A[i,j,i,j] += 1. |
A[i,j,i,j] += 1. |
85 |
self.__pde_prec.setValue(D=1./self.eta) |
self.__pde_prec.setValue(D=1/eta) #1./self.eta |
86 |
self.__pde_u.setValue(A=A*self.eta,q=fixed_u_mask,Y=f,y=surface_stress) |
self.__pde_u.setValue(A=A*self.eta,q=fixed_u_mask,Y=f,y=surface_stress) |
87 |
|
|
88 |
def B(self,arg): |
def B(self,arg): |
96 |
s1=util.interpolate(p1,Function(self.domain)) |
s1=util.interpolate(p1,Function(self.domain)) |
97 |
return util.integrate(s0*s1) |
return util.integrate(s0*s1) |
98 |
|
|
99 |
|
def inner_a(self,a0,a1): |
100 |
|
p0=util.interpolate(a0[1],Function(self.domain)) |
101 |
|
p1=util.interpolate(a1[1],Function(self.domain)) |
102 |
|
alfa=(1/self.vol)*util.integrate(p0) |
103 |
|
beta=(1/self.vol)*util.integrate(p1) |
104 |
|
v0=util.grad(a0[0]) |
105 |
|
v1=util.grad(a1[0]) |
106 |
|
#print "NORM",alfa,beta,util.integrate((p0-alfa)*(p1-beta))+util.integrate(util.inner(v0,v1)) |
107 |
|
return util.integrate((p0-alfa)*(p1-beta)+((1/self.eta)**2)*util.inner(v0,v1)) |
108 |
|
|
109 |
|
|
110 |
def getStress(self,u): |
def getStress(self,u): |
111 |
mg=util.grad(u) |
mg=util.grad(u) |
112 |
return 2.*self.eta*util.symmetric(mg) |
return 2.*self.eta*util.symmetric(mg) |
119 |
self.__pde_u.setValue(X=-self.getStress(u)-p*util.kronecker(self.domain)) |
self.__pde_u.setValue(X=-self.getStress(u)-p*util.kronecker(self.domain)) |
120 |
return self.__pde_u.getSolution(verbose=self.show_details) |
return self.__pde_u.getSolution(verbose=self.show_details) |
121 |
|
|
122 |
|
|
123 |
def solve_prec(self,p): |
def solve_prec(self,p): |
124 |
|
#proj=Projector(domain=self.domain, reduce = True, fast=False) |
125 |
self.__pde_prec.setTolerance(self.getSubProblemTolerance()) |
self.__pde_prec.setTolerance(self.getSubProblemTolerance()) |
126 |
self.__pde_prec.setValue(Y=p) |
self.__pde_prec.setValue(Y=p) |
127 |
q=self.__pde_prec.getSolution(verbose=self.show_details) |
q=self.__pde_prec.getSolution(verbose=self.show_details) |
128 |
|
q0=util.interpolate(q,Function(self.domain)) |
129 |
|
q-=(1/self.vol)*util.integrate(q0) |
130 |
return q |
return q |
131 |
|
|
132 |
def stoppingcriterium(self,Bv,v,p): |
def stoppingcriterium(self,Bv,v,p): |
133 |
n_r=util.sqrt(self.inner(Bv,Bv)) |
n_r=util.sqrt(self.inner(Bv,Bv)) |
134 |
n_v=util.Lsup(v) |
n_v=util.Lsup(v) |