|
# $Id$ |
|
1 |
|
|
2 |
## @file linearPDEs.py |
######################################################## |
3 |
|
# |
4 |
|
# Copyright (c) 2003-2008 by University of Queensland |
5 |
|
# Earth Systems Science Computational Center (ESSCC) |
6 |
|
# http://www.uq.edu.au/esscc |
7 |
|
# |
8 |
|
# Primary Business: Queensland, Australia |
9 |
|
# Licensed under the Open Software License version 3.0 |
10 |
|
# http://www.opensource.org/licenses/osl-3.0.php |
11 |
|
# |
12 |
|
######################################################## |
13 |
|
|
14 |
|
__copyright__="""Copyright (c) 2003-2008 by University of Queensland |
15 |
|
Earth Systems Science Computational Center (ESSCC) |
16 |
|
http://www.uq.edu.au/esscc |
17 |
|
Primary Business: Queensland, Australia""" |
18 |
|
__license__="""Licensed under the Open Software License version 3.0 |
19 |
|
http://www.opensource.org/licenses/osl-3.0.php""" |
20 |
|
__url__="http://www.uq.edu.au/esscc/escript-finley" |
21 |
|
|
22 |
""" |
""" |
23 |
@brief Functions and classes for linear PDEs |
The module provides an interface to define and solve linear partial |
24 |
|
differential equations (PDEs) and Transport problems within L{escript}. |
25 |
|
L{linearPDEs} does not provide any solver capabilities in itself but hands the |
26 |
|
PDE over to the PDE solver library defined through the L{Domain<escript.Domain>} |
27 |
|
of the PDE. The general interface is provided through the L{LinearPDE} class. |
28 |
|
L{TransportProblem} provides an interface to initial value problems dominated |
29 |
|
by its advective terms. |
30 |
|
|
31 |
|
@var __author__: name of author |
32 |
|
@var __copyright__: copyrights |
33 |
|
@var __license__: licence agreement |
34 |
|
@var __url__: url entry point on documentation |
35 |
|
@var __version__: version |
36 |
|
@var __date__: date of the version |
37 |
""" |
""" |
38 |
|
|
39 |
|
import math |
40 |
import escript |
import escript |
41 |
import util |
import util |
42 |
import numarray |
import numarray |
43 |
|
|
44 |
def identifyDomain(domain=None,data={}): |
__author__="Lutz Gross, l.gross@uq.edu.au" |
|
""" |
|
|
@brief Return the Domain which is equal to the input domain (if not None) |
|
|
and is the domain of all Data objects in the dictionary data. |
|
|
An exception is raised if this is not possible |
|
|
|
|
|
@param domain |
|
|
@param data |
|
|
""" |
|
|
# get the domain used by any Data object in the list data: |
|
|
data_domain=None |
|
|
for d in data.itervalues(): |
|
|
if isinstance(d,escript.Data): |
|
|
if not d.isEmpty(): data_domain=d.getDomain() |
|
|
# check if domain and data_domain are identical? |
|
|
if domain == None: |
|
|
if data_domain == None: |
|
|
raise ValueError,"Undefined PDE domain. Specify a domain or use a Data class object as coefficient" |
|
|
else: |
|
|
if data_domain == None: |
|
|
data_domain=domain |
|
|
else: |
|
|
if not data_domain == domain: |
|
|
raise ValueError,"Domain of coefficients doesnot match specified domain" |
|
|
# now we check if all Data class object coefficients are defined on data_domain: |
|
|
for i,d in data.iteritems(): |
|
|
if isinstance(d,escript.Data): |
|
|
if not d.isEmpty(): |
|
|
if not data_domain==d.getDomain(): |
|
|
raise ValueError,"Illegal domain for coefficient %s."%i |
|
|
# done: |
|
|
return data_domain |
|
|
|
|
|
def identifyNumEquationsAndSolutions(dim,coef={}): |
|
|
# get number of equations and number of unknowns: |
|
|
numEquations=0 |
|
|
numSolutions=0 |
|
|
for i in coef.iterkeys(): |
|
|
if not coef[i].isEmpty(): |
|
|
res=_PDECoefficientTypes[i].estimateNumEquationsAndNumSolutions(coef[i].getShape(),dim) |
|
|
if res==None: |
|
|
raise ValueError,"Illegal shape %s of coefficient %s"%(coef[i].getShape().__str__(),i) |
|
|
else: |
|
|
numEquations=max(numEquations,res[0]) |
|
|
numSolutions=max(numSolutions,res[1]) |
|
|
return numEquations,numSolutions |
|
45 |
|
|
46 |
|
|
47 |
def _CompTuple2(t1,t2): |
class IllegalCoefficient(ValueError): |
48 |
""" |
""" |
49 |
@brief |
Exception that is raised if an illegal coefficient of the general or |
50 |
|
particular PDE is requested. |
51 |
|
""" |
52 |
|
pass |
53 |
|
|
54 |
|
class IllegalCoefficientValue(ValueError): |
55 |
|
""" |
56 |
|
Exception that is raised if an incorrect value for a coefficient is used. |
57 |
|
""" |
58 |
|
pass |
59 |
|
|
60 |
@param t1 |
class IllegalCoefficientFunctionSpace(ValueError): |
61 |
@param t2 |
""" |
62 |
|
Exception that is raised if an incorrect function space for a coefficient |
63 |
|
is used. |
64 |
|
""" |
65 |
|
|
66 |
|
class UndefinedPDEError(ValueError): |
67 |
|
""" |
68 |
|
Exception that is raised if a PDE is not fully defined yet. |
69 |
""" |
""" |
70 |
dif=t1[0]+t1[1]-(t2[0]+t2[1]) |
pass |
|
if dif<0: return 1 |
|
|
elif dif>0: return -1 |
|
|
else: return 0 |
|
71 |
|
|
72 |
class PDECoefficientType: |
class PDECoef(object): |
73 |
""" |
""" |
74 |
@brief |
A class for describing a PDE coefficient. |
75 |
|
|
76 |
|
@cvar INTERIOR: indicator that coefficient is defined on the interior of |
77 |
|
the PDE domain |
78 |
|
@cvar BOUNDARY: indicator that coefficient is defined on the boundary of |
79 |
|
the PDE domain |
80 |
|
@cvar CONTACT: indicator that coefficient is defined on the contact region |
81 |
|
within the PDE domain |
82 |
|
@cvar INTERIOR_REDUCED: indicator that coefficient is defined on the |
83 |
|
interior of the PDE domain using a reduced |
84 |
|
integration order |
85 |
|
@cvar BOUNDARY_REDUCED: indicator that coefficient is defined on the |
86 |
|
boundary of the PDE domain using a reduced |
87 |
|
integration order |
88 |
|
@cvar CONTACT_REDUCED: indicator that coefficient is defined on the contact |
89 |
|
region within the PDE domain using a reduced |
90 |
|
integration order |
91 |
|
@cvar SOLUTION: indicator that coefficient is defined through a solution of |
92 |
|
the PDE |
93 |
|
@cvar REDUCED: indicator that coefficient is defined through a reduced |
94 |
|
solution of the PDE |
95 |
|
@cvar BY_EQUATION: indicator that the dimension of the coefficient shape is |
96 |
|
defined by the number of PDE equations |
97 |
|
@cvar BY_SOLUTION: indicator that the dimension of the coefficient shape is |
98 |
|
defined by the number of PDE solutions |
99 |
|
@cvar BY_DIM: indicator that the dimension of the coefficient shape is |
100 |
|
defined by the spatial dimension |
101 |
|
@cvar OPERATOR: indicator that the the coefficient alters the operator of |
102 |
|
the PDE |
103 |
|
@cvar RIGHTHANDSIDE: indicator that the the coefficient alters the right |
104 |
|
hand side of the PDE |
105 |
|
@cvar BOTH: indicator that the the coefficient alters the operator as well |
106 |
|
as the right hand side of the PDE |
107 |
|
|
108 |
""" |
""" |
|
# identifier for location of Data objects defining coefficients |
|
109 |
INTERIOR=0 |
INTERIOR=0 |
110 |
BOUNDARY=1 |
BOUNDARY=1 |
111 |
CONTACT=2 |
CONTACT=2 |
112 |
CONTINUOUS=3 |
SOLUTION=3 |
113 |
# identifier in the pattern of coefficients: |
REDUCED=4 |
114 |
# the pattern is a tuple of EQUATION,SOLUTION,DIM where DIM represents the spatial dimension, EQUATION the number of equations and SOLUTION the |
BY_EQUATION=5 |
115 |
# number of unknowns. |
BY_SOLUTION=6 |
116 |
EQUATION=3 |
BY_DIM=7 |
117 |
SOLUTION=4 |
OPERATOR=10 |
118 |
DIM=5 |
RIGHTHANDSIDE=11 |
119 |
# indicator for what is altered if the coefficient is altered: |
BOTH=12 |
120 |
OPERATOR=5 |
INTERIOR_REDUCED=13 |
121 |
RIGHTHANDSIDE=6 |
BOUNDARY_REDUCED=14 |
122 |
BOTH=7 |
CONTACT_REDUCED=15 |
123 |
def __init__(self,where,pattern,altering): |
|
124 |
""" |
def __init__(self, where, pattern, altering): |
125 |
@brief Initialise a PDE Coefficient type |
""" |
126 |
|
Initialises a PDE coefficient type. |
127 |
|
|
128 |
|
@param where: describes where the coefficient lives |
129 |
|
@type where: one of L{INTERIOR}, L{BOUNDARY}, L{CONTACT}, L{SOLUTION}, |
130 |
|
L{REDUCED}, L{INTERIOR_REDUCED}, L{BOUNDARY_REDUCED}, |
131 |
|
L{CONTACT_REDUCED} |
132 |
|
@param pattern: describes the shape of the coefficient and how the shape |
133 |
|
is built for a given spatial dimension and numbers of |
134 |
|
equations and solutions in then PDE. For instance, |
135 |
|
(L{BY_EQUATION},L{BY_SOLUTION},L{BY_DIM}) describes a |
136 |
|
rank 3 coefficient which is instantiated as shape (3,2,2) |
137 |
|
in case of three equations and two solution components |
138 |
|
on a 2-dimensional domain. In the case of single equation |
139 |
|
and a single solution component the shape components |
140 |
|
marked by L{BY_EQUATION} or L{BY_SOLUTION} are dropped. |
141 |
|
In this case the example would be read as (2,). |
142 |
|
@type pattern: C{tuple} of L{BY_EQUATION}, L{BY_SOLUTION}, L{BY_DIM} |
143 |
|
@param altering: indicates what part of the PDE is altered if the |
144 |
|
coefficient is altered |
145 |
|
@type altering: one of L{OPERATOR}, L{RIGHTHANDSIDE}, L{BOTH} |
146 |
""" |
""" |
147 |
|
super(PDECoef, self).__init__() |
148 |
self.what=where |
self.what=where |
149 |
self.pattern=pattern |
self.pattern=pattern |
150 |
self.altering=altering |
self.altering=altering |
151 |
|
self.resetValue() |
152 |
|
|
153 |
def getFunctionSpace(self,domain): |
def resetValue(self): |
154 |
""" |
""" |
155 |
@brief defines the FunctionSpace of the coefficient on the domain |
Resets the coefficient value to the default. |
|
|
|
|
@param domain |
|
156 |
""" |
""" |
157 |
if self.what==self.INTERIOR: return escript.Function(domain) |
self.value=escript.Data() |
158 |
elif self.what==self.BOUNDARY: return escript.FunctionOnBoundary(domain) |
|
159 |
elif self.what==self.CONTACT: return escript.FunctionOnContactZero(domain) |
def getFunctionSpace(self,domain,reducedEquationOrder=False,reducedSolutionOrder=False): |
160 |
elif self.what==self.CONTINUOUS: return escript.ContinuousFunction(domain) |
""" |
161 |
|
Returns the L{FunctionSpace<escript.FunctionSpace>} of the coefficient. |
162 |
|
|
163 |
|
@param domain: domain on which the PDE uses the coefficient |
164 |
|
@type domain: L{Domain<escript.Domain>} |
165 |
|
@param reducedEquationOrder: True to indicate that reduced order is used |
166 |
|
to represent the equation |
167 |
|
@type reducedEquationOrder: C{bool} |
168 |
|
@param reducedSolutionOrder: True to indicate that reduced order is used |
169 |
|
to represent the solution |
170 |
|
@type reducedSolutionOrder: C{bool} |
171 |
|
@return: L{FunctionSpace<escript.FunctionSpace>} of the coefficient |
172 |
|
@rtype: L{FunctionSpace<escript.FunctionSpace>} |
173 |
|
""" |
174 |
|
if self.what==self.INTERIOR: |
175 |
|
return escript.Function(domain) |
176 |
|
elif self.what==self.INTERIOR_REDUCED: |
177 |
|
return escript.ReducedFunction(domain) |
178 |
|
elif self.what==self.BOUNDARY: |
179 |
|
return escript.FunctionOnBoundary(domain) |
180 |
|
elif self.what==self.BOUNDARY_REDUCED: |
181 |
|
return escript.ReducedFunctionOnBoundary(domain) |
182 |
|
elif self.what==self.CONTACT: |
183 |
|
return escript.FunctionOnContactZero(domain) |
184 |
|
elif self.what==self.CONTACT_REDUCED: |
185 |
|
return escript.ReducedFunctionOnContactZero(domain) |
186 |
|
elif self.what==self.SOLUTION: |
187 |
|
if reducedEquationOrder and reducedSolutionOrder: |
188 |
|
return escript.ReducedSolution(domain) |
189 |
|
else: |
190 |
|
return escript.Solution(domain) |
191 |
|
elif self.what==self.REDUCED: |
192 |
|
return escript.ReducedSolution(domain) |
193 |
|
|
194 |
|
def getValue(self): |
195 |
|
""" |
196 |
|
Returns the value of the coefficient. |
197 |
|
|
198 |
|
@return: value of the coefficient |
199 |
|
@rtype: L{Data<escript.Data>} |
200 |
|
""" |
201 |
|
return self.value |
202 |
|
|
203 |
|
def setValue(self,domain,numEquations=1,numSolutions=1,reducedEquationOrder=False,reducedSolutionOrder=False,newValue=None): |
204 |
|
""" |
205 |
|
Sets the value of the coefficient to a new value. |
206 |
|
|
207 |
|
@param domain: domain on which the PDE uses the coefficient |
208 |
|
@type domain: L{Domain<escript.Domain>} |
209 |
|
@param numEquations: number of equations of the PDE |
210 |
|
@type numEquations: C{int} |
211 |
|
@param numSolutions: number of components of the PDE solution |
212 |
|
@type numSolutions: C{int} |
213 |
|
@param reducedEquationOrder: True to indicate that reduced order is used |
214 |
|
to represent the equation |
215 |
|
@type reducedEquationOrder: C{bool} |
216 |
|
@param reducedSolutionOrder: True to indicate that reduced order is used |
217 |
|
to represent the solution |
218 |
|
@type reducedSolutionOrder: C{bool} |
219 |
|
@param newValue: number of components of the PDE solution |
220 |
|
@type newValue: any object that can be converted into a |
221 |
|
L{Data<escript.Data>} object with the appropriate shape |
222 |
|
and L{FunctionSpace<escript.FunctionSpace>} |
223 |
|
@raise IllegalCoefficientValue: if the shape of the assigned value does |
224 |
|
not match the shape of the coefficient |
225 |
|
@raise IllegalCoefficientFunctionSpace: if unable to interpolate value |
226 |
|
to appropriate function space |
227 |
|
""" |
228 |
|
if newValue==None: |
229 |
|
newValue=escript.Data() |
230 |
|
elif isinstance(newValue,escript.Data): |
231 |
|
if not newValue.isEmpty(): |
232 |
|
if not newValue.getFunctionSpace() == self.getFunctionSpace(domain,reducedEquationOrder,reducedSolutionOrder): |
233 |
|
try: |
234 |
|
newValue=escript.Data(newValue,self.getFunctionSpace(domain,reducedEquationOrder,reducedSolutionOrder)) |
235 |
|
except: |
236 |
|
raise IllegalCoefficientFunctionSpace,"Unable to interpolate coefficient to function space %s"%self.getFunctionSpace(domain) |
237 |
|
else: |
238 |
|
newValue=escript.Data(newValue,self.getFunctionSpace(domain,reducedEquationOrder,reducedSolutionOrder)) |
239 |
|
if not newValue.isEmpty(): |
240 |
|
if not self.getShape(domain,numEquations,numSolutions)==newValue.getShape(): |
241 |
|
raise IllegalCoefficientValue,"Expected shape of coefficient is %s but actual shape is %s."%(self.getShape(domain,numEquations,numSolutions),newValue.getShape()) |
242 |
|
self.value=newValue |
243 |
|
|
244 |
def isAlteringOperator(self): |
def isAlteringOperator(self): |
245 |
""" |
""" |
246 |
@brief return true if the operator of the PDE is changed when the coefficient is changed |
Checks if the coefficient alters the operator of the PDE. |
247 |
""" |
|
248 |
|
@return: True if the operator of the PDE is changed when the |
249 |
|
coefficient is changed |
250 |
|
@rtype: C{bool} |
251 |
|
""" |
252 |
if self.altering==self.OPERATOR or self.altering==self.BOTH: |
if self.altering==self.OPERATOR or self.altering==self.BOTH: |
253 |
return not None |
return not None |
254 |
else: |
else: |
256 |
|
|
257 |
def isAlteringRightHandSide(self): |
def isAlteringRightHandSide(self): |
258 |
""" |
""" |
259 |
@brief return true if the right hand side of the PDE is changed when the coefficient is changed |
Checks if the coefficient alters the right hand side of the PDE. |
260 |
""" |
|
261 |
|
@rtype: C{bool} |
262 |
|
@return: True if the right hand side of the PDE is changed when the |
263 |
|
coefficient is changed, C{None} otherwise. |
264 |
|
""" |
265 |
if self.altering==self.RIGHTHANDSIDE or self.altering==self.BOTH: |
if self.altering==self.RIGHTHANDSIDE or self.altering==self.BOTH: |
266 |
return not None |
return not None |
267 |
else: |
else: |
268 |
return None |
return None |
269 |
|
|
270 |
def estimateNumEquationsAndNumSolutions(self,shape=(),dim=3): |
def estimateNumEquationsAndNumSolutions(self,domain,shape=()): |
271 |
""" |
""" |
272 |
@brief tries to estimate the number of equations in a given tensor shape for a given spatial dimension dim |
Tries to estimate the number of equations and number of solutions if |
273 |
|
the coefficient has the given shape. |
274 |
|
|
275 |
@param shape |
@param domain: domain on which the PDE uses the coefficient |
276 |
@param dim |
@type domain: L{Domain<escript.Domain>} |
277 |
|
@param shape: suggested shape of the coefficient |
278 |
|
@type shape: C{tuple} of C{int} values |
279 |
|
@return: the number of equations and number of solutions of the PDE if |
280 |
|
the coefficient has given shape. If no appropriate numbers |
281 |
|
could be identified, C{None} is returned |
282 |
|
@rtype: C{tuple} of two C{int} values or C{None} |
283 |
""" |
""" |
284 |
|
dim=domain.getDim() |
285 |
if len(shape)>0: |
if len(shape)>0: |
286 |
num=max(shape)+1 |
num=max(shape)+1 |
287 |
else: |
else: |
288 |
num=1 |
num=1 |
289 |
search=[] |
search=[] |
290 |
for u in range(num): |
if self.definesNumEquation() and self.definesNumSolutions(): |
291 |
for e in range(num): |
for u in range(num): |
292 |
search.append((e,u)) |
for e in range(num): |
293 |
search.sort(_CompTuple2) |
search.append((e,u)) |
294 |
for item in search: |
search.sort(self.__CompTuple2) |
295 |
s=self.buildShape(item[0],item[1],dim) |
for item in search: |
296 |
|
s=self.getShape(domain,item[0],item[1]) |
297 |
if len(s)==0 and len(shape)==0: |
if len(s)==0 and len(shape)==0: |
298 |
return (1,1) |
return (1,1) |
299 |
else: |
else: |
300 |
if s==shape: return item |
if s==shape: return item |
301 |
return None |
elif self.definesNumEquation(): |
302 |
|
for e in range(num,0,-1): |
303 |
def buildShape(self,e=1,u=1,dim=3): |
s=self.getShape(domain,e,0) |
304 |
""" |
if len(s)==0 and len(shape)==0: |
305 |
@brief builds the required shape for a given number of equations e, number of unknowns u and spatial dimension dim |
return (1,None) |
306 |
|
else: |
307 |
|
if s==shape: return (e,None) |
308 |
|
|
309 |
@param e |
elif self.definesNumSolutions(): |
310 |
@param u |
for u in range(num,0,-1): |
311 |
@param dim |
s=self.getShape(domain,0,u) |
312 |
""" |
if len(s)==0 and len(shape)==0: |
313 |
s=() |
return (None,1) |
|
for i in self.pattern: |
|
|
if i==self.EQUATION: |
|
|
if e>1: s=s+(e,) |
|
|
elif i==self.SOLUTION: |
|
|
if u>1: s=s+(u,) |
|
314 |
else: |
else: |
315 |
s=s+(dim,) |
if s==shape: return (None,u) |
316 |
return s |
return None |
317 |
|
|
318 |
_PDECoefficientTypes={ |
def definesNumSolutions(self): |
319 |
"A" : PDECoefficientType(PDECoefficientType.INTERIOR,(PDECoefficientType.EQUATION,PDECoefficientType.DIM,PDECoefficientType.SOLUTION,PDECoefficientType.DIM),PDECoefficientType.OPERATOR), |
""" |
320 |
"B" : PDECoefficientType(PDECoefficientType.INTERIOR,(PDECoefficientType.EQUATION,PDECoefficientType.DIM,PDECoefficientType.SOLUTION),PDECoefficientType.OPERATOR), |
Checks if the coefficient allows to estimate the number of solution |
321 |
"C" : PDECoefficientType(PDECoefficientType.INTERIOR,(PDECoefficientType.EQUATION,PDECoefficientType.SOLUTION,PDECoefficientType.DIM),PDECoefficientType.OPERATOR), |
components. |
|
"D" : PDECoefficientType(PDECoefficientType.INTERIOR,(PDECoefficientType.EQUATION,PDECoefficientType.SOLUTION),PDECoefficientType.OPERATOR), |
|
|
"X" : PDECoefficientType(PDECoefficientType.INTERIOR,(PDECoefficientType.EQUATION,PDECoefficientType.DIM),PDECoefficientType.RIGHTHANDSIDE), |
|
|
"Y" : PDECoefficientType(PDECoefficientType.INTERIOR,(PDECoefficientType.EQUATION,),PDECoefficientType.RIGHTHANDSIDE), |
|
|
"d" : PDECoefficientType(PDECoefficientType.BOUNDARY,(PDECoefficientType.EQUATION,PDECoefficientType.SOLUTION),PDECoefficientType.OPERATOR), |
|
|
"y" : PDECoefficientType(PDECoefficientType.BOUNDARY,(PDECoefficientType.EQUATION,),PDECoefficientType.RIGHTHANDSIDE), |
|
|
"d_contact" : PDECoefficientType(PDECoefficientType.CONTACT,(PDECoefficientType.EQUATION,PDECoefficientType.SOLUTION),PDECoefficientType.OPERATOR), |
|
|
"y_contact" : PDECoefficientType(PDECoefficientType.CONTACT,(PDECoefficientType.EQUATION,),PDECoefficientType.RIGHTHANDSIDE), |
|
|
"r" : PDECoefficientType(PDECoefficientType.CONTINUOUS,(PDECoefficientType.EQUATION,),PDECoefficientType.RIGHTHANDSIDE), |
|
|
"q" : PDECoefficientType(PDECoefficientType.CONTINUOUS,(PDECoefficientType.SOLUTION,),PDECoefficientType.BOTH), |
|
|
} |
|
322 |
|
|
323 |
class LinearPDE: |
@return: True if the coefficient allows an estimate of the number of |
324 |
""" |
solution components, False otherwise |
325 |
@brief Class to define a linear PDE |
@rtype: C{bool} |
326 |
|
""" |
327 |
class to define a linear PDE of the form |
for i in self.pattern: |
328 |
|
if i==self.BY_SOLUTION: return True |
329 |
|
return False |
330 |
|
|
331 |
-(A_{ijkl}u_{k,l})_{,j} -(B_{ijk}u_k)_{,j} + C_{ikl}u_{k,l} +D_{ik}u_k = - (X_{ij})_{,j} + Y_i |
def definesNumEquation(self): |
332 |
|
""" |
333 |
|
Checks if the coefficient allows to estimate the number of equations. |
334 |
|
|
335 |
with boundary conditons: |
@return: True if the coefficient allows an estimate of the number of |
336 |
|
equations, False otherwise |
337 |
|
@rtype: C{bool} |
338 |
|
""" |
339 |
|
for i in self.pattern: |
340 |
|
if i==self.BY_EQUATION: return True |
341 |
|
return False |
342 |
|
|
343 |
n_j*(A_{ijkl}u_{k,l}+B_{ijk}u_k)_{,j} + d_{ik}u_k = - n_j*X_{ij} + y_i |
def __CompTuple2(self,t1,t2): |
344 |
|
""" |
345 |
|
Compares two tuples of possible number of equations and number of |
346 |
|
solutions. |
347 |
|
|
348 |
and contact conditions |
@param t1: the first tuple |
349 |
|
@param t2: the second tuple |
350 |
|
@return: 0, 1, or -1 |
351 |
|
""" |
352 |
|
|
353 |
n_j*(A_{ijkl}u_{k,l}+B_{ijk}u_k)_{,j} + d_contact_{ik}[u_k] = - n_j*X_{ij} + y_contact_i |
dif=t1[0]+t1[1]-(t2[0]+t2[1]) |
354 |
|
if dif<0: return 1 |
355 |
|
elif dif>0: return -1 |
356 |
|
else: return 0 |
357 |
|
|
358 |
and constraints: |
def getShape(self,domain,numEquations=1,numSolutions=1): |
359 |
|
""" |
360 |
|
Builds the required shape of the coefficient. |
361 |
|
|
362 |
u_i=r_i where q_i>0 |
@param domain: domain on which the PDE uses the coefficient |
363 |
|
@type domain: L{Domain<escript.Domain>} |
364 |
|
@param numEquations: number of equations of the PDE |
365 |
|
@type numEquations: C{int} |
366 |
|
@param numSolutions: number of components of the PDE solution |
367 |
|
@type numSolutions: C{int} |
368 |
|
@return: shape of the coefficient |
369 |
|
@rtype: C{tuple} of C{int} values |
370 |
|
""" |
371 |
|
dim=domain.getDim() |
372 |
|
s=() |
373 |
|
for i in self.pattern: |
374 |
|
if i==self.BY_EQUATION: |
375 |
|
if numEquations>1: s=s+(numEquations,) |
376 |
|
elif i==self.BY_SOLUTION: |
377 |
|
if numSolutions>1: s=s+(numSolutions,) |
378 |
|
else: |
379 |
|
s=s+(dim,) |
380 |
|
return s |
381 |
|
|
382 |
|
class LinearProblem(object): |
383 |
""" |
""" |
384 |
DEFAULT_METHOD=util.DEFAULT_METHOD |
This is the base class to define a general linear PDE-type problem for |
385 |
DIRECT=util.DIRECT |
for an unknown function M{u} on a given domain defined through a |
386 |
CHOLEVSKY=util.CHOLEVSKY |
L{Domain<escript.Domain>} object. The problem can be given as a single |
387 |
PCG=util.PCG |
equation or as a system of equations. |
388 |
CR=util.CR |
|
389 |
CGS=util.CGS |
The class assumes that some sort of assembling process is required to form |
390 |
BICGSTAB=util.BICGSTAB |
a problem of the form |
391 |
SSOR=util.SSOR |
|
392 |
GMRES=util.GMRES |
M{L u=f} |
393 |
PRES20=util.PRES20 |
|
394 |
|
where M{L} is an operator and M{f} is the right hand side. This operator |
395 |
|
problem will be solved to get the unknown M{u}. |
396 |
|
|
397 |
|
@cvar DEFAULT: The default method used to solve the system of linear |
398 |
|
equations |
399 |
|
@cvar DIRECT: The direct solver based on LDU factorization |
400 |
|
@cvar CHOLEVSKY: The direct solver based on LDLt factorization (can only be |
401 |
|
applied for symmetric PDEs) |
402 |
|
@cvar PCG: The preconditioned conjugate gradient method (can only be applied |
403 |
|
for symmetric PDEs) |
404 |
|
@cvar CR: The conjugate residual method |
405 |
|
@cvar CGS: The conjugate gradient square method |
406 |
|
@cvar BICGSTAB: The stabilized BiConjugate Gradient method |
407 |
|
@cvar TFQMR: Transport Free Quasi Minimal Residual method |
408 |
|
@cvar MINRES: Minimum residual method |
409 |
|
@cvar SSOR: The symmetric overrelaxation method |
410 |
|
@cvar ILU0: The incomplete LU factorization preconditioner with no fill-in |
411 |
|
@cvar ILUT: The incomplete LU factorization preconditioner with fill-in |
412 |
|
@cvar JACOBI: The Jacobi preconditioner |
413 |
|
@cvar GMRES: The Gram-Schmidt minimum residual method |
414 |
|
@cvar PRES20: Special GMRES with restart after 20 steps and truncation after |
415 |
|
5 residuals |
416 |
|
@cvar LUMPING: Matrix lumping |
417 |
|
@cvar NO_REORDERING: No matrix reordering allowed |
418 |
|
@cvar MINIMUM_FILL_IN: Reorder matrix to reduce fill-in during factorization |
419 |
|
@cvar NESTED_DISSECTION: Reorder matrix to improve load balancing during |
420 |
|
factorization |
421 |
|
@cvar PASO: PASO solver package |
422 |
|
@cvar SCSL: SGI SCSL solver library |
423 |
|
@cvar MKL: Intel's MKL solver library |
424 |
|
@cvar UMFPACK: The UMFPACK library |
425 |
|
@cvar TRILINOS: The TRILINOS parallel solver class library from Sandia Natl |
426 |
|
Labs |
427 |
|
@cvar ITERATIVE: The default iterative solver |
428 |
|
@cvar AMG: Algebraic Multi Grid |
429 |
|
@cvar RILU: Recursive ILU |
430 |
|
@cvar GS: Gauss-Seidel solver |
431 |
|
|
432 |
def __init__(self,domain,numEquations=0,numSolutions=0): |
""" |
433 |
""" |
DEFAULT= 0 |
434 |
@brief initializes a new linear PDE. |
DIRECT= 1 |
435 |
|
CHOLEVSKY= 2 |
436 |
|
PCG= 3 |
437 |
|
CR= 4 |
438 |
|
CGS= 5 |
439 |
|
BICGSTAB= 6 |
440 |
|
SSOR= 7 |
441 |
|
ILU0= 8 |
442 |
|
ILUT= 9 |
443 |
|
JACOBI= 10 |
444 |
|
GMRES= 11 |
445 |
|
PRES20= 12 |
446 |
|
LUMPING= 13 |
447 |
|
NO_REORDERING= 17 |
448 |
|
MINIMUM_FILL_IN= 18 |
449 |
|
NESTED_DISSECTION= 19 |
450 |
|
SCSL= 14 |
451 |
|
MKL= 15 |
452 |
|
UMFPACK= 16 |
453 |
|
ITERATIVE= 20 |
454 |
|
PASO= 21 |
455 |
|
AMG= 22 |
456 |
|
RILU = 23 |
457 |
|
TRILINOS = 24 |
458 |
|
NONLINEAR_GMRES = 25 |
459 |
|
TFQMR = 26 |
460 |
|
MINRES = 27 |
461 |
|
GS=28 |
462 |
|
|
463 |
|
SMALL_TOLERANCE=1.e-13 |
464 |
|
PACKAGE_KEY="package" |
465 |
|
METHOD_KEY="method" |
466 |
|
SYMMETRY_KEY="symmetric" |
467 |
|
TOLERANCE_KEY="tolerance" |
468 |
|
PRECONDITIONER_KEY="preconditioner" |
469 |
|
|
470 |
|
|
471 |
|
def __init__(self,domain,numEquations=None,numSolutions=None,debug=False): |
472 |
|
""" |
473 |
|
Initializes a linear problem. |
474 |
|
|
475 |
|
@param domain: domain of the PDE |
476 |
|
@type domain: L{Domain<escript.Domain>} |
477 |
|
@param numEquations: number of equations. If C{None} the number of |
478 |
|
equations is extracted from the coefficients. |
479 |
|
@param numSolutions: number of solution components. If C{None} the number |
480 |
|
of solution components is extracted from the |
481 |
|
coefficients. |
482 |
|
@param debug: if True debug information is printed |
483 |
|
|
|
@param args |
|
484 |
""" |
""" |
485 |
|
super(LinearProblem, self).__init__() |
486 |
|
|
487 |
# initialize attributes |
self.__debug=debug |
|
self.__debug=None |
|
488 |
self.__domain=domain |
self.__domain=domain |
489 |
self.__numEquations=numEquations |
self.__numEquations=numEquations |
490 |
self.__numSolutions=numSolutions |
self.__numSolutions=numSolutions |
491 |
self.cleanCoefficients() |
self.__altered_coefficients=False |
492 |
|
self.__reduce_equation_order=False |
493 |
self.__operator=escript.Operator() |
self.__reduce_solution_order=False |
|
self.__operator_isValid=False |
|
|
self.__righthandside=escript.Data() |
|
|
self.__righthandside_isValid=False |
|
|
self.__solution=escript.Data() |
|
|
self.__solution_isValid=False |
|
|
|
|
|
# set some default values: |
|
|
self.__homogeneous_constraint=True |
|
|
self.__row_function_space=escript.Solution(self.__domain) |
|
|
self.__column_function_space=escript.Solution(self.__domain) |
|
494 |
self.__tolerance=1.e-8 |
self.__tolerance=1.e-8 |
495 |
self.__solver_method=util.DEFAULT_METHOD |
self.__solver_method=self.DEFAULT |
496 |
self.__matrix_type=self.__domain.getSystemMatrixTypeId(util.DEFAULT_METHOD,False) |
self.__solver_package=self.DEFAULT |
497 |
|
self.__preconditioner=self.DEFAULT |
498 |
|
self.__is_RHS_valid=False |
499 |
|
self.__is_operator_valid=False |
500 |
self.__sym=False |
self.__sym=False |
501 |
self.__lumping=False |
self.__COEFFICIENTS={} |
502 |
|
# initialize things: |
503 |
|
self.resetAllCoefficients() |
504 |
|
self.__system_type=None |
505 |
|
self.updateSystemType() |
506 |
|
# ========================================================================== |
507 |
|
# general stuff: |
508 |
|
# ========================================================================== |
509 |
|
def __str__(self): |
510 |
|
""" |
511 |
|
Returns a string representation of the PDE. |
512 |
|
|
513 |
|
@return: a simple representation of the PDE |
514 |
|
@rtype: C{str} |
515 |
|
""" |
516 |
|
return "<LinearProblem %d>"%id(self) |
517 |
|
# ========================================================================== |
518 |
|
# debug : |
519 |
|
# ========================================================================== |
520 |
|
def setDebugOn(self): |
521 |
|
""" |
522 |
|
Switches debug output on. |
523 |
|
""" |
524 |
|
self.__debug=not None |
525 |
|
|
526 |
def getCoefficient(self,name): |
def setDebugOff(self): |
527 |
|
""" |
528 |
|
Switches debug output off. |
529 |
""" |
""" |
530 |
@brief return the value of the coefficient name |
self.__debug=None |
531 |
|
|
532 |
@param name |
def setDebug(self, flag): |
533 |
""" |
""" |
534 |
return self.__coefficient[name] |
Switches debug output on if C{flag} is True otherwise it is switched off. |
535 |
|
|
536 |
def setValue(self,**coefficients): |
@param flag: desired debug status |
537 |
""" |
@type flag: C{bool} |
538 |
@brief sets new values to coefficients |
""" |
539 |
|
if flag: |
540 |
|
self.setDebugOn() |
541 |
|
else: |
542 |
|
self.setDebugOff() |
543 |
|
|
544 |
@param coefficients |
def trace(self,text): |
545 |
""" |
""" |
546 |
self._setValue(**coefficients) |
Prints the text message if debug mode is switched on. |
|
|
|
547 |
|
|
548 |
def _setValue(self,**coefficients): |
@param text: message to be printed |
549 |
""" |
@type text: C{string} |
550 |
@brief sets new values to coefficients |
""" |
551 |
|
if self.__debug: print "%s: %s"%(str(self),text) |
552 |
|
|
553 |
@param coefficients |
# ========================================================================== |
554 |
""" |
# some service functions: |
555 |
|
# ========================================================================== |
556 |
# get the dictionary of the coefficinets been altered: |
def introduceCoefficients(self,**coeff): |
557 |
alteredCoefficients={} |
""" |
558 |
for i,d in coefficients.iteritems(): |
Introduces new coefficients into the problem. |
|
if self.hasCoefficient(i): |
|
|
if d == None: |
|
|
alteredCoefficients[i]=escript.Data() |
|
|
elif isinstance(d,escript.Data): |
|
|
if d.isEmpty(): |
|
|
alteredCoefficients[i]=escript.Data() |
|
|
else: |
|
|
alteredCoefficients[i]=escript.Data(d,self.getFunctionSpaceOfCoefficient(i)) |
|
|
else: |
|
|
alteredCoefficients[i]=escript.Data(d,self.getFunctionSpaceOfCoefficient(i)) |
|
|
else: |
|
|
raise ValueError,"Attempt to set undefined coefficient %s"%i |
|
|
# if numEquations and numSolutions is undefined we try identify their values based on the coefficients: |
|
|
if self.__numEquations<1 or self.__numSolutions<1: |
|
|
numEquations,numSolutions=identifyNumEquationsAndSolutions(self.getDomain().getDim(),alteredCoefficients) |
|
|
if self.__numEquations<1 and numEquations>0: self.__numEquations=numEquations |
|
|
if self.__numSolutions<1 and numSolutions>0: self.__numSolutions=numSolutions |
|
|
if self.debug() and self.__numEquations>0: print "PDE Debug: identified number of equations is ",self.__numEquations |
|
|
if self.debug() and self.__numSolutions>0: print "PDE Debug: identified number of solutions is ",self.__numSolutions |
|
559 |
|
|
560 |
# now we check the shape of the coefficient if numEquations and numSolutions are set: |
Use:: |
|
if self.__numEquations>0 and self.__numSolutions>0: |
|
|
for i in self.__coefficient.iterkeys(): |
|
|
if alteredCoefficients.has_key(i) and not alteredCoefficients[i].isEmpty(): |
|
|
if not self.getShapeOfCoefficient(i)==alteredCoefficients[i].getShape(): |
|
|
raise ValueError,"Expected shape for coefficient %s is %s but actual shape is %s."%(i,self.getShapeOfCoefficient(i),alteredCoefficients[i].getShape()) |
|
|
else: |
|
|
if not self.__coefficient[i].isEmpty(): |
|
|
if not self.getShapeOfCoefficient(i)==self.__coefficient[i].getShape(): |
|
|
raise ValueError,"Expected shape for coefficient %s is %s but actual shape is %s."%(i,self.getShapeOfCoefficient(i),self.__coefficient[i].getShape()) |
|
|
# overwrite new values: |
|
|
for i,d in alteredCoefficients.iteritems(): |
|
|
if self.debug(): print "PDE Debug: Coefficient %s has been altered."%i |
|
|
self.__coefficient[i]=d |
|
|
self.alteredCoefficient(i) |
|
|
|
|
|
# reset the HomogeneousConstraintFlag: |
|
|
self.__setHomogeneousConstraintFlag() |
|
|
if not "q" in alteredCoefficients and not self.__homogeneous_constraint: self.__rebuildSystem() |
|
|
|
|
|
def cleanCoefficients(self): |
|
|
""" |
|
|
@brief resets all coefficients to default values. |
|
|
""" |
|
|
self.__coefficient={} |
|
|
for i in _PDECoefficientTypes.iterkeys(): |
|
|
self.__coefficient[i]=escript.Data() |
|
561 |
|
|
562 |
def getShapeOfCoefficient(self,name): |
p.introduceCoefficients(A=PDECoef(...), B=PDECoef(...)) |
563 |
|
|
564 |
|
to introduce the coefficients M{A} ans M{B}. |
565 |
|
""" |
566 |
|
for name, type in coeff.items(): |
567 |
|
if not isinstance(type,PDECoef): |
568 |
|
raise ValueError,"coefficient %s has no type."%name |
569 |
|
self.__COEFFICIENTS[name]=type |
570 |
|
self.__COEFFICIENTS[name].resetValue() |
571 |
|
self.trace("coefficient %s has been introduced."%name) |
572 |
|
|
573 |
|
def getDomain(self): |
574 |
""" |
""" |
575 |
@brief return the shape of the coefficient name |
Returns the domain of the PDE. |
576 |
|
|
577 |
@param name |
@return: the domain of the PDE |
578 |
|
@rtype: L{Domain<escript.Domain>} |
579 |
""" |
""" |
580 |
if self.hasCoefficient(name): |
return self.__domain |
|
return _PDECoefficientTypes[name].buildShape(self.getNumEquations(),self.getNumSolutions(),self.getDomain().getDim()) |
|
|
else: |
|
|
raise ValueError,"Solution coefficient %s requested"%name |
|
581 |
|
|
582 |
def getFunctionSpaceOfCoefficient(self,name): |
def getDim(self): |
583 |
""" |
""" |
584 |
@brief return the atoms of the coefficient name |
Returns the spatial dimension of the PDE. |
585 |
|
|
586 |
@param name |
@return: the spatial dimension of the PDE domain |
587 |
|
@rtype: C{int} |
588 |
""" |
""" |
589 |
if self.hasCoefficient(name): |
return self.getDomain().getDim() |
|
return _PDECoefficientTypes[name].getFunctionSpace(self.getDomain()) |
|
|
else: |
|
|
raise ValueError,"Solution coefficient %s requested"%name |
|
590 |
|
|
591 |
def alteredCoefficient(self,name): |
def getNumEquations(self): |
592 |
""" |
""" |
593 |
@brief annonced that coefficient name has been changed |
Returns the number of equations. |
594 |
|
|
595 |
|
@return: the number of equations |
596 |
|
@rtype: C{int} |
597 |
|
@raise UndefinedPDEError: if the number of equations is not specified yet |
598 |
|
""" |
599 |
|
if self.__numEquations==None: |
600 |
|
if self.__numSolutions==None: |
601 |
|
raise UndefinedPDEError,"Number of equations is undefined. Please specify argument numEquations." |
602 |
|
else: |
603 |
|
self.__numEquations=self.__numSolutions |
604 |
|
return self.__numEquations |
605 |
|
|
606 |
@param name |
def getNumSolutions(self): |
607 |
""" |
""" |
608 |
if self.hasCoefficient(name): |
Returns the number of unknowns. |
|
if _PDECoefficientTypes[name].isAlteringOperator(): self.__rebuildOperator() |
|
|
if _PDECoefficientTypes[name].isAlteringRightHandSide(): self.__rebuildRightHandSide() |
|
|
else: |
|
|
raise ValueError,"Solution coefficient %s requested"%name |
|
609 |
|
|
610 |
def __setHomogeneousConstraintFlag(self): |
@return: the number of unknowns |
611 |
""" |
@rtype: C{int} |
612 |
@brief checks if the constraints are homogeneous and sets self.__homogeneous_constraint accordingly. |
@raise UndefinedPDEError: if the number of unknowns is not specified yet |
613 |
""" |
""" |
614 |
self.__homogeneous_constraint=True |
if self.__numSolutions==None: |
615 |
q=self.getCoefficient("q") |
if self.__numEquations==None: |
616 |
r=self.getCoefficient("r") |
raise UndefinedPDEError,"Number of solution is undefined. Please specify argument numSolutions." |
617 |
if not q.isEmpty() and not r.isEmpty(): |
else: |
618 |
print (q*r).Lsup(), 1.e-13*r.Lsup() |
self.__numSolutions=self.__numEquations |
619 |
if (q*r).Lsup()>=1.e-13*r.Lsup(): self.__homogeneous_constraint=False |
return self.__numSolutions |
|
if self.debug(): |
|
|
if self.__homogeneous_constraint: |
|
|
print "PDE Debug: Constraints are homogeneous." |
|
|
else: |
|
|
print "PDE Debug: Constraints are inhomogeneous." |
|
|
|
|
620 |
|
|
621 |
def hasCoefficient(self,name): |
def reduceEquationOrder(self): |
622 |
""" |
""" |
623 |
@brief return true if name is the name of a coefficient |
Returns the status of order reduction for the equation. |
624 |
|
|
625 |
@param name |
@return: True if reduced interpolation order is used for the |
626 |
""" |
representation of the equation, False otherwise |
627 |
return self.__coefficient.has_key(name) |
@rtype: L{bool} |
628 |
|
""" |
629 |
|
return self.__reduce_equation_order |
630 |
|
|
631 |
|
def reduceSolutionOrder(self): |
632 |
|
""" |
633 |
|
Returns the status of order reduction for the solution. |
634 |
|
|
635 |
|
@return: True if reduced interpolation order is used for the |
636 |
|
representation of the solution, False otherwise |
637 |
|
@rtype: L{bool} |
638 |
|
""" |
639 |
|
return self.__reduce_solution_order |
640 |
|
|
641 |
def getFunctionSpaceForEquation(self): |
def getFunctionSpaceForEquation(self): |
642 |
""" |
""" |
643 |
@brief return true if the test functions should use reduced order |
Returns the L{FunctionSpace<escript.FunctionSpace>} used to discretize |
644 |
|
the equation. |
645 |
|
|
646 |
|
@return: representation space of equation |
647 |
|
@rtype: L{FunctionSpace<escript.FunctionSpace>} |
648 |
""" |
""" |
649 |
return self.__row_function_space |
if self.reduceEquationOrder(): |
650 |
|
return escript.ReducedSolution(self.getDomain()) |
651 |
|
else: |
652 |
|
return escript.Solution(self.getDomain()) |
653 |
|
|
654 |
def getFunctionSpaceForSolution(self): |
def getFunctionSpaceForSolution(self): |
655 |
""" |
""" |
656 |
@brief return true if the interpolation of the solution should use reduced order |
Returns the L{FunctionSpace<escript.FunctionSpace>} used to represent |
657 |
|
the solution. |
658 |
|
|
659 |
|
@return: representation space of solution |
660 |
|
@rtype: L{FunctionSpace<escript.FunctionSpace>} |
661 |
""" |
""" |
662 |
return self.__column_function_space |
if self.reduceSolutionOrder(): |
663 |
|
return escript.ReducedSolution(self.getDomain()) |
664 |
|
else: |
665 |
|
return escript.Solution(self.getDomain()) |
666 |
|
|
667 |
# ===== debug ============================================================== |
# ========================================================================== |
668 |
def setDebugOn(self): |
# solver settings: |
669 |
""" |
# ========================================================================== |
670 |
@brief |
def setSolverMethod(self,solver=None,preconditioner=None): |
671 |
""" |
""" |
672 |
self.__debug=not None |
Sets a new solver method and/or preconditioner. |
673 |
|
|
674 |
|
@param solver: new solver method to use |
675 |
|
@type solver: one of L{DEFAULT}, L{ITERATIVE} L{DIRECT}, L{CHOLEVSKY}, |
676 |
|
L{PCG}, L{CR}, L{CGS}, L{BICGSTAB}, L{SSOR}, L{GMRES}, |
677 |
|
L{TFQMR}, L{MINRES}, L{PRES20}, L{LUMPING}, L{AMG} |
678 |
|
@param preconditioner: new preconditioner to use |
679 |
|
@type preconditioner: one of L{DEFAULT}, L{JACOBI} L{ILU0}, L{ILUT}, |
680 |
|
L{SSOR}, L{RILU}, L{GS} |
681 |
|
|
682 |
|
@note: the solver method chosen may not be available by the selected |
683 |
|
package. |
684 |
|
""" |
685 |
|
if solver==None: solver=self.__solver_method |
686 |
|
if preconditioner==None: preconditioner=self.__preconditioner |
687 |
|
if solver==None: solver=self.DEFAULT |
688 |
|
if preconditioner==None: preconditioner=self.DEFAULT |
689 |
|
if not (solver,preconditioner)==self.getSolverMethod(): |
690 |
|
self.__solver_method=solver |
691 |
|
self.__preconditioner=preconditioner |
692 |
|
self.updateSystemType() |
693 |
|
self.trace("New solver is %s"%self.getSolverMethodName()) |
694 |
|
|
695 |
|
def getSolverMethodName(self): |
696 |
|
""" |
697 |
|
Returns the name of the solver currently used. |
698 |
|
|
699 |
|
@return: the name of the solver currently used |
700 |
|
@rtype: C{string} |
701 |
|
""" |
702 |
|
|
703 |
|
m=self.getSolverMethod() |
704 |
|
p=self.getSolverPackage() |
705 |
|
method="" |
706 |
|
if m[0]==self.DEFAULT: method="DEFAULT" |
707 |
|
elif m[0]==self.DIRECT: method= "DIRECT" |
708 |
|
elif m[0]==self.ITERATIVE: method= "ITERATIVE" |
709 |
|
elif m[0]==self.CHOLEVSKY: method= "CHOLEVSKY" |
710 |
|
elif m[0]==self.PCG: method= "PCG" |
711 |
|
elif m[0]==self.TFQMR: method= "TFQMR" |
712 |
|
elif m[0]==self.MINRES: method= "MINRES" |
713 |
|
elif m[0]==self.CR: method= "CR" |
714 |
|
elif m[0]==self.CGS: method= "CGS" |
715 |
|
elif m[0]==self.BICGSTAB: method= "BICGSTAB" |
716 |
|
elif m[0]==self.SSOR: method= "SSOR" |
717 |
|
elif m[0]==self.GMRES: method= "GMRES" |
718 |
|
elif m[0]==self.PRES20: method= "PRES20" |
719 |
|
elif m[0]==self.LUMPING: method= "LUMPING" |
720 |
|
elif m[0]==self.AMG: method= "AMG" |
721 |
|
if m[1]==self.DEFAULT: method+="+DEFAULT" |
722 |
|
elif m[1]==self.JACOBI: method+= "+JACOBI" |
723 |
|
elif m[1]==self.ILU0: method+= "+ILU0" |
724 |
|
elif m[1]==self.ILUT: method+= "+ILUT" |
725 |
|
elif m[1]==self.SSOR: method+= "+SSOR" |
726 |
|
elif m[1]==self.AMG: method+= "+AMG" |
727 |
|
elif m[1]==self.RILU: method+= "+RILU" |
728 |
|
elif m[1]==self.GS: method+= "+GS" |
729 |
|
if p==self.DEFAULT: package="DEFAULT" |
730 |
|
elif p==self.PASO: package= "PASO" |
731 |
|
elif p==self.MKL: package= "MKL" |
732 |
|
elif p==self.SCSL: package= "SCSL" |
733 |
|
elif p==self.UMFPACK: package= "UMFPACK" |
734 |
|
elif p==self.TRILINOS: package= "TRILINOS" |
735 |
|
else : method="unknown" |
736 |
|
return "%s solver of %s package"%(method,package) |
737 |
|
|
738 |
def setDebugOff(self): |
def getSolverMethod(self): |
|
""" |
|
|
@brief |
|
739 |
""" |
""" |
740 |
self.__debug=None |
Returns the solver method and preconditioner used. |
741 |
|
|
742 |
def debug(self): |
@return: the solver method currently used. |
743 |
|
@rtype: C{tuple} of C{int} |
744 |
""" |
""" |
745 |
@brief returns true if the PDE is in the debug mode |
return self.__solver_method,self.__preconditioner |
746 |
|
|
747 |
|
def setSolverPackage(self,package=None): |
748 |
""" |
""" |
749 |
return self.__debug |
Sets a new solver package. |
750 |
|
|
751 |
#===== Lumping =========================== |
@param package: the new solver package |
752 |
def setLumpingOn(self): |
@type package: one of L{DEFAULT}, L{PASO} L{SCSL}, L{MKL}, L{UMFPACK}, |
753 |
""" |
L{TRILINOS} |
754 |
@brief indicates to use matrix lumping |
""" |
755 |
""" |
if package==None: package=self.DEFAULT |
756 |
if not self.isUsingLumping(): |
if not package==self.getSolverPackage(): |
757 |
raise SystemError,"Lumping is not working yet! Talk to the experts" |
self.__solver_package=package |
758 |
if self.debug() : print "PDE Debug: lumping is set on" |
self.updateSystemType() |
759 |
self.__rebuildOperator() |
self.trace("New solver is %s"%self.getSolverMethodName()) |
|
self.__lumping=True |
|
760 |
|
|
761 |
def setLumpingOff(self): |
def getSolverPackage(self): |
762 |
""" |
""" |
763 |
@brief switches off matrix lumping |
Returns the package of the solver. |
|
""" |
|
|
if self.isUsingLumping(): |
|
|
if self.debug() : print "PDE Debug: lumping is set off" |
|
|
self.__rebuildOperator() |
|
|
self.__lumping=False |
|
764 |
|
|
765 |
def setLumping(self,flag=False): |
@return: the solver package currently being used |
766 |
""" |
@rtype: C{int} |
767 |
@brief set the matrix lumping flag to flag |
""" |
768 |
""" |
return self.__solver_package |
|
if flag: |
|
|
self.setLumpingOn() |
|
|
else: |
|
|
self.setLumpingOff() |
|
769 |
|
|
770 |
def isUsingLumping(self): |
def isUsingLumping(self): |
771 |
""" |
""" |
772 |
@brief |
Checks if matrix lumping is the current solver method. |
|
""" |
|
|
return self.__lumping |
|
|
|
|
|
#============ method business ========================================================= |
|
|
def setSolverMethod(self,solver=util.DEFAULT_METHOD): |
|
|
""" |
|
|
@brief sets a new solver |
|
|
""" |
|
|
if not solver==self.getSolverMethod(): |
|
|
self.__solver_method=solver |
|
|
if self.debug() : print "PDE Debug: New solver is %s"%solver |
|
|
self.__checkMatrixType() |
|
773 |
|
|
774 |
def getSolverMethod(self): |
@return: True if the current solver method is lumping |
775 |
""" |
@rtype: C{bool} |
776 |
@brief returns the solver method |
""" |
777 |
""" |
return self.getSolverMethod()[0]==self.LUMPING |
|
return self.__solver_method |
|
778 |
|
|
|
#============ tolerance business ========================================================= |
|
779 |
def setTolerance(self,tol=1.e-8): |
def setTolerance(self,tol=1.e-8): |
780 |
""" |
""" |
781 |
@brief resets the tolerance to tol. |
Resets the tolerance for the solver method to C{tol}. |
782 |
|
|
783 |
|
@param tol: new tolerance for the solver. If C{tol} is lower than the |
784 |
|
current tolerence the system will be resolved. |
785 |
|
@type tol: positive C{float} |
786 |
|
@raise ValueError: if tolerance is not positive |
787 |
""" |
""" |
788 |
if not tol>0: |
if not tol>0: |
789 |
raise ValueException,"Tolerance as to be positive" |
raise ValueError,"Tolerance has to be positive" |
790 |
if tol<self.getTolerance(): self.__rebuildSolution() |
if tol<self.getTolerance(): self.invalidateSolution() |
791 |
if self.debug() : print "PDE Debug: New tolerance %e",tol |
self.trace("New tolerance %e"%tol) |
792 |
self.__tolerance=tol |
self.__tolerance=tol |
793 |
return |
return |
794 |
|
|
795 |
def getTolerance(self): |
def getTolerance(self): |
796 |
""" |
""" |
797 |
@brief returns the tolerance set for the solution |
Returns the tolerance currently set for the solution. |
798 |
|
|
799 |
|
@return: tolerance currently used |
800 |
|
@rtype: C{float} |
801 |
""" |
""" |
802 |
return self.__tolerance |
return self.__tolerance |
803 |
|
|
804 |
#===== symmetry flag ========================== |
# ========================================================================== |
805 |
|
# symmetry flag: |
806 |
|
# ========================================================================== |
807 |
def isSymmetric(self): |
def isSymmetric(self): |
808 |
""" |
""" |
809 |
@brief returns true is the operator is considered to be symmetric |
Checks if symmetry is indicated. |
810 |
|
|
811 |
|
@return: True if a symmetric PDE is indicated, False otherwise |
812 |
|
@rtype: C{bool} |
813 |
""" |
""" |
814 |
return self.__sym |
return self.__sym |
815 |
|
|
816 |
def setSymmetryOn(self): |
def setSymmetryOn(self): |
817 |
""" |
""" |
818 |
@brief sets the symmetry flag to true |
Sets the symmetry flag. |
819 |
""" |
""" |
820 |
if not self.isSymmetric(): |
if not self.isSymmetric(): |
821 |
if self.debug() : print "PDE Debug: Operator is set to be symmetric" |
self.trace("PDE is set to be symmetric") |
822 |
self.__sym=True |
self.__sym=True |
823 |
self.__checkMatrixType() |
self.updateSystemType() |
824 |
|
|
825 |
def setSymmetryOff(self): |
def setSymmetryOff(self): |
826 |
""" |
""" |
827 |
@brief sets the symmetry flag to false |
Clears the symmetry flag. |
828 |
""" |
""" |
829 |
if self.isSymmetric(): |
if self.isSymmetric(): |
830 |
if self.debug() : print "PDE Debug: Operator is set to be unsymmetric" |
self.trace("PDE is set to be nonsymmetric") |
831 |
self.__sym=False |
self.__sym=False |
832 |
self.__checkMatrixType() |
self.updateSystemType() |
833 |
|
|
834 |
def setSymmetryTo(self,flag=False): |
def setSymmetryTo(self,flag=False): |
835 |
""" |
""" |
836 |
@brief sets the symmetry flag to flag |
Sets the symmetry flag to C{flag}. |
837 |
|
|
838 |
@param flag |
@param flag: If True, the symmetry flag is set otherwise reset. |
839 |
""" |
@type flag: C{bool} |
840 |
if flag: |
""" |
841 |
self.setSymmetryOn() |
if flag: |
842 |
else: |
self.setSymmetryOn() |
843 |
self.setSymmetryOff() |
else: |
844 |
|
self.setSymmetryOff() |
845 |
|
|
846 |
#===== order reduction ========================== |
# ========================================================================== |
847 |
|
# function space handling for the equation as well as the solution |
848 |
|
# ========================================================================== |
849 |
def setReducedOrderOn(self): |
def setReducedOrderOn(self): |
850 |
""" |
""" |
851 |
@brief switches to on reduced order |
Switches reduced order on for solution and equation representation. |
852 |
|
|
853 |
|
@raise RuntimeError: if order reduction is altered after a coefficient has |
854 |
|
been set |
855 |
""" |
""" |
856 |
self.setReducedOrderForSolutionOn() |
self.setReducedOrderForSolutionOn() |
857 |
self.setReducedOrderForEquationOn() |
self.setReducedOrderForEquationOn() |
858 |
|
|
859 |
def setReducedOrderOff(self): |
def setReducedOrderOff(self): |
860 |
""" |
""" |
861 |
@brief switches to full order |
Switches reduced order off for solution and equation representation |
862 |
|
|
863 |
|
@raise RuntimeError: if order reduction is altered after a coefficient has |
864 |
|
been set |
865 |
""" |
""" |
866 |
self.setReducedOrderForSolutionOff() |
self.setReducedOrderForSolutionOff() |
867 |
self.setReducedOrderForEquationOff() |
self.setReducedOrderForEquationOff() |
868 |
|
|
869 |
def setReducedOrderTo(self,flag=False): |
def setReducedOrderTo(self,flag=False): |
870 |
""" |
""" |
871 |
@brief sets order according to flag |
Sets order reduction state for both solution and equation representation |
872 |
|
according to flag. |
873 |
|
|
874 |
@param flag |
@param flag: if True, the order reduction is switched on for both solution |
875 |
|
and equation representation, otherwise or if flag is not |
876 |
|
present order reduction is switched off |
877 |
|
@type flag: C{bool} |
878 |
|
@raise RuntimeError: if order reduction is altered after a coefficient has |
879 |
|
been set |
880 |
""" |
""" |
881 |
self.setReducedOrderForSolutionTo(flag) |
self.setReducedOrderForSolutionTo(flag) |
882 |
self.setReducedOrderForEquationTo(flag) |
self.setReducedOrderForEquationTo(flag) |
|
|
|
883 |
|
|
884 |
#===== order reduction solution ========================== |
|
885 |
def setReducedOrderForSolutionOn(self): |
def setReducedOrderForSolutionOn(self): |
886 |
""" |
""" |
887 |
@brief switches to reduced order to interpolate solution |
Switches reduced order on for solution representation. |
888 |
|
|
889 |
|
@raise RuntimeError: if order reduction is altered after a coefficient has |
890 |
|
been set |
891 |
""" |
""" |
892 |
new_fs=escript.ReducedSolution(self.getDomain()) |
if not self.__reduce_solution_order: |
893 |
if self.getFunctionSpaceForSolution()!=new_fs: |
if self.__altered_coefficients: |
894 |
if self.debug() : print "PDE Debug: Reduced order is used to interpolate solution." |
raise RuntimeError,"order cannot be altered after coefficients have been defined." |
895 |
self.__column_function_space=new_fs |
self.trace("Reduced order is used for solution representation.") |
896 |
self.__rebuildSystem(deep=True) |
self.__reduce_solution_order=True |
897 |
|
self.initializeSystem() |
898 |
|
|
899 |
def setReducedOrderForSolutionOff(self): |
def setReducedOrderForSolutionOff(self): |
900 |
""" |
""" |
901 |
@brief switches to full order to interpolate solution |
Switches reduced order off for solution representation |
902 |
|
|
903 |
|
@raise RuntimeError: if order reduction is altered after a coefficient has |
904 |
|
been set. |
905 |
""" |
""" |
906 |
new_fs=escript.Solution(self.getDomain()) |
if self.__reduce_solution_order: |
907 |
if self.getFunctionSpaceForSolution()!=new_fs: |
if self.__altered_coefficients: |
908 |
if self.debug() : print "PDE Debug: Full order is used to interpolate solution." |
raise RuntimeError,"order cannot be altered after coefficients have been defined." |
909 |
self.__column_function_space=new_fs |
self.trace("Full order is used to interpolate solution.") |
910 |
self.__rebuildSystem(deep=True) |
self.__reduce_solution_order=False |
911 |
|
self.initializeSystem() |
912 |
|
|
913 |
def setReducedOrderForSolutionTo(self,flag=False): |
def setReducedOrderForSolutionTo(self,flag=False): |
914 |
""" |
""" |
915 |
@brief sets order for test functions according to flag |
Sets order reduction state for solution representation according to flag. |
916 |
|
|
917 |
@param flag |
@param flag: if flag is True, the order reduction is switched on for |
918 |
|
solution representation, otherwise or if flag is not present |
919 |
|
order reduction is switched off |
920 |
|
@type flag: C{bool} |
921 |
|
@raise RuntimeError: if order reduction is altered after a coefficient has |
922 |
|
been set |
923 |
""" |
""" |
924 |
if flag: |
if flag: |
925 |
self.setReducedOrderForSolutionOn() |
self.setReducedOrderForSolutionOn() |
926 |
else: |
else: |
927 |
self.setReducedOrderForSolutionOff() |
self.setReducedOrderForSolutionOff() |
928 |
|
|
|
#===== order reduction equation ========================== |
|
929 |
def setReducedOrderForEquationOn(self): |
def setReducedOrderForEquationOn(self): |
930 |
""" |
""" |
931 |
@brief switches to reduced order for test functions |
Switches reduced order on for equation representation. |
932 |
|
|
933 |
|
@raise RuntimeError: if order reduction is altered after a coefficient has |
934 |
|
been set |
935 |
""" |
""" |
936 |
new_fs=escript.ReducedSolution(self.getDomain()) |
if not self.__reduce_equation_order: |
937 |
if self.getFunctionSpaceForEquation()!=new_fs: |
if self.__altered_coefficients: |
938 |
if self.debug() : print "PDE Debug: Reduced order is used for test functions." |
raise RuntimeError,"order cannot be altered after coefficients have been defined." |
939 |
self.__row_function_space=new_fs |
self.trace("Reduced order is used for test functions.") |
940 |
self.__rebuildSystem(deep=True) |
self.__reduce_equation_order=True |
941 |
|
self.initializeSystem() |
942 |
|
|
943 |
def setReducedOrderForEquationOff(self): |
def setReducedOrderForEquationOff(self): |
944 |
""" |
""" |
945 |
@brief switches to full order for test functions |
Switches reduced order off for equation representation. |
946 |
|
|
947 |
|
@raise RuntimeError: if order reduction is altered after a coefficient has |
948 |
|
been set |
949 |
""" |
""" |
950 |
new_fs=escript.Solution(self.getDomain()) |
if self.__reduce_equation_order: |
951 |
if self.getFunctionSpaceForEquation()!=new_fs: |
if self.__altered_coefficients: |
952 |
if self.debug() : print "PDE Debug: Full order is used for test functions." |
raise RuntimeError,"order cannot be altered after coefficients have been defined." |
953 |
self.__row_function_space=new_fs |
self.trace("Full order is used for test functions.") |
954 |
self.__rebuildSystem(deep=True) |
self.__reduce_equation_order=False |
955 |
|
self.initializeSystem() |
956 |
|
|
957 |
def setReducedOrderForEquationTo(self,flag=False): |
def setReducedOrderForEquationTo(self,flag=False): |
958 |
""" |
""" |
959 |
@brief sets order for test functions according to flag |
Sets order reduction state for equation representation according to flag. |
960 |
|
|
961 |
@param flag |
@param flag: if flag is True, the order reduction is switched on for |
962 |
|
equation representation, otherwise or if flag is not present |
963 |
|
order reduction is switched off |
964 |
|
@type flag: C{bool} |
965 |
|
@raise RuntimeError: if order reduction is altered after a coefficient has |
966 |
|
been set |
967 |
""" |
""" |
968 |
if flag: |
if flag: |
969 |
self.setReducedOrderForEquationOn() |
self.setReducedOrderForEquationOn() |
970 |
else: |
else: |
971 |
self.setReducedOrderForEquationOff() |
self.setReducedOrderForEquationOff() |
|
|
|
972 |
|
|
973 |
# ==== initialization ===================================================================== |
def updateSystemType(self): |
974 |
def __makeNewOperator(self): |
""" |
975 |
|
Reassesses the system type and, if a new matrix is needed, resets the |
976 |
|
system. |
977 |
|
""" |
978 |
|
new_system_type=self.getRequiredSystemType() |
979 |
|
if not new_system_type==self.__system_type: |
980 |
|
self.trace("Matrix type is now %d."%new_system_type) |
981 |
|
self.__system_type=new_system_type |
982 |
|
self.initializeSystem() |
983 |
|
|
984 |
|
def getSystemType(self): |
985 |
|
""" |
986 |
|
Returns the current system type. |
987 |
|
""" |
988 |
|
return self.__system_type |
989 |
|
|
990 |
|
def checkSymmetricTensor(self,name,verbose=True): |
991 |
|
""" |
992 |
|
Tests a coefficient for symmetry. |
993 |
|
|
994 |
|
@param name: name of the coefficient |
995 |
|
@type name: C{str} |
996 |
|
@param verbose: if set to True or not present a report on coefficients |
997 |
|
which break the symmetry is printed. |
998 |
|
@type verbose: C{bool} |
999 |
|
@return: True if coefficient C{name} is symmetric |
1000 |
|
@rtype: C{bool} |
1001 |
|
""" |
1002 |
|
A=self.getCoefficient(name) |
1003 |
|
verbose=verbose or self.__debug |
1004 |
|
out=True |
1005 |
|
if not A.isEmpty(): |
1006 |
|
tol=util.Lsup(A)*self.SMALL_TOLERANCE |
1007 |
|
s=A.getShape() |
1008 |
|
if A.getRank() == 4: |
1009 |
|
if s[0]==s[2] and s[1] == s[3]: |
1010 |
|
for i in range(s[0]): |
1011 |
|
for j in range(s[1]): |
1012 |
|
for k in range(s[2]): |
1013 |
|
for l in range(s[3]): |
1014 |
|
if util.Lsup(A[i,j,k,l]-A[k,l,i,j])>tol: |
1015 |
|
if verbose: print "non-symmetric problem as %s[%d,%d,%d,%d]!=%s[%d,%d,%d,%d]"%(name,i,j,k,l,name,k,l,i,j) |
1016 |
|
out=False |
1017 |
|
else: |
1018 |
|
if verbose: print "non-symmetric problem because of inappropriate shape %s of coefficient %s."%(s,name) |
1019 |
|
out=False |
1020 |
|
elif A.getRank() == 2: |
1021 |
|
if s[0]==s[1]: |
1022 |
|
for j in range(s[0]): |
1023 |
|
for l in range(s[1]): |
1024 |
|
if util.Lsup(A[j,l]-A[l,j])>tol: |
1025 |
|
if verbose: print "non-symmetric problem because %s[%d,%d]!=%s[%d,%d]"%(name,j,l,name,l,j) |
1026 |
|
out=False |
1027 |
|
else: |
1028 |
|
if verbose: print "non-symmetric problem because of inappropriate shape %s of coefficient %s."%(s,name) |
1029 |
|
out=False |
1030 |
|
elif A.getRank() == 0: |
1031 |
|
pass |
1032 |
|
else: |
1033 |
|
raise ValueError,"Cannot check rank %s of %s."%(A.getRank(),name) |
1034 |
|
return out |
1035 |
|
|
1036 |
|
def checkReciprocalSymmetry(self,name0,name1,verbose=True): |
1037 |
|
""" |
1038 |
|
Tests two coefficients for reciprocal symmetry. |
1039 |
|
|
1040 |
|
@param name0: name of the first coefficient |
1041 |
|
@type name0: C{str} |
1042 |
|
@param name1: name of the second coefficient |
1043 |
|
@type name1: C{str} |
1044 |
|
@param verbose: if set to True or not present a report on coefficients |
1045 |
|
which break the symmetry is printed |
1046 |
|
@type verbose: C{bool} |
1047 |
|
@return: True if coefficients C{name0} and C{name1} are reciprocally |
1048 |
|
symmetric. |
1049 |
|
@rtype: C{bool} |
1050 |
|
""" |
1051 |
|
B=self.getCoefficient(name0) |
1052 |
|
C=self.getCoefficient(name1) |
1053 |
|
verbose=verbose or self.__debug |
1054 |
|
out=True |
1055 |
|
if B.isEmpty() and not C.isEmpty(): |
1056 |
|
if verbose: print "non-symmetric problem because %s is not present but %s is"%(name0,name1) |
1057 |
|
out=False |
1058 |
|
elif not B.isEmpty() and C.isEmpty(): |
1059 |
|
if verbose: print "non-symmetric problem because %s is not present but %s is"%(name0,name1) |
1060 |
|
out=False |
1061 |
|
elif not B.isEmpty() and not C.isEmpty(): |
1062 |
|
sB=B.getShape() |
1063 |
|
sC=C.getShape() |
1064 |
|
tol=(util.Lsup(B)+util.Lsup(C))*self.SMALL_TOLERANCE/2. |
1065 |
|
if len(sB) != len(sC): |
1066 |
|
if verbose: print "non-symmetric problem because ranks of %s (=%s) and %s (=%s) are different."%(name0,len(sB),name1,len(sC)) |
1067 |
|
out=False |
1068 |
|
else: |
1069 |
|
if len(sB)==0: |
1070 |
|
if util.Lsup(B-C)>tol: |
1071 |
|
if verbose: print "non-symmetric problem because %s!=%s"%(name0,name1) |
1072 |
|
out=False |
1073 |
|
elif len(sB)==1: |
1074 |
|
if sB[0]==sC[0]: |
1075 |
|
for j in range(sB[0]): |
1076 |
|
if util.Lsup(B[j]-C[j])>tol: |
1077 |
|
if verbose: print "non-symmetric PDE because %s[%d]!=%s[%d]"%(name0,j,name1,j) |
1078 |
|
out=False |
1079 |
|
else: |
1080 |
|
if verbose: print "non-symmetric problem because of inappropriate shapes %s and %s of coefficients %s and %s, respectively."%(sB,sC,name0,name1) |
1081 |
|
elif len(sB)==3: |
1082 |
|
if sB[0]==sC[1] and sB[1]==sC[2] and sB[2]==sC[0]: |
1083 |
|
for i in range(sB[0]): |
1084 |
|
for j in range(sB[1]): |
1085 |
|
for k in range(sB[2]): |
1086 |
|
if util.Lsup(B[i,j,k]-C[k,i,j])>tol: |
1087 |
|
if verbose: print "non-symmetric problem because %s[%d,%d,%d]!=%s[%d,%d,%d]"%(name0,i,j,k,name1,k,i,j) |
1088 |
|
out=False |
1089 |
|
else: |
1090 |
|
if verbose: print "non-symmetric problem because of inappropriate shapes %s and %s of coefficients %s and %s, respectively."%(sB,sC,name0,name1) |
1091 |
|
else: |
1092 |
|
raise ValueError,"Cannot check rank %s of %s and %s."%(len(sB),name0,name1) |
1093 |
|
return out |
1094 |
|
|
1095 |
|
def getCoefficient(self,name): |
1096 |
|
""" |
1097 |
|
Returns the value of the coefficient C{name}. |
1098 |
|
|
1099 |
|
@param name: name of the coefficient requested |
1100 |
|
@type name: C{string} |
1101 |
|
@return: the value of the coefficient |
1102 |
|
@rtype: L{Data<escript.Data>} |
1103 |
|
@raise IllegalCoefficient: if C{name} is not a coefficient of the PDE |
1104 |
|
""" |
1105 |
|
if self.hasCoefficient(name): |
1106 |
|
return self.__COEFFICIENTS[name].getValue() |
1107 |
|
else: |
1108 |
|
raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name |
1109 |
|
|
1110 |
|
def hasCoefficient(self,name): |
1111 |
|
""" |
1112 |
|
Returns True if C{name} is the name of a coefficient. |
1113 |
|
|
1114 |
|
@param name: name of the coefficient enquired |
1115 |
|
@type name: C{string} |
1116 |
|
@return: True if C{name} is the name of a coefficient of the general PDE, |
1117 |
|
False otherwise |
1118 |
|
@rtype: C{bool} |
1119 |
|
""" |
1120 |
|
return self.__COEFFICIENTS.has_key(name) |
1121 |
|
|
1122 |
|
def createCoefficient(self, name): |
1123 |
|
""" |
1124 |
|
Creates a L{Data<escript.Data>} object corresponding to coefficient |
1125 |
|
C{name}. |
1126 |
|
|
1127 |
|
@return: the coefficient C{name} initialized to 0 |
1128 |
|
@rtype: L{Data<escript.Data>} |
1129 |
|
@raise IllegalCoefficient: if C{name} is not a coefficient of the PDE |
1130 |
|
""" |
1131 |
|
if self.hasCoefficient(name): |
1132 |
|
return escript.Data(0.,self.getShapeOfCoefficient(name),self.getFunctionSpaceForCoefficient(name)) |
1133 |
|
else: |
1134 |
|
raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name |
1135 |
|
|
1136 |
|
def getFunctionSpaceForCoefficient(self,name): |
1137 |
|
""" |
1138 |
|
Returns the L{FunctionSpace<escript.FunctionSpace>} to be used for |
1139 |
|
coefficient C{name}. |
1140 |
|
|
1141 |
|
@param name: name of the coefficient enquired |
1142 |
|
@type name: C{string} |
1143 |
|
@return: the function space to be used for coefficient C{name} |
1144 |
|
@rtype: L{FunctionSpace<escript.FunctionSpace>} |
1145 |
|
@raise IllegalCoefficient: if C{name} is not a coefficient of the PDE |
1146 |
|
""" |
1147 |
|
if self.hasCoefficient(name): |
1148 |
|
return self.__COEFFICIENTS[name].getFunctionSpace(self.getDomain()) |
1149 |
|
else: |
1150 |
|
raise ValueError,"unknown coefficient %s requested"%name |
1151 |
|
|
1152 |
|
def getShapeOfCoefficient(self,name): |
1153 |
|
""" |
1154 |
|
Returns the shape of the coefficient C{name}. |
1155 |
|
|
1156 |
|
@param name: name of the coefficient enquired |
1157 |
|
@type name: C{string} |
1158 |
|
@return: the shape of the coefficient C{name} |
1159 |
|
@rtype: C{tuple} of C{int} |
1160 |
|
@raise IllegalCoefficient: if C{name} is not a coefficient of the PDE |
1161 |
|
""" |
1162 |
|
if self.hasCoefficient(name): |
1163 |
|
return self.__COEFFICIENTS[name].getShape(self.getDomain(),self.getNumEquations(),self.getNumSolutions()) |
1164 |
|
else: |
1165 |
|
raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name |
1166 |
|
|
1167 |
|
def resetAllCoefficients(self): |
1168 |
|
""" |
1169 |
|
Resets all coefficients to their default values. |
1170 |
|
""" |
1171 |
|
for i in self.__COEFFICIENTS.iterkeys(): |
1172 |
|
self.__COEFFICIENTS[i].resetValue() |
1173 |
|
|
1174 |
|
def alteredCoefficient(self,name): |
1175 |
|
""" |
1176 |
|
Announces that coefficient C{name} has been changed. |
1177 |
|
|
1178 |
|
@param name: name of the coefficient affected |
1179 |
|
@type name: C{string} |
1180 |
|
@raise IllegalCoefficient: if C{name} is not a coefficient of the PDE |
1181 |
|
@note: if C{name} is q or r, the method will not trigger a rebuild of the |
1182 |
|
system as constraints are applied to the solved system. |
1183 |
|
""" |
1184 |
|
if self.hasCoefficient(name): |
1185 |
|
self.trace("Coefficient %s has been altered."%name) |
1186 |
|
if not ((name=="q" or name=="r") and self.isUsingLumping()): |
1187 |
|
if self.__COEFFICIENTS[name].isAlteringOperator(): self.invalidateOperator() |
1188 |
|
if self.__COEFFICIENTS[name].isAlteringRightHandSide(): self.invalidateRightHandSide() |
1189 |
|
else: |
1190 |
|
raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name |
1191 |
|
|
1192 |
|
def validSolution(self): |
1193 |
""" |
""" |
1194 |
@brief |
Marks the solution as valid. |
1195 |
""" |
""" |
1196 |
return self.getDomain().newOperator( \ |
self.__is_solution_valid=True |
|
self.getNumEquations(), \ |
|
|
self.getFunctionSpaceForEquation(), \ |
|
|
self.getNumSolutions(), \ |
|
|
self.getFunctionSpaceForSolution(), \ |
|
|
self.__matrix_type) |
|
1197 |
|
|
1198 |
def __makeNewRightHandSide(self): |
def invalidateSolution(self): |
1199 |
""" |
""" |
1200 |
@brief |
Indicates the PDE has to be resolved if the solution is requested. |
1201 |
""" |
""" |
1202 |
return escript.Data(0.,(self.getNumEquations(),),self.getFunctionSpaceForEquation(),True) |
self.trace("System will be resolved.") |
1203 |
|
self.__is_solution_valid=False |
1204 |
|
|
1205 |
def __makeNewSolution(self): |
def isSolutionValid(self): |
1206 |
""" |
""" |
1207 |
@brief |
Returns True if the solution is still valid. |
1208 |
""" |
""" |
1209 |
return escript.Data(0.,(self.getNumSolutions(),),self.getFunctionSpaceForSolution(),True) |
return self.__is_solution_valid |
1210 |
|
|
1211 |
def __getFreshOperator(self): |
def validOperator(self): |
1212 |
""" |
""" |
1213 |
@brief |
Marks the operator as valid. |
1214 |
""" |
""" |
1215 |
if self.__operator.isEmpty(): |
self.__is_operator_valid=True |
|
self.__operator=self.__makeNewOperator() |
|
|
if self.debug() : print "PDE Debug: New operator allocated" |
|
|
else: |
|
|
self.__operator.setValue(0.) |
|
|
if self.debug() : print "PDE Debug: Operator reset to zero" |
|
|
return self.__operator |
|
1216 |
|
|
1217 |
def __getFreshRightHandSide(self): |
def invalidateOperator(self): |
1218 |
""" |
""" |
1219 |
@brief |
Indicates the operator has to be rebuilt next time it is used. |
1220 |
""" |
""" |
1221 |
if self.__righthandside.isEmpty(): |
self.trace("Operator will be rebuilt.") |
1222 |
self.__righthandside=self.__makeNewRightHandSide() |
self.invalidateSolution() |
1223 |
if self.debug() : print "PDE Debug: New right hand side allocated" |
self.__is_operator_valid=False |
1224 |
else: |
|
1225 |
print "fix self.__righthandside*=0" |
def isOperatorValid(self): |
1226 |
self.__righthandside*=0. |
""" |
1227 |
if self.debug() : print "PDE Debug: Right hand side reset to zero" |
Returns True if the operator is still valid. |
1228 |
return self.__righthandside |
""" |
1229 |
|
return self.__is_operator_valid |
1230 |
|
|
1231 |
|
def validRightHandSide(self): |
1232 |
|
""" |
1233 |
|
Marks the right hand side as valid. |
1234 |
|
""" |
1235 |
|
self.__is_RHS_valid=True |
1236 |
|
|
1237 |
# ==== rebuild switches ===================================================================== |
def invalidateRightHandSide(self): |
|
def __rebuildSolution(self,deep=False): |
|
1238 |
""" |
""" |
1239 |
@brief indicates the PDE has to be reolved if the solution is requested |
Indicates the right hand side has to be rebuilt next time it is used. |
1240 |
""" |
""" |
1241 |
if self.__solution_isValid and self.debug() : print "PDE Debug: PDE has to be resolved." |
if self.isRightHandSideValid(): self.trace("Right hand side has to be rebuilt.") |
1242 |
self.__solution_isValid=False |
self.invalidateSolution() |
1243 |
if deep: self.__solution=escript.Data(deep) |
self.__is_RHS_valid=False |
1244 |
|
|
1245 |
|
def isRightHandSideValid(self): |
1246 |
|
""" |
1247 |
|
Returns True if the operator is still valid. |
1248 |
|
""" |
1249 |
|
return self.__is_RHS_valid |
1250 |
|
|
1251 |
def __rebuildOperator(self,deep=False): |
def invalidateSystem(self): |
1252 |
""" |
""" |
1253 |
@brief indicates the operator has to be rebuilt next time it is used |
Announces that everything has to be rebuilt. |
1254 |
""" |
""" |
1255 |
if self.__operator_isValid and self.debug() : print "PDE Debug: Operator has to be rebuilt." |
if self.isRightHandSideValid(): self.trace("System has to be rebuilt.") |
1256 |
self.__rebuildSolution(deep) |
self.invalidateSolution() |
1257 |
self.__operator_isValid=False |
self.invalidateOperator() |
1258 |
if deep: self.__operator=escript.Operator() |
self.invalidateRightHandSide() |
1259 |
|
|
1260 |
def __rebuildRightHandSide(self,deep=False): |
def isSystemValid(self): |
1261 |
""" |
""" |
1262 |
@brief indicates the right hand side has to be rebuild next time it is used |
Returns True if the system (including solution) is still vaild. |
1263 |
""" |
""" |
1264 |
if self.__righthandside_isValid and self.debug() : print "PDE Debug: Right hand side has to be rebuilt." |
return self.isSolutionValid() and self.isOperatorValid() and self.isRightHandSideValid() |
|
self.__rebuildSolution(deep) |
|
|
self.__righthandside_isValid=False |
|
|
if not self.__homogeneous_constraint: self.__rebuildOperator() |
|
|
if deep: self.__righthandside=escript.Data() |
|
1265 |
|
|
1266 |
def __rebuildSystem(self,deep=False): |
def initializeSystem(self): |
1267 |
""" |
""" |
1268 |
@brief annonced that all coefficient name has been changed |
Resets the system clearing the operator, right hand side and solution. |
1269 |
""" |
""" |
1270 |
self.__rebuildSolution(deep) |
self.trace("New System has been created.") |
1271 |
self.__rebuildOperator(deep) |
self.__operator=escript.Operator() |
1272 |
self.__rebuildRightHandSide(deep) |
self.__righthandside=escript.Data() |
1273 |
|
self.__solution=escript.Data() |
1274 |
def __checkMatrixType(self): |
self.invalidateSystem() |
1275 |
|
|
1276 |
|
def getOperator(self): |
1277 |
""" |
""" |
1278 |
@brief reassess the matrix type and, if needed, initiates an operator rebuild |
Returns the operator of the linear problem. |
1279 |
|
|
1280 |
|
@return: the operator of the problem |
1281 |
""" |
""" |
1282 |
new_matrix_type=self.getDomain().getSystemMatrixTypeId(self.getSolverMethod(),self.isSymmetric()) |
return self.getSystem()[0] |
|
if not new_matrix_type==self.__matrix_type: |
|
|
if self.debug() : print "PDE Debug: Matrix type is now %d."%new_matrix_type |
|
|
self.__matrix_type=new_matrix_type |
|
|
self.__rebuildOperator(deep=True) |
|
1283 |
|
|
1284 |
#============ assembling ======================================================= |
def getRightHandSide(self): |
1285 |
def __copyConstraint(self,u): |
""" |
1286 |
""" |
Returns the right hand side of the linear problem. |
1287 |
@brief copies the constrint condition into u |
|
1288 |
""" |
@return: the right hand side of the problem |
1289 |
q=self.getCoefficient("q") |
@rtype: L{Data<escript.Data>} |
1290 |
r=self.getCoefficient("r") |
""" |
1291 |
if not q.isEmpty(): |
return self.getSystem()[1] |
|
if r.isEmpty(): |
|
|
r2=escript.Data(0,u.getShape(),u.getFunctionSpace()) |
|
|
else: |
|
|
r2=escript.Data(r,u.getFunctionSpace()) |
|
|
u.copyWithMask(r2,escript.Data(q,u.getFunctionSpace())) |
|
1292 |
|
|
1293 |
def __applyConstraint(self,rhs_update=True): |
def createRightHandSide(self): |
1294 |
""" |
""" |
1295 |
@brief applies the constraints defined by q and r to the system |
Returns an instance of a new right hand side. |
1296 |
""" |
""" |
1297 |
q=self.getCoefficient("q") |
self.trace("New right hand side is allocated.") |
1298 |
r=self.getCoefficient("r") |
if self.getNumEquations()>1: |
1299 |
if not q.isEmpty() and not self.__operator.isEmpty(): |
return escript.Data(0.,(self.getNumEquations(),),self.getFunctionSpaceForEquation(),True) |
1300 |
# q is the row and column mask to indicate where constraints are set: |
else: |
1301 |
row_q=escript.Data(q,self.getFunctionSpaceForEquation()) |
return escript.Data(0.,(),self.getFunctionSpaceForEquation(),True) |
|
col_q=escript.Data(q,self.getFunctionSpaceForSolution()) |
|
|
u=self.__makeNewSolution() |
|
|
if r.isEmpty(): |
|
|
r_s=self.__makeNewSolution() |
|
|
else: |
|
|
r_s=escript.Data(r,self.getFunctionSpaceForSolution()) |
|
|
u.copyWithMask(r_s,col_q) |
|
|
if not self.__righthandside.isEmpty() and rhs_update: self.__righthandside-=self.__operator*u |
|
|
self.__operator.nullifyRowsAndCols(row_q,col_q,1.) |
|
1302 |
|
|
1303 |
def getOperator(self): |
def createSolution(self): |
1304 |
|
""" |
1305 |
|
Returns an instance of a new solution. |
1306 |
|
""" |
1307 |
|
self.trace("New solution is allocated.") |
1308 |
|
if self.getNumSolutions()>1: |
1309 |
|
return escript.Data(0.,(self.getNumSolutions(),),self.getFunctionSpaceForSolution(),True) |
1310 |
|
else: |
1311 |
|
return escript.Data(0.,(),self.getFunctionSpaceForSolution(),True) |
1312 |
|
|
1313 |
|
def resetSolution(self): |
1314 |
""" |
""" |
1315 |
@brief returns the operator of the PDE |
Sets the solution to zero. |
1316 |
""" |
""" |
1317 |
if not self.__operator_isValid: |
if self.__solution.isEmpty(): |
1318 |
# some Constraints are applying for a lumpled stifness matrix: |
self.__solution=self.createSolution() |
1319 |
|
else: |
1320 |
|
self.__solution.setToZero() |
1321 |
|
self.trace("Solution is reset to zero.") |
1322 |
|
|
1323 |
|
def setSolution(self,u): |
1324 |
|
""" |
1325 |
|
Sets the solution assuming that makes the system valid. |
1326 |
|
""" |
1327 |
|
self.__solution=u |
1328 |
|
self.validSolution() |
1329 |
|
|
1330 |
|
def getCurrentSolution(self): |
1331 |
|
""" |
1332 |
|
Returns the solution in its current state. |
1333 |
|
""" |
1334 |
|
if self.__solution.isEmpty(): self.__solution=self.createSolution() |
1335 |
|
return self.__solution |
1336 |
|
|
1337 |
|
def resetRightHandSide(self): |
1338 |
|
""" |
1339 |
|
Sets the right hand side to zero. |
1340 |
|
""" |
1341 |
|
if self.__righthandside.isEmpty(): |
1342 |
|
self.__righthandside=self.createRightHandSide() |
1343 |
|
else: |
1344 |
|
self.__righthandside.setToZero() |
1345 |
|
self.trace("Right hand side is reset to zero.") |
1346 |
|
|
1347 |
|
def getCurrentRightHandSide(self): |
1348 |
|
""" |
1349 |
|
Returns the right hand side in its current state. |
1350 |
|
""" |
1351 |
|
if self.__righthandside.isEmpty(): self.__righthandside=self.createRightHandSide() |
1352 |
|
return self.__righthandside |
1353 |
|
|
1354 |
|
def resetOperator(self): |
1355 |
|
""" |
1356 |
|
Makes sure that the operator is instantiated and returns it initialized |
1357 |
|
with zeros. |
1358 |
|
""" |
1359 |
|
if self.__operator.isEmpty(): |
1360 |
if self.isUsingLumping(): |
if self.isUsingLumping(): |
1361 |
if self.getFunctionSpaceForEquation()==self.getFunctionSpaceForSolution(): |
self.__operator=self.createSolution() |
|
raise TypeError,"Lumped matrix requires same order for equations and unknowns" |
|
|
if not self.getCoefficient("A").isEmpty(): |
|
|
raise Warning,"Lumped matrix does not allow coefficient A" |
|
|
if not self.getCoefficient("B").isEmpty(): |
|
|
raise Warning,"Lumped matrix does not allow coefficient B" |
|
|
if not self.getCoefficient("C").isEmpty(): |
|
|
raise Warning,"Lumped matrix does not allow coefficient C" |
|
|
if not self.getCoefficient("D").isEmpty(): |
|
|
raise Warning,"Lumped matrix does not allow coefficient D" |
|
|
if self.debug() : print "PDE Debug: New lumped operator is built." |
|
|
mat=self.__makeNewOperator(self) |
|
1362 |
else: |
else: |
1363 |
if self.debug() : print "PDE Debug: New operator is built." |
self.__operator=self.createOperator() |
1364 |
mat=self.__getFreshOperator() |
else: |
|
|
|
|
self.getDomain().addPDEToSystem(mat,escript.Data(), \ |
|
|
self.getCoefficient("A"), \ |
|
|
self.getCoefficient("B"), \ |
|
|
self.getCoefficient("C"), \ |
|
|
self.getCoefficient("D"), \ |
|
|
escript.Data(), \ |
|
|
escript.Data(), \ |
|
|
self.getCoefficient("d"), \ |
|
|
escript.Data(),\ |
|
|
self.getCoefficient("d_contact"), \ |
|
|
escript.Data()) |
|
1365 |
if self.isUsingLumping(): |
if self.isUsingLumping(): |
1366 |
self.__operator=mat*escript.Data(1,(self.getNumSolutions(),),self.getFunctionSpaceOfSolution(),True) |
self.__operator.setToZero() |
1367 |
else: |
else: |
1368 |
self.__applyConstraint(rhs_update=False) |
self.__operator.resetValues() |
1369 |
self.__operator_isValid=True |
self.trace("Operator reset to zero") |
1370 |
|
|
1371 |
|
def getCurrentOperator(self): |
1372 |
|
""" |
1373 |
|
Returns the operator in its current state. |
1374 |
|
""" |
1375 |
|
if self.__operator.isEmpty(): self.resetOperator() |
1376 |
return self.__operator |
return self.__operator |
1377 |
|
|
1378 |
def getRightHandSide(self,ignoreConstraint=False): |
def setValue(self,**coefficients): |
1379 |
|
""" |
1380 |
|
Sets new values to coefficients. |
1381 |
|
|
1382 |
|
@raise IllegalCoefficient: if an unknown coefficient keyword is used |
1383 |
|
""" |
1384 |
|
# check if the coefficients are legal: |
1385 |
|
for i in coefficients.iterkeys(): |
1386 |
|
if not self.hasCoefficient(i): |
1387 |
|
raise IllegalCoefficient,"Attempt to set unknown coefficient %s"%i |
1388 |
|
# if the number of unknowns or equations is still unknown we try to estimate them: |
1389 |
|
if self.__numEquations==None or self.__numSolutions==None: |
1390 |
|
for i,d in coefficients.iteritems(): |
1391 |
|
if hasattr(d,"shape"): |
1392 |
|
s=d.shape |
1393 |
|
elif hasattr(d,"getShape"): |
1394 |
|
s=d.getShape() |
1395 |
|
else: |
1396 |
|
s=numarray.array(d).shape |
1397 |
|
if s!=None: |
1398 |
|
# get number of equations and number of unknowns: |
1399 |
|
res=self.__COEFFICIENTS[i].estimateNumEquationsAndNumSolutions(self.getDomain(),s) |
1400 |
|
if res==None: |
1401 |
|
raise IllegalCoefficientValue,"Illegal shape %s of coefficient %s"%(s,i) |
1402 |
|
else: |
1403 |
|
if self.__numEquations==None: self.__numEquations=res[0] |
1404 |
|
if self.__numSolutions==None: self.__numSolutions=res[1] |
1405 |
|
if self.__numEquations==None: raise UndefinedPDEError,"unidentified number of equations" |
1406 |
|
if self.__numSolutions==None: raise UndefinedPDEError,"unidentified number of solutions" |
1407 |
|
# now we check the shape of the coefficient if numEquations and numSolutions are set: |
1408 |
|
for i,d in coefficients.iteritems(): |
1409 |
|
try: |
1410 |
|
self.__COEFFICIENTS[i].setValue(self.getDomain(), |
1411 |
|
self.getNumEquations(),self.getNumSolutions(), |
1412 |
|
self.reduceEquationOrder(),self.reduceSolutionOrder(),d) |
1413 |
|
self.alteredCoefficient(i) |
1414 |
|
except IllegalCoefficientFunctionSpace,m: |
1415 |
|
# if the function space is wrong then we try the reduced version: |
1416 |
|
i_red=i+"_reduced" |
1417 |
|
if (not i_red in coefficients.keys()) and i_red in self.__COEFFICIENTS.keys(): |
1418 |
|
try: |
1419 |
|
self.__COEFFICIENTS[i_red].setValue(self.getDomain(), |
1420 |
|
self.getNumEquations(),self.getNumSolutions(), |
1421 |
|
self.reduceEquationOrder(),self.reduceSolutionOrder(),d) |
1422 |
|
self.alteredCoefficient(i_red) |
1423 |
|
except IllegalCoefficientValue,m: |
1424 |
|
raise IllegalCoefficientValue("Coefficient %s:%s"%(i,m)) |
1425 |
|
except IllegalCoefficientFunctionSpace,m: |
1426 |
|
raise IllegalCoefficientFunctionSpace("Coefficient %s:%s"%(i,m)) |
1427 |
|
else: |
1428 |
|
raise IllegalCoefficientFunctionSpace("Coefficient %s:%s"%(i,m)) |
1429 |
|
except IllegalCoefficientValue,m: |
1430 |
|
raise IllegalCoefficientValue("Coefficient %s:%s"%(i,m)) |
1431 |
|
self.__altered_coefficients=True |
1432 |
|
|
1433 |
|
# ========================================================================== |
1434 |
|
# methods that are typically overwritten when implementing a particular |
1435 |
|
# linear problem |
1436 |
|
# ========================================================================== |
1437 |
|
def getRequiredSystemType(self): |
1438 |
|
""" |
1439 |
|
Returns the system type which needs to be used by the current set up. |
1440 |
|
|
1441 |
|
@note: Typically this method is overwritten when implementing a |
1442 |
|
particular linear problem. |
1443 |
|
""" |
1444 |
|
return 0 |
1445 |
|
|
1446 |
|
def createOperator(self): |
1447 |
""" |
""" |
1448 |
@brief returns the right hand side of the PDE |
Returns an instance of a new operator. |
1449 |
|
|
1450 |
@param ignoreConstraint |
@note: This method is overwritten when implementing a particular |
1451 |
|
linear problem. |
1452 |
""" |
""" |
1453 |
if not self.__righthandside_isValid: |
return escript.Operator() |
1454 |
if self.debug() : print "PDE Debug: New right hand side is built." |
|
1455 |
self.getDomain().addPDEToRHS(self.__getFreshRightHandSide(), \ |
def checkSymmetry(self,verbose=True): |
1456 |
self.getCoefficient("X"), \ |
""" |
1457 |
self.getCoefficient("Y"),\ |
Tests the PDE for symmetry. |
1458 |
self.getCoefficient("y"),\ |
|
1459 |
self.getCoefficient("y_contact")) |
@param verbose: if set to True or not present a report on coefficients |
1460 |
self.__righthandside_isValid=True |
which break the symmetry is printed |
1461 |
if ignoreConstraint: self.__copyConstraint(self.__righthandside) |
@type verbose: C{bool} |
1462 |
return self.__righthandside |
@return: True if the problem is symmetric |
1463 |
|
@rtype: C{bool} |
1464 |
|
@note: Typically this method is overwritten when implementing a |
1465 |
|
particular linear problem. |
1466 |
|
""" |
1467 |
|
out=True |
1468 |
|
return out |
1469 |
|
|
1470 |
|
def getSolution(self,**options): |
1471 |
|
""" |
1472 |
|
Returns the solution of the problem. |
1473 |
|
|
1474 |
|
@return: the solution |
1475 |
|
@rtype: L{Data<escript.Data>} |
1476 |
|
|
1477 |
|
@note: This method is overwritten when implementing a particular |
1478 |
|
linear problem. |
1479 |
|
""" |
1480 |
|
return self.getCurrentSolution() |
1481 |
|
|
1482 |
|
def getSystem(self): |
1483 |
|
""" |
1484 |
|
Returns the operator and right hand side of the PDE. |
1485 |
|
|
1486 |
|
@return: the discrete version of the PDE |
1487 |
|
@rtype: C{tuple} of L{Operator,<escript.Operator>} and L{Data<escript.Data>}. |
1488 |
|
|
1489 |
|
@note: This method is overwritten when implementing a particular |
1490 |
|
linear problem. |
1491 |
|
""" |
1492 |
|
return (self.getCurrentOperator(), self.getCurrentRightHandSide()) |
1493 |
|
#===================== |
1494 |
|
|
1495 |
|
class LinearPDE(LinearProblem): |
1496 |
|
""" |
1497 |
|
This class is used to define a general linear, steady, second order PDE |
1498 |
|
for an unknown function M{u} on a given domain defined through a |
1499 |
|
L{Domain<escript.Domain>} object. |
1500 |
|
|
1501 |
|
For a single PDE having a solution with a single component the linear PDE |
1502 |
|
is defined in the following form: |
1503 |
|
|
1504 |
|
M{-(grad(A[j,l]+A_reduced[j,l])*grad(u)[l]+(B[j]+B_reduced[j])u)[j]+(C[l]+C_reduced[l])*grad(u)[l]+(D+D_reduced)=-grad(X+X_reduced)[j,j]+(Y+Y_reduced)} |
1505 |
|
|
1506 |
|
where M{grad(F)} denotes the spatial derivative of M{F}. Einstein's |
1507 |
|
summation convention, ie. summation over indexes appearing twice in a term |
1508 |
|
of a sum performed, is used. |
1509 |
|
The coefficients M{A}, M{B}, M{C}, M{D}, M{X} and M{Y} have to be specified |
1510 |
|
through L{Data<escript.Data>} objects in L{Function<escript.Function>} and |
1511 |
|
the coefficients M{A_reduced}, M{B_reduced}, M{C_reduced}, M{D_reduced}, |
1512 |
|
M{X_reduced} and M{Y_reduced} have to be specified through |
1513 |
|
L{Data<escript.Data>} objects in L{ReducedFunction<escript.ReducedFunction>}. |
1514 |
|
It is also allowed to use objects that can be converted into such |
1515 |
|
L{Data<escript.Data>} objects. M{A} and M{A_reduced} are rank two, M{B}, |
1516 |
|
M{C}, M{X}, M{B_reduced}, M{C_reduced} and M{X_reduced} are rank one and |
1517 |
|
M{D}, M{D_reduced}, M{Y} and M{Y_reduced} are scalar. |
1518 |
|
|
1519 |
|
The following natural boundary conditions are considered: |
1520 |
|
|
1521 |
|
M{n[j]*((A[i,j]+A_reduced[i,j])*grad(u)[l]+(B+B_reduced)[j]*u)+(d+d_reduced)*u=n[j]*(X[j]+X_reduced[j])+y} |
1522 |
|
|
1523 |
|
where M{n} is the outer normal field. Notice that the coefficients M{A}, |
1524 |
|
M{A_reduced}, M{B}, M{B_reduced}, M{X} and M{X_reduced} are defined in the |
1525 |
|
PDE. The coefficients M{d} and M{y} are each a scalar in |
1526 |
|
L{FunctionOnBoundary<escript.FunctionOnBoundary>} and the coefficients |
1527 |
|
M{d_reduced} and M{y_reduced} are each a scalar in |
1528 |
|
L{ReducedFunctionOnBoundary<escript.ReducedFunctionOnBoundary>}. |
1529 |
|
|
1530 |
|
Constraints for the solution prescribe the value of the solution at certain |
1531 |
|
locations in the domain. They have the form |
1532 |
|
|
1533 |
|
M{u=r} where M{q>0} |
1534 |
|
|
1535 |
|
M{r} and M{q} are each scalar where M{q} is the characteristic function |
1536 |
|
defining where the constraint is applied. The constraints override any |
1537 |
|
other condition set by the PDE or the boundary condition. |
1538 |
|
|
1539 |
|
The PDE is symmetrical if |
1540 |
|
|
1541 |
|
M{A[i,j]=A[j,i]} and M{B[j]=C[j]} and M{A_reduced[i,j]=A_reduced[j,i]} |
1542 |
|
and M{B_reduced[j]=C_reduced[j]} |
1543 |
|
|
1544 |
|
For a system of PDEs and a solution with several components the PDE has the |
1545 |
|
form |
1546 |
|
|
1547 |
|
M{-grad((A[i,j,k,l]+A_reduced[i,j,k,l])*grad(u[k])[l]+(B[i,j,k]+B_reduced[i,j,k])*u[k])[j]+(C[i,k,l]+C_reduced[i,k,l])*grad(u[k])[l]+(D[i,k]+D_reduced[i,k]*u[k] =-grad(X[i,j]+X_reduced[i,j])[j]+Y[i]+Y_reduced[i] } |
1548 |
|
|
1549 |
|
M{A} and M{A_reduced} are of rank four, M{B}, M{B_reduced}, M{C} and |
1550 |
|
M{C_reduced} are each of rank three, M{D}, M{D_reduced}, M{X_reduced} and |
1551 |
|
M{X} are each of rank two and M{Y} and M{Y_reduced} are of rank one. |
1552 |
|
The natural boundary conditions take the form: |
1553 |
|
|
1554 |
|
M{n[j]*((A[i,j,k,l]+A_reduced[i,j,k,l])*grad(u[k])[l]+(B[i,j,k]+B_reduced[i,j,k])*u[k])+(d[i,k]+d_reduced[i,k])*u[k]=n[j]*(X[i,j]+X_reduced[i,j])+y[i]+y_reduced[i]} |
1555 |
|
|
1556 |
|
The coefficient M{d} is of rank two and M{y} is of rank one both in |
1557 |
|
L{FunctionOnBoundary<escript.FunctionOnBoundary>}. The coefficients |
1558 |
|
M{d_reduced} is of rank two and M{y_reduced} is of rank one both in |
1559 |
|
L{ReducedFunctionOnBoundary<escript.ReducedFunctionOnBoundary>}. |
1560 |
|
|
1561 |
|
Constraints take the form |
1562 |
|
|
1563 |
|
M{u[i]=r[i]} where M{q[i]>0} |
1564 |
|
|
1565 |
|
M{r} and M{q} are each rank one. Notice that at some locations not |
1566 |
|
necessarily all components must have a constraint. |
1567 |
|
|
1568 |
|
The system of PDEs is symmetrical if |
1569 |
|
|
1570 |
|
- M{A[i,j,k,l]=A[k,l,i,j]} |
1571 |
|
- M{A_reduced[i,j,k,l]=A_reduced[k,l,i,j]} |
1572 |
|
- M{B[i,j,k]=C[k,i,j]} |
1573 |
|
- M{B_reduced[i,j,k]=C_reduced[k,i,j]} |
1574 |
|
- M{D[i,k]=D[i,k]} |
1575 |
|
- M{D_reduced[i,k]=D_reduced[i,k]} |
1576 |
|
- M{d[i,k]=d[k,i]} |
1577 |
|
- M{d_reduced[i,k]=d_reduced[k,i]} |
1578 |
|
|
1579 |
|
L{LinearPDE} also supports solution discontinuities over a contact region |
1580 |
|
in the domain. To specify the conditions across the discontinuity we are |
1581 |
|
using the generalised flux M{J} which, in the case of a system of PDEs |
1582 |
|
and several components of the solution, is defined as |
1583 |
|
|
1584 |
|
M{J[i,j]=(A[i,j,k,l]+A_reduced[[i,j,k,l])*grad(u[k])[l]+(B[i,j,k]+B_reduced[i,j,k])*u[k]-X[i,j]-X_reduced[i,j]} |
1585 |
|
|
1586 |
|
For the case of single solution component and single PDE M{J} is defined as |
1587 |
|
|
1588 |
|
M{J[j]=(A[i,j]+A_reduced[i,j])*grad(u)[j]+(B[i]+B_reduced[i])*u-X[i]-X_reduced[i]} |
1589 |
|
|
1590 |
|
In the context of discontinuities M{n} denotes the normal on the |
1591 |
|
discontinuity pointing from side 0 towards side 1 calculated from |
1592 |
|
L{getNormal<escript.FunctionSpace.getNormal>} of L{FunctionOnContactZero<escript.FunctionOnContactZero>}. |
1593 |
|
For a system of PDEs the contact condition takes the form |
1594 |
|
|
1595 |
|
M{n[j]*J0[i,j]=n[j]*J1[i,j]=(y_contact[i]+y_contact_reduced[i])- (d_contact[i,k]+d_contact_reduced[i,k])*jump(u)[k]} |
1596 |
|
|
1597 |
|
where M{J0} and M{J1} are the fluxes on side 0 and side 1 of the |
1598 |
|
discontinuity, respectively. M{jump(u)}, which is the difference of the |
1599 |
|
solution at side 1 and at side 0, denotes the jump of M{u} across |
1600 |
|
discontinuity along the normal calculated by L{jump<util.jump>}. |
1601 |
|
The coefficient M{d_contact} is of rank two and M{y_contact} is of rank one |
1602 |
|
both in L{FunctionOnContactZero<escript.FunctionOnContactZero>} or |
1603 |
|
L{FunctionOnContactOne<escript.FunctionOnContactOne>}. |
1604 |
|
The coefficient M{d_contact_reduced} is of rank two and M{y_contact_reduced} |
1605 |
|
is of rank one both in L{ReducedFunctionOnContactZero<escript.ReducedFunctionOnContactZero>} |
1606 |
|
or L{ReducedFunctionOnContactOne<escript.ReducedFunctionOnContactOne>}. |
1607 |
|
In case of a single PDE and a single component solution the contact |
1608 |
|
condition takes the form |
1609 |
|
|
1610 |
|
M{n[j]*J0_{j}=n[j]*J1_{j}=(y_contact+y_contact_reduced)-(d_contact+y_contact_reduced)*jump(u)} |
1611 |
|
|
1612 |
|
In this case the coefficient M{d_contact} and M{y_contact} are each scalar |
1613 |
|
both in L{FunctionOnContactZero<escript.FunctionOnContactZero>} or |
1614 |
|
L{FunctionOnContactOne<escript.FunctionOnContactOne>} and the coefficient |
1615 |
|
M{d_contact_reduced} and M{y_contact_reduced} are each scalar both in |
1616 |
|
L{ReducedFunctionOnContactZero<escript.ReducedFunctionOnContactZero>} or |
1617 |
|
L{ReducedFunctionOnContactOne<escript.ReducedFunctionOnContactOne>}. |
1618 |
|
|
1619 |
|
Typical usage:: |
1620 |
|
|
1621 |
|
p = LinearPDE(dom) |
1622 |
|
p.setValue(A=kronecker(dom), D=1, Y=0.5) |
1623 |
|
u = p.getSolution() |
1624 |
|
|
1625 |
|
""" |
1626 |
|
|
1627 |
|
def __init__(self,domain,numEquations=None,numSolutions=None,debug=False): |
1628 |
|
""" |
1629 |
|
Initializes a new linear PDE. |
1630 |
|
|
1631 |
|
@param domain: domain of the PDE |
1632 |
|
@type domain: L{Domain<escript.Domain>} |
1633 |
|
@param numEquations: number of equations. If C{None} the number of |
1634 |
|
equations is extracted from the PDE coefficients. |
1635 |
|
@param numSolutions: number of solution components. If C{None} the number |
1636 |
|
of solution components is extracted from the PDE |
1637 |
|
coefficients. |
1638 |
|
@param debug: if True debug information is printed |
1639 |
|
|
1640 |
|
""" |
1641 |
|
super(LinearPDE, self).__init__(domain,numEquations,numSolutions,debug) |
1642 |
|
# |
1643 |
|
# the coefficients of the PDE: |
1644 |
|
# |
1645 |
|
self.introduceCoefficients( |
1646 |
|
A=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,PDECoef.BY_DIM,PDECoef.BY_SOLUTION,PDECoef.BY_DIM),PDECoef.OPERATOR), |
1647 |
|
B=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,PDECoef.BY_DIM,PDECoef.BY_SOLUTION),PDECoef.OPERATOR), |
1648 |
|
C=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION,PDECoef.BY_DIM),PDECoef.OPERATOR), |
1649 |
|
D=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION),PDECoef.OPERATOR), |
1650 |
|
X=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,PDECoef.BY_DIM),PDECoef.RIGHTHANDSIDE), |
1651 |
|
Y=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE), |
1652 |
|
d=PDECoef(PDECoef.BOUNDARY,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION),PDECoef.OPERATOR), |
1653 |
|
y=PDECoef(PDECoef.BOUNDARY,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE), |
1654 |
|
d_contact=PDECoef(PDECoef.CONTACT,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION),PDECoef.OPERATOR), |
1655 |
|
y_contact=PDECoef(PDECoef.CONTACT,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE), |
1656 |
|
A_reduced=PDECoef(PDECoef.INTERIOR_REDUCED,(PDECoef.BY_EQUATION,PDECoef.BY_DIM,PDECoef.BY_SOLUTION,PDECoef.BY_DIM),PDECoef.OPERATOR), |
1657 |
|
B_reduced=PDECoef(PDECoef.INTERIOR_REDUCED,(PDECoef.BY_EQUATION,PDECoef.BY_DIM,PDECoef.BY_SOLUTION),PDECoef.OPERATOR), |
1658 |
|
C_reduced=PDECoef(PDECoef.INTERIOR_REDUCED,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION,PDECoef.BY_DIM),PDECoef.OPERATOR), |
1659 |
|
D_reduced=PDECoef(PDECoef.INTERIOR_REDUCED,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION),PDECoef.OPERATOR), |
1660 |
|
X_reduced=PDECoef(PDECoef.INTERIOR_REDUCED,(PDECoef.BY_EQUATION,PDECoef.BY_DIM),PDECoef.RIGHTHANDSIDE), |
1661 |
|
Y_reduced=PDECoef(PDECoef.INTERIOR_REDUCED,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE), |
1662 |
|
d_reduced=PDECoef(PDECoef.BOUNDARY_REDUCED,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION),PDECoef.OPERATOR), |
1663 |
|
y_reduced=PDECoef(PDECoef.BOUNDARY_REDUCED,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE), |
1664 |
|
d_contact_reduced=PDECoef(PDECoef.CONTACT_REDUCED,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION),PDECoef.OPERATOR), |
1665 |
|
y_contact_reduced=PDECoef(PDECoef.CONTACT_REDUCED,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE), |
1666 |
|
r=PDECoef(PDECoef.SOLUTION,(PDECoef.BY_SOLUTION,),PDECoef.RIGHTHANDSIDE), |
1667 |
|
q=PDECoef(PDECoef.SOLUTION,(PDECoef.BY_SOLUTION,),PDECoef.BOTH) ) |
1668 |
|
|
1669 |
|
def __str__(self): |
1670 |
|
""" |
1671 |
|
Returns the string representation of the PDE. |
1672 |
|
|
1673 |
|
@return: a simple representation of the PDE |
1674 |
|
@rtype: C{str} |
1675 |
|
""" |
1676 |
|
return "<LinearPDE %d>"%id(self) |
1677 |
|
|
1678 |
|
def getRequiredSystemType(self): |
1679 |
|
""" |
1680 |
|
Returns the system type which needs to be used by the current set up. |
1681 |
|
""" |
1682 |
|
return self.getDomain().getSystemMatrixTypeId(self.getSolverMethod()[0],self.getSolverMethod()[1],self.getSolverPackage(),self.isSymmetric()) |
1683 |
|
|
1684 |
|
def checkSymmetry(self,verbose=True): |
1685 |
|
""" |
1686 |
|
Tests the PDE for symmetry. |
1687 |
|
|
1688 |
|
@param verbose: if set to True or not present a report on coefficients |
1689 |
|
which break the symmetry is printed. |
1690 |
|
@type verbose: C{bool} |
1691 |
|
@return: True if the PDE is symmetric |
1692 |
|
@rtype: L{bool} |
1693 |
|
@note: This is a very expensive operation. It should be used for |
1694 |
|
degugging only! The symmetry flag is not altered. |
1695 |
|
""" |
1696 |
|
out=True |
1697 |
|
out=out and self.checkSymmetricTensor("A", verbose) |
1698 |
|
out=out and self.checkSymmetricTensor("A_reduced", verbose) |
1699 |
|
out=out and self.checkReciprocalSymmetry("B","C", verbose) |
1700 |
|
out=out and self.checkReciprocalSymmetry("B_reduced","C_reduced", verbose) |
1701 |
|
out=out and self.checkSymmetricTensor("D", verbose) |
1702 |
|
out=out and self.checkSymmetricTensor("D_reduced", verbose) |
1703 |
|
out=out and self.checkSymmetricTensor("d", verbose) |
1704 |
|
out=out and self.checkSymmetricTensor("d_reduced", verbose) |
1705 |
|
out=out and self.checkSymmetricTensor("d_contact", verbose) |
1706 |
|
out=out and self.checkSymmetricTensor("d_contact_reduced", verbose) |
1707 |
|
return out |
1708 |
|
|
1709 |
|
def createOperator(self): |
1710 |
|
""" |
1711 |
|
Returns an instance of a new operator. |
1712 |
|
""" |
1713 |
|
self.trace("New operator is allocated.") |
1714 |
|
return self.getDomain().newOperator( \ |
1715 |
|
self.getNumEquations(), \ |
1716 |
|
self.getFunctionSpaceForEquation(), \ |
1717 |
|
self.getNumSolutions(), \ |
1718 |
|
self.getFunctionSpaceForSolution(), \ |
1719 |
|
self.getSystemType()) |
1720 |
|
|
1721 |
|
def getSolution(self,**options): |
1722 |
|
""" |
1723 |
|
Returns the solution of the PDE. |
1724 |
|
|
1725 |
|
@return: the solution |
1726 |
|
@rtype: L{Data<escript.Data>} |
1727 |
|
@param options: solver options |
1728 |
|
@keyword verbose: True to get some information during PDE solution |
1729 |
|
@type verbose: C{bool} |
1730 |
|
@keyword reordering: reordering scheme to be used during elimination. |
1731 |
|
Allowed values are L{NO_REORDERING}, |
1732 |
|
L{MINIMUM_FILL_IN} and L{NESTED_DISSECTION} |
1733 |
|
@keyword iter_max: maximum number of iteration steps allowed |
1734 |
|
@keyword drop_tolerance: threshold for dropping in L{ILUT} |
1735 |
|
@keyword drop_storage: maximum of allowed memory in L{ILUT} |
1736 |
|
@keyword truncation: maximum number of residuals in L{GMRES} |
1737 |
|
@keyword restart: restart cycle length in L{GMRES} |
1738 |
|
""" |
1739 |
|
if not self.isSolutionValid(): |
1740 |
|
mat,f=self.getSystem() |
1741 |
|
if self.isUsingLumping(): |
1742 |
|
self.setSolution(f*1/mat) |
1743 |
|
else: |
1744 |
|
options[self.TOLERANCE_KEY]=self.getTolerance() |
1745 |
|
options[self.METHOD_KEY]=self.getSolverMethod()[0] |
1746 |
|
options[self.PRECONDITIONER_KEY]=self.getSolverMethod()[1] |
1747 |
|
options[self.PACKAGE_KEY]=self.getSolverPackage() |
1748 |
|
options[self.SYMMETRY_KEY]=self.isSymmetric() |
1749 |
|
self.trace("PDE is resolved.") |
1750 |
|
self.trace("solver options: %s"%str(options)) |
1751 |
|
self.setSolution(mat.solve(f,options)) |
1752 |
|
return self.getCurrentSolution() |
1753 |
|
|
1754 |
def getSystem(self): |
def getSystem(self): |
1755 |
""" |
""" |
1756 |
@brief |
Returns the operator and right hand side of the PDE. |
1757 |
|
|
1758 |
|
@return: the discrete version of the PDE |
1759 |
|
@rtype: C{tuple} of L{Operator,<escript.Operator>} and |
1760 |
|
L{Data<escript.Data>} |
1761 |
""" |
""" |
1762 |
if not self.__operator_isValid and not self.__righthandside_isValid: |
if not self.isOperatorValid() or not self.isRightHandSideValid(): |
|
if self.debug() : print "PDE Debug: New PDE is built." |
|
1763 |
if self.isUsingLumping(): |
if self.isUsingLumping(): |
1764 |
self.getRightHandSide(ignoreConstraint=True) |
if not self.isOperatorValid(): |
1765 |
self.getOperator() |
if not self.getFunctionSpaceForEquation()==self.getFunctionSpaceForSolution(): |
1766 |
|
raise TypeError,"Lumped matrix requires same order for equations and unknowns" |
1767 |
|
if not self.getCoefficient("A").isEmpty(): |
1768 |
|
raise ValueError,"coefficient A in lumped matrix may not be present." |
1769 |
|
if not self.getCoefficient("B").isEmpty(): |
1770 |
|
raise ValueError,"coefficient B in lumped matrix may not be present." |
1771 |
|
if not self.getCoefficient("C").isEmpty(): |
1772 |
|
raise ValueError,"coefficient C in lumped matrix may not be present." |
1773 |
|
if not self.getCoefficient("d_contact").isEmpty(): |
1774 |
|
raise ValueError,"coefficient d_contact in lumped matrix may not be present." |
1775 |
|
if not self.getCoefficient("A_reduced").isEmpty(): |
1776 |
|
raise ValueError,"coefficient A_reduced in lumped matrix may not be present." |
1777 |
|
if not self.getCoefficient("B_reduced").isEmpty(): |
1778 |
|
raise ValueError,"coefficient B_reduced in lumped matrix may not be present." |
1779 |
|
if not self.getCoefficient("C_reduced").isEmpty(): |
1780 |
|
raise ValueError,"coefficient C_reduced in lumped matrix may not be present." |
1781 |
|
if not self.getCoefficient("d_contact_reduced").isEmpty(): |
1782 |
|
raise ValueError,"coefficient d_contact_reduced in lumped matrix may not be present." |
1783 |
|
D=self.getCoefficient("D") |
1784 |
|
d=self.getCoefficient("d") |
1785 |
|
D_reduced=self.getCoefficient("D_reduced") |
1786 |
|
d_reduced=self.getCoefficient("d_reduced") |
1787 |
|
if not D.isEmpty(): |
1788 |
|
if self.getNumSolutions()>1: |
1789 |
|
D_times_e=util.matrix_mult(D,numarray.ones((self.getNumSolutions(),))) |
1790 |
|
else: |
1791 |
|
D_times_e=D |
1792 |
|
else: |
1793 |
|
D_times_e=escript.Data() |
1794 |
|
if not d.isEmpty(): |
1795 |
|
if self.getNumSolutions()>1: |
1796 |
|
d_times_e=util.matrix_mult(d,numarray.ones((self.getNumSolutions(),))) |
1797 |
|
else: |
1798 |
|
d_times_e=d |
1799 |
|
else: |
1800 |
|
d_times_e=escript.Data() |
1801 |
|
|
1802 |
|
if not D_reduced.isEmpty(): |
1803 |
|
if self.getNumSolutions()>1: |
1804 |
|
D_reduced_times_e=util.matrix_mult(D_reduced,numarray.ones((self.getNumSolutions(),))) |
1805 |
|
else: |
1806 |
|
D_reduced_times_e=D_reduced |
1807 |
|
else: |
1808 |
|
D_reduced_times_e=escript.Data() |
1809 |
|
if not d_reduced.isEmpty(): |
1810 |
|
if self.getNumSolutions()>1: |
1811 |
|
d_reduced_times_e=util.matrix_mult(d_reduced,numarray.ones((self.getNumSolutions(),))) |
1812 |
|
else: |
1813 |
|
d_reduced_times_e=d_reduced |
1814 |
|
else: |
1815 |
|
d_reduced_times_e=escript.Data() |
1816 |
|
|
1817 |
|
self.resetOperator() |
1818 |
|
operator=self.getCurrentOperator() |
1819 |
|
if False and hasattr(self.getDomain(), "addPDEToLumpedSystem") : |
1820 |
|
self.getDomain().addPDEToLumpedSystem(operator, D_times_e, d_times_e) |
1821 |
|
self.getDomain().addPDEToLumpedSystem(operator, D_reduced_times_e, d_reduced_times_e) |
1822 |
|
else: |
1823 |
|
self.getDomain().addPDEToRHS(operator, \ |
1824 |
|
escript.Data(), \ |
1825 |
|
D_times_e, \ |
1826 |
|
d_times_e,\ |
1827 |
|
escript.Data()) |
1828 |
|
self.getDomain().addPDEToRHS(operator, \ |
1829 |
|
escript.Data(), \ |
1830 |
|
D_reduced_times_e, \ |
1831 |
|
d_reduced_times_e,\ |
1832 |
|
escript.Data()) |
1833 |
|
self.trace("New lumped operator has been built.") |
1834 |
|
if not self.isRightHandSideValid(): |
1835 |
|
self.resetRightHandSide() |
1836 |
|
righthandside=self.getCurrentRightHandSide() |
1837 |
|
self.getDomain().addPDEToRHS(righthandside, \ |
1838 |
|
self.getCoefficient("X"), \ |
1839 |
|
self.getCoefficient("Y"),\ |
1840 |
|
self.getCoefficient("y"),\ |
1841 |
|
self.getCoefficient("y_contact")) |
1842 |
|
self.getDomain().addPDEToRHS(righthandside, \ |
1843 |
|
self.getCoefficient("X_reduced"), \ |
1844 |
|
self.getCoefficient("Y_reduced"),\ |
1845 |
|
self.getCoefficient("y_reduced"),\ |
1846 |
|
self.getCoefficient("y_contact_reduced")) |
1847 |
|
self.trace("New right hand side has been built.") |
1848 |
|
self.validRightHandSide() |
1849 |
|
self.insertConstraint(rhs_only=False) |
1850 |
|
self.validOperator() |
1851 |
else: |
else: |
1852 |
self.getDomain().addPDEToSystem(self.__getFreshOperator(),self.__getFreshRightHandSide(), \ |
if not self.isOperatorValid() and not self.isRightHandSideValid(): |
1853 |
|
self.resetRightHandSide() |
1854 |
|
righthandside=self.getCurrentRightHandSide() |
1855 |
|
self.resetOperator() |
1856 |
|
operator=self.getCurrentOperator() |
1857 |
|
self.getDomain().addPDEToSystem(operator,righthandside, \ |
1858 |
|
self.getCoefficient("A"), \ |
1859 |
|
self.getCoefficient("B"), \ |
1860 |
|
self.getCoefficient("C"), \ |
1861 |
|
self.getCoefficient("D"), \ |
1862 |
|
self.getCoefficient("X"), \ |
1863 |
|
self.getCoefficient("Y"), \ |
1864 |
|
self.getCoefficient("d"), \ |
1865 |
|
self.getCoefficient("y"), \ |
1866 |
|
self.getCoefficient("d_contact"), \ |
1867 |
|
self.getCoefficient("y_contact")) |
1868 |
|
self.getDomain().addPDEToSystem(operator,righthandside, \ |
1869 |
|
self.getCoefficient("A_reduced"), \ |
1870 |
|
self.getCoefficient("B_reduced"), \ |
1871 |
|
self.getCoefficient("C_reduced"), \ |
1872 |
|
self.getCoefficient("D_reduced"), \ |
1873 |
|
self.getCoefficient("X_reduced"), \ |
1874 |
|
self.getCoefficient("Y_reduced"), \ |
1875 |
|
self.getCoefficient("d_reduced"), \ |
1876 |
|
self.getCoefficient("y_reduced"), \ |
1877 |
|
self.getCoefficient("d_contact_reduced"), \ |
1878 |
|
self.getCoefficient("y_contact_reduced")) |
1879 |
|
self.insertConstraint(rhs_only=False) |
1880 |
|
self.trace("New system has been built.") |
1881 |
|
self.validOperator() |
1882 |
|
self.validRightHandSide() |
1883 |
|
elif not self.isRightHandSideValid(): |
1884 |
|
self.resetRightHandSide() |
1885 |
|
righthandside=self.getCurrentRightHandSide() |
1886 |
|
self.getDomain().addPDEToRHS(righthandside, |
1887 |
|
self.getCoefficient("X"), \ |
1888 |
|
self.getCoefficient("Y"),\ |
1889 |
|
self.getCoefficient("y"),\ |
1890 |
|
self.getCoefficient("y_contact")) |
1891 |
|
self.getDomain().addPDEToRHS(righthandside, |
1892 |
|
self.getCoefficient("X_reduced"), \ |
1893 |
|
self.getCoefficient("Y_reduced"),\ |
1894 |
|
self.getCoefficient("y_reduced"),\ |
1895 |
|
self.getCoefficient("y_contact_reduced")) |
1896 |
|
self.insertConstraint(rhs_only=True) |
1897 |
|
self.trace("New right hand side has been built.") |
1898 |
|
self.validRightHandSide() |
1899 |
|
elif not self.isOperatorValid(): |
1900 |
|
self.resetOperator() |
1901 |
|
operator=self.getCurrentOperator() |
1902 |
|
self.getDomain().addPDEToSystem(operator,escript.Data(), \ |
1903 |
self.getCoefficient("A"), \ |
self.getCoefficient("A"), \ |
1904 |
self.getCoefficient("B"), \ |
self.getCoefficient("B"), \ |
1905 |
self.getCoefficient("C"), \ |
self.getCoefficient("C"), \ |
1906 |
self.getCoefficient("D"), \ |
self.getCoefficient("D"), \ |
1907 |
self.getCoefficient("X"), \ |
escript.Data(), \ |
1908 |
self.getCoefficient("Y"), \ |
escript.Data(), \ |
1909 |
self.getCoefficient("d"), \ |
self.getCoefficient("d"), \ |
1910 |
self.getCoefficient("y"), \ |
escript.Data(),\ |
1911 |
self.getCoefficient("d_contact"), \ |
self.getCoefficient("d_contact"), \ |
1912 |
self.getCoefficient("y_contact")) |
escript.Data()) |
1913 |
self.__operator_isValid=True |
self.getDomain().addPDEToSystem(operator,escript.Data(), \ |
1914 |
self.__righthandside_isValid=True |
self.getCoefficient("A_reduced"), \ |
1915 |
self.__applyConstraint() |
self.getCoefficient("B_reduced"), \ |
1916 |
self.__copyConstraint(self.__righthandside) |
self.getCoefficient("C_reduced"), \ |
1917 |
elif not self.__operator_isValid: |
self.getCoefficient("D_reduced"), \ |
1918 |
self.getOperator() |
escript.Data(), \ |
1919 |
elif not self.__righthandside_isValid: |
escript.Data(), \ |
1920 |
self.getRightHandSide() |
self.getCoefficient("d_reduced"), \ |
1921 |
return (self.__operator,self.__righthandside) |
escript.Data(),\ |
1922 |
|
self.getCoefficient("d_contact_reduced"), \ |
1923 |
def solve(self,**options): |
escript.Data()) |
1924 |
""" |
self.insertConstraint(rhs_only=False) |
1925 |
@brief solve the PDE |
self.trace("New operator has been built.") |
1926 |
|
self.validOperator() |
1927 |
@param options |
return (self.getCurrentOperator(), self.getCurrentRightHandSide()) |
1928 |
""" |
|
1929 |
mat,f=self.getSystem() |
def insertConstraint(self, rhs_only=False): |
1930 |
if self.isUsingLumping(): |
""" |
1931 |
out=f/mat |
Applies the constraints defined by q and r to the PDE. |
1932 |
self.__copyConstraint(out) |
|
1933 |
else: |
@param rhs_only: if True only the right hand side is altered by the |
1934 |
options[util.TOLERANCE_KEY]=self.getTolerance() |
constraint |
1935 |
options[util.METHOD_KEY]=self.getSolverMethod() |
@type rhs_only: C{bool} |
1936 |
options[util.SYMMETRY_KEY]=self.isSymmetric() |
""" |
1937 |
if self.debug() : print "PDE Debug: solver options: ",options |
q=self.getCoefficient("q") |
1938 |
out=mat.solve(f,options) |
r=self.getCoefficient("r") |
1939 |
return out |
righthandside=self.getCurrentRightHandSide() |
1940 |
|
operator=self.getCurrentOperator() |
1941 |
|
|
1942 |
def getSolution(self,**options): |
if not q.isEmpty(): |
1943 |
""" |
if r.isEmpty(): |
1944 |
@brief returns the solution of the PDE |
r_s=self.createSolution() |
1945 |
|
else: |
1946 |
|
r_s=r |
1947 |
|
if not rhs_only and not operator.isEmpty(): |
1948 |
|
if self.isUsingLumping(): |
1949 |
|
operator.copyWithMask(escript.Data(1.,q.getShape(),q.getFunctionSpace()),q) |
1950 |
|
else: |
1951 |
|
row_q=escript.Data(q,self.getFunctionSpaceForEquation()) |
1952 |
|
col_q=escript.Data(q,self.getFunctionSpaceForSolution()) |
1953 |
|
u=self.createSolution() |
1954 |
|
u.copyWithMask(r_s,col_q) |
1955 |
|
righthandside-=operator*u |
1956 |
|
operator.nullifyRowsAndCols(row_q,col_q,1.) |
1957 |
|
righthandside.copyWithMask(r_s,q) |
1958 |
|
|
1959 |
@param options |
def setValue(self,**coefficients): |
1960 |
""" |
""" |
1961 |
if not self.__solution_isValid: |
Sets new values to coefficients. |
1962 |
if self.debug() : print "PDE Debug: PDE is resolved." |
|
1963 |
self.__solution=self.solve(**options) |
@param coefficients: new values assigned to coefficients |
1964 |
self.__solution_isValid=True |
@keyword A: value for coefficient C{A} |
1965 |
return self.__solution |
@type A: any type that can be cast to a L{Data<escript.Data>} object on |
1966 |
#============ some serivice functions ===================================================== |
L{Function<escript.Function>} |
1967 |
def getDomain(self): |
@keyword A_reduced: value for coefficient C{A_reduced} |
1968 |
|
@type A_reduced: any type that can be cast to a L{Data<escript.Data>} |
1969 |
|
object on L{ReducedFunction<escript.ReducedFunction>} |
1970 |
|
@keyword B: value for coefficient C{B} |
1971 |
|
@type B: any type that can be cast to a L{Data<escript.Data>} object on |
1972 |
|
L{Function<escript.Function>} |
1973 |
|
@keyword B_reduced: value for coefficient C{B_reduced} |
1974 |
|
@type B_reduced: any type that can be cast to a L{Data<escript.Data>} |
1975 |
|
object on L{ReducedFunction<escript.ReducedFunction>} |
1976 |
|
@keyword C: value for coefficient C{C} |
1977 |
|
@type C: any type that can be cast to a L{Data<escript.Data>} object on |
1978 |
|
L{Function<escript.Function>} |
1979 |
|
@keyword C_reduced: value for coefficient C{C_reduced} |
1980 |
|
@type C_reduced: any type that can be cast to a L{Data<escript.Data>} |
1981 |
|
object on L{ReducedFunction<escript.ReducedFunction>} |
1982 |
|
@keyword D: value for coefficient C{D} |
1983 |
|
@type D: any type that can be cast to a L{Data<escript.Data>} object on |
1984 |
|
L{Function<escript.Function>} |
1985 |
|
@keyword D_reduced: value for coefficient C{D_reduced} |
1986 |
|
@type D_reduced: any type that can be cast to a L{Data<escript.Data>} |
1987 |
|
object on L{ReducedFunction<escript.ReducedFunction>} |
1988 |
|
@keyword X: value for coefficient C{X} |
1989 |
|
@type X: any type that can be cast to a L{Data<escript.Data>} object on |
1990 |
|
L{Function<escript.Function>} |
1991 |
|
@keyword X_reduced: value for coefficient C{X_reduced} |
1992 |
|
@type X_reduced: any type that can be cast to a L{Data<escript.Data>} |
1993 |
|
object on L{ReducedFunction<escript.ReducedFunction>} |
1994 |
|
@keyword Y: value for coefficient C{Y} |
1995 |
|
@type Y: any type that can be cast to a L{Data<escript.Data>} object on |
1996 |
|
L{Function<escript.Function>} |
1997 |
|
@keyword Y_reduced: value for coefficient C{Y_reduced} |
1998 |
|
@type Y_reduced: any type that can be cast to a L{Data<escript.Data>} |
1999 |
|
object on L{ReducedFunction<escript.Function>} |
2000 |
|
@keyword d: value for coefficient C{d} |
2001 |
|
@type d: any type that can be cast to a L{Data<escript.Data>} object on |
2002 |
|
L{FunctionOnBoundary<escript.FunctionOnBoundary>} |
2003 |
|
@keyword d_reduced: value for coefficient C{d_reduced} |
2004 |
|
@type d_reduced: any type that can be cast to a L{Data<escript.Data>} |
2005 |
|
object on L{ReducedFunctionOnBoundary<escript.ReducedFunctionOnBoundary>} |
2006 |
|
@keyword y: value for coefficient C{y} |
2007 |
|
@type y: any type that can be cast to a L{Data<escript.Data>} object on |
2008 |
|
L{FunctionOnBoundary<escript.FunctionOnBoundary>} |
2009 |
|
@keyword d_contact: value for coefficient C{d_contact} |
2010 |
|
@type d_contact: any type that can be cast to a L{Data<escript.Data>} |
2011 |
|
object on L{FunctionOnContactOne<escript.FunctionOnContactOne>} |
2012 |
|
or L{FunctionOnContactZero<escript.FunctionOnContactZero>} |
2013 |
|
@keyword d_contact_reduced: value for coefficient C{d_contact_reduced} |
2014 |
|
@type d_contact_reduced: any type that can be cast to a L{Data<escript.Data>} |
2015 |
|
object on L{ReducedFunctionOnContactOne<escript.ReducedFunctionOnContactOne>} |
2016 |
|
or L{ReducedFunctionOnContactZero<escript.ReducedFunctionOnContactZero>} |
2017 |
|
@keyword y_contact: value for coefficient C{y_contact} |
2018 |
|
@type y_contact: any type that can be cast to a L{Data<escript.Data>} |
2019 |
|
object on L{FunctionOnContactOne<escript.FunctionOnContactOne>} |
2020 |
|
or L{FunctionOnContactZero<escript.FunctionOnContactZero>} |
2021 |
|
@keyword y_contact_reduced: value for coefficient C{y_contact_reduced} |
2022 |
|
@type y_contact_reduced: any type that can be cast to a L{Data<escript.Data>} |
2023 |
|
object on L{ReducedFunctionOnContactOne<escript.FunctionOnContactOne>} |
2024 |
|
or L{ReducedFunctionOnContactZero<escript.FunctionOnContactZero>} |
2025 |
|
@keyword r: values prescribed to the solution at the locations of |
2026 |
|
constraints |
2027 |
|
@type r: any type that can be cast to a L{Data<escript.Data>} object on |
2028 |
|
L{Solution<escript.Solution>} or L{ReducedSolution<escript.ReducedSolution>} |
2029 |
|
depending on whether reduced order is used for the solution |
2030 |
|
@keyword q: mask for location of constraints |
2031 |
|
@type q: any type that can be cast to a L{Data<escript.Data>} object on |
2032 |
|
L{Solution<escript.Solution>} or L{ReducedSolution<escript.ReducedSolution>} |
2033 |
|
depending on whether reduced order is used for the |
2034 |
|
representation of the equation |
2035 |
|
@raise IllegalCoefficient: if an unknown coefficient keyword is used |
2036 |
|
""" |
2037 |
|
super(LinearPDE,self).setValue(**coefficients) |
2038 |
|
# check if the systrem is inhomogeneous: |
2039 |
|
if len(coefficients)>0 and not self.isUsingLumping(): |
2040 |
|
q=self.getCoefficient("q") |
2041 |
|
r=self.getCoefficient("r") |
2042 |
|
if not q.isEmpty() and not r.isEmpty(): |
2043 |
|
if util.Lsup(q*r)>0.: |
2044 |
|
self.trace("Inhomogeneous constraint detected.") |
2045 |
|
self.invalidateSystem() |
2046 |
|
|
2047 |
|
|
2048 |
|
def getResidual(self,u=None): |
2049 |
|
""" |
2050 |
|
Returns the residual of u or the current solution if u is not present. |
2051 |
|
|
2052 |
|
@param u: argument in the residual calculation. It must be representable |
2053 |
|
in L{self.getFunctionSpaceForSolution()}. If u is not present |
2054 |
|
or equals C{None} the current solution is used. |
2055 |
|
@type u: L{Data<escript.Data>} or None |
2056 |
|
@return: residual of u |
2057 |
|
@rtype: L{Data<escript.Data>} |
2058 |
""" |
""" |
2059 |
@brief returns the domain of the PDE |
if u==None: |
2060 |
|
return self.getOperator()*self.getSolution()-self.getRightHandSide() |
2061 |
|
else: |
2062 |
|
return self.getOperator()*escript.Data(u,self.getFunctionSpaceForSolution())-self.getRightHandSide() |
2063 |
|
|
2064 |
|
def getFlux(self,u=None): |
2065 |
""" |
""" |
2066 |
return self.__domain |
Returns the flux M{J} for a given M{u}. |
2067 |
|
|
2068 |
def getDim(self): |
M{J[i,j]=(A[i,j,k,l]+A_reduced[A[i,j,k,l]]*grad(u[k])[l]+(B[i,j,k]+B_reduced[i,j,k])u[k]-X[i,j]-X_reduced[i,j]} |
2069 |
|
|
2070 |
|
or |
2071 |
|
|
2072 |
|
M{J[j]=(A[i,j]+A_reduced[i,j])*grad(u)[l]+(B[j]+B_reduced[j])u-X[j]-X_reduced[j]} |
2073 |
|
|
2074 |
|
@param u: argument in the flux. If u is not present or equals L{None} the |
2075 |
|
current solution is used. |
2076 |
|
@type u: L{Data<escript.Data>} or None |
2077 |
|
@return: flux |
2078 |
|
@rtype: L{Data<escript.Data>} |
2079 |
|
""" |
2080 |
|
if u==None: u=self.getSolution() |
2081 |
|
return util.tensormult(self.getCoefficient("A"),util.grad(u,Funtion(self.getDomain))) \ |
2082 |
|
+util.matrixmult(self.getCoefficient("B"),u) \ |
2083 |
|
-util.self.getCoefficient("X") \ |
2084 |
|
+util.tensormult(self.getCoefficient("A_reduced"),util.grad(u,ReducedFuntion(self.getDomain))) \ |
2085 |
|
+util.matrixmult(self.getCoefficient("B_reduced"),u) \ |
2086 |
|
-util.self.getCoefficient("X_reduced") |
2087 |
|
|
2088 |
|
|
2089 |
|
class Poisson(LinearPDE): |
2090 |
|
""" |
2091 |
|
Class to define a Poisson equation problem. This is generally a |
2092 |
|
L{LinearPDE} of the form |
2093 |
|
|
2094 |
|
M{-grad(grad(u)[j])[j] = f} |
2095 |
|
|
2096 |
|
with natural boundary conditions |
2097 |
|
|
2098 |
|
M{n[j]*grad(u)[j] = 0 } |
2099 |
|
|
2100 |
|
and constraints: |
2101 |
|
|
2102 |
|
M{u=0} where M{q>0} |
2103 |
|
|
2104 |
|
""" |
2105 |
|
|
2106 |
|
def __init__(self,domain,debug=False): |
2107 |
""" |
""" |
2108 |
@brief returns the spatial dimension of the PDE |
Initializes a new Poisson equation. |
2109 |
|
|
2110 |
|
@param domain: domain of the PDE |
2111 |
|
@type domain: L{Domain<escript.Domain>} |
2112 |
|
@param debug: if True debug information is printed |
2113 |
|
|
2114 |
""" |
""" |
2115 |
return self.getDomain().getDim() |
super(Poisson, self).__init__(domain,1,1,debug) |
2116 |
|
self.introduceCoefficients( |
2117 |
|
f=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE), |
2118 |
|
f_reduced=PDECoef(PDECoef.INTERIOR_REDUCED,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE)) |
2119 |
|
self.setSymmetryOn() |
2120 |
|
|
2121 |
def getNumEquations(self): |
def setValue(self,**coefficients): |
2122 |
""" |
""" |
2123 |
@brief returns the number of equations |
Sets new values to coefficients. |
2124 |
|
|
2125 |
|
@param coefficients: new values assigned to coefficients |
2126 |
|
@keyword f: value for right hand side M{f} |
2127 |
|
@type f: any type that can be cast to a L{Scalar<escript.Scalar>} object |
2128 |
|
on L{Function<escript.Function>} |
2129 |
|
@keyword q: mask for location of constraints |
2130 |
|
@type q: any type that can be cast to a rank zero L{Data<escript.Data>} |
2131 |
|
object on L{Solution<escript.Solution>} or |
2132 |
|
L{ReducedSolution<escript.ReducedSolution>} depending on whether |
2133 |
|
reduced order is used for the representation of the equation |
2134 |
|
@raise IllegalCoefficient: if an unknown coefficient keyword is used |
2135 |
""" |
""" |
2136 |
if self.__numEquations>0: |
super(Poisson, self).setValue(**coefficients) |
2137 |
return self.__numEquations |
|
2138 |
|
|
2139 |
|
def getCoefficient(self,name): |
2140 |
|
""" |
2141 |
|
Returns the value of the coefficient C{name} of the general PDE. |
2142 |
|
|
2143 |
|
@param name: name of the coefficient requested |
2144 |
|
@type name: C{string} |
2145 |
|
@return: the value of the coefficient C{name} |
2146 |
|
@rtype: L{Data<escript.Data>} |
2147 |
|
@raise IllegalCoefficient: invalid coefficient name |
2148 |
|
@note: This method is called by the assembling routine to map the Poisson |
2149 |
|
equation onto the general PDE. |
2150 |
|
""" |
2151 |
|
if name == "A" : |
2152 |
|
return escript.Data(util.kronecker(self.getDim()),escript.Function(self.getDomain())) |
2153 |
|
elif name == "Y" : |
2154 |
|
return self.getCoefficient("f") |
2155 |
|
elif name == "Y_reduced" : |
2156 |
|
return self.getCoefficient("f_reduced") |
2157 |
else: |
else: |
2158 |
raise ValueError,"Number of equations is undefined. Please specify argument numEquations." |
return super(Poisson, self).getCoefficient(name) |
2159 |
|
|
2160 |
def getNumSolutions(self): |
class Helmholtz(LinearPDE): |
2161 |
|
""" |
2162 |
|
Class to define a Helmholtz equation problem. This is generally a |
2163 |
|
L{LinearPDE} of the form |
2164 |
|
|
2165 |
|
M{S{omega}*u - grad(k*grad(u)[j])[j] = f} |
2166 |
|
|
2167 |
|
with natural boundary conditions |
2168 |
|
|
2169 |
|
M{k*n[j]*grad(u)[j] = g- S{alpha}u } |
2170 |
|
|
2171 |
|
and constraints: |
2172 |
|
|
2173 |
|
M{u=r} where M{q>0} |
2174 |
|
|
2175 |
|
""" |
2176 |
|
|
2177 |
|
def __init__(self,domain,debug=False): |
2178 |
""" |
""" |
2179 |
@brief returns the number of unknowns |
Initializes a new Helmholtz equation. |
2180 |
|
|
2181 |
|
@param domain: domain of the PDE |
2182 |
|
@type domain: L{Domain<escript.Domain>} |
2183 |
|
@param debug: if True debug information is printed |
2184 |
|
|
2185 |
|
""" |
2186 |
|
super(Helmholtz, self).__init__(domain,1,1,debug) |
2187 |
|
self.introduceCoefficients( |
2188 |
|
omega=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,),PDECoef.OPERATOR), |
2189 |
|
k=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,),PDECoef.OPERATOR), |
2190 |
|
f=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE), |
2191 |
|
f_reduced=PDECoef(PDECoef.INTERIOR_REDUCED,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE), |
2192 |
|
alpha=PDECoef(PDECoef.BOUNDARY,(PDECoef.BY_EQUATION,),PDECoef.OPERATOR), |
2193 |
|
g=PDECoef(PDECoef.BOUNDARY,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE), |
2194 |
|
g_reduced=PDECoef(PDECoef.BOUNDARY_REDUCED,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE)) |
2195 |
|
self.setSymmetryOn() |
2196 |
|
|
2197 |
|
def setValue(self,**coefficients): |
2198 |
|
""" |
2199 |
|
Sets new values to coefficients. |
2200 |
|
|
2201 |
|
@param coefficients: new values assigned to coefficients |
2202 |
|
@keyword omega: value for coefficient M{S{omega}} |
2203 |
|
@type omega: any type that can be cast to a L{Scalar<escript.Scalar>} |
2204 |
|
object on L{Function<escript.Function>} |
2205 |
|
@keyword k: value for coefficient M{k} |
2206 |
|
@type k: any type that can be cast to a L{Scalar<escript.Scalar>} object |
2207 |
|
on L{Function<escript.Function>} |
2208 |
|
@keyword f: value for right hand side M{f} |
2209 |
|
@type f: any type that can be cast to a L{Scalar<escript.Scalar>} object |
2210 |
|
on L{Function<escript.Function>} |
2211 |
|
@keyword alpha: value for right hand side M{S{alpha}} |
2212 |
|
@type alpha: any type that can be cast to a L{Scalar<escript.Scalar>} |
2213 |
|
object on L{FunctionOnBoundary<escript.FunctionOnBoundary>} |
2214 |
|
@keyword g: value for right hand side M{g} |
2215 |
|
@type g: any type that can be cast to a L{Scalar<escript.Scalar>} object |
2216 |
|
on L{FunctionOnBoundary<escript.FunctionOnBoundary>} |
2217 |
|
@keyword r: prescribed values M{r} for the solution in constraints |
2218 |
|
@type r: any type that can be cast to a L{Scalar<escript.Scalar>} object |
2219 |
|
on L{Solution<escript.Solution>} or |
2220 |
|
L{ReducedSolution<escript.ReducedSolution>} depending on whether |
2221 |
|
reduced order is used for the representation of the equation |
2222 |
|
@keyword q: mask for the location of constraints |
2223 |
|
@type q: any type that can be cast to a L{Scalar<escript.Scalar>} object |
2224 |
|
on L{Solution<escript.Solution>} or |
2225 |
|
L{ReducedSolution<escript.ReducedSolution>} depending on whether |
2226 |
|
reduced order is used for the representation of the equation |
2227 |
|
@raise IllegalCoefficient: if an unknown coefficient keyword is used |
2228 |
|
""" |
2229 |
|
super(Helmholtz, self).setValue(**coefficients) |
2230 |
|
|
2231 |
|
def getCoefficient(self,name): |
2232 |
""" |
""" |
2233 |
if self.__numSolutions>0: |
Returns the value of the coefficient C{name} of the general PDE. |
2234 |
return self.__numSolutions |
|
2235 |
|
@param name: name of the coefficient requested |
2236 |
|
@type name: C{string} |
2237 |
|
@return: the value of the coefficient C{name} |
2238 |
|
@rtype: L{Data<escript.Data>} |
2239 |
|
@raise IllegalCoefficient: invalid name |
2240 |
|
""" |
2241 |
|
if name == "A" : |
2242 |
|
return escript.Data(numarray.identity(self.getDim()),escript.Function(self.getDomain()))*self.getCoefficient("k") |
2243 |
|
elif name == "D" : |
2244 |
|
return self.getCoefficient("omega") |
2245 |
|
elif name == "Y" : |
2246 |
|
return self.getCoefficient("f") |
2247 |
|
elif name == "d" : |
2248 |
|
return self.getCoefficient("alpha") |
2249 |
|
elif name == "y" : |
2250 |
|
return self.getCoefficient("g") |
2251 |
|
elif name == "Y_reduced" : |
2252 |
|
return self.getCoefficient("f_reduced") |
2253 |
|
elif name == "y_reduced" : |
2254 |
|
return self.getCoefficient("g_reduced") |
2255 |
else: |
else: |
2256 |
raise ValueError,"Number of solution is undefined. Please specify argument numSolutions." |
return super(Helmholtz, self).getCoefficient(name) |
2257 |
|
|
2258 |
|
class LameEquation(LinearPDE): |
2259 |
|
""" |
2260 |
|
Class to define a Lame equation problem. This problem is defined as: |
2261 |
|
|
2262 |
|
M{-grad(S{mu}*(grad(u[i])[j]+grad(u[j])[i]))[j] - grad(S{lambda}*grad(u[k])[k])[j] = F_i -grad(S{sigma}[ij])[j] } |
2263 |
|
|
2264 |
|
with natural boundary conditions: |
2265 |
|
|
2266 |
|
M{n[j]*(S{mu}*(grad(u[i])[j]+grad(u[j])[i]) + S{lambda}*grad(u[k])[k]) = f_i +n[j]*S{sigma}[ij] } |
2267 |
|
|
2268 |
|
and constraints: |
2269 |
|
|
2270 |
|
M{u[i]=r[i]} where M{q[i]>0} |
2271 |
|
|
2272 |
|
""" |
2273 |
|
|
2274 |
def checkSymmetry(self): |
def __init__(self,domain,debug=False): |
2275 |
""" |
""" |
2276 |
@brief returns if the Operator is symmetric. This is a very expensive operation!!! The symmetry flag is not altered. |
Initializes a new Lame equation. |
2277 |
|
|
2278 |
|
@param domain: domain of the PDE |
2279 |
|
@type domain: L{Domain<escript.Domain>} |
2280 |
|
@param debug: if True debug information is printed |
2281 |
|
|
2282 |
|
""" |
2283 |
|
super(LameEquation, self).__init__(domain,\ |
2284 |
|
domain.getDim(),domain.getDim(),debug) |
2285 |
|
self.introduceCoefficients(lame_lambda=PDECoef(PDECoef.INTERIOR,(),PDECoef.OPERATOR), |
2286 |
|
lame_mu=PDECoef(PDECoef.INTERIOR,(),PDECoef.OPERATOR), |
2287 |
|
F=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE), |
2288 |
|
sigma=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,PDECoef.BY_DIM),PDECoef.RIGHTHANDSIDE), |
2289 |
|
f=PDECoef(PDECoef.BOUNDARY,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE)) |
2290 |
|
self.setSymmetryOn() |
2291 |
|
|
2292 |
|
def setValues(self,**coefficients): |
2293 |
|
""" |
2294 |
|
Sets new values to coefficients. |
2295 |
|
|
2296 |
|
@param coefficients: new values assigned to coefficients |
2297 |
|
@keyword lame_mu: value for coefficient M{S{mu}} |
2298 |
|
@type lame_mu: any type that can be cast to a L{Scalar<escript.Scalar>} |
2299 |
|
object on L{Function<escript.Function>} |
2300 |
|
@keyword lame_lambda: value for coefficient M{S{lambda}} |
2301 |
|
@type lame_lambda: any type that can be cast to a L{Scalar<escript.Scalar>} |
2302 |
|
object on L{Function<escript.Function>} |
2303 |
|
@keyword F: value for internal force M{F} |
2304 |
|
@type F: any type that can be cast to a L{Vector<escript.Vector>} object |
2305 |
|
on L{Function<escript.Function>} |
2306 |
|
@keyword sigma: value for initial stress M{S{sigma}} |
2307 |
|
@type sigma: any type that can be cast to a L{Tensor<escript.Tensor>} |
2308 |
|
object on L{Function<escript.Function>} |
2309 |
|
@keyword f: value for external force M{f} |
2310 |
|
@type f: any type that can be cast to a L{Vector<escript.Vector>} object |
2311 |
|
on L{FunctionOnBoundary<escript.FunctionOnBoundary>} |
2312 |
|
@keyword r: prescribed values M{r} for the solution in constraints |
2313 |
|
@type r: any type that can be cast to a L{Vector<escript.Vector>} object |
2314 |
|
on L{Solution<escript.Solution>} or |
2315 |
|
L{ReducedSolution<escript.ReducedSolution>} depending on whether |
2316 |
|
reduced order is used for the representation of the equation |
2317 |
|
@keyword q: mask for the location of constraints |
2318 |
|
@type q: any type that can be cast to a L{Vector<escript.Vector>} object |
2319 |
|
on L{Solution<escript.Solution>} or |
2320 |
|
L{ReducedSolution<escript.ReducedSolution>} depending on whether |
2321 |
|
reduced order is used for the representation of the equation |
2322 |
|
@raise IllegalCoefficient: if an unknown coefficient keyword is used |
2323 |
|
""" |
2324 |
|
super(LameEquation, self).setValues(**coefficients) |
2325 |
|
|
2326 |
|
def getCoefficient(self,name): |
2327 |
|
""" |
2328 |
|
Returns the value of the coefficient C{name} of the general PDE. |
2329 |
|
|
2330 |
|
@param name: name of the coefficient requested |
2331 |
|
@type name: C{string} |
2332 |
|
@return: the value of the coefficient C{name} |
2333 |
|
@rtype: L{Data<escript.Data>} |
2334 |
|
@raise IllegalCoefficient: invalid coefficient name |
2335 |
|
""" |
2336 |
|
if name == "A" : |
2337 |
|
out =self.createCoefficient("A") |
2338 |
|
for i in range(self.getDim()): |
2339 |
|
for j in range(self.getDim()): |
2340 |
|
out[i,i,j,j] += self.getCoefficient("lame_lambda") |
2341 |
|
out[i,j,j,i] += self.getCoefficient("lame_mu") |
2342 |
|
out[i,j,i,j] += self.getCoefficient("lame_mu") |
2343 |
|
return out |
2344 |
|
elif name == "X" : |
2345 |
|
return self.getCoefficient("sigma") |
2346 |
|
elif name == "Y" : |
2347 |
|
return self.getCoefficient("F") |
2348 |
|
elif name == "y" : |
2349 |
|
return self.getCoefficient("f") |
2350 |
|
else: |
2351 |
|
return super(LameEquation, self).getCoefficient(name) |
2352 |
|
|
2353 |
|
def LinearSinglePDE(domain,debug=False): |
2354 |
|
""" |
2355 |
|
Defines a single linear PDE. |
2356 |
|
|
2357 |
|
@param domain: domain of the PDE |
2358 |
|
@type domain: L{Domain<escript.Domain>} |
2359 |
|
@param debug: if True debug information is printed |
2360 |
|
@rtype: L{LinearPDE} |
2361 |
|
""" |
2362 |
|
return LinearPDE(domain,numEquations=1,numSolutions=1,debug=debug) |
2363 |
|
|
2364 |
|
def LinearPDESystem(domain,debug=False): |
2365 |
|
""" |
2366 |
|
Defines a system of linear PDEs. |
2367 |
|
|
2368 |
|
@param domain: domain of the PDEs |
2369 |
|
@type domain: L{Domain<escript.Domain>} |
2370 |
|
@param debug: if True debug information is printed |
2371 |
|
@rtype: L{LinearPDE} |
2372 |
|
""" |
2373 |
|
return LinearPDE(domain,numEquations=domain.getDim(),numSolutions=domain.getDim(),debug=debug) |
2374 |
|
|
2375 |
|
|
2376 |
|
class TransportPDE(LinearProblem): |
2377 |
|
""" |
2378 |
|
This class is used to define a transport problem given by a general linear, |
2379 |
|
time dependent, second order PDE for an unknown, non-negative, |
2380 |
|
time-dependent function M{u} on a given domain defined through a |
2381 |
|
L{Domain<escript.Domain>} object. |
2382 |
|
|
2383 |
|
For a single equation with a solution with a single component the transport |
2384 |
|
problem is defined in the following form: |
2385 |
|
|
2386 |
|
M{(M+M_reduced)*u_t=-(grad(A[j,l]+A_reduced[j,l])*grad(u)[l]+(B[j]+B_reduced[j])u)[j]+(C[l]+C_reduced[l])*grad(u)[l]+(D+D_reduced)-grad(X+X_reduced)[j,j]+(Y+Y_reduced)} |
2387 |
|
|
2388 |
|
where M{u_t} denotes the time derivative of M{u} and M{grad(F)} denotes the |
2389 |
|
spatial derivative of M{F}. Einstein's summation convention, ie. summation |
2390 |
|
over indexes appearing twice in a term of a sum performed, is used. |
2391 |
|
The coefficients M{M}, M{A}, M{B}, M{C}, M{D}, M{X} and M{Y} have to be |
2392 |
|
specified through L{Data<escript.Data>} objects in L{Function<escript.Function>} |
2393 |
|
and the coefficients M{M_reduced}, M{A_reduced}, M{B_reduced}, M{C_reduced}, |
2394 |
|
M{D_reduced}, M{X_reduced} and M{Y_reduced} have to be specified through |
2395 |
|
L{Data<escript.Data>} objects in L{ReducedFunction<escript.ReducedFunction>}. |
2396 |
|
It is also allowed to use objects that can be converted into such |
2397 |
|
L{Data<escript.Data>} objects. M{M} and M{M_reduced} are scalar, M{A} and |
2398 |
|
M{A_reduced} are rank two, M{B}, M{C}, M{X}, M{B_reduced}, M{C_reduced} and |
2399 |
|
M{X_reduced} are rank one and M{D}, M{D_reduced}, M{Y} and M{Y_reduced} are |
2400 |
|
scalar. |
2401 |
|
|
2402 |
|
The following natural boundary conditions are considered: |
2403 |
|
|
2404 |
|
M{n[j]*((A[i,j]+A_reduced[i,j])*grad(u)[l]+(B+B_reduced)[j]*u+X[j]+X_reduced[j])+(d+d_reduced)*u+y+y_reduced=(m+m_reduced)*u_t} |
2405 |
|
|
2406 |
|
where M{n} is the outer normal field. Notice that the coefficients M{A}, |
2407 |
|
M{A_reduced}, M{B}, M{B_reduced}, M{X} and M{X_reduced} are defined in the |
2408 |
|
transport problem. The coefficients M{m}, M{d} and M{y} are each a scalar in |
2409 |
|
L{FunctionOnBoundary<escript.FunctionOnBoundary>} and the coefficients |
2410 |
|
M{m_reduced}, M{d_reduced} and M{y_reduced} are each a scalar in |
2411 |
|
L{ReducedFunctionOnBoundary<escript.ReducedFunctionOnBoundary>}. |
2412 |
|
|
2413 |
|
Constraints for the solution prescribing the value of the solution at |
2414 |
|
certain locations in the domain have the form |
2415 |
|
|
2416 |
|
M{u_t=r} where M{q>0} |
2417 |
|
|
2418 |
|
M{r} and M{q} are each scalar where M{q} is the characteristic function |
2419 |
|
defining where the constraint is applied. The constraints override any other |
2420 |
|
condition set by the transport problem or the boundary condition. |
2421 |
|
|
2422 |
|
The transport problem is symmetrical if |
2423 |
|
|
2424 |
|
M{A[i,j]=A[j,i]} and M{B[j]=C[j]} and M{A_reduced[i,j]=A_reduced[j,i]} and |
2425 |
|
M{B_reduced[j]=C_reduced[j]} |
2426 |
|
|
2427 |
|
For a system and a solution with several components the transport problem |
2428 |
|
has the form |
2429 |
|
|
2430 |
|
M{(M[i,k]+M_reduced[i,k])*u[k]_t=-grad((A[i,j,k,l]+A_reduced[i,j,k,l])*grad(u[k])[l]+(B[i,j,k]+B_reduced[i,j,k])*u[k])[j]+(C[i,k,l]+C_reduced[i,k,l])*grad(u[k])[l]+(D[i,k]+D_reduced[i,k]*u[k]-grad(X[i,j]+X_reduced[i,j])[j]+Y[i]+Y_reduced[i] } |
2431 |
|
|
2432 |
|
M{A} and M{A_reduced} are of rank four, M{B}, M{B_reduced}, M{C} and |
2433 |
|
M{C_reduced} are each of rank three, M{M}, M{M_reduced}, M{D}, M{D_reduced}, |
2434 |
|
M{X_reduced} and M{X} are each of rank two and M{Y} and M{Y_reduced} are of |
2435 |
|
rank one. The natural boundary conditions take the form: |
2436 |
|
|
2437 |
|
M{n[j]*((A[i,j,k,l]+A_reduced[i,j,k,l])*grad(u[k])[l]+(B[i,j,k]+B_reduced[i,j,k])*u[k]+X[i,j]+X_reduced[i,j])+(d[i,k]+d_reduced[i,k])*u[k]+y[i]+y_reduced[i]= (m[i,k]+m_reduced[i,k])*u[k]_t} |
2438 |
|
|
2439 |
|
The coefficient M{d} and M{m} are of rank two and M{y} is of rank one with |
2440 |
|
all in L{FunctionOnBoundary<escript.FunctionOnBoundary>}. The coefficients |
2441 |
|
M{d_reduced} and M{m_reduced} are of rank two and M{y_reduced} is of rank |
2442 |
|
one all in L{ReducedFunctionOnBoundary<escript.ReducedFunctionOnBoundary>}. |
2443 |
|
|
2444 |
|
Constraints take the form |
2445 |
|
|
2446 |
|
M{u[i]_t=r[i]} where M{q[i]>0} |
2447 |
|
|
2448 |
|
M{r} and M{q} are each rank one. Notice that at some locations not |
2449 |
|
necessarily all components must have a constraint. |
2450 |
|
|
2451 |
|
The transport problem is symmetrical if |
2452 |
|
|
2453 |
|
- M{M[i,k]=M[i,k]} |
2454 |
|
- M{M_reduced[i,k]=M_reduced[i,k]} |
2455 |
|
- M{A[i,j,k,l]=A[k,l,i,j]} |
2456 |
|
- M{A_reduced[i,j,k,l]=A_reduced[k,l,i,j]} |
2457 |
|
- M{B[i,j,k]=C[k,i,j]} |
2458 |
|
- M{B_reduced[i,j,k]=C_reduced[k,i,j]} |
2459 |
|
- M{D[i,k]=D[i,k]} |
2460 |
|
- M{D_reduced[i,k]=D_reduced[i,k]} |
2461 |
|
- M{m[i,k]=m[k,i]} |
2462 |
|
- M{m_reduced[i,k]=m_reduced[k,i]} |
2463 |
|
- M{d[i,k]=d[k,i]} |
2464 |
|
- M{d_reduced[i,k]=d_reduced[k,i]} |
2465 |
|
|
2466 |
|
L{TransportPDE} also supports solution discontinuities over a contact region |
2467 |
|
in the domain. To specify the conditions across the discontinuity we are |
2468 |
|
using the generalised flux M{J} which, in the case of a system of PDEs and |
2469 |
|
several components of the solution, is defined as |
2470 |
|
|
2471 |
|
M{J[i,j]=(A[i,j,k,l]+A_reduced[[i,j,k,l])*grad(u[k])[l]+(B[i,j,k]+B_reduced[i,j,k])*u[k]+X[i,j]+X_reduced[i,j]} |
2472 |
|
|
2473 |
|
For the case of single solution component and single PDE M{J} is defined as |
2474 |
|
|
2475 |
|
M{J[j]=(A[i,j]+A_reduced[i,j])*grad(u)[j]+(B[i]+B_reduced[i])*u+X[i]+X_reduced[i]} |
2476 |
|
|
2477 |
|
In the context of discontinuities M{n} denotes the normal on the |
2478 |
|
discontinuity pointing from side 0 towards side 1 calculated from |
2479 |
|
L{getNormal<escript.FunctionSpace.getNormal>} of L{FunctionOnContactZero<escript.FunctionOnContactZero>}. |
2480 |
|
For a system of transport problems the contact condition takes the form |
2481 |
|
|
2482 |
|
M{n[j]*J0[i,j]=n[j]*J1[i,j]=(y_contact[i]+y_contact_reduced[i])- (d_contact[i,k]+d_contact_reduced[i,k])*jump(u)[k]} |
2483 |
|
|
2484 |
|
where M{J0} and M{J1} are the fluxes on side 0 and side 1 of the |
2485 |
|
discontinuity, respectively. M{jump(u)}, which is the difference of the |
2486 |
|
solution at side 1 and at side 0, denotes the jump of M{u} across |
2487 |
|
discontinuity along the normal calculated by L{jump<util.jump>}. |
2488 |
|
The coefficient M{d_contact} is of rank two and M{y_contact} is of rank one |
2489 |
|
both in L{FunctionOnContactZero<escript.FunctionOnContactZero>} or L{FunctionOnContactOne<escript.FunctionOnContactOne>}. |
2490 |
|
The coefficient M{d_contact_reduced} is of rank two and M{y_contact_reduced} |
2491 |
|
is of rank one both in L{ReducedFunctionOnContactZero<escript.ReducedFunctionOnContactZero>} or L{ReducedFunctionOnContactOne<escript.ReducedFunctionOnContactOne>}. |
2492 |
|
In case of a single PDE and a single component solution the contact |
2493 |
|
condition takes the form |
2494 |
|
|
2495 |
|
M{n[j]*J0_{j}=n[j]*J1_{j}=(y_contact+y_contact_reduced)-(d_contact+y_contact_reduced)*jump(u)} |
2496 |
|
|
2497 |
|
In this case the coefficient M{d_contact} and M{y_contact} are each scalar |
2498 |
|
both in L{FunctionOnContactZero<escript.FunctionOnContactZero>} or |
2499 |
|
L{FunctionOnContactOne<escript.FunctionOnContactOne>} and the coefficient |
2500 |
|
M{d_contact_reduced} and M{y_contact_reduced} are each scalar both in |
2501 |
|
L{ReducedFunctionOnContactZero<escript.ReducedFunctionOnContactZero>} or |
2502 |
|
L{ReducedFunctionOnContactOne<escript.ReducedFunctionOnContactOne>}. |
2503 |
|
|
2504 |
|
Typical usage:: |
2505 |
|
|
2506 |
|
p = TransportPDE(dom) |
2507 |
|
p.setValue(M=1., C=[-1.,0.]) |
2508 |
|
p.setInitialSolution(u=exp(-length(dom.getX()-[0.1,0.1])**2) |
2509 |
|
t = 0 |
2510 |
|
dt = 0.1 |
2511 |
|
while (t < 1.): |
2512 |
|
u = p.solve(dt) |
2513 |
|
|
2514 |
|
""" |
2515 |
|
def __init__(self,domain,numEquations=None,numSolutions=None, theta=0.5,debug=False): |
2516 |
|
""" |
2517 |
|
Initializes a transport problem. |
2518 |
|
|
2519 |
|
@param domain: domain of the PDE |
2520 |
|
@type domain: L{Domain<escript.Domain>} |
2521 |
|
@param numEquations: number of equations. If C{None} the number of |
2522 |
|
equations is extracted from the coefficients. |
2523 |
|
@param numSolutions: number of solution components. If C{None} the number |
2524 |
|
of solution components is extracted from the |
2525 |
|
coefficients. |
2526 |
|
@param debug: if True debug information is printed |
2527 |
|
@param theta: time stepping control: |
2528 |
|
- 1.0: backward Euler, |
2529 |
|
- 0.0: forward Euler, |
2530 |
|
- 0.5: Crank-Nicholson scheme, U{see here<http://en.wikipedia.org/wiki/Crank-Nicolson_method>} |
2531 |
|
|
2532 |
|
""" |
2533 |
|
if theta<0 or theta>1: |
2534 |
|
raise ValueError,"theta (=%s) needs to be between 0 and 1 (inclusive)."%theta |
2535 |
|
super(TransportPDE, self).__init__(domain,numEquations,numSolutions,debug) |
2536 |
|
|
2537 |
|
self.__theta=theta |
2538 |
|
self.setConstraintWeightingFactor() |
2539 |
|
# |
2540 |
|
# the coefficients of the transport problem |
2541 |
|
# |
2542 |
|
self.introduceCoefficients( |
2543 |
|
M=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION),PDECoef.OPERATOR), |
2544 |
|
A=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,PDECoef.BY_DIM,PDECoef.BY_SOLUTION,PDECoef.BY_DIM),PDECoef.OPERATOR), |
2545 |
|
B=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,PDECoef.BY_DIM,PDECoef.BY_SOLUTION),PDECoef.OPERATOR), |
2546 |
|
C=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION,PDECoef.BY_DIM),PDECoef.OPERATOR), |
2547 |
|
D=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION),PDECoef.OPERATOR), |
2548 |
|
X=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,PDECoef.BY_DIM),PDECoef.RIGHTHANDSIDE), |
2549 |
|
Y=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE), |
2550 |
|
m=PDECoef(PDECoef.BOUNDARY,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION),PDECoef.OPERATOR), |
2551 |
|
d=PDECoef(PDECoef.BOUNDARY,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION),PDECoef.OPERATOR), |
2552 |
|
y=PDECoef(PDECoef.BOUNDARY,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE), |
2553 |
|
d_contact=PDECoef(PDECoef.CONTACT,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION),PDECoef.OPERATOR), |
2554 |
|
y_contact=PDECoef(PDECoef.CONTACT,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE), |
2555 |
|
M_reduced=PDECoef(PDECoef.INTERIOR_REDUCED,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION),PDECoef.OPERATOR), |
2556 |
|
A_reduced=PDECoef(PDECoef.INTERIOR_REDUCED,(PDECoef.BY_EQUATION,PDECoef.BY_DIM,PDECoef.BY_SOLUTION,PDECoef.BY_DIM),PDECoef.OPERATOR), |
2557 |
|
B_reduced=PDECoef(PDECoef.INTERIOR_REDUCED,(PDECoef.BY_EQUATION,PDECoef.BY_DIM,PDECoef.BY_SOLUTION),PDECoef.OPERATOR), |
2558 |
|
C_reduced=PDECoef(PDECoef.INTERIOR_REDUCED,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION,PDECoef.BY_DIM),PDECoef.OPERATOR), |
2559 |
|
D_reduced=PDECoef(PDECoef.INTERIOR_REDUCED,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION),PDECoef.OPERATOR), |
2560 |
|
X_reduced=PDECoef(PDECoef.INTERIOR_REDUCED,(PDECoef.BY_EQUATION,PDECoef.BY_DIM),PDECoef.RIGHTHANDSIDE), |
2561 |
|
Y_reduced=PDECoef(PDECoef.INTERIOR_REDUCED,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE), |
2562 |
|
m_reduced=PDECoef(PDECoef.BOUNDARY_REDUCED,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION),PDECoef.OPERATOR), |
2563 |
|
d_reduced=PDECoef(PDECoef.BOUNDARY_REDUCED,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION),PDECoef.OPERATOR), |
2564 |
|
y_reduced=PDECoef(PDECoef.BOUNDARY_REDUCED,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE), |
2565 |
|
d_contact_reduced=PDECoef(PDECoef.CONTACT_REDUCED,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION),PDECoef.OPERATOR), |
2566 |
|
y_contact_reduced=PDECoef(PDECoef.CONTACT_REDUCED,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE), |
2567 |
|
r=PDECoef(PDECoef.SOLUTION,(PDECoef.BY_SOLUTION,),PDECoef.RIGHTHANDSIDE), |
2568 |
|
q=PDECoef(PDECoef.SOLUTION,(PDECoef.BY_SOLUTION,),PDECoef.BOTH) ) |
2569 |
|
|
2570 |
|
def __str__(self): |
2571 |
|
""" |
2572 |
|
Returns the string representation of the transport problem. |
2573 |
|
|
2574 |
|
@return: a simple representation of the transport problem |
2575 |
|
@rtype: C{str} |
2576 |
|
""" |
2577 |
|
return "<TransportPDE %d>"%id(self) |
2578 |
|
|
2579 |
|
def getTheta(self): |
2580 |
|
""" |
2581 |
|
Returns the time stepping control parameter. |
2582 |
|
@rtype: float |
2583 |
|
""" |
2584 |
|
return self.__theta |
2585 |
|
|
2586 |
|
|
2587 |
|
def checkSymmetry(self,verbose=True): |
2588 |
|
""" |
2589 |
|
Tests the transport problem for symmetry. |
2590 |
|
|
2591 |
|
@param verbose: if set to True or not present a report on coefficients |
2592 |
|
which break the symmetry is printed. |
2593 |
|
@type verbose: C{bool} |
2594 |
|
@return: True if the PDE is symmetric |
2595 |
|
@rtype: C{bool} |
2596 |
|
@note: This is a very expensive operation. It should be used for |
2597 |
|
degugging only! The symmetry flag is not altered. |
2598 |
|
""" |
2599 |
|
out=True |
2600 |
|
out=out and self.checkSymmetricTensor("M", verbose) |
2601 |
|
out=out and self.checkSymmetricTensor("M_reduced", verbose) |
2602 |
|
out=out and self.checkSymmetricTensor("A", verbose) |
2603 |
|
out=out and self.checkSymmetricTensor("A_reduced", verbose) |
2604 |
|
out=out and self.checkReciprocalSymmetry("B","C", verbose) |
2605 |
|
out=out and self.checkReciprocalSymmetry("B_reduced","C_reduced", verbose) |
2606 |
|
out=out and self.checkSymmetricTensor("D", verbose) |
2607 |
|
out=out and self.checkSymmetricTensor("D_reduced", verbose) |
2608 |
|
out=out and self.checkSymmetricTensor("m", verbose) |
2609 |
|
out=out and self.checkSymmetricTensor("m_reduced", verbose) |
2610 |
|
out=out and self.checkSymmetricTensor("d", verbose) |
2611 |
|
out=out and self.checkSymmetricTensor("d_reduced", verbose) |
2612 |
|
out=out and self.checkSymmetricTensor("d_contact", verbose) |
2613 |
|
out=out and self.checkSymmetricTensor("d_contact_reduced", verbose) |
2614 |
|
return out |
2615 |
|
|
2616 |
|
def setValue(self,**coefficients): |
2617 |
""" |
""" |
2618 |
raise SystemError,"checkSymmetry is not implemented yet" |
Sets new values to coefficients. |
2619 |
|
|
2620 |
return None |
@param coefficients: new values assigned to coefficients |
2621 |
|
@keyword M: value for coefficient C{M} |
2622 |
|
@type M: any type that can be cast to a L{Data<escript.Data>} object on |
2623 |
|
L{Function<escript.Function>} |
2624 |
|
@keyword M_reduced: value for coefficient C{M_reduced} |
2625 |
|
@type M_reduced: any type that can be cast to a L{Data<escript.Data>} |
2626 |
|
object on L{Function<escript.ReducedFunction>} |
2627 |
|
@keyword A: value for coefficient C{A} |
2628 |
|
@type A: any type that can be cast to a L{Data<escript.Data>} object on |
2629 |
|
L{Function<escript.Function>} |
2630 |
|
@keyword A_reduced: value for coefficient C{A_reduced} |
2631 |
|
@type A_reduced: any type that can be cast to a L{Data<escript.Data>} |
2632 |
|
object on L{ReducedFunction<escript.ReducedFunction>} |
2633 |
|
@keyword B: value for coefficient C{B} |
2634 |
|
@type B: any type that can be cast to a L{Data<escript.Data>} object on |
2635 |
|
L{Function<escript.Function>} |
2636 |
|
@keyword B_reduced: value for coefficient C{B_reduced} |
2637 |
|
@type B_reduced: any type that can be cast to a L{Data<escript.Data>} |
2638 |
|
object on L{ReducedFunction<escript.ReducedFunction>} |
2639 |
|
@keyword C: value for coefficient C{C} |
2640 |
|
@type C: any type that can be cast to a L{Data<escript.Data>} object on |
2641 |
|
L{Function<escript.Function>} |
2642 |
|
@keyword C_reduced: value for coefficient C{C_reduced} |
2643 |
|
@type C_reduced: any type that can be cast to a L{Data<escript.Data>} |
2644 |
|
object on L{ReducedFunction<escript.ReducedFunction>} |
2645 |
|
@keyword D: value for coefficient C{D} |
2646 |
|
@type D: any type that can be cast to a L{Data<escript.Data>} object on |
2647 |
|
L{Function<escript.Function>} |
2648 |
|
@keyword D_reduced: value for coefficient C{D_reduced} |
2649 |
|
@type D_reduced: any type that can be cast to a L{Data<escript.Data>} |
2650 |
|
object on L{ReducedFunction<escript.ReducedFunction>} |
2651 |
|
@keyword X: value for coefficient C{X} |
2652 |
|
@type X: any type that can be cast to a L{Data<escript.Data>} object on |
2653 |
|
L{Function<escript.Function>} |
2654 |
|
@keyword X_reduced: value for coefficient C{X_reduced} |
2655 |
|
@type X_reduced: any type that can be cast to a L{Data<escript.Data>} |
2656 |
|
object on L{ReducedFunction<escript.ReducedFunction>} |
2657 |
|
@keyword Y: value for coefficient C{Y} |
2658 |
|
@type Y: any type that can be cast to a L{Data<escript.Data>} object on |
2659 |
|
L{Function<escript.Function>} |
2660 |
|
@keyword Y_reduced: value for coefficient C{Y_reduced} |
2661 |
|
@type Y_reduced: any type that can be cast to a L{Data<escript.Data>} |
2662 |
|
object on L{ReducedFunction<escript.Function>} |
2663 |
|
@keyword m: value for coefficient C{m} |
2664 |
|
@type m: any type that can be cast to a L{Data<escript.Data>} object on |
2665 |
|
L{FunctionOnBoundary<escript.FunctionOnBoundary>} |
2666 |
|
@keyword m_reduced: value for coefficient C{m_reduced} |
2667 |
|
@type m_reduced: any type that can be cast to a L{Data<escript.Data>} |
2668 |
|
object on L{FunctionOnBoundary<escript.ReducedFunctionOnBoundary>} |
2669 |
|
@keyword d: value for coefficient C{d} |
2670 |
|
@type d: any type that can be cast to a L{Data<escript.Data>} object on |
2671 |
|
L{FunctionOnBoundary<escript.FunctionOnBoundary>} |
2672 |
|
@keyword d_reduced: value for coefficient C{d_reduced} |
2673 |
|
@type d_reduced: any type that can be cast to a L{Data<escript.Data>} |
2674 |
|
object on L{ReducedFunctionOnBoundary<escript.ReducedFunctionOnBoundary>} |
2675 |
|
@keyword y: value for coefficient C{y} |
2676 |
|
@type y: any type that can be cast to a L{Data<escript.Data>} object on |
2677 |
|
L{FunctionOnBoundary<escript.FunctionOnBoundary>} |
2678 |
|
@keyword d_contact: value for coefficient C{d_contact} |
2679 |
|
@type d_contact: any type that can be cast to a L{Data<escript.Data>} |
2680 |
|
object on L{FunctionOnContactOne<escript.FunctionOnContactOne>} or L{FunctionOnContactZero<escript.FunctionOnContactZero>} |
2681 |
|
@keyword d_contact_reduced: value for coefficient C{d_contact_reduced} |
2682 |
|
@type d_contact_reduced: any type that can be cast to a L{Data<escript.Data>} object on L{ReducedFunctionOnContactOne<escript.ReducedFunctionOnContactOne>} or L{ReducedFunctionOnContactZero<escript.ReducedFunctionOnContactZero>} |
2683 |
|
@keyword y_contact: value for coefficient C{y_contact} |
2684 |
|
@type y_contact: any type that can be cast to a L{Data<escript.Data>} |
2685 |
|
object on L{FunctionOnContactOne<escript.FunctionOnContactOne>} or L{FunctionOnContactZero<escript.FunctionOnContactZero>} |
2686 |
|
@keyword y_contact_reduced: value for coefficient C{y_contact_reduced} |
2687 |
|
@type y_contact_reduced: any type that can be cast to a L{Data<escript.Data>} object on L{ReducedFunctionOnContactOne<escript.FunctionOnContactOne>} or L{ReducedFunctionOnContactZero<escript.FunctionOnContactZero>} |
2688 |
|
@keyword r: values prescribed to the solution at the locations of constraints |
2689 |
|
@type r: any type that can be cast to a L{Data<escript.Data>} object on |
2690 |
|
L{Solution<escript.Solution>} or L{ReducedSolution<escript.ReducedSolution>} |
2691 |
|
depending on whether reduced order is used for the solution |
2692 |
|
@keyword q: mask for the location of constraints |
2693 |
|
@type q: any type that can be cast to a L{Data<escript.Data>} object on |
2694 |
|
L{Solution<escript.Solution>} or |
2695 |
|
L{ReducedSolution<escript.ReducedSolution>} depending on whether |
2696 |
|
reduced order is used for the representation of the equation |
2697 |
|
@raise IllegalCoefficient: if an unknown coefficient keyword is used |
2698 |
|
""" |
2699 |
|
super(TransportPDE,self).setValue(**coefficients) |
2700 |
|
|
2701 |
|
def createOperator(self): |
2702 |
|
""" |
2703 |
|
Returns an instance of a new transport operator. |
2704 |
|
""" |
2705 |
|
self.trace("New Transport problem is allocated.") |
2706 |
|
return self.getDomain().newTransportProblem( \ |
2707 |
|
self.getTheta(), |
2708 |
|
self.getNumEquations(), \ |
2709 |
|
self.getFunctionSpaceForSolution(), \ |
2710 |
|
self.getSystemType()) |
2711 |
|
|
2712 |
|
def setInitialSolution(self,u): |
2713 |
|
""" |
2714 |
|
Sets the initial solution. |
2715 |
|
|
2716 |
|
@param u: new initial solution |
2717 |
|
@type u: any object that can be interpolated to a L{Data<escript.Data>} |
2718 |
|
object on L{Solution<escript.Solution>} or L{ReducedSolution<escript.ReducedSolution>} |
2719 |
|
@note: C{u} must be non-negative |
2720 |
|
""" |
2721 |
|
u2=util.interpolate(u,self.getFunctionSpaceForSolution()) |
2722 |
|
if self.getNumSolutions() == 1: |
2723 |
|
if u2.getShape()!=(): |
2724 |
|
raise ValueError,"Illegal shape %s of initial solution."%(u2.getShape(),) |
2725 |
|
else: |
2726 |
|
if u2.getShape()!=(self.getNumSolutions(),): |
2727 |
|
raise ValueError,"Illegal shape %s of initial solution."%(u2.getShape(),) |
2728 |
|
self.getOperator().setInitialValue(u2) |
2729 |
|
|
2730 |
def getFlux(self,u): |
def getRequiredSystemType(self): |
2731 |
""" |
""" |
2732 |
@brief returns the flux J_ij for a given u |
Returns the system type which needs to be used by the current set up. |
2733 |
|
|
2734 |
|
@return: a code to indicate the type of transport problem scheme used |
2735 |
|
@rtype: C{float} |
2736 |
|
""" |
2737 |
|
return self.getDomain().getTransportTypeId(self.getSolverMethod()[0],self.getSolverMethod()[1],self.getSolverPackage(),self.isSymmetric()) |
2738 |
|
|
2739 |
|
def getUnlimitedTimeStepSize(self): |
2740 |
|
""" |
2741 |
|
Returns the value returned by the C{getSafeTimeStepSize} method to |
2742 |
|
indicate no limit on the safe time step size. |
2743 |
|
|
2744 |
J_ij=A_{ijkl}u_{k,l}+B_{ijk}u_k-X_{ij} |
@return: the value used to indicate that no limit is set to the time |
2745 |
|
step size |
2746 |
|
@rtype: C{float} |
2747 |
|
@note: Typically the biggest positive float is returned |
2748 |
|
""" |
2749 |
|
return self.getOperator().getUnlimitedTimeStepSize() |
2750 |
|
|
2751 |
@param u argument of the operator |
def getSafeTimeStepSize(self): |
2752 |
|
""" |
2753 |
|
Returns a safe time step size to do the next time step. |
2754 |
|
|
2755 |
|
@return: safe time step size |
2756 |
|
@rtype: C{float} |
2757 |
|
@note: If not C{getSafeTimeStepSize()} < C{getUnlimitedTimeStepSize()} |
2758 |
|
any time step size can be used. |
2759 |
""" |
""" |
2760 |
raise SystemError,"getFlux is not implemented yet" |
return self.getOperator().getSafeTimeStepSize() |
|
return None |
|
2761 |
|
|
2762 |
def applyOperator(self,u): |
def setConstraintWeightingFactor(self,value=1./util.sqrt(util.EPSILON)): |
2763 |
""" |
""" |
2764 |
@brief applies the operator of the PDE to a given solution u in weak from |
Sets the weighting factor used to insert the constraints into the problem |
2765 |
|
|
2766 |
@param u argument of the operator |
@param value: value for the weighting factor |
2767 |
|
@type value: large positive C{float} |
2768 |
|
""" |
2769 |
|
if not value>0: |
2770 |
|
raise ValueError,"weighting factor needs to be positive." |
2771 |
|
self.__constraint_factor=value |
2772 |
|
self.trace("Weighting factor for constraints is set to %e."%value) |
2773 |
|
|
2774 |
|
def getConstraintWeightingFactor(self): |
2775 |
|
""" |
2776 |
|
returns the weighting factor used to insert the constraints into the problem |
2777 |
|
@return: value for the weighting factor |
2778 |
|
@rtype: C{float} |
2779 |
""" |
""" |
2780 |
return self.getOperator()*escript.Data(u,self.getFunctionSpaceForSolution()) |
return self.__constraint_factor |
2781 |
|
#==================================================================== |
2782 |
def getResidual(self,u): |
def getSolution(self,dt,**options): |
2783 |
""" |
""" |
2784 |
@brief return the residual of u in the weak from |
Returns the solution of the problem. |
2785 |
|
|
2786 |
|
@return: the solution |
2787 |
|
@rtype: L{Data<escript.Data>} |
2788 |
|
|
|
@param u |
|
2789 |
""" |
""" |
2790 |
return self.applyOperator(u)-self.getRightHandSide() |
if dt<=0: |
2791 |
|
raise ValueError,"step size needs to be positive." |
2792 |
|
options[self.TOLERANCE_KEY]=self.getTolerance() |
2793 |
|
self.setSolution(self.getOperator().solve(self.getRightHandSide(),dt,options)) |
2794 |
|
self.validSolution() |
2795 |
|
return self.getCurrentSolution() |
2796 |
|
|
2797 |
|
def getSystem(self): |
2798 |
|
""" |
2799 |
|
Returns the operator and right hand side of the PDE. |
2800 |
|
|
2801 |
|
@return: the discrete version of the PDE |
2802 |
|
@rtype: C{tuple} of L{Operator,<escript.Operator>} and |
2803 |
|
L{Data<escript.Data>} |
2804 |
|
|
2805 |
|
""" |
2806 |
|
if not self.isOperatorValid() or not self.isRightHandSideValid(): |
2807 |
|
self.resetRightHandSide() |
2808 |
|
righthandside=self.getCurrentRightHandSide() |
2809 |
|
self.resetOperator() |
2810 |
|
operator=self.getCurrentOperator() |
2811 |
|
self.getDomain().addPDEToTransportProblem( |
2812 |
|
operator, |
2813 |
|
righthandside, |
2814 |
|
self.getCoefficient("M"), |
2815 |
|
self.getCoefficient("A"), |
2816 |
|
self.getCoefficient("B"), |
2817 |
|
self.getCoefficient("C"), |
2818 |
|
self.getCoefficient("D"), |
2819 |
|
self.getCoefficient("X"), |
2820 |
|
self.getCoefficient("Y"), |
2821 |
|
self.getCoefficient("d"), |
2822 |
|
self.getCoefficient("y"), |
2823 |
|
self.getCoefficient("d_contact"), |
2824 |
|
self.getCoefficient("y_contact")) |
2825 |
|
self.getDomain().addPDEToTransportProblem( |
2826 |
|
operator, |
2827 |
|
righthandside, |
2828 |
|
self.getCoefficient("M_reduced"), |
2829 |
|
self.getCoefficient("A_reduced"), |
2830 |
|
self.getCoefficient("B_reduced"), |
2831 |
|
self.getCoefficient("C_reduced"), |
2832 |
|
self.getCoefficient("D_reduced"), |
2833 |
|
self.getCoefficient("X_reduced"), |
2834 |
|
self.getCoefficient("Y_reduced"), |
2835 |
|
self.getCoefficient("d_reduced"), |
2836 |
|
self.getCoefficient("y_reduced"), |
2837 |
|
self.getCoefficient("d_contact_reduced"), |
2838 |
|
self.getCoefficient("y_contact_reduced")) |
2839 |
|
operator.insertConstraint(righthandside,self.getCoefficient("q"),self.getCoefficient("r"),self.getConstraintWeightingFactor()) |
2840 |
|
self.trace("New system has been built.") |
2841 |
|
self.validOperator() |
2842 |
|
self.validRightHandSide() |
2843 |
|
return (self.getCurrentOperator(), self.getCurrentRightHandSide()) |
2844 |
|
|
|
class Poisson(LinearPDE): |
|
|
""" |
|
|
@brief Class to define a Poisson equstion problem: |
|
|
|
|
|
class to define a linear PDE of the form |
|
|
|
|
|
-u_{,jj} = f |
|
|
|
|
|
with boundary conditons: |
|
|
|
|
|
n_j*u_{,j} = 0 |
|
|
|
|
|
and constraints: |
|
|
|
|
|
u=0 where q>0 |
|
|
|
|
|
""" |
|
|
|
|
|
def __init__(self,domain=None,f=escript.Data(),q=escript.Data()): |
|
|
LinearPDE.__init__(self,domain=identifyDomain(domain,{"f":f, "q":q})) |
|
|
self._setValue(A=numarray.identity(self.getDomain().getDim())) |
|
|
self.setSymmetryOn() |
|
|
self.setValue(f,q) |
|
|
|
|
|
def setValue(self,f=escript.Data(),q=escript.Data()): |
|
|
self._setValue(Y=f,q=q) |
|
|
|
|
|
|
|
|
# $Log$ |
|
|
# Revision 1.3 2004/12/17 07:43:10 jgs |
|
|
# *** empty log message *** |
|
|
# |
|
|
# Revision 1.1.2.3 2004/12/16 00:12:34 gross |
|
|
# __init__ of LinearPDE does not accept any coefficients anymore |
|
|
# |
|
|
# Revision 1.1.2.2 2004/12/14 03:55:01 jgs |
|
|
# *** empty log message *** |
|
|
# |
|
|
# Revision 1.1.2.1 2004/12/12 22:53:47 gross |
|
|
# linearPDE has been renamed LinearPDE |
|
|
# |
|
|
# Revision 1.1.1.1.2.7 2004/12/07 10:13:08 gross |
|
|
# GMRES added |
|
|
# |
|
|
# Revision 1.1.1.1.2.6 2004/12/07 03:19:50 gross |
|
|
# options for GMRES and PRES20 added |
|
|
# |
|
|
# Revision 1.1.1.1.2.5 2004/12/01 06:25:15 gross |
|
|
# some small changes |
|
|
# |
|
|
# Revision 1.1.1.1.2.4 2004/11/24 01:50:21 gross |
|
|
# Finley solves 4M unknowns now |
|
|
# |
|
|
# Revision 1.1.1.1.2.3 2004/11/15 06:05:26 gross |
|
|
# poisson solver added |
|
|
# |
|
|
# Revision 1.1.1.1.2.2 2004/11/12 06:58:15 gross |
|
|
# a lot of changes to get the linearPDE class running: most important change is that there is no matrix format exposed to the user anymore. the format is chosen by the Domain according to the solver and symmetry |
|
|
# |
|
|
# Revision 1.1.1.1.2.1 2004/10/28 22:59:22 gross |
|
|
# finley's RecTest.py is running now: problem in SystemMatrixAdapater fixed |
|
|
# |
|
|
# Revision 1.1.1.1 2004/10/26 06:53:56 jgs |
|
|
# initial import of project esys2 |
|
|
# |
|
|
# Revision 1.3.2.3 2004/10/26 06:43:48 jgs |
|
|
# committing Lutz's and Paul's changes to brach jgs |
|
|
# |
|
|
# Revision 1.3.4.1 2004/10/20 05:32:51 cochrane |
|
|
# Added incomplete Doxygen comments to files, or merely put the docstrings that already exist into Doxygen form. |
|
|
# |
|
|
# Revision 1.3 2004/09/23 00:53:23 jgs |
|
|
# minor fixes |
|
|
# |
|
|
# Revision 1.1 2004/08/28 12:58:06 gross |
|
|
# SimpleSolve is not running yet: problem with == of functionsspace |
|
|
# |
|
|
# |
|