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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 104 - (show annotations)
Fri Dec 17 07:43:12 2004 UTC (15 years ago) by jgs
Original Path: trunk/esys2/escript/py_src/linearPDEs.py
File MIME type: text/x-python
File size: 36737 byte(s)
*** empty log message ***

1 # $Id$
2
3 ## @file linearPDEs.py
4
5 """
6 @brief Functions and classes for linear PDEs
7 """
8
9 import escript
10 import util
11 import numarray
12
13 def identifyDomain(domain=None,data={}):
14 """
15 @brief Return the Domain which is equal to the input domain (if not None)
16 and is the domain of all Data objects in the dictionary data.
17 An exception is raised if this is not possible
18
19 @param domain
20 @param data
21 """
22 # get the domain used by any Data object in the list data:
23 data_domain=None
24 for d in data.itervalues():
25 if isinstance(d,escript.Data):
26 if not d.isEmpty(): data_domain=d.getDomain()
27 # check if domain and data_domain are identical?
28 if domain == None:
29 if data_domain == None:
30 raise ValueError,"Undefined PDE domain. Specify a domain or use a Data class object as coefficient"
31 else:
32 if data_domain == None:
33 data_domain=domain
34 else:
35 if not data_domain == domain:
36 raise ValueError,"Domain of coefficients doesnot match specified domain"
37 # now we check if all Data class object coefficients are defined on data_domain:
38 for i,d in data.iteritems():
39 if isinstance(d,escript.Data):
40 if not d.isEmpty():
41 if not data_domain==d.getDomain():
42 raise ValueError,"Illegal domain for coefficient %s."%i
43 # done:
44 return data_domain
45
46 def identifyNumEquationsAndSolutions(dim,coef={}):
47 # get number of equations and number of unknowns:
48 numEquations=0
49 numSolutions=0
50 for i in coef.iterkeys():
51 if not coef[i].isEmpty():
52 res=_PDECoefficientTypes[i].estimateNumEquationsAndNumSolutions(coef[i].getShape(),dim)
53 if res==None:
54 raise ValueError,"Illegal shape %s of coefficient %s"%(coef[i].getShape().__str__(),i)
55 else:
56 numEquations=max(numEquations,res[0])
57 numSolutions=max(numSolutions,res[1])
58 return numEquations,numSolutions
59
60
61 def _CompTuple2(t1,t2):
62 """
63 @brief
64
65 @param t1
66 @param t2
67 """
68 dif=t1[0]+t1[1]-(t2[0]+t2[1])
69 if dif<0: return 1
70 elif dif>0: return -1
71 else: return 0
72
73 class PDECoefficientType:
74 """
75 @brief
76 """
77 # identifier for location of Data objects defining coefficients
78 INTERIOR=0
79 BOUNDARY=1
80 CONTACT=2
81 CONTINUOUS=3
82 # identifier in the pattern of coefficients:
83 # the pattern is a tuple of EQUATION,SOLUTION,DIM where DIM represents the spatial dimension, EQUATION the number of equations and SOLUTION the
84 # number of unknowns.
85 EQUATION=3
86 SOLUTION=4
87 DIM=5
88 # indicator for what is altered if the coefficient is altered:
89 OPERATOR=5
90 RIGHTHANDSIDE=6
91 BOTH=7
92 def __init__(self,where,pattern,altering):
93 """
94 @brief Initialise a PDE Coefficient type
95 """
96 self.what=where
97 self.pattern=pattern
98 self.altering=altering
99
100 def getFunctionSpace(self,domain):
101 """
102 @brief defines the FunctionSpace of the coefficient on the domain
103
104 @param domain
105 """
106 if self.what==self.INTERIOR: return escript.Function(domain)
107 elif self.what==self.BOUNDARY: return escript.FunctionOnBoundary(domain)
108 elif self.what==self.CONTACT: return escript.FunctionOnContactZero(domain)
109 elif self.what==self.CONTINUOUS: return escript.ContinuousFunction(domain)
110
111 def isAlteringOperator(self):
112 """
113 @brief return true if the operator of the PDE is changed when the coefficient is changed
114 """
115 if self.altering==self.OPERATOR or self.altering==self.BOTH:
116 return not None
117 else:
118 return None
119
120 def isAlteringRightHandSide(self):
121 """
122 @brief return true if the right hand side of the PDE is changed when the coefficient is changed
123 """
124 if self.altering==self.RIGHTHANDSIDE or self.altering==self.BOTH:
125 return not None
126 else:
127 return None
128
129 def estimateNumEquationsAndNumSolutions(self,shape=(),dim=3):
130 """
131 @brief tries to estimate the number of equations in a given tensor shape for a given spatial dimension dim
132
133 @param shape
134 @param dim
135 """
136 if len(shape)>0:
137 num=max(shape)+1
138 else:
139 num=1
140 search=[]
141 for u in range(num):
142 for e in range(num):
143 search.append((e,u))
144 search.sort(_CompTuple2)
145 for item in search:
146 s=self.buildShape(item[0],item[1],dim)
147 if len(s)==0 and len(shape)==0:
148 return (1,1)
149 else:
150 if s==shape: return item
151 return None
152
153 def buildShape(self,e=1,u=1,dim=3):
154 """
155 @brief builds the required shape for a given number of equations e, number of unknowns u and spatial dimension dim
156
157 @param e
158 @param u
159 @param dim
160 """
161 s=()
162 for i in self.pattern:
163 if i==self.EQUATION:
164 if e>1: s=s+(e,)
165 elif i==self.SOLUTION:
166 if u>1: s=s+(u,)
167 else:
168 s=s+(dim,)
169 return s
170
171 _PDECoefficientTypes={
172 "A" : PDECoefficientType(PDECoefficientType.INTERIOR,(PDECoefficientType.EQUATION,PDECoefficientType.DIM,PDECoefficientType.SOLUTION,PDECoefficientType.DIM),PDECoefficientType.OPERATOR),
173 "B" : PDECoefficientType(PDECoefficientType.INTERIOR,(PDECoefficientType.EQUATION,PDECoefficientType.DIM,PDECoefficientType.SOLUTION),PDECoefficientType.OPERATOR),
174 "C" : PDECoefficientType(PDECoefficientType.INTERIOR,(PDECoefficientType.EQUATION,PDECoefficientType.SOLUTION,PDECoefficientType.DIM),PDECoefficientType.OPERATOR),
175 "D" : PDECoefficientType(PDECoefficientType.INTERIOR,(PDECoefficientType.EQUATION,PDECoefficientType.SOLUTION),PDECoefficientType.OPERATOR),
176 "X" : PDECoefficientType(PDECoefficientType.INTERIOR,(PDECoefficientType.EQUATION,PDECoefficientType.DIM),PDECoefficientType.RIGHTHANDSIDE),
177 "Y" : PDECoefficientType(PDECoefficientType.INTERIOR,(PDECoefficientType.EQUATION,),PDECoefficientType.RIGHTHANDSIDE),
178 "d" : PDECoefficientType(PDECoefficientType.BOUNDARY,(PDECoefficientType.EQUATION,PDECoefficientType.SOLUTION),PDECoefficientType.OPERATOR),
179 "y" : PDECoefficientType(PDECoefficientType.BOUNDARY,(PDECoefficientType.EQUATION,),PDECoefficientType.RIGHTHANDSIDE),
180 "d_contact" : PDECoefficientType(PDECoefficientType.CONTACT,(PDECoefficientType.EQUATION,PDECoefficientType.SOLUTION),PDECoefficientType.OPERATOR),
181 "y_contact" : PDECoefficientType(PDECoefficientType.CONTACT,(PDECoefficientType.EQUATION,),PDECoefficientType.RIGHTHANDSIDE),
182 "r" : PDECoefficientType(PDECoefficientType.CONTINUOUS,(PDECoefficientType.EQUATION,),PDECoefficientType.RIGHTHANDSIDE),
183 "q" : PDECoefficientType(PDECoefficientType.CONTINUOUS,(PDECoefficientType.SOLUTION,),PDECoefficientType.BOTH),
184 }
185
186 class LinearPDE:
187 """
188 @brief Class to define a linear PDE
189
190 class to define a linear PDE of the form
191
192 -(A_{ijkl}u_{k,l})_{,j} -(B_{ijk}u_k)_{,j} + C_{ikl}u_{k,l} +D_{ik}u_k = - (X_{ij})_{,j} + Y_i
193
194 with boundary conditons:
195
196 n_j*(A_{ijkl}u_{k,l}+B_{ijk}u_k)_{,j} + d_{ik}u_k = - n_j*X_{ij} + y_i
197
198 and contact conditions
199
200 n_j*(A_{ijkl}u_{k,l}+B_{ijk}u_k)_{,j} + d_contact_{ik}[u_k] = - n_j*X_{ij} + y_contact_i
201
202 and constraints:
203
204 u_i=r_i where q_i>0
205
206 """
207 DEFAULT_METHOD=util.DEFAULT_METHOD
208 DIRECT=util.DIRECT
209 CHOLEVSKY=util.CHOLEVSKY
210 PCG=util.PCG
211 CR=util.CR
212 CGS=util.CGS
213 BICGSTAB=util.BICGSTAB
214 SSOR=util.SSOR
215 GMRES=util.GMRES
216 PRES20=util.PRES20
217
218 def __init__(self,domain,numEquations=0,numSolutions=0):
219 """
220 @brief initializes a new linear PDE.
221
222 @param args
223 """
224
225 # initialize attributes
226 self.__debug=None
227 self.__domain=domain
228 self.__numEquations=numEquations
229 self.__numSolutions=numSolutions
230 self.cleanCoefficients()
231
232 self.__operator=escript.Operator()
233 self.__operator_isValid=False
234 self.__righthandside=escript.Data()
235 self.__righthandside_isValid=False
236 self.__solution=escript.Data()
237 self.__solution_isValid=False
238
239 # set some default values:
240 self.__homogeneous_constraint=True
241 self.__row_function_space=escript.Solution(self.__domain)
242 self.__column_function_space=escript.Solution(self.__domain)
243 self.__tolerance=1.e-8
244 self.__solver_method=util.DEFAULT_METHOD
245 self.__matrix_type=self.__domain.getSystemMatrixTypeId(util.DEFAULT_METHOD,False)
246 self.__sym=False
247 self.__lumping=False
248
249 def getCoefficient(self,name):
250 """
251 @brief return the value of the coefficient name
252
253 @param name
254 """
255 return self.__coefficient[name]
256
257 def setValue(self,**coefficients):
258 """
259 @brief sets new values to coefficients
260
261 @param coefficients
262 """
263 self._setValue(**coefficients)
264
265
266 def _setValue(self,**coefficients):
267 """
268 @brief sets new values to coefficients
269
270 @param coefficients
271 """
272
273 # get the dictionary of the coefficinets been altered:
274 alteredCoefficients={}
275 for i,d in coefficients.iteritems():
276 if self.hasCoefficient(i):
277 if d == None:
278 alteredCoefficients[i]=escript.Data()
279 elif isinstance(d,escript.Data):
280 if d.isEmpty():
281 alteredCoefficients[i]=escript.Data()
282 else:
283 alteredCoefficients[i]=escript.Data(d,self.getFunctionSpaceOfCoefficient(i))
284 else:
285 alteredCoefficients[i]=escript.Data(d,self.getFunctionSpaceOfCoefficient(i))
286 else:
287 raise ValueError,"Attempt to set undefined coefficient %s"%i
288 # if numEquations and numSolutions is undefined we try identify their values based on the coefficients:
289 if self.__numEquations<1 or self.__numSolutions<1:
290 numEquations,numSolutions=identifyNumEquationsAndSolutions(self.getDomain().getDim(),alteredCoefficients)
291 if self.__numEquations<1 and numEquations>0: self.__numEquations=numEquations
292 if self.__numSolutions<1 and numSolutions>0: self.__numSolutions=numSolutions
293 if self.debug() and self.__numEquations>0: print "PDE Debug: identified number of equations is ",self.__numEquations
294 if self.debug() and self.__numSolutions>0: print "PDE Debug: identified number of solutions is ",self.__numSolutions
295
296 # now we check the shape of the coefficient if numEquations and numSolutions are set:
297 if self.__numEquations>0 and self.__numSolutions>0:
298 for i in self.__coefficient.iterkeys():
299 if alteredCoefficients.has_key(i) and not alteredCoefficients[i].isEmpty():
300 if not self.getShapeOfCoefficient(i)==alteredCoefficients[i].getShape():
301 raise ValueError,"Expected shape for coefficient %s is %s but actual shape is %s."%(i,self.getShapeOfCoefficient(i),alteredCoefficients[i].getShape())
302 else:
303 if not self.__coefficient[i].isEmpty():
304 if not self.getShapeOfCoefficient(i)==self.__coefficient[i].getShape():
305 raise ValueError,"Expected shape for coefficient %s is %s but actual shape is %s."%(i,self.getShapeOfCoefficient(i),self.__coefficient[i].getShape())
306 # overwrite new values:
307 for i,d in alteredCoefficients.iteritems():
308 if self.debug(): print "PDE Debug: Coefficient %s has been altered."%i
309 self.__coefficient[i]=d
310 self.alteredCoefficient(i)
311
312 # reset the HomogeneousConstraintFlag:
313 self.__setHomogeneousConstraintFlag()
314 if not "q" in alteredCoefficients and not self.__homogeneous_constraint: self.__rebuildSystem()
315
316 def cleanCoefficients(self):
317 """
318 @brief resets all coefficients to default values.
319 """
320 self.__coefficient={}
321 for i in _PDECoefficientTypes.iterkeys():
322 self.__coefficient[i]=escript.Data()
323
324 def getShapeOfCoefficient(self,name):
325 """
326 @brief return the shape of the coefficient name
327
328 @param name
329 """
330 if self.hasCoefficient(name):
331 return _PDECoefficientTypes[name].buildShape(self.getNumEquations(),self.getNumSolutions(),self.getDomain().getDim())
332 else:
333 raise ValueError,"Solution coefficient %s requested"%name
334
335 def getFunctionSpaceOfCoefficient(self,name):
336 """
337 @brief return the atoms of the coefficient name
338
339 @param name
340 """
341 if self.hasCoefficient(name):
342 return _PDECoefficientTypes[name].getFunctionSpace(self.getDomain())
343 else:
344 raise ValueError,"Solution coefficient %s requested"%name
345
346 def alteredCoefficient(self,name):
347 """
348 @brief annonced that coefficient name has been changed
349
350 @param name
351 """
352 if self.hasCoefficient(name):
353 if _PDECoefficientTypes[name].isAlteringOperator(): self.__rebuildOperator()
354 if _PDECoefficientTypes[name].isAlteringRightHandSide(): self.__rebuildRightHandSide()
355 else:
356 raise ValueError,"Solution coefficient %s requested"%name
357
358 def __setHomogeneousConstraintFlag(self):
359 """
360 @brief checks if the constraints are homogeneous and sets self.__homogeneous_constraint accordingly.
361 """
362 self.__homogeneous_constraint=True
363 q=self.getCoefficient("q")
364 r=self.getCoefficient("r")
365 if not q.isEmpty() and not r.isEmpty():
366 print (q*r).Lsup(), 1.e-13*r.Lsup()
367 if (q*r).Lsup()>=1.e-13*r.Lsup(): self.__homogeneous_constraint=False
368 if self.debug():
369 if self.__homogeneous_constraint:
370 print "PDE Debug: Constraints are homogeneous."
371 else:
372 print "PDE Debug: Constraints are inhomogeneous."
373
374
375 def hasCoefficient(self,name):
376 """
377 @brief return true if name is the name of a coefficient
378
379 @param name
380 """
381 return self.__coefficient.has_key(name)
382
383 def getFunctionSpaceForEquation(self):
384 """
385 @brief return true if the test functions should use reduced order
386 """
387 return self.__row_function_space
388
389 def getFunctionSpaceForSolution(self):
390 """
391 @brief return true if the interpolation of the solution should use reduced order
392 """
393 return self.__column_function_space
394
395 # ===== debug ==============================================================
396 def setDebugOn(self):
397 """
398 @brief
399 """
400 self.__debug=not None
401
402 def setDebugOff(self):
403 """
404 @brief
405 """
406 self.__debug=None
407
408 def debug(self):
409 """
410 @brief returns true if the PDE is in the debug mode
411 """
412 return self.__debug
413
414 #===== Lumping ===========================
415 def setLumpingOn(self):
416 """
417 @brief indicates to use matrix lumping
418 """
419 if not self.isUsingLumping():
420 raise SystemError,"Lumping is not working yet! Talk to the experts"
421 if self.debug() : print "PDE Debug: lumping is set on"
422 self.__rebuildOperator()
423 self.__lumping=True
424
425 def setLumpingOff(self):
426 """
427 @brief switches off matrix lumping
428 """
429 if self.isUsingLumping():
430 if self.debug() : print "PDE Debug: lumping is set off"
431 self.__rebuildOperator()
432 self.__lumping=False
433
434 def setLumping(self,flag=False):
435 """
436 @brief set the matrix lumping flag to flag
437 """
438 if flag:
439 self.setLumpingOn()
440 else:
441 self.setLumpingOff()
442
443 def isUsingLumping(self):
444 """
445 @brief
446 """
447 return self.__lumping
448
449 #============ method business =========================================================
450 def setSolverMethod(self,solver=util.DEFAULT_METHOD):
451 """
452 @brief sets a new solver
453 """
454 if not solver==self.getSolverMethod():
455 self.__solver_method=solver
456 if self.debug() : print "PDE Debug: New solver is %s"%solver
457 self.__checkMatrixType()
458
459 def getSolverMethod(self):
460 """
461 @brief returns the solver method
462 """
463 return self.__solver_method
464
465 #============ tolerance business =========================================================
466 def setTolerance(self,tol=1.e-8):
467 """
468 @brief resets the tolerance to tol.
469 """
470 if not tol>0:
471 raise ValueException,"Tolerance as to be positive"
472 if tol<self.getTolerance(): self.__rebuildSolution()
473 if self.debug() : print "PDE Debug: New tolerance %e",tol
474 self.__tolerance=tol
475 return
476 def getTolerance(self):
477 """
478 @brief returns the tolerance set for the solution
479 """
480 return self.__tolerance
481
482 #===== symmetry flag ==========================
483 def isSymmetric(self):
484 """
485 @brief returns true is the operator is considered to be symmetric
486 """
487 return self.__sym
488
489 def setSymmetryOn(self):
490 """
491 @brief sets the symmetry flag to true
492 """
493 if not self.isSymmetric():
494 if self.debug() : print "PDE Debug: Operator is set to be symmetric"
495 self.__sym=True
496 self.__checkMatrixType()
497
498 def setSymmetryOff(self):
499 """
500 @brief sets the symmetry flag to false
501 """
502 if self.isSymmetric():
503 if self.debug() : print "PDE Debug: Operator is set to be unsymmetric"
504 self.__sym=False
505 self.__checkMatrixType()
506
507 def setSymmetryTo(self,flag=False):
508 """
509 @brief sets the symmetry flag to flag
510
511 @param flag
512 """
513 if flag:
514 self.setSymmetryOn()
515 else:
516 self.setSymmetryOff()
517
518 #===== order reduction ==========================
519 def setReducedOrderOn(self):
520 """
521 @brief switches to on reduced order
522 """
523 self.setReducedOrderForSolutionOn()
524 self.setReducedOrderForEquationOn()
525
526 def setReducedOrderOff(self):
527 """
528 @brief switches to full order
529 """
530 self.setReducedOrderForSolutionOff()
531 self.setReducedOrderForEquationOff()
532
533 def setReducedOrderTo(self,flag=False):
534 """
535 @brief sets order according to flag
536
537 @param flag
538 """
539 self.setReducedOrderForSolutionTo(flag)
540 self.setReducedOrderForEquationTo(flag)
541
542
543 #===== order reduction solution ==========================
544 def setReducedOrderForSolutionOn(self):
545 """
546 @brief switches to reduced order to interpolate solution
547 """
548 new_fs=escript.ReducedSolution(self.getDomain())
549 if self.getFunctionSpaceForSolution()!=new_fs:
550 if self.debug() : print "PDE Debug: Reduced order is used to interpolate solution."
551 self.__column_function_space=new_fs
552 self.__rebuildSystem(deep=True)
553
554 def setReducedOrderForSolutionOff(self):
555 """
556 @brief switches to full order to interpolate solution
557 """
558 new_fs=escript.Solution(self.getDomain())
559 if self.getFunctionSpaceForSolution()!=new_fs:
560 if self.debug() : print "PDE Debug: Full order is used to interpolate solution."
561 self.__column_function_space=new_fs
562 self.__rebuildSystem(deep=True)
563
564 def setReducedOrderForSolutionTo(self,flag=False):
565 """
566 @brief sets order for test functions according to flag
567
568 @param flag
569 """
570 if flag:
571 self.setReducedOrderForSolutionOn()
572 else:
573 self.setReducedOrderForSolutionOff()
574
575 #===== order reduction equation ==========================
576 def setReducedOrderForEquationOn(self):
577 """
578 @brief switches to reduced order for test functions
579 """
580 new_fs=escript.ReducedSolution(self.getDomain())
581 if self.getFunctionSpaceForEquation()!=new_fs:
582 if self.debug() : print "PDE Debug: Reduced order is used for test functions."
583 self.__row_function_space=new_fs
584 self.__rebuildSystem(deep=True)
585
586 def setReducedOrderForEquationOff(self):
587 """
588 @brief switches to full order for test functions
589 """
590 new_fs=escript.Solution(self.getDomain())
591 if self.getFunctionSpaceForEquation()!=new_fs:
592 if self.debug() : print "PDE Debug: Full order is used for test functions."
593 self.__row_function_space=new_fs
594 self.__rebuildSystem(deep=True)
595
596 def setReducedOrderForEquationTo(self,flag=False):
597 """
598 @brief sets order for test functions according to flag
599
600 @param flag
601 """
602 if flag:
603 self.setReducedOrderForEquationOn()
604 else:
605 self.setReducedOrderForEquationOff()
606
607
608 # ==== initialization =====================================================================
609 def __makeNewOperator(self):
610 """
611 @brief
612 """
613 return self.getDomain().newOperator( \
614 self.getNumEquations(), \
615 self.getFunctionSpaceForEquation(), \
616 self.getNumSolutions(), \
617 self.getFunctionSpaceForSolution(), \
618 self.__matrix_type)
619
620 def __makeNewRightHandSide(self):
621 """
622 @brief
623 """
624 return escript.Data(0.,(self.getNumEquations(),),self.getFunctionSpaceForEquation(),True)
625
626 def __makeNewSolution(self):
627 """
628 @brief
629 """
630 return escript.Data(0.,(self.getNumSolutions(),),self.getFunctionSpaceForSolution(),True)
631
632 def __getFreshOperator(self):
633 """
634 @brief
635 """
636 if self.__operator.isEmpty():
637 self.__operator=self.__makeNewOperator()
638 if self.debug() : print "PDE Debug: New operator allocated"
639 else:
640 self.__operator.setValue(0.)
641 if self.debug() : print "PDE Debug: Operator reset to zero"
642 return self.__operator
643
644 def __getFreshRightHandSide(self):
645 """
646 @brief
647 """
648 if self.__righthandside.isEmpty():
649 self.__righthandside=self.__makeNewRightHandSide()
650 if self.debug() : print "PDE Debug: New right hand side allocated"
651 else:
652 print "fix self.__righthandside*=0"
653 self.__righthandside*=0.
654 if self.debug() : print "PDE Debug: Right hand side reset to zero"
655 return self.__righthandside
656
657 # ==== rebuild switches =====================================================================
658 def __rebuildSolution(self,deep=False):
659 """
660 @brief indicates the PDE has to be reolved if the solution is requested
661 """
662 if self.__solution_isValid and self.debug() : print "PDE Debug: PDE has to be resolved."
663 self.__solution_isValid=False
664 if deep: self.__solution=escript.Data(deep)
665
666
667 def __rebuildOperator(self,deep=False):
668 """
669 @brief indicates the operator has to be rebuilt next time it is used
670 """
671 if self.__operator_isValid and self.debug() : print "PDE Debug: Operator has to be rebuilt."
672 self.__rebuildSolution(deep)
673 self.__operator_isValid=False
674 if deep: self.__operator=escript.Operator()
675
676 def __rebuildRightHandSide(self,deep=False):
677 """
678 @brief indicates the right hand side has to be rebuild next time it is used
679 """
680 if self.__righthandside_isValid and self.debug() : print "PDE Debug: Right hand side has to be rebuilt."
681 self.__rebuildSolution(deep)
682 self.__righthandside_isValid=False
683 if not self.__homogeneous_constraint: self.__rebuildOperator()
684 if deep: self.__righthandside=escript.Data()
685
686 def __rebuildSystem(self,deep=False):
687 """
688 @brief annonced that all coefficient name has been changed
689 """
690 self.__rebuildSolution(deep)
691 self.__rebuildOperator(deep)
692 self.__rebuildRightHandSide(deep)
693
694 def __checkMatrixType(self):
695 """
696 @brief reassess the matrix type and, if needed, initiates an operator rebuild
697 """
698 new_matrix_type=self.getDomain().getSystemMatrixTypeId(self.getSolverMethod(),self.isSymmetric())
699 if not new_matrix_type==self.__matrix_type:
700 if self.debug() : print "PDE Debug: Matrix type is now %d."%new_matrix_type
701 self.__matrix_type=new_matrix_type
702 self.__rebuildOperator(deep=True)
703
704 #============ assembling =======================================================
705 def __copyConstraint(self,u):
706 """
707 @brief copies the constrint condition into u
708 """
709 q=self.getCoefficient("q")
710 r=self.getCoefficient("r")
711 if not q.isEmpty():
712 if r.isEmpty():
713 r2=escript.Data(0,u.getShape(),u.getFunctionSpace())
714 else:
715 r2=escript.Data(r,u.getFunctionSpace())
716 u.copyWithMask(r2,escript.Data(q,u.getFunctionSpace()))
717
718 def __applyConstraint(self,rhs_update=True):
719 """
720 @brief applies the constraints defined by q and r to the system
721 """
722 q=self.getCoefficient("q")
723 r=self.getCoefficient("r")
724 if not q.isEmpty() and not self.__operator.isEmpty():
725 # q is the row and column mask to indicate where constraints are set:
726 row_q=escript.Data(q,self.getFunctionSpaceForEquation())
727 col_q=escript.Data(q,self.getFunctionSpaceForSolution())
728 u=self.__makeNewSolution()
729 if r.isEmpty():
730 r_s=self.__makeNewSolution()
731 else:
732 r_s=escript.Data(r,self.getFunctionSpaceForSolution())
733 u.copyWithMask(r_s,col_q)
734 if not self.__righthandside.isEmpty() and rhs_update: self.__righthandside-=self.__operator*u
735 self.__operator.nullifyRowsAndCols(row_q,col_q,1.)
736
737 def getOperator(self):
738 """
739 @brief returns the operator of the PDE
740 """
741 if not self.__operator_isValid:
742 # some Constraints are applying for a lumpled stifness matrix:
743 if self.isUsingLumping():
744 if self.getFunctionSpaceForEquation()==self.getFunctionSpaceForSolution():
745 raise TypeError,"Lumped matrix requires same order for equations and unknowns"
746 if not self.getCoefficient("A").isEmpty():
747 raise Warning,"Lumped matrix does not allow coefficient A"
748 if not self.getCoefficient("B").isEmpty():
749 raise Warning,"Lumped matrix does not allow coefficient B"
750 if not self.getCoefficient("C").isEmpty():
751 raise Warning,"Lumped matrix does not allow coefficient C"
752 if not self.getCoefficient("D").isEmpty():
753 raise Warning,"Lumped matrix does not allow coefficient D"
754 if self.debug() : print "PDE Debug: New lumped operator is built."
755 mat=self.__makeNewOperator(self)
756 else:
757 if self.debug() : print "PDE Debug: New operator is built."
758 mat=self.__getFreshOperator()
759
760 self.getDomain().addPDEToSystem(mat,escript.Data(), \
761 self.getCoefficient("A"), \
762 self.getCoefficient("B"), \
763 self.getCoefficient("C"), \
764 self.getCoefficient("D"), \
765 escript.Data(), \
766 escript.Data(), \
767 self.getCoefficient("d"), \
768 escript.Data(),\
769 self.getCoefficient("d_contact"), \
770 escript.Data())
771 if self.isUsingLumping():
772 self.__operator=mat*escript.Data(1,(self.getNumSolutions(),),self.getFunctionSpaceOfSolution(),True)
773 else:
774 self.__applyConstraint(rhs_update=False)
775 self.__operator_isValid=True
776 return self.__operator
777
778 def getRightHandSide(self,ignoreConstraint=False):
779 """
780 @brief returns the right hand side of the PDE
781
782 @param ignoreConstraint
783 """
784 if not self.__righthandside_isValid:
785 if self.debug() : print "PDE Debug: New right hand side is built."
786 self.getDomain().addPDEToRHS(self.__getFreshRightHandSide(), \
787 self.getCoefficient("X"), \
788 self.getCoefficient("Y"),\
789 self.getCoefficient("y"),\
790 self.getCoefficient("y_contact"))
791 self.__righthandside_isValid=True
792 if ignoreConstraint: self.__copyConstraint(self.__righthandside)
793 return self.__righthandside
794
795 def getSystem(self):
796 """
797 @brief
798 """
799 if not self.__operator_isValid and not self.__righthandside_isValid:
800 if self.debug() : print "PDE Debug: New PDE is built."
801 if self.isUsingLumping():
802 self.getRightHandSide(ignoreConstraint=True)
803 self.getOperator()
804 else:
805 self.getDomain().addPDEToSystem(self.__getFreshOperator(),self.__getFreshRightHandSide(), \
806 self.getCoefficient("A"), \
807 self.getCoefficient("B"), \
808 self.getCoefficient("C"), \
809 self.getCoefficient("D"), \
810 self.getCoefficient("X"), \
811 self.getCoefficient("Y"), \
812 self.getCoefficient("d"), \
813 self.getCoefficient("y"), \
814 self.getCoefficient("d_contact"), \
815 self.getCoefficient("y_contact"))
816 self.__operator_isValid=True
817 self.__righthandside_isValid=True
818 self.__applyConstraint()
819 self.__copyConstraint(self.__righthandside)
820 elif not self.__operator_isValid:
821 self.getOperator()
822 elif not self.__righthandside_isValid:
823 self.getRightHandSide()
824 return (self.__operator,self.__righthandside)
825
826 def solve(self,**options):
827 """
828 @brief solve the PDE
829
830 @param options
831 """
832 mat,f=self.getSystem()
833 if self.isUsingLumping():
834 out=f/mat
835 self.__copyConstraint(out)
836 else:
837 options[util.TOLERANCE_KEY]=self.getTolerance()
838 options[util.METHOD_KEY]=self.getSolverMethod()
839 options[util.SYMMETRY_KEY]=self.isSymmetric()
840 if self.debug() : print "PDE Debug: solver options: ",options
841 out=mat.solve(f,options)
842 return out
843
844 def getSolution(self,**options):
845 """
846 @brief returns the solution of the PDE
847
848 @param options
849 """
850 if not self.__solution_isValid:
851 if self.debug() : print "PDE Debug: PDE is resolved."
852 self.__solution=self.solve(**options)
853 self.__solution_isValid=True
854 return self.__solution
855 #============ some serivice functions =====================================================
856 def getDomain(self):
857 """
858 @brief returns the domain of the PDE
859 """
860 return self.__domain
861
862 def getDim(self):
863 """
864 @brief returns the spatial dimension of the PDE
865 """
866 return self.getDomain().getDim()
867
868 def getNumEquations(self):
869 """
870 @brief returns the number of equations
871 """
872 if self.__numEquations>0:
873 return self.__numEquations
874 else:
875 raise ValueError,"Number of equations is undefined. Please specify argument numEquations."
876
877 def getNumSolutions(self):
878 """
879 @brief returns the number of unknowns
880 """
881 if self.__numSolutions>0:
882 return self.__numSolutions
883 else:
884 raise ValueError,"Number of solution is undefined. Please specify argument numSolutions."
885
886
887 def checkSymmetry(self):
888 """
889 @brief returns if the Operator is symmetric. This is a very expensive operation!!! The symmetry flag is not altered.
890 """
891 raise SystemError,"checkSymmetry is not implemented yet"
892
893 return None
894
895 def getFlux(self,u):
896 """
897 @brief returns the flux J_ij for a given u
898
899 J_ij=A_{ijkl}u_{k,l}+B_{ijk}u_k-X_{ij}
900
901 @param u argument of the operator
902
903 """
904 raise SystemError,"getFlux is not implemented yet"
905 return None
906
907 def applyOperator(self,u):
908 """
909 @brief applies the operator of the PDE to a given solution u in weak from
910
911 @param u argument of the operator
912
913 """
914 return self.getOperator()*escript.Data(u,self.getFunctionSpaceForSolution())
915
916 def getResidual(self,u):
917 """
918 @brief return the residual of u in the weak from
919
920 @param u
921 """
922 return self.applyOperator(u)-self.getRightHandSide()
923
924 class Poisson(LinearPDE):
925 """
926 @brief Class to define a Poisson equstion problem:
927
928 class to define a linear PDE of the form
929
930 -u_{,jj} = f
931
932 with boundary conditons:
933
934 n_j*u_{,j} = 0
935
936 and constraints:
937
938 u=0 where q>0
939
940 """
941
942 def __init__(self,domain=None,f=escript.Data(),q=escript.Data()):
943 LinearPDE.__init__(self,domain=identifyDomain(domain,{"f":f, "q":q}))
944 self._setValue(A=numarray.identity(self.getDomain().getDim()))
945 self.setSymmetryOn()
946 self.setValue(f,q)
947
948 def setValue(self,f=escript.Data(),q=escript.Data()):
949 self._setValue(Y=f,q=q)
950
951
952 # $Log$
953 # Revision 1.3 2004/12/17 07:43:10 jgs
954 # *** empty log message ***
955 #
956 # Revision 1.1.2.3 2004/12/16 00:12:34 gross
957 # __init__ of LinearPDE does not accept any coefficients anymore
958 #
959 # Revision 1.1.2.2 2004/12/14 03:55:01 jgs
960 # *** empty log message ***
961 #
962 # Revision 1.1.2.1 2004/12/12 22:53:47 gross
963 # linearPDE has been renamed LinearPDE
964 #
965 # Revision 1.1.1.1.2.7 2004/12/07 10:13:08 gross
966 # GMRES added
967 #
968 # Revision 1.1.1.1.2.6 2004/12/07 03:19:50 gross
969 # options for GMRES and PRES20 added
970 #
971 # Revision 1.1.1.1.2.5 2004/12/01 06:25:15 gross
972 # some small changes
973 #
974 # Revision 1.1.1.1.2.4 2004/11/24 01:50:21 gross
975 # Finley solves 4M unknowns now
976 #
977 # Revision 1.1.1.1.2.3 2004/11/15 06:05:26 gross
978 # poisson solver added
979 #
980 # Revision 1.1.1.1.2.2 2004/11/12 06:58:15 gross
981 # a lot of changes to get the linearPDE class running: most important change is that there is no matrix format exposed to the user anymore. the format is chosen by the Domain according to the solver and symmetry
982 #
983 # Revision 1.1.1.1.2.1 2004/10/28 22:59:22 gross
984 # finley's RecTest.py is running now: problem in SystemMatrixAdapater fixed
985 #
986 # Revision 1.1.1.1 2004/10/26 06:53:56 jgs
987 # initial import of project esys2
988 #
989 # Revision 1.3.2.3 2004/10/26 06:43:48 jgs
990 # committing Lutz's and Paul's changes to brach jgs
991 #
992 # Revision 1.3.4.1 2004/10/20 05:32:51 cochrane
993 # Added incomplete Doxygen comments to files, or merely put the docstrings that already exist into Doxygen form.
994 #
995 # Revision 1.3 2004/09/23 00:53:23 jgs
996 # minor fixes
997 #
998 # Revision 1.1 2004/08/28 12:58:06 gross
999 # SimpleSolve is not running yet: problem with == of functionsspace
1000 #
1001 #

Properties

Name Value
svn:eol-style native
svn:keywords Author Date Id Revision

  ViewVC Help
Powered by ViewVC 1.1.26