/[escript]/trunk/escript/py_src/linearPDEs.py
ViewVC logotype

Annotation of /trunk/escript/py_src/linearPDEs.py

Parent Directory Parent Directory | Revision Log Revision Log


Revision 431 - (hide annotations)
Fri Jan 13 05:07:10 2006 UTC (13 years, 8 months ago) by gross
File MIME type: text/x-python
File size: 110838 byte(s)
new implementation of ILU0
1 jgs 102 # $Id$
2    
3 jgs 149 #
4     # COPYRIGHT ACcESS 2004 - All Rights Reserved
5 jgs 151 #
6 jgs 149 # This software is the property of ACcESS. No part of this code
7     # may be copied in any form or by any means without the expressed written
8     # consent of ACcESS. Copying, use or modification of this software
9     # by any unauthorised person is illegal unless that
10     # person has a software license agreement with ACcESS.
11     #
12     """
13 jgs 151 The module provides an interface to define and solve linear partial
14     differential equations (PDEs) within L{escript}. L{linearPDEs} does not provide any
15     solver capabilities in itself but hands the PDE over to
16 jgs 149 the PDE solver library defined through the L{Domain<escript.Domain>} of the PDE.
17     The general interface is provided through the L{LinearPDE} class. The
18     L{AdvectivePDE} which is derived from the L{LinearPDE} class
19     provides an interface to PDE dominated by its advective terms. The L{Poisson},
20 gross 345 L{Helmholtz}, L{LameEquation}, L{AdvectionDiffusion}
21 jgs 149 classs which are also derived form the L{LinearPDE} class should be used
22     to define of solve these sepecial PDEs.
23 jgs 102
24 jgs 149 @var __author__: name of author
25     @var __licence__: licence agreement
26     @var __url__: url entry point on documentation
27     @var __version__: version
28     @var __date__: date of the version
29 jgs 102 """
30    
31     import escript
32     import util
33     import numarray
34    
35 jgs 149 __author__="Lutz Gross, l.gross@uq.edu.au"
36     __licence__="contact: esys@access.uq.edu.au"
37     __url__="http://www.iservo.edu.au/esys/escript"
38     __version__="$Revision$"
39     __date__="$Date$"
40 jgs 102
41 jgs 149
42 jgs 148 class IllegalCoefficient(ValueError):
43 jgs 102 """
44 jgs 148 raised if an illegal coefficient of the general ar particular PDE is requested.
45     """
46 jgs 102
47 jgs 148 class IllegalCoefficientValue(ValueError):
48 jgs 102 """
49 jgs 148 raised if an incorrect value for a coefficient is used.
50     """
51 jgs 122
52 jgs 148 class UndefinedPDEError(ValueError):
53     """
54     raised if a PDE is not fully defined yet.
55     """
56 jgs 102
57 jgs 151 class PDECoefficient(object):
58 jgs 102 """
59 jgs 149 A class for describing a PDE coefficient
60    
61     @cvar INTERIOR: indicator that coefficient is defined on the interior of the PDE domain
62     @cvar BOUNDARY: indicator that coefficient is defined on the boundary of the PDE domain
63     @cvar CONTACT: indicator that coefficient is defined on the contact region within the PDE domain
64     @cvar SOLUTION: indicator that coefficient is defined trough a solution of the PDE
65 jgs 150 @cvar REDUCED: indicator that coefficient is defined trough a reduced solution of the PDE
66 jgs 149 @cvar BY_EQUATION: indicator that the dimension of the coefficient shape is defined by the number PDE equations
67     @cvar BY_SOLUTION: indicator that the dimension of the coefficient shape is defined by the number PDE solutions
68     @cvar BY_DIM: indicator that the dimension of the coefficient shape is defined by the spatial dimension
69     @cvar OPERATOR: indicator that the the coefficient alters the operator of the PDE
70     @cvar RIGHTHANDSIDE: indicator that the the coefficient alters the right hand side of the PDE
71     @cvar BOTH: indicator that the the coefficient alters the operator as well as the right hand side of the PDE
72    
73 jgs 102 """
74     INTERIOR=0
75     BOUNDARY=1
76     CONTACT=2
77 jgs 149 SOLUTION=3
78 jgs 150 REDUCED=4
79 jgs 149 BY_EQUATION=5
80     BY_SOLUTION=6
81     BY_DIM=7
82     OPERATOR=10
83     RIGHTHANDSIDE=11
84     BOTH=12
85    
86 jgs 102 def __init__(self,where,pattern,altering):
87     """
88 jgs 122 Initialise a PDE Coefficient type
89 jgs 151
90 jgs 149 @param where: describes where the coefficient lives
91 jgs 150 @type where: one of L{INTERIOR}, L{BOUNDARY}, L{CONTACT}, L{SOLUTION}, L{REDUCED}
92 jgs 151 @param pattern: describes the shape of the coefficient and how the shape is build for a given
93 jgs 149 spatial dimension and numbers of equation and solution in then PDE. For instance,
94     (L{BY_EQUATION},L{BY_SOLUTION},L{BY_DIM}) descrbes a rank 3 coefficient which
95     is instanciated as shape (3,2,2) in case of a three equations and two solution components
96     on a 2-dimensional domain. In the case of single equation and a single solution component
97     the shape compoments marked by L{BY_EQUATION} or L{BY_SOLUTION} are dropped. In this case
98     the example would be read as (2,).
99     @type pattern: C{tuple} of L{BY_EQUATION}, L{BY_SOLUTION}, L{BY_DIM}
100     @param altering: indicates what part of the PDE is altered if the coefficiennt is altered
101     @type altering: one of L{OPERATOR}, L{RIGHTHANDSIDE}, L{BOTH}
102    
103 jgs 102 """
104 jgs 151 super(PDECoefficient, self).__init__()
105 jgs 102 self.what=where
106     self.pattern=pattern
107     self.altering=altering
108 jgs 108 self.resetValue()
109 jgs 102
110 jgs 108 def resetValue(self):
111     """
112 jgs 122 resets coefficient value to default
113 jgs 108 """
114     self.value=escript.Data()
115    
116 jgs 149 def getFunctionSpace(self,domain,reducedEquationOrder=False,reducedSolutionOrder=False):
117 jgs 102 """
118 jgs 149 defines the L{FunctionSpace<escript.FunctionSpace>} of the coefficient
119 jgs 102
120 jgs 149 @param domain: domain on which the PDE uses the coefficient
121     @type domain: L{Domain<escript.Domain>}
122     @param reducedEquationOrder: True to indicate that reduced order is used to represent the equation
123     @type domain: C{bool}
124     @param reducedSolutionOrder: True to indicate that reduced order is used to represent the solution
125     @type domain: C{bool}
126     @return: L{FunctionSpace<escript.FunctionSpace>} of the coefficient
127     @rtype: L{FunctionSpace<escript.FunctionSpace>}
128 jgs 102 """
129 jgs 151 if self.what==self.INTERIOR:
130 jgs 149 return escript.Function(domain)
131 jgs 151 elif self.what==self.BOUNDARY:
132 jgs 149 return escript.FunctionOnBoundary(domain)
133 jgs 151 elif self.what==self.CONTACT:
134 jgs 149 return escript.FunctionOnContactZero(domain)
135 jgs 151 elif self.what==self.SOLUTION:
136 jgs 149 if reducedEquationOrder and reducedSolutionOrder:
137     return escript.ReducedSolution(domain)
138     else:
139     return escript.Solution(domain)
140 jgs 151 elif self.what==self.REDUCED:
141     return escript.ReducedSolution(domain)
142 jgs 102
143 jgs 108 def getValue(self):
144     """
145 jgs 149 returns the value of the coefficient
146    
147     @return: value of the coefficient
148     @rtype: L{Data<escript.Data>}
149 jgs 108 """
150     return self.value
151 jgs 148
152 jgs 149 def setValue(self,domain,numEquations=1,numSolutions=1,reducedEquationOrder=False,reducedSolutionOrder=False,newValue=None):
153 jgs 108 """
154 jgs 149 set the value of the coefficient to a new value
155    
156     @param domain: domain on which the PDE uses the coefficient
157     @type domain: L{Domain<escript.Domain>}
158     @param numEquations: number of equations of the PDE
159     @type numEquations: C{int}
160     @param numSolutions: number of components of the PDE solution
161     @type numSolutions: C{int}
162     @param reducedEquationOrder: True to indicate that reduced order is used to represent the equation
163     @type domain: C{bool}
164     @param reducedSolutionOrder: True to indicate that reduced order is used to represent the solution
165     @type domain: C{bool}
166     @param newValue: number of components of the PDE solution
167     @type newValue: any object that can be converted into a L{Data<escript.Data>} object with the appropriate shape and L{FunctionSpace<escript.FunctionSpace>}
168     @raise IllegalCoefficientValue: if the shape of the assigned value does not match the shape of the coefficient
169 jgs 108 """
170 jgs 148 if newValue==None:
171     newValue=escript.Data()
172     elif isinstance(newValue,escript.Data):
173     if not newValue.isEmpty():
174 jgs 149 try:
175     newValue=escript.Data(newValue,self.getFunctionSpace(domain,reducedEquationOrder,reducedSolutionOrder))
176     except:
177     raise IllegalCoefficientValue,"Unable to interpolate coefficient to function space %s"%self.getFunctionSpace(domain)
178 jgs 148 else:
179 jgs 149 newValue=escript.Data(newValue,self.getFunctionSpace(domain,reducedEquationOrder,reducedSolutionOrder))
180 jgs 148 if not newValue.isEmpty():
181     if not self.getShape(domain,numEquations,numSolutions)==newValue.getShape():
182 jgs 149 raise IllegalCoefficientValue,"Expected shape of coefficient is %s but actual shape is %s."%(self.getShape(domain,numEquations,numSolutions),newValue.getShape())
183 jgs 108 self.value=newValue
184 jgs 148
185 jgs 102 def isAlteringOperator(self):
186     """
187 jgs 149 checks if the coefficient alters the operator of the PDE
188    
189     @return: True if the operator of the PDE is changed when the coefficient is changed
190     @rtype: C{bool}
191 jgs 102 """
192     if self.altering==self.OPERATOR or self.altering==self.BOTH:
193     return not None
194     else:
195     return None
196    
197     def isAlteringRightHandSide(self):
198     """
199 jgs 149 checks if the coefficeint alters the right hand side of the PDE
200    
201     @rtype: C{bool}
202     @return: True if the right hand side of the PDE is changed when the coefficient is changed
203 jgs 102 """
204     if self.altering==self.RIGHTHANDSIDE or self.altering==self.BOTH:
205     return not None
206     else:
207     return None
208    
209 jgs 148 def estimateNumEquationsAndNumSolutions(self,domain,shape=()):
210 jgs 102 """
211 jgs 149 tries to estimate the number of equations and number of solutions if the coefficient has the given shape
212 jgs 102
213 jgs 149 @param domain: domain on which the PDE uses the coefficient
214     @type domain: L{Domain<escript.Domain>}
215     @param shape: suggested shape of the coefficient
216     @type shape: C{tuple} of C{int} values
217     @return: the number of equations and number of solutions of the PDE is the coefficient has shape s.
218     If no appropriate numbers could be identified, C{None} is returned
219     @rtype: C{tuple} of two C{int} values or C{None}
220 jgs 102 """
221 jgs 148 dim=domain.getDim()
222 jgs 102 if len(shape)>0:
223     num=max(shape)+1
224     else:
225     num=1
226     search=[]
227 jgs 149 if self.definesNumEquation() and self.definesNumSolutions():
228     for u in range(num):
229     for e in range(num):
230     search.append((e,u))
231     search.sort(self.__CompTuple2)
232     for item in search:
233 jgs 148 s=self.getShape(domain,item[0],item[1])
234 jgs 102 if len(s)==0 and len(shape)==0:
235     return (1,1)
236     else:
237     if s==shape: return item
238 jgs 149 elif self.definesNumEquation():
239     for e in range(num,0,-1):
240     s=self.getShape(domain,e,0)
241     if len(s)==0 and len(shape)==0:
242     return (1,None)
243     else:
244     if s==shape: return (e,None)
245    
246     elif self.definesNumSolutions():
247     for u in range(num,0,-1):
248     s=self.getShape(domain,0,u)
249     if len(s)==0 and len(shape)==0:
250     return (None,1)
251     else:
252     if s==shape: return (None,u)
253 jgs 102 return None
254 jgs 149 def definesNumSolutions(self):
255     """
256     checks if the coefficient allows to estimate the number of solution components
257 jgs 102
258 jgs 149 @return: True if the coefficient allows an estimate of the number of solution components
259     @rtype: C{bool}
260     """
261     for i in self.pattern:
262     if i==self.BY_SOLUTION: return True
263     return False
264    
265     def definesNumEquation(self):
266     """
267     checks if the coefficient allows to estimate the number of equations
268    
269     @return: True if the coefficient allows an estimate of the number of equations
270     @rtype: C{bool}
271     """
272     for i in self.pattern:
273     if i==self.BY_EQUATION: return True
274     return False
275    
276     def __CompTuple2(self,t1,t2):
277     """
278     Compare two tuples of possible number of equations and number of solutions
279    
280     @param t1: The first tuple
281     @param t2: The second tuple
282    
283     """
284    
285     dif=t1[0]+t1[1]-(t2[0]+t2[1])
286     if dif<0: return 1
287     elif dif>0: return -1
288     else: return 0
289    
290 jgs 148 def getShape(self,domain,numEquations=1,numSolutions=1):
291 jgs 149 """
292     builds the required shape of the coefficient
293 jgs 102
294 jgs 149 @param domain: domain on which the PDE uses the coefficient
295     @type domain: L{Domain<escript.Domain>}
296     @param numEquations: number of equations of the PDE
297     @type numEquations: C{int}
298     @param numSolutions: number of components of the PDE solution
299     @type numSolutions: C{int}
300     @return: shape of the coefficient
301     @rtype: C{tuple} of C{int} values
302     """
303     dim=domain.getDim()
304     s=()
305     for i in self.pattern:
306     if i==self.BY_EQUATION:
307 jgs 148 if numEquations>1: s=s+(numEquations,)
308 jgs 149 elif i==self.BY_SOLUTION:
309 jgs 148 if numSolutions>1: s=s+(numSolutions,)
310 jgs 102 else:
311     s=s+(dim,)
312 jgs 149 return s
313 jgs 102
314 jgs 151 class LinearPDE(object):
315 jgs 102 """
316 jgs 149 This class is used to define a general linear, steady, second order PDE
317     for an unknown function M{u} on a given domain defined through a L{Domain<escript.Domain>} object.
318 jgs 102
319 jgs 149 For a single PDE with a solution with a single component the linear PDE is defined in the following form:
320 jgs 151
321 jgs 149 M{-grad(A[j,l]*grad(u)[l]+B[j]u)[j]+C[l]*grad(u)[l]+D*u =-grad(X)[j,j]+Y}
322 jgs 102
323 jgs 151 where M{grad(F)} denotes the spatial derivative of M{F}. Einstein's summation convention,
324     ie. summation over indexes appearing twice in a term of a sum is performed, is used.
325     The coefficients M{A}, M{B}, M{C}, M{D}, M{X} and M{Y} have to be specified through L{Data<escript.Data>} objects in the
326     L{Function<escript.Function>} on the PDE or objects that can be converted into such L{Data<escript.Data>} objects.
327     M{A} is a rank two, M{B}, M{C} and M{X} are rank one and M{D} and M{Y} are scalar.
328 jgs 102
329 jgs 149 The following natural boundary conditions are considered:
330 jgs 102
331 jgs 149 M{n[j]*(A[i,j]*grad(u)[l]+B[j]*u)+d*u=n[j]*X[j]+y}
332 jgs 102
333 jgs 149 where M{n} is the outer normal field calculated by L{getNormal<escript.FunctionSpace.getNormal>} of L{FunctionOnBoundary<escript.FunctionOnBoundary>}.
334 jgs 151 Notice that the coefficients M{A}, M{B} and M{X} are defined in the PDE. The coefficients M{d} and M{y} are
335     each a scalar in the L{FunctionOnBoundary<escript.FunctionOnBoundary>}.
336 jgs 102
337    
338 jgs 149 Constraints for the solution prescribing the value of the solution at certain locations in the domain. They have the form
339 jgs 102
340 jgs 149 M{u=r} where M{q>0}
341    
342     M{r} and M{q} are each scalar where M{q} is the characteristic function defining where the constraint is applied.
343     The constraints override any other condition set by the PDE or the boundary condition.
344 jgs 151
345 jgs 149 The PDE is symmetrical if
346 jgs 151
347 jgs 149 M{A[i,j]=A[j,i]} and M{B[j]=C[j]}
348    
349     For a system of PDEs and a solution with several components the PDE has the form
350    
351     M{-grad(A[i,j,k,l]*grad(u[k])[l]+B[i,j,k]*u[k])[j]+C[i,k,l]*grad(u[k])[l]+D[i,k]*u[k] =-grad(X[i,j])[j]+Y[i] }
352    
353 jgs 151 M{A} is a ramk four, M{B} and M{C} are each a rank three, M{D} and M{X} are each a rank two and M{Y} is a rank one.
354 jgs 149 The natural boundary conditions take the form:
355    
356     M{n[j]*(A[i,j,k,l]*grad(u[k])[l]+B[i,j,k]*u[k])+d[i,k]*u[k]=n[j]*X[i,j]+y[i]}
357    
358    
359     The coefficient M{d} is a rank two and M{y} is a rank one both in the L{FunctionOnBoundary<escript.FunctionOnBoundary>}. Constraints take the form
360    
361    
362     M{u[i]=r[i]} where M{q[i]>0}
363    
364 jgs 151 M{r} and M{q} are each rank one. Notice that at some locations not necessarily all components must have a constraint.
365 jgs 149
366     The system of PDEs is symmetrical if
367    
368     - M{A[i,j,k,l]=A[k,l,i,j]}
369     - M{B[i,j,k]=C[k,i,j]}
370     - M{D[i,k]=D[i,k]}
371     - M{d[i,k]=d[k,i]}
372    
373 jgs 151 L{LinearPDE} also supports solution discontinuities over a contact region in the domain. To specify the conditions across the
374 jgs 149 discontinuity we are using the generalised flux M{J} which is in the case of a systems of PDEs and several components of the solution
375 jgs 151 defined as
376 jgs 149
377     M{J[i,j]=A[i,j,k,l]*grad(u[k])[l]+B[i,j,k]*u[k]-X[i,j]}
378    
379     For the case of single solution component and single PDE M{J} is defined
380    
381     M{J_{j}=A[i,j]*grad(u)[j]+B[i]*u-X[i]}
382    
383     In the context of discontinuities M{n} denotes the normal on the discontinuity pointing from side 0 towards side 1
384     calculated from L{getNormal<escript.FunctionSpace.getNormal>} of L{FunctionOnContactZero<escript.FunctionOnContactZero>}. For a system of PDEs
385     the contact condition takes the form
386    
387     M{n[j]*J0[i,j]=n[j]*J1[i,j]=y_contact[i]- d_contact[i,k]*jump(u)[k]}
388 jgs 151
389 jgs 149 where M{J0} and M{J1} are the fluxes on side 0 and side 1 of the discontinuity, respectively. M{jump(u)}, which is the difference
390 jgs 151 of the solution at side 1 and at side 0, denotes the jump of M{u} across discontinuity along the normal calcualted by
391 jgs 149 L{jump<util.jump>}.
392     The coefficient M{d_contact} is a rank two and M{y_contact} is a rank one both in the L{FunctionOnContactZero<escript.FunctionOnContactZero>} or L{FunctionOnContactOne<escript.FunctionOnContactOne>}.
393     In case of a single PDE and a single component solution the contact condition takes the form
394    
395     M{n[j]*J0_{j}=n[j]*J1_{j}=y_contact-d_contact*jump(u)}
396    
397     In this case the the coefficient M{d_contact} and M{y_contact} are eaach scalar
398     both in the L{FunctionOnContactZero<escript.FunctionOnContactZero>} or L{FunctionOnContactOne<escript.FunctionOnContactOne>}.
399    
400 jgs 150 @cvar DEFAULT: The default method used to solve the system of linear equations
401 jgs 149 @cvar DIRECT: The direct solver based on LDU factorization
402     @cvar CHOLEVSKY: The direct solver based on LDLt factorization (can only be applied for symmetric PDEs)
403     @cvar PCG: The preconditioned conjugate gradient method (can only be applied for symmetric PDEs)
404 jgs 151 @cvar CR: The conjugate residual method
405 jgs 149 @cvar CGS: The conjugate gardient square method
406     @cvar BICGSTAB: The stabilized BiConjugate Gradient method.
407     @cvar SSOR: The symmetric overrealaxtion method
408     @cvar ILU0: The incomplete LU factorization preconditioner with no fill in
409     @cvar ILUT: The incomplete LU factorization preconditioner with will in
410     @cvar JACOBI: The Jacobi preconditioner
411     @cvar GMRES: The Gram-Schmidt minimum residual method
412     @cvar PRES20: Special GMRES with restart after 20 steps and truncation after 5 residuals
413     @cvar LUMPING: Matrix lumping.
414     @cvar NO_REORDERING: No matrix reordering allowed
415     @cvar MINIMUM_FILL_IN: Reorder matrix to reduce fill-in during factorization
416     @cvar NESTED_DISSECTION: Reorder matrix to improve load balancing during factorization
417 jgs 150 @cvar PASO: PASO solver package
418     @cvar SCSL: SGI SCSL solver library
419     @cvar MKL: Intel's MKL solver library
420 jgs 151 @cvar UMFPACK: the UMFPACK library
421 jgs 150 @cvar ITERATIVE: The default iterative solver
422 gross 430 @cvar AMG: algebraic multi grid
423     @cvar RILU: recursive ILU
424 jgs 149
425 jgs 102 """
426 jgs 150 DEFAULT= 0
427     DIRECT= 1
428     CHOLEVSKY= 2
429     PCG= 3
430     CR= 4
431     CGS= 5
432     BICGSTAB= 6
433     SSOR= 7
434     ILU0= 8
435     ILUT= 9
436     JACOBI= 10
437     GMRES= 11
438     PRES20= 12
439     LUMPING= 13
440     NO_REORDERING= 17
441     MINIMUM_FILL_IN= 18
442     NESTED_DISSECTION= 19
443     SCSL= 14
444     MKL= 15
445     UMFPACK= 16
446     ITERATIVE= 20
447     PASO= 21
448 gross 430 AMG= 22
449     RILU = 23
450 jgs 150
451 jgs 149 __TOL=1.e-13
452 jgs 150 __PACKAGE_KEY="package"
453 jgs 149 __METHOD_KEY="method"
454     __SYMMETRY_KEY="symmetric"
455     __TOLERANCE_KEY="tolerance"
456 gross 387 __PRECONDITIONER_KEY="preconditioner"
457 jgs 102
458 jgs 148
459     def __init__(self,domain,numEquations=None,numSolutions=None,debug=False):
460 jgs 102 """
461 jgs 148 initializes a new linear PDE
462 jgs 102
463 jgs 148 @param domain: domain of the PDE
464 jgs 149 @type domain: L{Domain<escript.Domain>}
465 jgs 148 @param numEquations: number of equations. If numEquations==None the number of equations
466     is exracted from the PDE coefficients.
467     @param numSolutions: number of solution components. If numSolutions==None the number of solution components
468     is exracted from the PDE coefficients.
469     @param debug: if True debug informations are printed.
470    
471 jgs 102 """
472 jgs 151 super(LinearPDE, self).__init__()
473 jgs 148 #
474     # the coefficients of the general PDE:
475     #
476     self.__COEFFICIENTS_OF_GENEARL_PDE={
477 jgs 149 "A" : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.BY_EQUATION,PDECoefficient.BY_DIM,PDECoefficient.BY_SOLUTION,PDECoefficient.BY_DIM),PDECoefficient.OPERATOR),
478     "B" : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.BY_EQUATION,PDECoefficient.BY_DIM,PDECoefficient.BY_SOLUTION),PDECoefficient.OPERATOR),
479     "C" : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.BY_EQUATION,PDECoefficient.BY_SOLUTION,PDECoefficient.BY_DIM),PDECoefficient.OPERATOR),
480     "D" : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.BY_EQUATION,PDECoefficient.BY_SOLUTION),PDECoefficient.OPERATOR),
481     "X" : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.BY_EQUATION,PDECoefficient.BY_DIM),PDECoefficient.RIGHTHANDSIDE),
482     "Y" : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.BY_EQUATION,),PDECoefficient.RIGHTHANDSIDE),
483     "d" : PDECoefficient(PDECoefficient.BOUNDARY,(PDECoefficient.BY_EQUATION,PDECoefficient.BY_SOLUTION),PDECoefficient.OPERATOR),
484     "y" : PDECoefficient(PDECoefficient.BOUNDARY,(PDECoefficient.BY_EQUATION,),PDECoefficient.RIGHTHANDSIDE),
485     "d_contact" : PDECoefficient(PDECoefficient.CONTACT,(PDECoefficient.BY_EQUATION,PDECoefficient.BY_SOLUTION),PDECoefficient.OPERATOR),
486     "y_contact" : PDECoefficient(PDECoefficient.CONTACT,(PDECoefficient.BY_EQUATION,),PDECoefficient.RIGHTHANDSIDE),
487     "r" : PDECoefficient(PDECoefficient.SOLUTION,(PDECoefficient.BY_SOLUTION,),PDECoefficient.RIGHTHANDSIDE),
488     "q" : PDECoefficient(PDECoefficient.SOLUTION,(PDECoefficient.BY_SOLUTION,),PDECoefficient.BOTH)}
489 jgs 102
490 jgs 148 # COEFFICIENTS can be overwritten by subclasses:
491     self.COEFFICIENTS=self.__COEFFICIENTS_OF_GENEARL_PDE
492 jgs 149 self.__altered_coefficients=False
493 jgs 102 # initialize attributes
494 jgs 148 self.__debug=debug
495 jgs 104 self.__domain=domain
496     self.__numEquations=numEquations
497     self.__numSolutions=numSolutions
498 jgs 148 self.__resetSystem()
499 jgs 102
500     # set some default values:
501 jgs 149 self.__reduce_equation_order=False
502     self.__reduce_solution_order=False
503 jgs 102 self.__tolerance=1.e-8
504 jgs 150 self.__solver_method=self.DEFAULT
505     self.__solver_package=self.DEFAULT
506 gross 387 self.__preconditioner=self.DEFAULT
507 jgs 150 self.__matrix_type=self.__domain.getSystemMatrixTypeId(self.DEFAULT,self.DEFAULT,False)
508 jgs 102 self.__sym=False
509    
510 jgs 148 self.resetCoefficients()
511     self.trace("PDE Coeffients are %s"%str(self.COEFFICIENTS.keys()))
512     # =============================================================================
513     # general stuff:
514     # =============================================================================
515     def __str__(self):
516 jgs 149 """
517     returns string representation of the PDE
518    
519     @return: a simple representation of the PDE
520     @rtype: C{str}
521     """
522     return "<LinearPDE %d>"%id(self)
523 jgs 148 # =============================================================================
524     # debug :
525     # =============================================================================
526     def setDebugOn(self):
527 jgs 149 """
528 jgs 148 switches on debugging
529 jgs 108 """
530 jgs 148 self.__debug=not None
531    
532     def setDebugOff(self):
533 jgs 108 """
534 jgs 148 switches off debugging
535     """
536     self.__debug=None
537 jgs 108
538 jgs 148 def trace(self,text):
539     """
540 jgs 149 print the text message if debugging is swiched on.
541     @param text: message
542     @type text: C{string}
543 jgs 102 """
544 jgs 148 if self.__debug: print "%s: %s"%(str(self),text)
545 jgs 102
546 jgs 148 # =============================================================================
547     # some service functions:
548     # =============================================================================
549     def getDomain(self):
550 jgs 102 """
551 jgs 148 returns the domain of the PDE
552 jgs 102
553 jgs 149 @return: the domain of the PDE
554     @rtype: L{Domain<escript.Domain>}
555 jgs 108 """
556 jgs 148 return self.__domain
557 jgs 122
558 jgs 148 def getDim(self):
559 jgs 108 """
560 jgs 148 returns the spatial dimension of the PDE
561 jgs 108
562 jgs 149 @return: the spatial dimension of the PDE domain
563     @rtype: C{int}
564 jgs 148 """
565     return self.getDomain().getDim()
566 jgs 102
567 jgs 148 def getNumEquations(self):
568     """
569     returns the number of equations
570 jgs 102
571 jgs 149 @return: the number of equations
572     @rtype: C{int}
573 jgs 148 @raise UndefinedPDEError: if the number of equations is not be specified yet.
574     """
575     if self.__numEquations==None:
576     raise UndefinedPDEError,"Number of equations is undefined. Please specify argument numEquations."
577     else:
578     return self.__numEquations
579 jgs 147
580 jgs 148 def getNumSolutions(self):
581     """
582     returns the number of unknowns
583 jgs 147
584 jgs 149 @return: the number of unknowns
585     @rtype: C{int}
586 jgs 148 @raise UndefinedPDEError: if the number of unknowns is not be specified yet.
587     """
588     if self.__numSolutions==None:
589     raise UndefinedPDEError,"Number of solution is undefined. Please specify argument numSolutions."
590     else:
591     return self.__numSolutions
592    
593 jgs 149 def reduceEquationOrder(self):
594     """
595     return status for order reduction for equation
596    
597     @return: return True is reduced interpolation order is used for the represenation of the equation
598     @rtype: L{bool}
599     """
600     return self.__reduce_equation_order
601    
602     def reduceSolutionOrder(self):
603     """
604     return status for order reduction for the solution
605    
606     @return: return True is reduced interpolation order is used for the represenation of the solution
607     @rtype: L{bool}
608     """
609     return self.__reduce_solution_order
610 jgs 151
611 jgs 108 def getFunctionSpaceForEquation(self):
612     """
613 jgs 149 returns the L{FunctionSpace<escript.FunctionSpace>} used to discretize the equation
614 jgs 148
615 jgs 149 @return: representation space of equation
616     @rtype: L{FunctionSpace<escript.FunctionSpace>}
617 jgs 108 """
618 jgs 149 if self.reduceEquationOrder():
619     return escript.ReducedSolution(self.getDomain())
620     else:
621     return escript.Solution(self.getDomain())
622 jgs 108
623     def getFunctionSpaceForSolution(self):
624     """
625 jgs 149 returns the L{FunctionSpace<escript.FunctionSpace>} used to represent the solution
626 jgs 148
627 jgs 149 @return: representation space of solution
628     @rtype: L{FunctionSpace<escript.FunctionSpace>}
629 jgs 108 """
630 jgs 149 if self.reduceSolutionOrder():
631     return escript.ReducedSolution(self.getDomain())
632     else:
633     return escript.Solution(self.getDomain())
634 jgs 108
635 jgs 102
636 jgs 148 def getOperator(self):
637 jgs 102 """
638 jgs 148 provides access to the operator of the PDE
639 jgs 102
640 jgs 149 @return: the operator of the PDE
641     @rtype: L{Operator<escript.Operator>}
642 jgs 108 """
643 jgs 148 m=self.getSystem()[0]
644     if self.isUsingLumping():
645     return self.copyConstraint(1./m)
646 jgs 147 else:
647 jgs 148 return m
648 jgs 147
649 jgs 148 def getRightHandSide(self):
650 jgs 147 """
651 jgs 148 provides access to the right hand side of the PDE
652 jgs 149 @return: the right hand side of the PDE
653     @rtype: L{Data<escript.Data>}
654 jgs 147 """
655 jgs 148 r=self.getSystem()[1]
656     if self.isUsingLumping():
657     return self.copyConstraint(r)
658 jgs 147 else:
659 jgs 148 return r
660 jgs 147
661 jgs 148 def applyOperator(self,u=None):
662 jgs 102 """
663 jgs 148 applies the operator of the PDE to a given u or the solution of PDE if u is not present.
664 jgs 102
665 jgs 148 @param u: argument of the operator. It must be representable in C{elf.getFunctionSpaceForSolution()}. If u is not present or equals L{None}
666     the current solution is used.
667 jgs 149 @type u: L{Data<escript.Data>} or None
668     @return: image of u
669     @rtype: L{Data<escript.Data>}
670 jgs 102 """
671 jgs 149 if u==None:
672 jgs 148 return self.getOperator()*self.getSolution()
673 jgs 102 else:
674 jgs 148 self.getOperator()*escript.Data(u,self.getFunctionSpaceForSolution())
675 jgs 102
676 jgs 148 def getResidual(self,u=None):
677 jgs 102 """
678 jgs 148 return the residual of u or the current solution if u is not present.
679 jgs 102
680 jgs 148 @param u: argument in the residual calculation. It must be representable in C{elf.getFunctionSpaceForSolution()}. If u is not present or equals L{None}
681     the current solution is used.
682 jgs 149 @type u: L{Data<escript.Data>} or None
683     @return: residual of u
684     @rtype: L{Data<escript.Data>}
685 jgs 102 """
686 jgs 148 return self.applyOperator(u)-self.getRightHandSide()
687 jgs 102
688 jgs 148 def checkSymmetry(self,verbose=True):
689     """
690     test the PDE for symmetry.
691 jgs 102
692 jgs 149 @param verbose: if equal to True or not present a report on coefficients which are breaking the symmetry is printed.
693     @type verbose: C{bool}
694     @return: True if the PDE is symmetric.
695     @rtype: L{Data<escript.Data>}
696 jgs 148 @note: This is a very expensive operation. It should be used for degugging only! The symmetry flag is not altered.
697     """
698 jgs 149 verbose=verbose or self.__debug
699 jgs 148 out=True
700     if self.getNumSolutions()!=self.getNumEquations():
701     if verbose: print "non-symmetric PDE because of different number of equations and solutions"
702     out=False
703     else:
704     A=self.getCoefficientOfGeneralPDE("A")
705     if not A.isEmpty():
706 jgs 149 tol=util.Lsup(A)*self.__TOL
707 jgs 148 if self.getNumSolutions()>1:
708     for i in range(self.getNumEquations()):
709     for j in range(self.getDim()):
710     for k in range(self.getNumSolutions()):
711     for l in range(self.getDim()):
712     if util.Lsup(A[i,j,k,l]-A[k,l,i,j])>tol:
713     if verbose: print "non-symmetric PDE because A[%d,%d,%d,%d]!=A[%d,%d,%d,%d]"%(i,j,k,l,k,l,i,j)
714     out=False
715     else:
716     for j in range(self.getDim()):
717     for l in range(self.getDim()):
718     if util.Lsup(A[j,l]-A[l,j])>tol:
719     if verbose: print "non-symmetric PDE because A[%d,%d]!=A[%d,%d]"%(j,l,l,j)
720     out=False
721     B=self.getCoefficientOfGeneralPDE("B")
722     C=self.getCoefficientOfGeneralPDE("C")
723     if B.isEmpty() and not C.isEmpty():
724     if verbose: print "non-symmetric PDE because B is not present but C is"
725     out=False
726     elif not B.isEmpty() and C.isEmpty():
727     if verbose: print "non-symmetric PDE because C is not present but B is"
728     out=False
729     elif not B.isEmpty() and not C.isEmpty():
730 jgs 149 tol=(util.Lsup(B)+util.Lsup(C))*self.__TOL/2.
731 jgs 148 if self.getNumSolutions()>1:
732     for i in range(self.getNumEquations()):
733     for j in range(self.getDim()):
734     for k in range(self.getNumSolutions()):
735     if util.Lsup(B[i,j,k]-C[k,i,j])>tol:
736     if verbose: print "non-symmetric PDE because B[%d,%d,%d]!=C[%d,%d,%d]"%(i,j,k,k,i,j)
737     out=False
738     else:
739     for j in range(self.getDim()):
740     if util.Lsup(B[j]-C[j])>tol:
741     if verbose: print "non-symmetric PDE because B[%d]!=C[%d]"%(j,j)
742     out=False
743     if self.getNumSolutions()>1:
744     D=self.getCoefficientOfGeneralPDE("D")
745     if not D.isEmpty():
746 jgs 149 tol=util.Lsup(D)*self.__TOL
747 jgs 148 for i in range(self.getNumEquations()):
748     for k in range(self.getNumSolutions()):
749     if util.Lsup(D[i,k]-D[k,i])>tol:
750     if verbose: print "non-symmetric PDE because D[%d,%d]!=D[%d,%d]"%(i,k,k,i)
751     out=False
752 jgs 149 d=self.getCoefficientOfGeneralPDE("d")
753     if not d.isEmpty():
754     tol=util.Lsup(d)*self.__TOL
755     for i in range(self.getNumEquations()):
756     for k in range(self.getNumSolutions()):
757     if util.Lsup(d[i,k]-d[k,i])>tol:
758     if verbose: print "non-symmetric PDE because d[%d,%d]!=d[%d,%d]"%(i,k,k,i)
759     out=False
760     d_contact=self.getCoefficientOfGeneralPDE("d_contact")
761     if not d_contact.isEmpty():
762     tol=util.Lsup(d_contact)*self.__TOL
763     for i in range(self.getNumEquations()):
764     for k in range(self.getNumSolutions()):
765     if util.Lsup(d_contact[i,k]-d_contact[k,i])>tol:
766     if verbose: print "non-symmetric PDE because d_contact[%d,%d]!=d_contact[%d,%d]"%(i,k,k,i)
767     out=False
768 jgs 148 return out
769    
770     def getSolution(self,**options):
771 jgs 102 """
772 jgs 149 returns the solution of the PDE. If the solution is not valid the PDE is solved.
773 jgs 102
774 jgs 148 @return: the solution
775 jgs 149 @rtype: L{Data<escript.Data>}
776 jgs 148 @param options: solver options
777 jgs 149 @keyword verbose: True to get some information during PDE solution
778     @type verbose: C{bool}
779     @keyword reordering: reordering scheme to be used during elimination. Allowed values are
780     L{NO_REORDERING}, L{MINIMUM_FILL_IN}, L{NESTED_DISSECTION}
781 jgs 148 @keyword iter_max: maximum number of iteration steps allowed.
782 jgs 149 @keyword drop_tolerance: threshold for drupping in L{ILUT}
783     @keyword drop_storage: maximum of allowed memory in L{ILUT}
784     @keyword truncation: maximum number of residuals in L{GMRES}
785     @keyword restart: restart cycle length in L{GMRES}
786 jgs 102 """
787 jgs 148 if not self.__solution_isValid:
788     mat,f=self.getSystem()
789     if self.isUsingLumping():
790     self.__solution=self.copyConstraint(f*mat)
791     else:
792 jgs 149 options[self.__TOLERANCE_KEY]=self.getTolerance()
793 gross 387 options[self.__METHOD_KEY]=self.getSolverMethod()[0]
794     options[self.__PRECONDITIONER_KEY]=self.getSolverMethod()[1]
795 jgs 150 options[self.__PACKAGE_KEY]=self.getSolverPackage()
796 jgs 149 options[self.__SYMMETRY_KEY]=self.isSymmetric()
797 jgs 148 self.trace("PDE is resolved.")
798     self.trace("solver options: %s"%str(options))
799     self.__solution=mat.solve(f,options)
800     self.__solution_isValid=True
801     return self.__solution
802 jgs 102
803 jgs 148 def getFlux(self,u=None):
804     """
805 jgs 149 returns the flux M{J} for a given M{u}
806 jgs 102
807 jgs 149 M{J[i,j]=A[i,j,k,l]*grad(u[k])[l]+B[i,j,k]u[k]-X[i,j]}
808 jgs 102
809 jgs 151 or
810 jgs 149
811     M{J[j]=A[i,j]*grad(u)[l]+B[j]u-X[j]}
812    
813 jgs 148 @param u: argument in the flux. If u is not present or equals L{None} the current solution is used.
814 jgs 149 @type u: L{Data<escript.Data>} or None
815     @return: flux
816     @rtype: L{Data<escript.Data>}
817 jgs 148 """
818     if u==None: u=self.getSolution()
819     return util.tensormult(self.getCoefficientOfGeneralPDE("A"),util.grad(u))+util.matrixmult(self.getCoefficientOfGeneralPDE("B"),u)-util.self.getCoefficientOfGeneralPDE("X")
820     # =============================================================================
821     # solver settings:
822     # =============================================================================
823 gross 387 def setSolverMethod(self,solver=None,preconditioner=None):
824 jgs 102 """
825 jgs 122 sets a new solver
826 jgs 148
827 jgs 149 @param solver: sets a new solver method.
828 gross 431 @type solver: one of L{DEFAULT}, L{ITERATIVE} L{DIRECT}, L{CHOLEVSKY}, L{PCG}, L{CR}, L{CGS}, L{BICGSTAB}, L{SSOR}, L{GMRES}, L{PRES20}, L{LUMPING}, L{AMG}
829 gross 387 @param preconditioner: sets a new solver method.
830 gross 431 @type preconditioner: one of L{DEFAULT}, L{JACOBI} L{ILU0}, L{ILUT},L{SSOR}, L{RILU}
831 jgs 102 """
832 jgs 150 if solver==None: solve=self.DEFAULT
833 gross 387 if preconditioner==None: preconditioner=self.DEFAULT
834     if not (solver,preconditioner)==self.getSolverMethod():
835 jgs 102 self.__solver_method=solver
836 gross 387 self.__preconditioner=preconditioner
837 jgs 102 self.__checkMatrixType()
838 jgs 148 self.trace("New solver is %s"%self.getSolverMethodName())
839 jgs 102
840 jgs 148 def getSolverMethodName(self):
841     """
842     returns the name of the solver currently used
843    
844 jgs 149 @return: the name of the solver currently used.
845 jgs 148 @rtype: C{string}
846     """
847    
848     m=self.getSolverMethod()
849 jgs 150 p=self.getSolverPackage()
850 gross 425 method=""
851 gross 387 if m[0]==self.DEFAULT: method="DEFAULT"
852     elif m[0]==self.DIRECT: method= "DIRECT"
853     elif m[0]==self.ITERATIVE: method= "ITERATIVE"
854     elif m[0]==self.CHOLEVSKY: method= "CHOLEVSKY"
855     elif m[0]==self.PCG: method= "PCG"
856     elif m[0]==self.CR: method= "CR"
857     elif m[0]==self.CGS: method= "CGS"
858     elif m[0]==self.BICGSTAB: method= "BICGSTAB"
859     elif m[0]==self.SSOR: method= "SSOR"
860     elif m[0]==self.GMRES: method= "GMRES"
861     elif m[0]==self.PRES20: method= "PRES20"
862     elif m[0]==self.LUMPING: method= "LUMPING"
863 gross 431 elif m[0]==self.AMG: method= "AMG"
864 gross 425 if m[1]==self.DEFAULT: method+="+DEFAULT"
865     elif m[1]==self.JACOBI: method+= "+JACOBI"
866     elif m[1]==self.ILU0: method+= "+ILU0"
867     elif m[1]==self.ILUT: method+= "+ILUT"
868     elif m[1]==self.SSOR: method+= "+SSOR"
869 gross 431 elif m[1]==self.AMG: method+= "+AMG"
870     elif m[1]==self.RILU: method+= "+RILU"
871 jgs 150 if p==self.DEFAULT: package="DEFAULT"
872     elif p==self.PASO: package= "PASO"
873     elif p==self.MKL: package= "MKL"
874     elif p==self.SCSL: package= "SCSL"
875     elif p==self.UMFPACK: package= "UMFPACK"
876     else : method="unknown"
877     return "%s solver of %s package"%(method,package)
878 jgs 148
879 jgs 149
880 jgs 102 def getSolverMethod(self):
881     """
882 jgs 122 returns the solver method
883 jgs 149
884 jgs 151 @return: the solver method currently be used.
885 jgs 149 @rtype: C{int}
886 jgs 102 """
887 gross 387 return self.__solver_method,self.__preconditioner
888 jgs 102
889 jgs 150 def setSolverPackage(self,package=None):
890     """
891     sets a new solver package
892    
893     @param solver: sets a new solver method.
894     @type solver: one of L{DEFAULT}, L{PASO} L{SCSL}, L{MKL}, L{UMLPACK}
895     """
896     if package==None: package=self.DEFAULT
897     if not package==self.getSolverPackage():
898     self.__solver_method=solver
899     self.__checkMatrixType()
900     self.trace("New solver is %s"%self.getSolverMethodName())
901    
902     def getSolverPackage(self):
903     """
904     returns the package of the solver
905    
906 jgs 151 @return: the solver package currently being used.
907 jgs 150 @rtype: C{int}
908     """
909     return self.__solver_package
910    
911 jgs 148 def isUsingLumping(self):
912     """
913     checks if matrix lumping is used a solver method
914    
915 jgs 149 @return: True is lumping is currently used a solver method.
916 jgs 148 @rtype: C{bool}
917     """
918 gross 387 return self.getSolverMethod()[0]==self.LUMPING
919 jgs 148
920 jgs 102 def setTolerance(self,tol=1.e-8):
921     """
922 jgs 149 resets the tolerance for the solver method to tol where for an appropriate norm M{|.|}
923 jgs 148
924 jgs 149 M{|L{getResidual}()|<tol*|L{getRightHandSide}()|}
925 jgs 148
926 jgs 149 defines the stopping criterion.
927 jgs 148
928     @param tol: new tolerance for the solver. If the tol is lower then the current tolerence
929     the system will be resolved.
930 jgs 149 @type tol: positive C{float}
931 jgs 148 @raise ValueException: if tolerance is not positive.
932 jgs 102 """
933     if not tol>0:
934     raise ValueException,"Tolerance as to be positive"
935 jgs 148 if tol<self.getTolerance(): self.__invalidateSolution()
936     self.trace("New tolerance %e"%tol)
937 jgs 102 self.__tolerance=tol
938     return
939 jgs 148
940 jgs 102 def getTolerance(self):
941     """
942 jgs 122 returns the tolerance set for the solution
943 jgs 148
944     @return: tolerance currently used.
945     @rtype: C{float}
946 jgs 102 """
947     return self.__tolerance
948    
949 jgs 148 # =============================================================================
950     # symmetry flag:
951     # =============================================================================
952 jgs 102 def isSymmetric(self):
953     """
954 jgs 148 checks if symmetry is indicated.
955 jgs 149
956     @return: True is a symmetric PDE is indicated, otherwise False is returned
957     @rtype: C{bool}
958 jgs 102 """
959     return self.__sym
960    
961     def setSymmetryOn(self):
962     """
963 jgs 148 sets the symmetry flag.
964 jgs 102 """
965     if not self.isSymmetric():
966 jgs 148 self.trace("PDE is set to be symmetric")
967 jgs 102 self.__sym=True
968     self.__checkMatrixType()
969    
970     def setSymmetryOff(self):
971     """
972 jgs 148 removes the symmetry flag.
973 jgs 102 """
974     if self.isSymmetric():
975 jgs 148 self.trace("PDE is set to be unsymmetric")
976 jgs 102 self.__sym=False
977     self.__checkMatrixType()
978    
979     def setSymmetryTo(self,flag=False):
980 jgs 148 """
981     sets the symmetry flag to flag
982 jgs 149
983 jgs 148 @param flag: If flag, the symmetry flag is set otherwise the symmetry flag is released.
984     @type flag: C{bool}
985     """
986     if flag:
987     self.setSymmetryOn()
988     else:
989     self.setSymmetryOff()
990 jgs 102
991 jgs 148 # =============================================================================
992     # function space handling for the equation as well as the solution
993     # =============================================================================
994 jgs 102 def setReducedOrderOn(self):
995     """
996 jgs 148 switches on reduced order for solution and equation representation
997 jgs 149
998     @raise RuntimeError: if order reduction is altered after a coefficient has been set.
999 jgs 102 """
1000     self.setReducedOrderForSolutionOn()
1001     self.setReducedOrderForEquationOn()
1002    
1003     def setReducedOrderOff(self):
1004     """
1005 jgs 148 switches off reduced order for solution and equation representation
1006 jgs 149
1007     @raise RuntimeError: if order reduction is altered after a coefficient has been set.
1008 jgs 102 """
1009     self.setReducedOrderForSolutionOff()
1010     self.setReducedOrderForEquationOff()
1011    
1012     def setReducedOrderTo(self,flag=False):
1013     """
1014 jgs 149 sets order reduction for both solution and equation representation according to flag.
1015     @param flag: if flag is True, the order reduction is switched on for both solution and equation representation, otherwise or
1016 jgs 148 if flag is not present order reduction is switched off
1017     @type flag: C{bool}
1018 jgs 149 @raise RuntimeError: if order reduction is altered after a coefficient has been set.
1019 jgs 102 """
1020     self.setReducedOrderForSolutionTo(flag)
1021     self.setReducedOrderForEquationTo(flag)
1022    
1023 jgs 148
1024 jgs 102 def setReducedOrderForSolutionOn(self):
1025     """
1026 jgs 148 switches on reduced order for solution representation
1027 jgs 149
1028     @raise RuntimeError: if order reduction is altered after a coefficient has been set.
1029 jgs 102 """
1030 jgs 149 if not self.__reduce_solution_order:
1031     if self.__altered_coefficients:
1032     raise RuntimeError,"order cannot be altered after coefficients have been defined."
1033 jgs 148 self.trace("Reduced order is used to solution representation.")
1034 jgs 149 self.__reduce_solution_order=True
1035 jgs 148 self.__resetSystem()
1036 jgs 102
1037     def setReducedOrderForSolutionOff(self):
1038     """
1039 jgs 148 switches off reduced order for solution representation
1040 jgs 149
1041     @raise RuntimeError: if order reduction is altered after a coefficient has been set.
1042 jgs 102 """
1043 jgs 149 if self.__reduce_solution_order:
1044     if self.__altered_coefficients:
1045     raise RuntimeError,"order cannot be altered after coefficients have been defined."
1046 jgs 148 self.trace("Full order is used to interpolate solution.")
1047 jgs 149 self.__reduce_solution_order=False
1048 jgs 148 self.__resetSystem()
1049 jgs 102
1050     def setReducedOrderForSolutionTo(self,flag=False):
1051     """
1052 jgs 122 sets order for test functions according to flag
1053 jgs 102
1054 jgs 149 @param flag: if flag is True, the order reduction is switched on for solution representation, otherwise or
1055 jgs 148 if flag is not present order reduction is switched off
1056     @type flag: C{bool}
1057 jgs 149 @raise RuntimeError: if order reduction is altered after a coefficient has been set.
1058 jgs 102 """
1059     if flag:
1060     self.setReducedOrderForSolutionOn()
1061     else:
1062     self.setReducedOrderForSolutionOff()
1063 jgs 148
1064 jgs 102 def setReducedOrderForEquationOn(self):
1065     """
1066 jgs 148 switches on reduced order for equation representation
1067 jgs 149
1068     @raise RuntimeError: if order reduction is altered after a coefficient has been set.
1069 jgs 102 """
1070 jgs 149 if not self.__reduce_equation_order:
1071     if self.__altered_coefficients:
1072     raise RuntimeError,"order cannot be altered after coefficients have been defined."
1073 jgs 148 self.trace("Reduced order is used for test functions.")
1074 jgs 149 self.__reduce_equation_order=True
1075 jgs 148 self.__resetSystem()
1076 jgs 102
1077     def setReducedOrderForEquationOff(self):
1078     """
1079 jgs 148 switches off reduced order for equation representation
1080 jgs 149
1081     @raise RuntimeError: if order reduction is altered after a coefficient has been set.
1082 jgs 102 """
1083 jgs 149 if self.__reduce_equation_order:
1084     if self.__altered_coefficients:
1085     raise RuntimeError,"order cannot be altered after coefficients have been defined."
1086 jgs 148 self.trace("Full order is used for test functions.")
1087 jgs 149 self.__reduce_equation_order=False
1088 jgs 148 self.__resetSystem()
1089 jgs 102
1090     def setReducedOrderForEquationTo(self,flag=False):
1091     """
1092 jgs 122 sets order for test functions according to flag
1093 jgs 102
1094 jgs 149 @param flag: if flag is True, the order reduction is switched on for equation representation, otherwise or
1095 jgs 148 if flag is not present order reduction is switched off
1096     @type flag: C{bool}
1097 jgs 149 @raise RuntimeError: if order reduction is altered after a coefficient has been set.
1098 jgs 102 """
1099     if flag:
1100     self.setReducedOrderForEquationOn()
1101     else:
1102     self.setReducedOrderForEquationOff()
1103 jgs 148
1104     # =============================================================================
1105     # private method:
1106     # =============================================================================
1107     def __checkMatrixType(self):
1108     """
1109     reassess the matrix type and, if a new matrix is needed, resets the system.
1110     """
1111 gross 387 new_matrix_type=self.getDomain().getSystemMatrixTypeId(self.getSolverMethod()[0],self.getSolverPackage(),self.isSymmetric())
1112 jgs 148 if not new_matrix_type==self.__matrix_type:
1113     self.trace("Matrix type is now %d."%new_matrix_type)
1114     self.__matrix_type=new_matrix_type
1115     self.__resetSystem()
1116 jgs 149 #
1117 jgs 148 # rebuild switches :
1118 jgs 149 #
1119 jgs 148 def __invalidateSolution(self):
1120     """
1121     indicates the PDE has to be resolved if the solution is requested
1122     """
1123     if self.__solution_isValid: self.trace("PDE has to be resolved.")
1124     self.__solution_isValid=False
1125    
1126     def __invalidateOperator(self):
1127     """
1128     indicates the operator has to be rebuilt next time it is used
1129     """
1130 jgs 149 if self.__operator_is_Valid: self.trace("Operator has to be rebuilt.")
1131 jgs 148 self.__invalidateSolution()
1132 jgs 149 self.__operator_is_Valid=False
1133 jgs 148
1134     def __invalidateRightHandSide(self):
1135     """
1136     indicates the right hand side has to be rebuild next time it is used
1137     """
1138     if self.__righthandside_isValid: self.trace("Right hand side has to be rebuilt.")
1139     self.__invalidateSolution()
1140     self.__righthandside_isValid=False
1141    
1142     def __invalidateSystem(self):
1143     """
1144     annonced that everthing has to be rebuild:
1145     """
1146     if self.__righthandside_isValid: self.trace("System has to be rebuilt.")
1147     self.__invalidateSolution()
1148     self.__invalidateOperator()
1149     self.__invalidateRightHandSide()
1150    
1151     def __resetSystem(self):
1152     """
1153     annonced that everthing has to be rebuild:
1154     """
1155     self.trace("New System is built from scratch.")
1156     self.__operator=escript.Operator()
1157 jgs 149 self.__operator_is_Valid=False
1158 jgs 148 self.__righthandside=escript.Data()
1159     self.__righthandside_isValid=False
1160     self.__solution=escript.Data()
1161     self.__solution_isValid=False
1162 jgs 149 #
1163 jgs 148 # system initialization:
1164 jgs 149 #
1165 jgs 121 def __getNewOperator(self):
1166 jgs 102 """
1167 jgs 148 returns an instance of a new operator
1168 jgs 102 """
1169 jgs 148 self.trace("New operator is allocated.")
1170 jgs 102 return self.getDomain().newOperator( \
1171     self.getNumEquations(), \
1172     self.getFunctionSpaceForEquation(), \
1173     self.getNumSolutions(), \
1174     self.getFunctionSpaceForSolution(), \
1175     self.__matrix_type)
1176    
1177 jgs 148 def __getNewRightHandSide(self):
1178 jgs 102 """
1179 jgs 148 returns an instance of a new right hand side
1180 jgs 102 """
1181 jgs 148 self.trace("New right hand side is allocated.")
1182 jgs 121 if self.getNumEquations()>1:
1183 jgs 148 return escript.Data(0.,(self.getNumEquations(),),self.getFunctionSpaceForEquation(),True)
1184 jgs 121 else:
1185 jgs 148 return escript.Data(0.,(),self.getFunctionSpaceForEquation(),True)
1186 jgs 102
1187 jgs 121 def __getNewSolution(self):
1188 jgs 102 """
1189 jgs 148 returns an instance of a new solution
1190 jgs 102 """
1191 jgs 148 self.trace("New solution is allocated.")
1192 jgs 121 if self.getNumSolutions()>1:
1193     return escript.Data(0.,(self.getNumSolutions(),),self.getFunctionSpaceForSolution(),True)
1194     else:
1195     return escript.Data(0.,(),self.getFunctionSpaceForSolution(),True)
1196 jgs 102
1197 jgs 148 def __makeFreshSolution(self):
1198     """
1199     makes sure that the solution is instantiated and returns it initialized by zeros
1200     """
1201     if self.__solution.isEmpty():
1202     self.__solution=self.__getNewSolution()
1203     else:
1204     self.__solution*=0
1205     self.trace("Solution is reset to zero.")
1206     return self.__solution
1207    
1208     def __makeFreshRightHandSide(self):
1209     """
1210     makes sure that the right hand side is instantiated and returns it initialized by zeros
1211     """
1212     if self.__righthandside.isEmpty():
1213     self.__righthandside=self.__getNewRightHandSide()
1214     else:
1215     self.__righthandside*=0
1216     self.trace("Right hand side is reset to zero.")
1217     return self.__righthandside
1218    
1219 jgs 121 def __makeFreshOperator(self):
1220 jgs 102 """
1221 jgs 148 makes sure that the operator is instantiated and returns it initialized by zeros
1222 jgs 102 """
1223     if self.__operator.isEmpty():
1224 jgs 121 self.__operator=self.__getNewOperator()
1225 jgs 102 else:
1226 jgs 149 self.__operator.resetValues()
1227 jgs 148 self.trace("Operator reset to zero")
1228 jgs 102 return self.__operator
1229    
1230 jgs 148 def __applyConstraint(self):
1231     """
1232     applies the constraints defined by q and r to the system
1233     """
1234     if not self.isUsingLumping():
1235     q=self.getCoefficientOfGeneralPDE("q")
1236     r=self.getCoefficientOfGeneralPDE("r")
1237     if not q.isEmpty() and not self.__operator.isEmpty():
1238     # q is the row and column mask to indicate where constraints are set:
1239     row_q=escript.Data(q,self.getFunctionSpaceForEquation())
1240     col_q=escript.Data(q,self.getFunctionSpaceForSolution())
1241     u=self.__getNewSolution()
1242     if r.isEmpty():
1243     r_s=self.__getNewSolution()
1244     else:
1245     r_s=escript.Data(r,self.getFunctionSpaceForSolution())
1246     u.copyWithMask(r_s,col_q)
1247 jgs 149 if not self.__righthandside.isEmpty():
1248 jgs 148 self.__righthandside-=self.__operator*u
1249     self.__righthandside=self.copyConstraint(self.__righthandside)
1250     self.__operator.nullifyRowsAndCols(row_q,col_q,1.)
1251     # =============================================================================
1252     # function giving access to coefficients of the general PDE:
1253     # =============================================================================
1254     def getCoefficientOfGeneralPDE(self,name):
1255 jgs 108 """
1256 jgs 148 return the value of the coefficient name of the general PDE.
1257    
1258 jgs 149 @note: This method is called by the assembling routine it can be overwritten
1259 jgs 148 to map coefficients of a particular PDE to the general PDE.
1260 jgs 149 @param name: name of the coefficient requested.
1261 jgs 148 @type name: C{string}
1262 jgs 149 @return: the value of the coefficient name
1263     @rtype: L{Data<escript.Data>}
1264 jgs 148 @raise IllegalCoefficient: if name is not one of coefficients
1265 jgs 149 M{A}, M{B}, M{C}, M{D}, M{X}, M{Y}, M{d}, M{y}, M{d_contact}, M{y_contact}, M{r} or M{q}.
1266 jgs 108 """
1267 jgs 148 if self.hasCoefficientOfGeneralPDE(name):
1268     return self.getCoefficient(name)
1269     else:
1270     raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name
1271 jgs 108
1272 jgs 148 def hasCoefficientOfGeneralPDE(self,name):
1273 jgs 108 """
1274 jgs 148 checks if name is a the name of a coefficient of the general PDE.
1275 jgs 149
1276 jgs 148 @param name: name of the coefficient enquired.
1277     @type name: C{string}
1278 jgs 149 @return: True if name is the name of a coefficient of the general PDE. Otherwise False.
1279     @rtype: C{bool}
1280    
1281 jgs 108 """
1282 jgs 148 return self.__COEFFICIENTS_OF_GENEARL_PDE.has_key(name)
1283 jgs 108
1284 jgs 148 def createCoefficientOfGeneralPDE(self,name):
1285 jgs 108 """
1286 jgs 148 returns a new instance of a coefficient for coefficient name of the general PDE
1287    
1288     @param name: name of the coefficient requested.
1289     @type name: C{string}
1290 jgs 149 @return: a coefficient name initialized to 0.
1291     @rtype: L{Data<escript.Data>}
1292 jgs 148 @raise IllegalCoefficient: if name is not one of coefficients
1293 jgs 149 M{A}, M{B}, M{C}, M{D}, M{X}, M{Y}, M{d}, M{y}, M{d_contact}, M{y_contact}, M{r} or M{q}.
1294 jgs 108 """
1295 jgs 148 if self.hasCoefficientOfGeneralPDE(name):
1296     return escript.Data(0,self.getShapeOfCoefficientOfGeneralPDE(name),self.getFunctionSpaceForCoefficientOfGeneralPDE(name))
1297 jgs 108 else:
1298 jgs 148 raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name
1299 jgs 108
1300 jgs 148 def getFunctionSpaceForCoefficientOfGeneralPDE(self,name):
1301 jgs 108 """
1302 jgs 149 return the L{FunctionSpace<escript.FunctionSpace>} to be used for coefficient name of the general PDE
1303 jgs 148
1304     @param name: name of the coefficient enquired.
1305     @type name: C{string}
1306 jgs 149 @return: the function space to be used for coefficient name
1307     @rtype: L{FunctionSpace<escript.FunctionSpace>}
1308 jgs 148 @raise IllegalCoefficient: if name is not one of coefficients
1309 jgs 149 M{A}, M{B}, M{C}, M{D}, M{X}, M{Y}, M{d}, M{y}, M{d_contact}, M{y_contact}, M{r} or M{q}.
1310 jgs 108 """
1311 jgs 148 if self.hasCoefficientOfGeneralPDE(name):
1312     return self.__COEFFICIENTS_OF_GENEARL_PDE[name].getFunctionSpace(self.getDomain())
1313 jgs 108 else:
1314 jgs 148 raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name
1315 jgs 108
1316 jgs 148 def getShapeOfCoefficientOfGeneralPDE(self,name):
1317     """
1318     return the shape of the coefficient name of the general PDE
1319 jgs 108
1320 jgs 148 @param name: name of the coefficient enquired.
1321     @type name: C{string}
1322 jgs 149 @return: the shape of the coefficient name
1323     @rtype: C{tuple} of C{int}
1324 jgs 148 @raise IllegalCoefficient: if name is not one of coefficients
1325 jgs 149 M{A}, M{B}, M{C}, M{D}, M{X}, M{Y}, M{d}, M{y}, M{d_contact}, M{y_contact}, M{r} or M{q}.
1326 jgs 148 """
1327     if self.hasCoefficientOfGeneralPDE(name):
1328     return self.__COEFFICIENTS_OF_GENEARL_PDE[name].getShape(self.getDomain(),self.getNumEquations(),self.getNumSolutions())
1329     else:
1330     raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name
1331 jgs 108
1332 jgs 148 # =============================================================================
1333     # functions giving access to coefficients of a particular PDE implementation:
1334     # =============================================================================
1335     def getCoefficient(self,name):
1336     """
1337     returns the value of the coefficient name
1338 jgs 108
1339 jgs 149 @param name: name of the coefficient requested.
1340 jgs 148 @type name: C{string}
1341 jgs 149 @return: the value of the coefficient name
1342     @rtype: L{Data<escript.Data>}
1343 jgs 148 @raise IllegalCoefficient: if name is not a coefficient of the PDE.
1344     """
1345     if self.hasCoefficient(name):
1346     return self.COEFFICIENTS[name].getValue()
1347     else:
1348     raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name
1349 jgs 108
1350 jgs 148 def hasCoefficient(self,name):
1351     """
1352     return True if name is the name of a coefficient
1353 jgs 108
1354 jgs 148 @param name: name of the coefficient enquired.
1355     @type name: C{string}
1356 jgs 149 @return: True if name is the name of a coefficient of the general PDE. Otherwise False.
1357     @rtype: C{bool}
1358 jgs 148 """
1359     return self.COEFFICIENTS.has_key(name)
1360 jgs 108
1361 jgs 148 def createCoefficient(self, name):
1362     """
1363 jgs 149 create a L{Data<escript.Data>} object corresponding to coefficient name
1364 jgs 108
1365 jgs 149 @return: a coefficient name initialized to 0.
1366     @rtype: L{Data<escript.Data>}
1367 jgs 148 @raise IllegalCoefficient: if name is not a coefficient of the PDE.
1368     """
1369     if self.hasCoefficient(name):
1370     return escript.Data(0.,self.getShapeOfCoefficient(name),self.getFunctionSpaceForCoefficient(name))
1371     else:
1372     raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name
1373    
1374     def getFunctionSpaceForCoefficient(self,name):
1375     """
1376 jgs 149 return the L{FunctionSpace<escript.FunctionSpace>} to be used for coefficient name
1377 jgs 148
1378     @param name: name of the coefficient enquired.
1379     @type name: C{string}
1380 jgs 149 @return: the function space to be used for coefficient name
1381     @rtype: L{FunctionSpace<escript.FunctionSpace>}
1382 jgs 148 @raise IllegalCoefficient: if name is not a coefficient of the PDE.
1383     """
1384     if self.hasCoefficient(name):
1385     return self.COEFFICIENTS[name].getFunctionSpace(self.getDomain())
1386     else:
1387     raise ValueError,"unknown coefficient %s requested"%name
1388     def getShapeOfCoefficient(self,name):
1389     """
1390 jgs 149 return the shape of the coefficient name
1391 jgs 148
1392     @param name: name of the coefficient enquired.
1393     @type name: C{string}
1394 jgs 149 @return: the shape of the coefficient name
1395     @rtype: C{tuple} of C{int}
1396 jgs 148 @raise IllegalCoefficient: if name is not a coefficient of the PDE.
1397     """
1398     if self.hasCoefficient(name):
1399     return self.COEFFICIENTS[name].getShape(self.getDomain(),self.getNumEquations(),self.getNumSolutions())
1400     else:
1401     raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name
1402    
1403     def resetCoefficients(self):
1404     """
1405     resets all coefficients to there default values.
1406     """
1407     for i in self.COEFFICIENTS.iterkeys():
1408     self.COEFFICIENTS[i].resetValue()
1409    
1410     def alteredCoefficient(self,name):
1411     """
1412     announce that coefficient name has been changed
1413    
1414     @param name: name of the coefficient enquired.
1415     @type name: C{string}
1416     @raise IllegalCoefficient: if name is not a coefficient of the PDE.
1417 jgs 149 @note: if name is q or r, the method will not trigger a rebuilt of the system as constraints are applied to the solved system.
1418 jgs 148 """
1419     if self.hasCoefficient(name):
1420     self.trace("Coefficient %s has been altered."%name)
1421 jgs 149 if not ((name=="q" or name=="r") and self.isUsingLumping()):
1422     if self.COEFFICIENTS[name].isAlteringOperator(): self.__invalidateOperator()
1423     if self.COEFFICIENTS[name].isAlteringRightHandSide(): self.__invalidateRightHandSide()
1424 jgs 148 else:
1425     raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name
1426    
1427     def copyConstraint(self,u):
1428 jgs 108 """
1429 jgs 149 copies the constraint into u and returns u.
1430    
1431 jgs 148 @param u: a function of rank 0 is a single PDE is solved and of shape (numSolution,) for a system of PDEs
1432 jgs 149 @type u: L{Data<escript.Data>}
1433     @return: the input u modified by the constraints.
1434     @rtype: L{Data<escript.Data>}
1435     @warning: u is altered if it has the appropriate L{FunctionSpace<escript.FunctionSpace>}
1436 jgs 108 """
1437 jgs 148 q=self.getCoefficientOfGeneralPDE("q")
1438     r=self.getCoefficientOfGeneralPDE("r")
1439     if not q.isEmpty():
1440     if u.isEmpty(): u=escript.Data(0.,q.getShape(),q.getFunctionSpace())
1441     if r.isEmpty():
1442     r=escript.Data(0,u.getShape(),u.getFunctionSpace())
1443     else:
1444     r=escript.Data(r,u.getFunctionSpace())
1445     u.copyWithMask(r,escript.Data(q,u.getFunctionSpace()))
1446     return u
1447    
1448     def setValue(self,**coefficients):
1449     """
1450     sets new values to coefficients
1451    
1452 jgs 149 @param coefficients: new values assigned to coefficients
1453     @keyword A: value for coefficient A.
1454     @type A: any type that can be casted to L{Data<escript.Data>} object on L{Function<escript.Function>}.
1455 jgs 148 @keyword B: value for coefficient B
1456 jgs 149 @type B: any type that can be casted to L{Data<escript.Data>} object on L{Function<escript.Function>}.
1457 jgs 148 @keyword C: value for coefficient C
1458 jgs 149 @type C: any type that can be casted to L{Data<escript.Data>} object on L{Function<escript.Function>}.
1459 jgs 148 @keyword D: value for coefficient D
1460 jgs 149 @type D: any type that can be casted to L{Data<escript.Data>} object on L{Function<escript.Function>}.
1461 jgs 148 @keyword X: value for coefficient X
1462 jgs 149 @type X: any type that can be casted to L{Data<escript.Data>} object on L{Function<escript.Function>}.
1463 jgs 148 @keyword Y: value for coefficient Y
1464 jgs 149 @type Y: any type that can be casted to L{Data<escript.Data>} object on L{Function<escript.Function>}.
1465 jgs 148 @keyword d: value for coefficient d
1466 jgs 149 @type d: any type that can be casted to L{Data<escript.Data>} object on L{FunctionOnBoundary<escript.FunctionOnBoundary>}.
1467 jgs 148 @keyword y: value for coefficient y
1468 jgs 149 @type y: any type that can be casted to L{Data<escript.Data>} object on L{FunctionOnBoundary<escript.FunctionOnBoundary>}.
1469 jgs 148 @keyword d_contact: value for coefficient d_contact
1470 jgs 149 @type d_contact: any type that can be casted to L{Data<escript.Data>} object on L{FunctionOnContactOne<escript.FunctionOnContactOne>}.
1471     or L{FunctionOnContactZero<escript.FunctionOnContactZero>}.
1472 jgs 148 @keyword y_contact: value for coefficient y_contact
1473 jgs 149 @type y_contact: any type that can be casted to L{Data<escript.Data>} object on L{FunctionOnContactOne<escript.FunctionOnContactOne>}.
1474     or L{FunctionOnContactZero<escript.FunctionOnContactZero>}.
1475 jgs 148 @keyword r: values prescribed to the solution at the locations of constraints
1476 jgs 149 @type r: any type that can be casted to L{Data<escript.Data>} object on L{Solution<escript.Solution>} or L{ReducedSolution<escript.ReducedSolution>}
1477 jgs 148 depending of reduced order is used for the solution.
1478     @keyword q: mask for location of constraints
1479 jgs 149 @type q: any type that can be casted to L{Data<escript.Data>} object on L{Solution<escript.Solution>} or L{ReducedSolution<escript.ReducedSolution>}
1480 jgs 148 depending of reduced order is used for the representation of the equation.
1481     @raise IllegalCoefficient: if an unknown coefficient keyword is used.
1482     """
1483 jgs 108 # check if the coefficients are legal:
1484     for i in coefficients.iterkeys():
1485     if not self.hasCoefficient(i):
1486 jgs 148 raise IllegalCoefficient,"Attempt to set unknown coefficient %s"%i
1487 jgs 108 # if the number of unknowns or equations is still unknown we try to estimate them:
1488 jgs 148 if self.__numEquations==None or self.__numSolutions==None:
1489 jgs 108 for i,d in coefficients.iteritems():
1490     if hasattr(d,"shape"):
1491     s=d.shape
1492     elif hasattr(d,"getShape"):
1493     s=d.getShape()
1494     else:
1495     s=numarray.array(d).shape
1496     if s!=None:
1497     # get number of equations and number of unknowns:
1498 jgs 148 res=self.COEFFICIENTS[i].estimateNumEquationsAndNumSolutions(self.getDomain(),s)
1499 jgs 108 if res==None:
1500 jgs 148 raise IllegalCoefficientValue,"Illegal shape %s of coefficient %s"%(s,i)
1501 jgs 108 else:
1502 jgs 148 if self.__numEquations==None: self.__numEquations=res[0]
1503     if self.__numSolutions==None: self.__numSolutions=res[1]
1504     if self.__numEquations==None: raise UndefinedPDEError,"unidententified number of equations"
1505     if self.__numSolutions==None: raise UndefinedPDEError,"unidententified number of solutions"
1506 jgs 108 # now we check the shape of the coefficient if numEquations and numSolutions are set:
1507     for i,d in coefficients.iteritems():
1508 jgs 148 try:
1509 jgs 149 self.COEFFICIENTS[i].setValue(self.getDomain(),self.getNumEquations(),self.getNumSolutions(),self.reduceEquationOrder(),self.reduceSolutionOrder(),d)
1510 jgs 148 except IllegalCoefficientValue,m:
1511     raise IllegalCoefficientValue("Coefficient %s:%s"%(i,m))
1512 jgs 108 self.alteredCoefficient(i)
1513    
1514 jgs 149 self.__altered_coefficients=True
1515 jgs 148 # check if the systrem is inhomogeneous:
1516     if len(coefficients)>0 and not self.isUsingLumping():
1517     q=self.getCoefficientOfGeneralPDE("q")
1518     r=self.getCoefficientOfGeneralPDE("r")
1519     homogeneous_constraint=True
1520     if not q.isEmpty() and not r.isEmpty():
1521 jgs 149 if util.Lsup(q*r)>=1.e-13*util.Lsup(r):
1522 jgs 148 self.trace("Inhomogeneous constraint detected.")
1523     self.__invalidateSystem()
1524 jgs 108
1525     def getSystem(self):
1526     """
1527 jgs 122 return the operator and right hand side of the PDE
1528 jgs 149
1529     @return: the discrete version of the PDE
1530     @rtype: C{tuple} of L{Operator,<escript.Operator>} and L{Data<escript.Data>}.
1531 jgs 108 """
1532 jgs 149 if not self.__operator_is_Valid or not self.__righthandside_isValid:
1533 jgs 108 if self.isUsingLumping():
1534 jgs 149 if not self.__operator_is_Valid:
1535 jgs 148 if not self.getFunctionSpaceForEquation()==self.getFunctionSpaceForSolution(): raise TypeError,"Lumped matrix requires same order for equations and unknowns"
1536 jgs 149 if not self.getCoefficientOfGeneralPDE("A").isEmpty(): raise Warning,"Using coefficient A in lumped matrix can produce wrong results"
1537     if not self.getCoefficientOfGeneralPDE("B").isEmpty(): raise Warning,"Using coefficient B in lumped matrix can produce wrong results"
1538     if not self.getCoefficientOfGeneralPDE("C").isEmpty(): raise Warning,"Using coefficient C in lumped matrix can produce wrong results"
1539 jgs 121 mat=self.__getNewOperator()
1540 jgs 108 self.getDomain().addPDEToSystem(mat,escript.Data(), \
1541 jgs 148 self.getCoefficientOfGeneralPDE("A"), \
1542     self.getCoefficientOfGeneralPDE("B"), \
1543     self.getCoefficientOfGeneralPDE("C"), \
1544     self.getCoefficientOfGeneralPDE("D"), \
1545 jgs 108 escript.Data(), \
1546     escript.Data(), \
1547 jgs 148 self.getCoefficientOfGeneralPDE("d"), \
1548 jgs 108 escript.Data(),\
1549 jgs 148 self.getCoefficientOfGeneralPDE("d_contact"), \
1550 jgs 108 escript.Data())
1551 jgs 148 self.__operator=1./(mat*escript.Data(1,(self.getNumSolutions(),),self.getFunctionSpaceForSolution(),True))
1552     del mat
1553     self.trace("New lumped operator has been built.")
1554 jgs 149 self.__operator_is_Valid=True
1555 jgs 108 if not self.__righthandside_isValid:
1556 jgs 121 self.getDomain().addPDEToRHS(self.__makeFreshRightHandSide(), \
1557 jgs 148 self.getCoefficientOfGeneralPDE("X"), \
1558     self.getCoefficientOfGeneralPDE("Y"),\
1559     self.getCoefficientOfGeneralPDE("y"),\
1560     self.getCoefficientOfGeneralPDE("y_contact"))
1561     self.trace("New right hand side as been built.")
1562 jgs 108 self.__righthandside_isValid=True
1563     else:
1564 jgs 149 if not self.__operator_is_Valid and not self.__righthandside_isValid:
1565 jgs 121 self.getDomain().addPDEToSystem(self.__makeFreshOperator(),self.__makeFreshRightHandSide(), \
1566 jgs 148 self.getCoefficientOfGeneralPDE("A"), \
1567     self.getCoefficientOfGeneralPDE("B"), \
1568     self.getCoefficientOfGeneralPDE("C"), \
1569     self.getCoefficientOfGeneralPDE("D"), \
1570     self.getCoefficientOfGeneralPDE("X"), \
1571     self.getCoefficientOfGeneralPDE("Y"), \
1572     self.getCoefficientOfGeneralPDE("d"), \
1573     self.getCoefficientOfGeneralPDE("y"), \
1574     self.getCoefficientOfGeneralPDE("d_contact"), \
1575     self.getCoefficientOfGeneralPDE("y_contact"))
1576 jgs 108 self.__applyConstraint()
1577 jgs 148 self.__righthandside=self.copyConstraint(self.__righthandside)
1578     self.trace("New system has been built.")
1579 jgs 149 self.__operator_is_Valid=True
1580 jgs 108 self.__righthandside_isValid=True
1581     elif not self.__righthandside_isValid:
1582 jgs 121 self.getDomain().addPDEToRHS(self.__makeFreshRightHandSide(), \
1583 jgs 148 self.getCoefficientOfGeneralPDE("X"), \
1584     self.getCoefficientOfGeneralPDE("Y"),\
1585     self.getCoefficientOfGeneralPDE("y"),\
1586     self.getCoefficientOfGeneralPDE("y_contact"))
1587     self.__righthandside=self.copyConstraint(self.__righthandside)
1588     self.trace("New right hand side has been built.")
1589 jgs 108 self.__righthandside_isValid=True
1590 jgs 149 elif not self.__operator_is_Valid:
1591 jgs 121 self.getDomain().addPDEToSystem(self.__makeFreshOperator(),escript.Data(), \
1592 jgs 148 self.getCoefficientOfGeneralPDE("A"), \
1593     self.getCoefficientOfGeneralPDE("B"), \
1594     self.getCoefficientOfGeneralPDE("C"), \
1595     self.getCoefficientOfGeneralPDE("D"), \
1596 jgs 108 escript.Data(), \
1597     escript.Data(), \
1598 jgs 148 self.getCoefficientOfGeneralPDE("d"), \
1599 jgs 108 escript.Data(),\
1600 jgs 148 self.getCoefficientOfGeneralPDE("d_contact"), \
1601 jgs 108 escript.Data())
1602     self.__applyConstraint()
1603 jgs 148 self.trace("New operator has been built.")
1604 jgs 149 self.__operator_is_Valid=True
1605 jgs 108 return (self.__operator,self.__righthandside)
1606 jgs 102
1607    
1608 jgs 149 class Poisson(LinearPDE):
1609     """
1610     Class to define a Poisson equation problem, which is genear L{LinearPDE} of the form
1611 jgs 102
1612 jgs 149 M{-grad(grad(u)[j])[j] = f}
1613 jgs 102
1614 jgs 149 with natural boundary conditons
1615    
1616     M{n[j]*grad(u)[j] = 0 }
1617    
1618     and constraints:
1619    
1620     M{u=0} where M{q>0}
1621    
1622 jgs 108 """
1623 jgs 148
1624 gross 328 def __init__(self,domain,debug=False):
1625 jgs 149 """
1626     initializes a new Poisson equation
1627 jgs 104
1628 jgs 149 @param domain: domain of the PDE
1629     @type domain: L{Domain<escript.Domain>}
1630     @param debug: if True debug informations are printed.
1631 jgs 102
1632 jgs 149 """
1633 jgs 151 super(Poisson, self).__init__(domain,1,1,debug)
1634 jgs 149 self.COEFFICIENTS={"f": PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.BY_EQUATION,),PDECoefficient.RIGHTHANDSIDE),
1635     "q": PDECoefficient(PDECoefficient.SOLUTION,(PDECoefficient.BY_EQUATION,),PDECoefficient.BOTH)}
1636     self.setSymmetryOn()
1637 jgs 102
1638 jgs 149 def setValue(self,**coefficients):
1639     """
1640     sets new values to coefficients
1641 jgs 102
1642 jgs 149 @param coefficients: new values assigned to coefficients
1643     @keyword f: value for right hand side M{f}
1644     @type f: any type that can be casted to L{Scalar<escript.Scalar>} object on L{Function<escript.Function>}.
1645     @keyword q: mask for location of constraints
1646     @type q: any type that can be casted to rank zeo L{Data<escript.Data>} object on L{Solution<escript.Solution>} or L{ReducedSolution<escript.ReducedSolution>}
1647     depending of reduced order is used for the representation of the equation.
1648     @raise IllegalCoefficient: if an unknown coefficient keyword is used.
1649     """
1650 jgs 151 super(Poisson, self).setValue(**coefficients)
1651 jgs 102
1652 jgs 149 def getCoefficientOfGeneralPDE(self,name):
1653     """
1654     return the value of the coefficient name of the general PDE
1655     @param name: name of the coefficient requested.
1656     @type name: C{string}
1657     @return: the value of the coefficient name
1658     @rtype: L{Data<escript.Data>}
1659     @raise IllegalCoefficient: if name is not one of coefficients
1660     M{A}, M{B}, M{C}, M{D}, M{X}, M{Y}, M{d}, M{y}, M{d_contact}, M{y_contact}, M{r} or M{q}.
1661     @note: This method is called by the assembling routine to map the Possion equation onto the general PDE.
1662     """
1663     if name == "A" :
1664 gross 328 return escript.Data(util.kronecker(self.getDim()),escript.Function(self.getDomain()))
1665 jgs 149 elif name == "B" :
1666     return escript.Data()
1667     elif name == "C" :
1668     return escript.Data()
1669     elif name == "D" :
1670     return escript.Data()
1671     elif name == "X" :
1672     return escript.Data()
1673     elif name == "Y" :
1674     return self.getCoefficient("f")
1675     elif name == "d" :
1676     return escript.Data()
1677     elif name == "y" :
1678     return escript.Data()
1679     elif name == "d_contact" :
1680     return escript.Data()
1681     elif name == "y_contact" :
1682     return escript.Data()
1683     elif name == "r" :
1684     return escript.Data()
1685     elif name == "q" :
1686     return self.getCoefficient("q")
1687     else:
1688     raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name
1689 jgs 102
1690 jgs 149 class Helmholtz(LinearPDE):
1691     """
1692     Class to define a Helmhotz equation problem, which is genear L{LinearPDE} of the form
1693    
1694     M{S{omega}*u - grad(k*grad(u)[j])[j] = f}
1695    
1696     with natural boundary conditons
1697    
1698     M{k*n[j]*grad(u)[j] = g- S{alpha}u }
1699    
1700 jgs 122 and constraints:
1701 jgs 102
1702 jgs 149 M{u=r} where M{q>0}
1703    
1704 jgs 110 """
1705 jgs 149
1706     def __init__(self,domain,debug=False):
1707     """
1708     initializes a new Poisson equation
1709    
1710     @param domain: domain of the PDE
1711     @type domain: L{Domain<escript.Domain>}
1712     @param debug: if True debug informations are printed.
1713    
1714     """
1715 jgs 151 super(Helmholtz, self).__init__(domain,1,1,debug)
1716 jgs 149 self.COEFFICIENTS={"omega": PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.BY_EQUATION,),PDECoefficient.OPERATOR),
1717     "k": PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.BY_EQUATION,),PDECoefficient.OPERATOR),
1718     "f": PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.BY_EQUATION,),PDECoefficient.RIGHTHANDSIDE),
1719     "alpha": PDECoefficient(PDECoefficient.BOUNDARY,(PDECoefficient.BY_EQUATION,),PDECoefficient.OPERATOR),
1720     "g": PDECoefficient(PDECoefficient.BOUNDARY,(PDECoefficient.BY_EQUATION,),PDECoefficient.RIGHTHANDSIDE),
1721     "r": PDECoefficient(PDECoefficient.SOLUTION,(PDECoefficient.BY_EQUATION,),PDECoefficient.BOTH),
1722     "q": PDECoefficient(PDECoefficient.SOLUTION,(PDECoefficient.BY_EQUATION,),PDECoefficient.BOTH)}
1723     self.setSymmetryOn()
1724    
1725     def setValue(self,**coefficients):
1726     """
1727     sets new values to coefficients
1728    
1729     @param coefficients: new values assigned to coefficients
1730     @keyword omega: value for coefficient M{S{omega}}
1731     @type omega: any type that can be casted to L{Scalar<escript.Scalar>} object on L{Function<escript.Function>}.
1732     @keyword k: value for coefficeint M{k}
1733     @type k: any type that can be casted to L{Scalar<escript.Scalar>} object on L{Function<escript.Function>}.
1734     @keyword f: value for right hand side M{f}
1735     @type f: any type that can be casted to L{Scalar<escript.Scalar>} object on L{Function<escript.Function>}.
1736     @keyword alpha: value for right hand side M{S{alpha}}
1737     @type alpha: any type that can be casted to L{Scalar<escript.Scalar>} object on L{FunctionOnBoundary<escript.FunctionOnBoundary>}.
1738     @keyword g: value for right hand side M{g}
1739     @type g: any type that can be casted to L{Scalar<escript.Scalar>} object on L{FunctionOnBoundary<escript.FunctionOnBoundary>}.
1740     @keyword r: prescribed values M{r} for the solution in constraints.
1741     @type r: any type that can be casted to L{Scalar<escript.Scalar>} object on L{Solution<escript.Solution>} or L{ReducedSolution<escript.ReducedSolution>}
1742     depending of reduced order is used for the representation of the equation.
1743     @keyword q: mask for location of constraints
1744     @type q: any type that can be casted to L{Scalar<escript.Scalar>} object on L{Solution<escript.Solution>} or L{ReducedSolution<escript.ReducedSolution>}
1745     depending of reduced order is used for the representation of the equation.
1746     @raise IllegalCoefficient: if an unknown coefficient keyword is used.
1747     """
1748 jgs 151 super(Helmholtz, self).setValue(**coefficients)
1749 jgs 149
1750     def getCoefficientOfGeneralPDE(self,name):
1751     """
1752     return the value of the coefficient name of the general PDE
1753    
1754     @param name: name of the coefficient requested.
1755     @type name: C{string}
1756     @return: the value of the coefficient name
1757     @rtype: L{Data<escript.Data>}
1758     @raise IllegalCoefficient: if name is not one of coefficients
1759     "A", M{B}, M{C}, M{D}, M{X}, M{Y}, M{d}, M{y}, M{d_contact}, M{y_contact}, M{r} or M{q}.
1760     @note: This method is called by the assembling routine to map the Possion equation onto the general PDE.
1761     """
1762     if name == "A" :
1763     return escript.Data(numarray.identity(self.getDim()),escript.Function(self.getDomain()))*self.getCoefficient("k")
1764     elif name == "B" :
1765     return escript.Data()
1766     elif name == "C" :
1767     return escript.Data()
1768     elif name == "D" :
1769     return self.getCoefficient("omega")
1770     elif name == "X" :
1771     return escript.Data()
1772     elif name == "Y" :
1773     return self.getCoefficient("f")
1774     elif name == "d" :
1775     return self.getCoefficient("alpha")
1776     elif name == "y" :
1777     return self.getCoefficient("g")
1778     elif name == "d_contact" :
1779     return escript.Data()
1780     elif name == "y_contact" :
1781     return escript.Data()
1782     elif name == "r" :
1783     return self.getCoefficient("r")
1784     elif name == "q" :
1785     return self.getCoefficient("q")
1786     else:
1787     raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name
1788    
1789     class LameEquation(LinearPDE):
1790     """
1791     Class to define a Lame equation problem:
1792    
1793     M{-grad(S{mu}*(grad(u[i])[j]+grad(u[j])[i]))[j] - grad(S{lambda}*grad(u[j])[i])[j] = F_i -grad(S{sigma}[i,j])[j] }
1794    
1795     with natural boundary conditons:
1796    
1797     M{n[j]*(S{mu}*(grad(u[i])[j]+grad(u[j])[i]) - S{lambda}*grad(u[j])[i]) = f_i -n[j]*S{sigma}[i,j] }
1798    
1799     and constraints:
1800    
1801     M{u[i]=r[i]} where M{q[i]>0}
1802    
1803     """
1804    
1805     def __init__(self,domain,debug=False):
1806 jgs 151 super(LameEquation, self).__init__(domain,\
1807     domain.getDim(),domain.getDim(),debug)
1808     self.COEFFICIENTS={ "lame_lambda" : PDECoefficient(PDECoefficient.INTERIOR,(),PDECoefficient.OPERATOR),
1809 jgs 149 "lame_mu" : PDECoefficient(PDECoefficient.INTERIOR,(),PDECoefficient.OPERATOR),
1810     "F" : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.BY_EQUATION,),PDECoefficient.RIGHTHANDSIDE),
1811     "sigma" : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.BY_EQUATION,PDECoefficient.BY_DIM),PDECoefficient.RIGHTHANDSIDE),
1812     "f" : PDECoefficient(PDECoefficient.BOUNDARY,(PDECoefficient.BY_EQUATION,),PDECoefficient.RIGHTHANDSIDE),
1813     "r" : PDECoefficient(PDECoefficient.SOLUTION,(PDECoefficient.BY_EQUATION,),PDECoefficient.BOTH),
1814     "q" : PDECoefficient(PDECoefficient.SOLUTION,(PDECoefficient.BY_EQUATION,),PDECoefficient.BOTH)}
1815 jgs 151 self.setSymmetryOn()
1816 jgs 149
1817     def setValue(self,**coefficients):
1818     """
1819     sets new values to coefficients
1820    
1821     @param coefficients: new values assigned to coefficients
1822     @keyword lame_mu: value for coefficient M{S{mu}}
1823     @type lame_mu: any type that can be casted to L{Scalar<escript.Scalar>} object on L{Function<escript.Function>}.
1824     @keyword lame_lambda: value for coefficient M{S{lambda}}
1825     @type lame_lambda: any type that can be casted to L{Scalar<escript.Scalar>} object on L{Function<escript.Function>}.
1826     @keyword F: value for internal force M{F}
1827     @type F: any type that can be casted to L{Vector<escript.Vector>} object on L{Function<escript.Function>}
1828     @keyword sigma: value for initial stress M{S{sigma}}
1829     @type sigma: any type that can be casted to L{Tensor<escript.Tensor>} object on L{Function<escript.Function>}
1830     @keyword f: value for extrenal force M{f}
1831     @type f: any type that can be casted to L{Vector<escript.Vector>} object on L{FunctionOnBoundary<escript.FunctionOnBoundary>}
1832     @keyword r: prescribed values M{r} for the solution in constraints.
1833     @type r: any type that can be casted to L{Vector<escript.Vector>} object on L{Solution<escript.Solution>} or L{ReducedSolution<escript.ReducedSolution>}
1834     depending of reduced order is used for the representation of the equation.
1835     @keyword q: mask for location of constraints
1836     @type q: any type that can be casted to L{Vector<escript.Vector>} object on L{Solution<escript.Solution>} or L{ReducedSolution<escript.ReducedSolution>}
1837     depending of reduced order is used for the representation of the equation.
1838     @raise IllegalCoefficient: if an unknown coefficient keyword is used.
1839     """
1840 jgs 151 super(LameEquation, self).setValue(**coefficients)
1841 jgs 149
1842     def getCoefficientOfGeneralPDE(self,name):
1843     """
1844     return the value of the coefficient name of the general PDE
1845    
1846     @param name: name of the coefficient requested.
1847     @type name: C{string}
1848     @return: the value of the coefficient name
1849     @rtype: L{Data<escript.Data>}
1850     @raise IllegalCoefficient: if name is not one of coefficients
1851     "A", M{B}, M{C}, M{D}, M{X}, M{Y}, M{d}, M{y}, M{d_contact}, M{y_contact}, M{r} or M{q}.
1852     @note: This method is called by the assembling routine to map the Possion equation onto the general PDE.
1853     """
1854     if name == "A" :
1855     out =self.createCoefficientOfGeneralPDE("A")
1856     for i in range(self.getDim()):
1857     for j in range(self.getDim()):
1858     out[i,i,j,j] += self.getCoefficient("lame_lambda")
1859     out[i,j,j,i] += self.getCoefficient("lame_mu")
1860     out[i,j,i,j] += self.getCoefficient("lame_mu")
1861     return out
1862     elif name == "B" :
1863     return escript.Data()
1864     elif name == "C" :
1865     return escript.Data()
1866     elif name == "D" :
1867     return escript.Data()
1868     elif name == "X" :
1869     return self.getCoefficient("sigma")
1870     elif name == "Y" :
1871     return self.getCoefficient("F")
1872     elif name == "d" :
1873     return escript.Data()
1874     elif name == "y" :
1875     return self.getCoefficient("f")
1876     elif name == "d_contact" :
1877     return escript.Data()
1878     elif name == "y_contact" :
1879     return escript.Data()
1880     elif name == "r" :
1881     return self.getCoefficient("r")
1882     elif name == "q" :
1883     return self.getCoefficient("q")
1884     else:
1885     raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name
1886    
1887     class AdvectivePDE(LinearPDE):
1888     """
1889     In cases of PDEs dominated by the advection terms M{B} and M{C} against the adevctive terms M{A}
1890     up-winding has been used. The L{AdvectivePDE} class applies SUPG upwinding to the advective terms.
1891    
1892     In the following we set
1893    
1894     M{Z[j]=C[j]-B[j]}
1895    
1896 jgs 151 or
1897 jgs 149
1898     M{Z[i,k,l]=C[i,k,l]-B[i,l,k]}
1899    
1900 jgs 151 To measure the dominance of the advective terms over the diffusive term M{A} the
1901 jgs 149 X{Pelclet number} M{P} is used. It is defined as
1902 jgs 151
1903 jgs 149 M{P=h|Z|/(2|A|)}
1904    
1905 jgs 151 where M{|.|} denotes the L{length<util.length>} of the arument and M{h} is the local cell size
1906 jgs 149 from L{getSize<escript.Domain.getSize>}. Where M{|A|==0} M{P} is M{S{infinity}}.
1907    
1908     From the X{Pelclet number} the stabilization parameters M{S{Xi}} and M{S{Xi}} are calculated:
1909    
1910     M{S{Xi}=S{xi}(P) h/|Z|}
1911    
1912     where M{S{xi}} is a suitable function of the Peclet number.
1913    
1914     In the case of a single PDE the coefficient are up-dated in the following way:
1915     - M{A[i,j] S{<-} A[i,j] + S{Xi} * Z[j] * Z[l]}
1916     - M{B[j] S{<-} B[j] + S{Xi} * C[j] * D}
1917     - M{C[j] S{<-} C[j] + S{Xi} * B[j] * D}
1918     - M{X[j] S{<-} X[j] + S{Xi} * Z[j] * Y}
1919    
1920     Similar for the case of a systems of PDEs:
1921     - M{A[i,j,k,l] S{<-} A[i,j,k,l]+ S{delta}[p,m] * S{Xi} * Z[p,i,j] * Z[m,k,l]}
1922     - M{B[i,j,k] S{<-} B[i,j,k] + S{delta}[p,m] * S{Xi} * D[p,k] * C[m,i,j]}
1923     - M{C[i,k,l] S{<-} C[i,k,l] + S{delta}[p,m] * S{Xi} * D[p,k] * B[m,l,i]}
1924     - M{X[i,j] S{<-} X[i,j] + S{delta}[p,m] * S{Xi} * Y[p] * Z[m,i,j]}
1925    
1926     where M{S{delta}} is L{kronecker}.
1927     Using upwinding in this form, introduces an additonal error which is proprtional to the cell size M{h}
1928     but with the intension to stabilize the solution.
1929    
1930     """
1931     def __init__(self,domain,numEquations=None,numSolutions=None,xi=None,debug=False):
1932     """
1933 jgs 151 creates a linear, steady, second order PDE on a L{Domain<escript.Domain>}
1934 jgs 149
1935     @param domain: domain of the PDE
1936     @type domain: L{Domain<escript.Domain>}
1937     @param numEquations: number of equations. If numEquations==None the number of equations
1938     is exracted from the PDE coefficients.
1939     @param numSolutions: number of solution components. If numSolutions==None the number of solution components
1940     is exracted from the PDE coefficients.
1941 jgs 151 @param xi: defines a function which returns for any given Preclet number as L{Scalar<escript.Scalar>} object the
1942 jgs 149 M{S{xi}}-value used to define the stabilization parameters. If equal to None, L{ELMAN_RAMAGE} is used.
1943     @type xi: callable object which returns a L{Scalar<escript.Scalar>} object.
1944     @param debug: if True debug informations are printed.
1945     """
1946 jgs 151 super(AdvectivePDE, self).__init__(domain,\
1947     numEquations,numSolutions,debug)
1948 jgs 148 if xi==None:
1949     self.__xi=AdvectivePDE.ELMAN_RAMAGE
1950     else:
1951     self.__xi=xi
1952 jgs 110 self.__Xi=escript.Data()
1953 jgs 102
1954 jgs 151 def setValue(**coefficients):
1955 jgs 149 """
1956     sets new values to coefficients
1957 jgs 102
1958 jgs 149 @param coefficients: new values assigned to coefficients
1959     @keyword A: value for coefficient A.
1960     @type A: any type that can be casted to L{Data<escript.Data>} object on L{Function<escript.Function>}.
1961     @keyword B: value for coefficient B
1962     @type B: any type that can be casted to L{Data<escript.Data>} object on L{Function<escript.Function>}.
1963     @keyword C: value for coefficient C
1964     @type C: any type that can be casted to L{Data<escript.Data>} object on L{Function<escript.Function>}.
1965     @keyword D: value for coefficient D
1966     @type D: any type that can be casted to L{Data<escript.Data>} object on L{Function<escript.Function>}.
1967     @keyword X: value for coefficient X
1968     @type X: any type that can be casted to L{Data<escript.Data>} object on L{Function<escript.Function>}.
1969     @keyword Y: value for coefficient Y
1970     @type Y: any type that can be casted to L{Data<escript.Data>} object on L{Function<escript.Function>}.
1971     @keyword d: value for coefficient d
1972     @type d: any type that can be casted to L{Data<escript.Data>} object on L{FunctionOnBoundary<escript.FunctionOnBoundary>}.
1973     @keyword y: value for coefficient y
1974     @type y: any type that can be casted to L{Data<escript.Data>} object on L{FunctionOnBoundary<escript.FunctionOnBoundary>}.
1975     @keyword d_contact: value for coefficient d_contact
1976     @type d_contact: any type that can be casted to L{Data<escript.Data>} object on L{FunctionOnContactOne<escript.FunctionOnContactOne>}.
1977     or L{FunctionOnContactZero<escript.FunctionOnContactZero>}.
1978     @keyword y_contact: value for coefficient y_contact
1979     @type y_contact: any type that can be casted to L{Data<escript.Data>} object on L{FunctionOnContactOne<escript.FunctionOnContactOne>}.
1980     or L{FunctionOnContactZero<escript.FunctionOnContactZero>}.
1981     @keyword r: values prescribed to the solution at the locations of constraints
1982     @type r: any type that can be casted to L{Data<escript.Data>} object on L{Solution<escript.Solution>} or L{ReducedSolution<escript.ReducedSolution>}
1983     depending of reduced order is used for the solution.
1984     @keyword q: mask for location of constraints
1985     @type q: any type that can be casted to L{Data<escript.Data>} object on L{Solution<escript.Solution>} or L{ReducedSolution<escript.ReducedSolution>}
1986     depending of reduced order is used for the representation of the equation.
1987     @raise IllegalCoefficient: if an unknown coefficient keyword is used.
1988 jgs 148
1989 jgs 149 """
1990     if "A" in coefficients.keys() or "B" in coefficients.keys() or "C" in coefficients.keys(): self.__Xi=escript.Data()
1991 jgs 151 super(AdvectivePDE, self).setValue(**coefficients)
1992 jgs 149
1993     def ELMAN_RAMAGE(self,P):
1994 jgs 151 """
1995     Predefined function to set a values for M{S{xi}} from a Preclet number M{P}.
1996 jgs 149 This function uses the method suggested by H.C. Elman and A. Ramage, I{SIAM J. Numer. Anal.}, B{40} (2002)
1997     - M{S{xi}(P)=0} for M{P<1}
1998     - M{S{xi}(P)=(1-1/P)/2} otherwise
1999    
2000 jgs 151 @param P: Preclet number
2001 jgs 149 @type P: L{Scalar<escript.Scalar>}
2002     @return: up-wind weightimg factor
2003     @rtype: L{Scalar<escript.Scalar>}
2004 jgs 151 """
2005 gross 304 return util.wherePositive(P-1.)*0.5*(1.-1./(P+1.e-15))
2006 jgs 149
2007     def SIMPLIFIED_BROOK_HUGHES(self,P):
2008 jgs 151 """
2009     Predefined function to set a values for M{S{xi}} from a Preclet number M{P}.
2010 jgs 149 The original methods is
2011 jgs 151
2012 jgs 149 M{S{xi}(P)=coth(P)-1/P}
2013    
2014     As the evaluation of M{coth} is expensive we are using the approximation:
2015 jgs 151
2016 jgs 149 - M{S{xi}(P)=P/3} where M{P<3}
2017     - M{S{xi}(P)=1/2} otherwise
2018    
2019 jgs 151 @param P: Preclet number
2020 jgs 149 @type P: L{Scalar<escript.Scalar>}
2021     @return: up-wind weightimg factor
2022     @rtype: L{Scalar<escript.Scalar>}
2023 jgs 151 """
2024 gross 304 c=util.whereNegative(P-3.)
2025 jgs 148 return P/6.*c+1./2.*(1.-c)
2026    
2027 jgs 149 def HALF(self,P):
2028 jgs 151 """
2029     Predefined function to set value M{1/2} for M{S{xi}}
2030    
2031     @param P: Preclet number
2032 jgs 149 @type P: L{Scalar<escript.Scalar>}
2033     @return: up-wind weightimg factor
2034     @rtype: L{Scalar<escript.Scalar>}
2035 jgs 151 """
2036 jgs 149 return escript.Scalar(0.5,P.getFunctionSpace())
2037 jgs 148
2038 gross 328 def __calculateXi(self,peclet_factor,flux,h):
2039     flux=util.Lsup(flux)
2040     if flux_max>0.:
2041     return h*self.__xi(flux*peclet_factor)/(flux+flux_max*self.__TOL)
2042 jgs 149 else:
2043     return 0.
2044    
2045     def __getXi(self):
2046 jgs 110 if self.__Xi.isEmpty():
2047     B=self.getCoefficient("B")
2048     C=self.getCoefficient("C")
2049     A=self.getCoefficient("A")
2050     h=self.getDomain().getSize()
2051     self.__Xi=escript.Scalar(0.,self.getFunctionSpaceForCoefficient("A"))
2052     if not C.isEmpty() or not B.isEmpty():
2053     if not C.isEmpty() and not B.isEmpty():
2054 gross 328 flux2=escript.Scalar(0,self.getFunctionSpaceForCoefficient("A"))
2055 jgs 110 if self.getNumEquations()>1:
2056     if self.getNumSolutions()>1:
2057     for i in range(self.getNumEquations()):
2058     for k in range(self.getNumSolutions()):
2059 gross 328 for l in range(self.getDim()): flux2+=(C[i,k,l]-B[i,l,k])**2
2060     # flux=C-util.reorderComponents(B,[0,2,1])
2061 jgs 110 else:
2062     for i in range(self.getNumEquations()):
2063 gross 328 for l in range(self.getDim()): flux2+=(C[i,l]-B[i,l])**2
2064     # flux=C-B
2065 jgs 110 else:
2066     if self.getNumSolutions()>1:
2067     for k in range(self.getNumSolutions()):
2068 gross 328 for l in range(self.getDim()): flux2+=(C[k,l]-B[l,k])**2
2069     # flux=C-util.reorderComponents(B,[1,0])
2070 jgs 110 else:
2071 gross 328 for l in range(self.getDim()): flux2+=(C[l]-B[l])**2
2072     #flux=C-B
2073     length_of_flux=util.sqrt(flux2)
2074 jgs 110 elif C.isEmpty():
2075 gross 328 length_of_flux=util.length(B)
2076     #flux=B
2077 jgs 110 else:
2078 gross 328 length_of_flux=util.length(C)
2079     #flux=C
2080 jgs 102
2081 gross 328 #length_of_flux=util.length(flux)
2082     flux_max=util.Lsup(length_of_flux)
2083     if flux_max>0.:
2084     # length_of_A=util.inner(flux,util.tensormutiply(A,flux))
2085 jgs 110 length_of_A=util.length(A)
2086     A_max=util.Lsup(length_of_A)
2087     if A_max>0:
2088 jgs 149 inv_A=1./(length_of_A+A_max*self.__TOL)
2089 jgs 110 else:
2090 jgs 149 inv_A=1./self.__TOL
2091 gross 328 peclet_number=length_of_flux*h/2*inv_A
2092 jgs 110 xi=self.__xi(peclet_number)
2093 gross 328 self.__Xi=h*xi/(length_of_flux+flux_max*self.__TOL)
2094 jgs 149 self.trace("preclet number = %e"%util.Lsup(peclet_number))
2095 jgs 110 return self.__Xi
2096 jgs 102
2097 jgs 148
2098     def getCoefficientOfGeneralPDE(self,name):
2099 jgs 108 """
2100 jgs 122 return the value of the coefficient name of the general PDE
2101    
2102 jgs 149 @param name: name of the coefficient requested.
2103     @type name: C{string}
2104     @return: the value of the coefficient name
2105     @rtype: L{Data<escript.Data>}
2106     @raise IllegalCoefficient: if name is not one of coefficients
2107     M{A}, M{B}, M{C}, M{D}, M{X}, M{Y}, M{d}, M{y}, M{d_contact}, M{y_contact}, M{r} or M{q}.
2108     @note: This method is called by the assembling routine to map the Possion equation onto the general PDE.
2109 jgs 108 """
2110 jgs 110 if not self.getNumEquations() == self.getNumSolutions():
2111     raise ValueError,"AdvectivePDE expects the number of solution componets and the number of equations to be equal."
2112    
2113 jgs 148 if name == "A" :
2114 jgs 108 A=self.getCoefficient("A")
2115     B=self.getCoefficient("B")
2116     C=self.getCoefficient("C")
2117 jgs 110 if B.isEmpty() and C.isEmpty():
2118     Aout=A
2119     else:
2120 jgs 148 if A.isEmpty():
2121 jgs 110 Aout=self.createNewCoefficient("A")
2122     else:
2123     Aout=A[:]
2124 jgs 149 Xi=self.__getXi()
2125 jgs 110 if self.getNumEquations()>1:
2126     for i in range(self.getNumEquations()):
2127     for j in range(self.getDim()):
2128 jgs 108 for k in range(self.getNumSolutions()):
2129 jgs 110 for l in range(self.getDim()):
2130     if not C.isEmpty() and not B.isEmpty():
2131 gross 328 # tmp=C-util.reorderComponents(B,[0,2,1])
2132     # Aout=Aout+Xi*util.generalTensorProduct(util.reorder(tmp,[1,2,0]),tmp,offset=1)
2133 jgs 110 for p in range(self.getNumEquations()): Aout[i,j,k,l]+=Xi*(C[p,i,j]-B[p,j,i])*(C[p,k,l]-B[p,l,k])
2134     elif C.isEmpty():
2135     for p in range(self.getNumEquations()): Aout[i,j,k,l]+=Xi*B[p,j,i]*B[p,l,k]
2136 gross 328 # Aout=Aout+Xi*util.generalTensorProduct(util.reorder(B,[2,1,0]),util.reorder(B,[0,2,1]),offset=1)
2137 jgs 110 else:
2138     for p in range(self.getNumEquations()): Aout[i,j,k,l]+=Xi*C[p,i,j]*C[p,k,l]
2139 gross 328 # Aout=Aout+Xi*util.generalTensorProduct(util.reorder(C,[1,2,0]),C,offset=1)
2140 jgs 110 else:
2141     for j in range(self.getDim()):
2142     for l in range(self.getDim()):
2143     if not C.isEmpty() and not B.isEmpty():
2144     Aout[j,l]+=Xi*(C[j]-B[j])*(C[l]-B[l])
2145     elif C.isEmpty():
2146     Aout[j,l]+=Xi*B[j]*B[l]
2147     else:
2148     Aout[j,l]+=Xi*C[j]*C[l]
2149 gross 328 # if not C.isEmpty() and not B.isEmpty():
2150     # tmp=C-B
2151     # Aout=Aout+Xi*util.outer(tmp,tmp)
2152     # elif C.isEmpty():
2153     # Aout=Aout+Xi*util.outer(B,B)
2154     # else:
2155     # Aout=Aout+Xi*util.outer(C,C)
2156 jgs 110 return Aout
2157 jgs 148 elif name == "B" :
2158 jgs 110 B=self.getCoefficient("B")
2159     C=self.getCoefficient("C")
2160     D=self.getCoefficient("D")
2161     if C.isEmpty() or D.isEmpty():
2162     Bout=B
2163     else:
2164 jgs 149 Xi=self.__getXi()
2165 jgs 148 if B.isEmpty():
2166 jgs 110 Bout=self.createNewCoefficient("B")
2167     else:
2168     Bout=B[:]
2169     if self.getNumEquations()>1:
2170     for k in range(self.getNumSolutions()):
2171 jgs 148 for p in range(self.getNumEquations()):
2172 jgs 110 tmp=Xi*D[p,k]
2173     for i in range(self.getNumEquations()):
2174     for j in range(self.getDim()):
2175     Bout[i,j,k]+=tmp*C[p,i,j]
2176 gross 328 # Bout=Bout+Xi*util.generalTensorProduct(util.reorder(C,[1,2,0]),D,offset=1)
2177 jgs 110 else:
2178     tmp=Xi*D
2179     for j in range(self.getDim()): Bout[j]+=tmp*C[j]
2180 gross 328 # Bout=Bout+Xi*D*C
2181 jgs 110 return Bout
2182 jgs 148 elif name == "C" :
2183 jgs 110 B=self.getCoefficient("B")
2184     C=self.getCoefficient("C")
2185     D=self.getCoefficient("D")
2186     if B.isEmpty() or D.isEmpty():
2187     Cout=C
2188     else:
2189 jgs 149 Xi=self.__getXi()
2190 jgs 148 if C.isEmpty():
2191 jgs 110 Cout=self.createNewCoefficient("C")
2192     else:
2193     Cout=C[:]
2194     if self.getNumEquations()>1:
2195     for k in range(self.getNumSolutions()):
2196     for p in range(self.getNumEquations()):
2197     tmp=Xi*D[p,k]
2198     for i in range(self.getNumEquations()):
2199     for l in range(self.getDim()):
2200     Cout[i,k,l]+=tmp*B[p,l,i]
2201 gross 328 # Cout=Cout+Xi*B[p,l,i]*D[p,k]
2202 jgs 110 else:
2203     tmp=Xi*D
2204     for j in range(self.getDim()): Cout[j]+=tmp*B[j]
2205 gross 328 # Cout=Cout+tmp*D*B
2206 jgs 110 return Cout
2207 jgs 148 elif name == "D" :
2208 jgs 108 return self.getCoefficient("D")
2209 jgs 148 elif name == "X" :
2210 jgs 110 X=self.getCoefficient("X")
2211     Y=self.getCoefficient("Y")
2212     B=self.getCoefficient("B")
2213     C=self.getCoefficient("C")
2214     if Y.isEmpty() or (B.isEmpty() and C.isEmpty()):
2215     Xout=X
2216     else:
2217     if X.isEmpty():
2218     Xout=self.createNewCoefficient("X")
2219     else:
2220     Xout=X[:]
2221 jgs 149 Xi=self.__getXi()
2222 jgs 110 if self.getNumEquations()>1:
2223 jgs 148 for p in range(self.getNumEquations()):
2224 jgs 110 tmp=Xi*Y[p]
2225     for i in range(self.getNumEquations()):
2226     for j in range(self.getDim()):
2227     if not C.isEmpty() and not B.isEmpty():
2228     Xout[i,j]+=tmp*(C[p,i,j]-B[p,j,i])
2229 gross 328 # Xout=X_out+Xi*util.inner(Y,C-util.reorderComponents(B,[0,2,1]),offset=1)
2230 jgs 110 elif C.isEmpty():
2231     Xout[i,j]-=tmp*B[p,j,i]
2232 gross 328 # Xout=X_out-Xi*util.inner(Y,util.reorderComponents(B,[0,2,1]),offset=1)
2233 jgs 110 else:
2234     Xout[i,j]+=tmp*C[p,i,j]
2235 gross 328 # Xout=X_out+Xi*util.inner(Y,C,offset=1)
2236 jgs 110 else:
2237     tmp=Xi*Y
2238     for j in range(self.getDim()):
2239     if not C.isEmpty() and not B.isEmpty():
2240     Xout[j]+=tmp*(C[j]-B[j])
2241 gross 328 # Xout=Xout+Xi*Y*(C-B)
2242 jgs 110 elif C.isEmpty():
2243     Xout[j]-=tmp*B[j]
2244 gross 328 # Xout=Xout-Xi*Y*B
2245 jgs 110 else:
2246     Xout[j]+=tmp*C[j]
2247 gross 328 # Xout=Xout+Xi*Y*C
2248 jgs 110 return Xout
2249 jgs 148 elif name == "Y" :
2250 jgs 108 return self.getCoefficient("Y")
2251 jgs 148 elif name == "d" :
2252 jgs 108 return self.getCoefficient("d")
2253 jgs 148 elif name == "y" :
2254 jgs 108 return self.getCoefficient("y")
2255 jgs 148 elif name == "d_contact" :
2256 jgs 108 return self.getCoefficient("d_contact")
2257     elif name == "y_contact" :
2258     return self.getCoefficient("y_contact")
2259 jgs 148 elif name == "r" :
2260 jgs 108 return self.getCoefficient("r")
2261 jgs 148 elif name == "q" :
2262 jgs 108 return self.getCoefficient("q")
2263     else:
2264 jgs 149 raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name
2265 jgs 108
2266 gross 345 class AdvectionDiffusion(LinearPDE):
2267     """
2268     Class to define PDE equation of the unisotropic advection-diffusion problem, which is genear L{LinearPDE} of the form
2269 jgs 108
2270 gross 348 M{S{omega}*u + inner(v,grad(u))- grad(matrixmult(k_bar,grad(u))[j])[j] = f}
2271 gross 345
2272     with natural boundary conditons
2273    
2274     M{n[j]*matrixmult(k,grad(u))[j] = g- S{alpha}u }
2275    
2276     and constraints:
2277    
2278     M{u=r} where M{q>0}
2279    
2280 gross 348 and
2281 gross 345
2282 gross 348 M{k_bar[i,j]=k[i,j]+upwind[i]*upwind[j]}
2283    
2284 gross 345 """
2285    
2286     def __init__(self,domain,debug=False):
2287     """
2288     initializes a new Poisson equation
2289    
2290     @param domain: domain of the PDE
2291     @type domain: L{Domain<escript.Domain>}
2292     @param debug: if True debug informations are printed.
2293    
2294     """
2295 gross 349 super(AdvectionDiffusion, self).__init__(domain,1,1,debug)
2296 gross 345 self.COEFFICIENTS={"omega": PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.BY_EQUATION,),PDECoefficient.OPERATOR),
2297     "k": PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.BY_DIM,PDECoefficient.BY_DIM),PDECoefficient.OPERATOR),
2298     "f": PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.BY_EQUATION,),PDECoefficient.RIGHTHANDSIDE),
2299     "v": PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.BY_DIM,),PDECoefficient.OPERATOR),
2300 gross 348 "upwind": PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.BY_DIM,),PDECoefficient.OPERATOR),
2301 gross 345 "alpha": PDECoefficient(PDECoefficient.BOUNDARY,(PDECoefficient.BY_EQUATION,),PDECoefficient.OPERATOR),
2302     "g": PDECoefficient(PDECoefficient.BOUNDARY,(PDECoefficient.BY_EQUATION,),PDECoefficient.RIGHTHANDSIDE),
2303     "r": PDECoefficient(PDECoefficient.SOLUTION,(PDECoefficient.BY_EQUATION,),PDECoefficient.BOTH),
2304     "q": PDECoefficient(PDECoefficient.SOLUTION,(PDECoefficient.BY_EQUATION,),PDECoefficient.BOTH)}
2305    
2306     def setValue(self,**coefficients):
2307     """
2308     sets new values to coefficients
2309    
2310     @param coefficients: new values assigned to coefficients
2311     @keyword omega: value for coefficient M{S{omega}}
2312     @type omega: any type that can be casted to L{Scalar<escript.Scalar>} object on L{Function<escript.Function>}.
2313     @keyword k: value for coefficient M{k}
2314     @type k: any type that can be casted to L{Tensor<escript.Tensor>} object on L{Function<escript.Function>}.
2315     @keyword v: value for coefficient M{v}
2316     @type v: any type that can be casted to L{Vector<escript.Vector>} object on L{Function<escript.Function>}.
2317 gross 348 @keyword upwind: value for upwind term M{upwind}
2318     @type upwind: any type that can be casted to L{Vector<escript.Vector>} object on L{Function<escript.Function>}.
2319 gross 345 @keyword f: value for right hand side M{f}
2320     @type f: any type that can be casted to L{Scalar<escript.Scalar>} object on L{Function<escript.Function>}.
2321     @keyword alpha: value for right hand side M{S{alpha}}
2322     @type alpha: any type that can be casted to L{Scalar<escript.Scalar>} object on L{FunctionOnBoundary<escript.FunctionOnBoundary>}.
2323     @keyword g: value for right hand side M{g}
2324     @type g: any type that can be casted to L{Scalar<escript.Scalar>} object on L{FunctionOnBoundary<escript.FunctionOnBoundary>}.
2325     @keyword r: prescribed values M{r} for the solution in constraints.
2326     @type r: any type that can be casted to L{Scalar<escript.Scalar>} object on L{Solution<escript.Solution>} or L{ReducedSolution<escript.ReducedSolution>}
2327     depending of reduced order is used for the representation of the equation.
2328     @keyword q: mask for location of constraints
2329     @type q: any type that can be casted to L{Scalar<escript.Scalar>} object on L{Solution<escript.Solution>} or L{ReducedSolution<escript.ReducedSolution>}
2330     depending of reduced order is used for the representation of the equation.
2331     @raise IllegalCoefficient: if an unknown coefficient keyword is used.
2332     """
2333 gross 349 super(AdvectionDiffusion, self).setValue(**coefficients)
2334 gross 345
2335     def getCoefficientOfGeneralPDE(self,name):
2336     """
2337     return the value of the coefficient name of the general PDE
2338    
2339     @param name: name of the coefficient requested.
2340     @type name: C{string}
2341     @return: the value of the coefficient name
2342     @rtype: L{Data<escript.Data>}
2343     @raise IllegalCoefficient: if name is not one of coefficients
2344     "A", M{B}, M{C}, M{D}, M{X}, M{Y}, M{d}, M{y}, M{d_contact}, M{y_contact}, M{r} or M{q}.
2345     @note: This method is called by the assembling routine to map the Possion equation onto the general PDE.
2346     """
2347     if name == "A" :
2348 gross 410 return self.getCoefficient("k")+util.outer(self.getCoefficient("upwind"),self.getCoefficient("upwind"))
2349 gross 345 elif name == "B" :
2350     return escript.Data()
2351     elif name == "C" :
2352 gross 410 return self.getCoefficient("v")
2353 gross 345 elif name == "D" :
2354     return self.getCoefficient(