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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 102 - (show annotations)
Wed Dec 15 07:08:39 2004 UTC (14 years, 8 months ago) by jgs
File MIME type: text/x-python
File size: 37355 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,**args):
219 """
220 @brief initializes a new linear PDE.
221
222 @param args
223 """
224
225 # initialize attributes
226 self.__debug=None
227 self.__domain=None
228 self.__numEquations=0
229 self.__numSolutions=0
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 # check the arguments
240 coef={}
241 for arg in args.iterkeys():
242 if arg=="domain":
243 self.__domain=args[arg]
244 elif arg=="numEquations":
245 self.__numEquations=args[arg]
246 elif arg=="numSolutions":
247 self.__numSolutions=args[arg]
248 elif _PDECoefficientTypes.has_key(arg):
249 coef[arg]=args[arg]
250 else:
251 raise ValueError,"Illegal argument %s"%arg
252
253 # get the domain of the PDE
254 self.__domain=identifyDomain(self.__domain,coef)
255
256 # set some default values:
257 self.__homogeneous_constraint=True
258 self.__row_function_space=escript.Solution(self.__domain)
259 self.__column_function_space=escript.Solution(self.__domain)
260 self.__tolerance=1.e-8
261 self.__solver_method=util.DEFAULT_METHOD
262 self.__matrix_type=self.__domain.getSystemMatrixTypeId(util.DEFAULT_METHOD,False)
263 self.__sym=False
264 self.__lumping=False
265 self.__numEquations=0
266 self.__numSolutions=0
267 # now we can set the ceofficients:
268 self._setCoefficient(**coef)
269
270 def getCoefficient(self,name):
271 """
272 @brief return the value of the coefficient name
273
274 @param name
275 """
276 return self.__coefficient[name]
277
278 def setValue(self,**coefficients):
279 """
280 @brief sets new values to coefficients
281
282 @param coefficients
283 """
284 self._setCoefficient(**coefficients)
285
286
287 def _setCoefficient(self,**coefficients):
288 """
289 @brief sets new values to coefficients
290
291 @param coefficients
292 """
293
294 # get the dictionary of the coefficinets been altered:
295 alteredCoefficients={}
296 for i,d in coefficients.iteritems():
297 if self.hasCoefficient(i):
298 if d == None:
299 alteredCoefficients[i]=escript.Data()
300 elif isinstance(d,escript.Data):
301 if d.isEmpty():
302 alteredCoefficients[i]=escript.Data()
303 else:
304 alteredCoefficients[i]=escript.Data(d,self.getFunctionSpaceOfCoefficient(i))
305 else:
306 if self.__numEquations>0 and self.__numSolutions>0:
307 alteredCoefficients[i]=escript.Data(d,self.getShapeOfCoefficient(i),self.getFunctionSpaceOfCoefficient(i))
308 else:
309 alteredCoefficients[i]=escript.Data(d,self.getFunctionSpaceOfCoefficient(i))
310 else:
311 raise ValueError,"Attempt to set undefined coefficient %s"%i
312 # if numEquations and numSolutions is undefined we try identify their values based on the coefficients:
313 if self.__numEquations<1 or self.__numSolutions<1:
314 numEquations,numSolutions=identifyNumEquationsAndSolutions(self.getDomain().getDim(),alteredCoefficients)
315 if self.__numEquations<1 and numEquations>0: self.__numEquations=numEquations
316 if self.__numSolutions<1 and numSolutions>0: self.__numSolutions=numSolutions
317 if self.debug() and self.__numEquations>0: print "PDE Debug: identified number of equations is ",self.__numEquations
318 if self.debug() and self.__numSolutions>0: print "PDE Debug: identified number of solutions is ",self.__numSolutions
319
320 # now we check the shape of the coefficient if numEquations and numSolutions are set:
321 if self.__numEquations>0 and self.__numSolutions>0:
322 for i in self.__coefficient.iterkeys():
323 if alteredCoefficients.has_key(i) and not alteredCoefficients[i].isEmpty():
324 if not self.getShapeOfCoefficient(i)==alteredCoefficients[i].getShape():
325 raise ValueError,"Expected shape for coefficient %s is %s but actual shape is %s."%(i,self.getShapeOfCoefficient(i),alteredCoefficients[i].getShape())
326 else:
327 if not self.__coefficient[i].isEmpty():
328 if not self.getShapeOfCoefficient(i)==self.__coefficient[i].getShape():
329 raise ValueError,"Expected shape for coefficient %s is %s but actual shape is %s."%(i,self.getShapeOfCoefficient(i),self.__coefficient[i].getShape())
330 # overwrite new values:
331 for i,d in alteredCoefficients.iteritems():
332 if self.debug(): print "PDE Debug: Coefficient %s has been altered."%i
333 self.__coefficient[i]=d
334 self.alteredCoefficient(i)
335
336 # reset the HomogeneousConstraintFlag:
337 self.__setHomogeneousConstraintFlag()
338 if not "q" in alteredCoefficients and not self.__homogeneous_constraint: self.__rebuildSystem()
339
340 def cleanCoefficients(self):
341 """
342 @brief resets all coefficients to default values.
343 """
344 self.__coefficient={}
345 for i in _PDECoefficientTypes.iterkeys():
346 self.__coefficient[i]=escript.Data()
347
348 def getShapeOfCoefficient(self,name):
349 """
350 @brief return the shape of the coefficient name
351
352 @param name
353 """
354 if self.hasCoefficient(name):
355 return _PDECoefficientTypes[name].buildShape(self.getNumEquations(),self.getNumSolutions(),self.getDomain().getDim())
356 else:
357 raise ValueError,"Solution coefficient %s requested"%name
358
359 def getFunctionSpaceOfCoefficient(self,name):
360 """
361 @brief return the atoms of the coefficient name
362
363 @param name
364 """
365 if self.hasCoefficient(name):
366 return _PDECoefficientTypes[name].getFunctionSpace(self.getDomain())
367 else:
368 raise ValueError,"Solution coefficient %s requested"%name
369
370 def alteredCoefficient(self,name):
371 """
372 @brief annonced that coefficient name has been changed
373
374 @param name
375 """
376 if self.hasCoefficient(name):
377 if _PDECoefficientTypes[name].isAlteringOperator(): self.__rebuildOperator()
378 if _PDECoefficientTypes[name].isAlteringRightHandSide(): self.__rebuildRightHandSide()
379 else:
380 raise ValueError,"Solution coefficient %s requested"%name
381
382 def __setHomogeneousConstraintFlag(self):
383 """
384 @brief checks if the constraints are homogeneous and sets self.__homogeneous_constraint accordingly.
385 """
386 self.__homogeneous_constraint=True
387 q=self.getCoefficient("q")
388 r=self.getCoefficient("r")
389 if not q.isEmpty() and not r.isEmpty():
390 print (q*r).Lsup(), 1.e-13*r.Lsup()
391 if (q*r).Lsup()>=1.e-13*r.Lsup(): self.__homogeneous_constraint=False
392 if self.debug():
393 if self.__homogeneous_constraint:
394 print "PDE Debug: Constraints are homogeneous."
395 else:
396 print "PDE Debug: Constraints are inhomogeneous."
397
398
399 def hasCoefficient(self,name):
400 """
401 @brief return true if name is the name of a coefficient
402
403 @param name
404 """
405 return self.__coefficient.has_key(name)
406
407 def getFunctionSpaceForEquation(self):
408 """
409 @brief return true if the test functions should use reduced order
410 """
411 return self.__row_function_space
412
413 def getFunctionSpaceForSolution(self):
414 """
415 @brief return true if the interpolation of the solution should use reduced order
416 """
417 return self.__column_function_space
418
419 # ===== debug ==============================================================
420 def setDebugOn(self):
421 """
422 @brief
423 """
424 self.__debug=not None
425
426 def setDebugOff(self):
427 """
428 @brief
429 """
430 self.__debug=None
431
432 def debug(self):
433 """
434 @brief returns true if the PDE is in the debug mode
435 """
436 return self.__debug
437
438 #===== Lumping ===========================
439 def setLumpingOn(self):
440 """
441 @brief indicates to use matrix lumping
442 """
443 if not self.isUsingLumping():
444 raise SystemError,"Lumping is not working yet! Talk to the experts"
445 if self.debug() : print "PDE Debug: lumping is set on"
446 self.__rebuildOperator()
447 self.__lumping=True
448
449 def setLumpingOff(self):
450 """
451 @brief switches off matrix lumping
452 """
453 if self.isUsingLumping():
454 if self.debug() : print "PDE Debug: lumping is set off"
455 self.__rebuildOperator()
456 self.__lumping=False
457
458 def setLumping(self,flag=False):
459 """
460 @brief set the matrix lumping flag to flag
461 """
462 if flag:
463 self.setLumpingOn()
464 else:
465 self.setLumpingOff()
466
467 def isUsingLumping(self):
468 """
469 @brief
470 """
471 return self.__lumping
472
473 #============ method business =========================================================
474 def setSolverMethod(self,solver=util.DEFAULT_METHOD):
475 """
476 @brief sets a new solver
477 """
478 if not solver==self.getSolverMethod():
479 self.__solver_method=solver
480 if self.debug() : print "PDE Debug: New solver is %s"%solver
481 self.__checkMatrixType()
482
483 def getSolverMethod(self):
484 """
485 @brief returns the solver method
486 """
487 return self.__solver_method
488
489 #============ tolerance business =========================================================
490 def setTolerance(self,tol=1.e-8):
491 """
492 @brief resets the tolerance to tol.
493 """
494 if not tol>0:
495 raise ValueException,"Tolerance as to be positive"
496 if tol<self.getTolerance(): self.__rebuildSolution()
497 if self.debug() : print "PDE Debug: New tolerance %e",tol
498 self.__tolerance=tol
499 return
500 def getTolerance(self):
501 """
502 @brief returns the tolerance set for the solution
503 """
504 return self.__tolerance
505
506 #===== symmetry flag ==========================
507 def isSymmetric(self):
508 """
509 @brief returns true is the operator is considered to be symmetric
510 """
511 return self.__sym
512
513 def setSymmetryOn(self):
514 """
515 @brief sets the symmetry flag to true
516 """
517 if not self.isSymmetric():
518 if self.debug() : print "PDE Debug: Operator is set to be symmetric"
519 self.__sym=True
520 self.__checkMatrixType()
521
522 def setSymmetryOff(self):
523 """
524 @brief sets the symmetry flag to false
525 """
526 if self.isSymmetric():
527 if self.debug() : print "PDE Debug: Operator is set to be unsymmetric"
528 self.__sym=False
529 self.__checkMatrixType()
530
531 def setSymmetryTo(self,flag=False):
532 """
533 @brief sets the symmetry flag to flag
534
535 @param flag
536 """
537 if flag:
538 self.setSymmetryOn()
539 else:
540 self.setSymmetryOff()
541
542 #===== order reduction ==========================
543 def setReducedOrderOn(self):
544 """
545 @brief switches to on reduced order
546 """
547 self.setReducedOrderForSolutionOn()
548 self.setReducedOrderForEquationOn()
549
550 def setReducedOrderOff(self):
551 """
552 @brief switches to full order
553 """
554 self.setReducedOrderForSolutionOff()
555 self.setReducedOrderForEquationOff()
556
557 def setReducedOrderTo(self,flag=False):
558 """
559 @brief sets order according to flag
560
561 @param flag
562 """
563 self.setReducedOrderForSolutionTo(flag)
564 self.setReducedOrderForEquationTo(flag)
565
566
567 #===== order reduction solution ==========================
568 def setReducedOrderForSolutionOn(self):
569 """
570 @brief switches to reduced order to interpolate solution
571 """
572 new_fs=escript.ReducedSolution(self.getDomain())
573 if self.getFunctionSpaceForSolution()!=new_fs:
574 if self.debug() : print "PDE Debug: Reduced order is used to interpolate solution."
575 self.__column_function_space=new_fs
576 self.__rebuildSystem(deep=True)
577
578 def setReducedOrderForSolutionOff(self):
579 """
580 @brief switches to full order to interpolate solution
581 """
582 new_fs=escript.Solution(self.getDomain())
583 if self.getFunctionSpaceForSolution()!=new_fs:
584 if self.debug() : print "PDE Debug: Full order is used to interpolate solution."
585 self.__column_function_space=new_fs
586 self.__rebuildSystem(deep=True)
587
588 def setReducedOrderForSolutionTo(self,flag=False):
589 """
590 @brief sets order for test functions according to flag
591
592 @param flag
593 """
594 if flag:
595 self.setReducedOrderForSolutionOn()
596 else:
597 self.setReducedOrderForSolutionOff()
598
599 #===== order reduction equation ==========================
600 def setReducedOrderForEquationOn(self):
601 """
602 @brief switches to reduced order for test functions
603 """
604 new_fs=escript.ReducedSolution(self.getDomain())
605 if self.getFunctionSpaceForEquation()!=new_fs:
606 if self.debug() : print "PDE Debug: Reduced order is used for test functions."
607 self.__row_function_space=new_fs
608 self.__rebuildSystem(deep=True)
609
610 def setReducedOrderForEquationOff(self):
611 """
612 @brief switches to full order for test functions
613 """
614 new_fs=escript.Solution(self.getDomain())
615 if self.getFunctionSpaceForEquation()!=new_fs:
616 if self.debug() : print "PDE Debug: Full order is used for test functions."
617 self.__row_function_space=new_fs
618 self.__rebuildSystem(deep=True)
619
620 def setReducedOrderForEquationTo(self,flag=False):
621 """
622 @brief sets order for test functions according to flag
623
624 @param flag
625 """
626 if flag:
627 self.setReducedOrderForEquationOn()
628 else:
629 self.setReducedOrderForEquationOff()
630
631
632 # ==== initialization =====================================================================
633 def __makeNewOperator(self):
634 """
635 @brief
636 """
637 return self.getDomain().newOperator( \
638 self.getNumEquations(), \
639 self.getFunctionSpaceForEquation(), \
640 self.getNumSolutions(), \
641 self.getFunctionSpaceForSolution(), \
642 self.__matrix_type)
643
644 def __makeNewRightHandSide(self):
645 """
646 @brief
647 """
648 return escript.Data(0.,(self.getNumEquations(),),self.getFunctionSpaceForEquation(),True)
649
650 def __makeNewSolution(self):
651 """
652 @brief
653 """
654 return escript.Data(0.,(self.getNumSolutions(),),self.getFunctionSpaceForSolution(),True)
655
656 def __getFreshOperator(self):
657 """
658 @brief
659 """
660 if self.__operator.isEmpty():
661 self.__operator=self.__makeNewOperator()
662 if self.debug() : print "PDE Debug: New operator allocated"
663 else:
664 self.__operator.setValue(0.)
665 if self.debug() : print "PDE Debug: Operator reset to zero"
666 return self.__operator
667
668 def __getFreshRightHandSide(self):
669 """
670 @brief
671 """
672 if self.__righthandside.isEmpty():
673 self.__righthandside=self.__makeNewRightHandSide()
674 if self.debug() : print "PDE Debug: New right hand side allocated"
675 else:
676 print "fix self.__righthandside*=0"
677 self.__righthandside*=0.
678 if self.debug() : print "PDE Debug: Right hand side reset to zero"
679 return self.__righthandside
680
681 # ==== rebuild switches =====================================================================
682 def __rebuildSolution(self,deep=False):
683 """
684 @brief indicates the PDE has to be reolved if the solution is requested
685 """
686 if self.__solution_isValid and self.debug() : print "PDE Debug: PDE has to be resolved."
687 self.__solution_isValid=False
688 if deep: self.__solution=escript.Data(deep)
689
690
691 def __rebuildOperator(self,deep=False):
692 """
693 @brief indicates the operator has to be rebuilt next time it is used
694 """
695 if self.__operator_isValid and self.debug() : print "PDE Debug: Operator has to be rebuilt."
696 self.__rebuildSolution(deep)
697 self.__operator_isValid=False
698 if deep: self.__operator=escript.Operator()
699
700 def __rebuildRightHandSide(self,deep=False):
701 """
702 @brief indicates the right hand side has to be rebuild next time it is used
703 """
704 if self.__righthandside_isValid and self.debug() : print "PDE Debug: Right hand side has to be rebuilt."
705 self.__rebuildSolution(deep)
706 self.__righthandside_isValid=False
707 if not self.__homogeneous_constraint: self.__rebuildOperator()
708 if deep: self.__righthandside=escript.Data()
709
710 def __rebuildSystem(self,deep=False):
711 """
712 @brief annonced that all coefficient name has been changed
713 """
714 self.__rebuildSolution(deep)
715 self.__rebuildOperator(deep)
716 self.__rebuildRightHandSide(deep)
717
718 def __checkMatrixType(self):
719 """
720 @brief reassess the matrix type and, if needed, initiates an operator rebuild
721 """
722 new_matrix_type=self.getDomain().getSystemMatrixTypeId(self.getSolverMethod(),self.isSymmetric())
723 if not new_matrix_type==self.__matrix_type:
724 if self.debug() : print "PDE Debug: Matrix type is now %d."%new_matrix_type
725 self.__matrix_type=new_matrix_type
726 self.__rebuildOperator(deep=True)
727
728 #============ assembling =======================================================
729 def __copyConstraint(self,u):
730 """
731 @brief copies the constrint condition into u
732 """
733 q=self.getCoefficient("q")
734 r=self.getCoefficient("r")
735 if not q.isEmpty():
736 if r.isEmpty():
737 r2=escript.Data(0,u.getShape(),u.getFunctionSpace())
738 else:
739 r2=escript.Data(r,u.getFunctionSpace())
740 u.copyWithMask(r2,escript.Data(q,u.getFunctionSpace()))
741
742 def __applyConstraint(self,rhs_update=True):
743 """
744 @brief applies the constraints defined by q and r to the system
745 """
746 q=self.getCoefficient("q")
747 r=self.getCoefficient("r")
748 if not q.isEmpty() and not self.__operator.isEmpty():
749 # q is the row and column mask to indicate where constraints are set:
750 row_q=escript.Data(q,self.getFunctionSpaceForEquation())
751 col_q=escript.Data(q,self.getFunctionSpaceForSolution())
752 u=self.__makeNewSolution()
753 if r.isEmpty():
754 r_s=self.__makeNewSolution()
755 else:
756 r_s=escript.Data(r,self.getFunctionSpaceForSolution())
757 u.copyWithMask(r_s,col_q)
758 if not self.__righthandside.isEmpty() and rhs_update: self.__righthandside-=self.__operator*u
759 self.__operator.nullifyRowsAndCols(row_q,col_q,1.)
760
761 def getOperator(self):
762 """
763 @brief returns the operator of the PDE
764 """
765 if not self.__operator_isValid:
766 # some Constraints are applying for a lumpled stifness matrix:
767 if self.isUsingLumping():
768 if self.getFunctionSpaceForEquation()==self.getFunctionSpaceForSolution():
769 raise TypeError,"Lumped matrix requires same order for equations and unknowns"
770 if not self.getCoefficient("A").isEmpty():
771 raise Warning,"Lumped matrix does not allow coefficient A"
772 if not self.getCoefficient("B").isEmpty():
773 raise Warning,"Lumped matrix does not allow coefficient B"
774 if not self.getCoefficient("C").isEmpty():
775 raise Warning,"Lumped matrix does not allow coefficient C"
776 if not self.getCoefficient("D").isEmpty():
777 raise Warning,"Lumped matrix does not allow coefficient D"
778 if self.debug() : print "PDE Debug: New lumped operator is built."
779 mat=self.__makeNewOperator(self)
780 else:
781 if self.debug() : print "PDE Debug: New operator is built."
782 mat=self.__getFreshOperator()
783
784 self.getDomain().addPDEToSystem(mat,escript.Data(), \
785 self.getCoefficient("A"), \
786 self.getCoefficient("B"), \
787 self.getCoefficient("C"), \
788 self.getCoefficient("D"), \
789 escript.Data(), \
790 escript.Data(), \
791 self.getCoefficient("d"), \
792 escript.Data(),\
793 self.getCoefficient("d_contact"), \
794 escript.Data())
795 if self.isUsingLumping():
796 self.__operator=mat*escript.Data(1,(self.getNumSolutions(),),self.getFunctionSpaceOfSolution(),True)
797 else:
798 self.__applyConstraint(rhs_update=False)
799 self.__operator_isValid=True
800 return self.__operator
801
802 def getRightHandSide(self,ignoreConstraint=False):
803 """
804 @brief returns the right hand side of the PDE
805
806 @param ignoreConstraint
807 """
808 if not self.__righthandside_isValid:
809 if self.debug() : print "PDE Debug: New right hand side is built."
810 self.getDomain().addPDEToRHS(self.__getFreshRightHandSide(), \
811 self.getCoefficient("X"), \
812 self.getCoefficient("Y"),\
813 self.getCoefficient("y"),\
814 self.getCoefficient("y_contact"))
815 self.__righthandside_isValid=True
816 if ignoreConstraint: self.__copyConstraint(self.__righthandside)
817 return self.__righthandside
818
819 def getSystem(self):
820 """
821 @brief
822 """
823 if not self.__operator_isValid and not self.__righthandside_isValid:
824 if self.debug() : print "PDE Debug: New PDE is built."
825 if self.isUsingLumping():
826 self.getRightHandSide(ignoreConstraint=True)
827 self.getOperator()
828 else:
829 self.getDomain().addPDEToSystem(self.__getFreshOperator(),self.__getFreshRightHandSide(), \
830 self.getCoefficient("A"), \
831 self.getCoefficient("B"), \
832 self.getCoefficient("C"), \
833 self.getCoefficient("D"), \
834 self.getCoefficient("X"), \
835 self.getCoefficient("Y"), \
836 self.getCoefficient("d"), \
837 self.getCoefficient("y"), \
838 self.getCoefficient("d_contact"), \
839 self.getCoefficient("y_contact"))
840 self.__operator_isValid=True
841 self.__righthandside_isValid=True
842 self.__applyConstraint()
843 self.__copyConstraint(self.__righthandside)
844 elif not self.__operator_isValid:
845 self.getOperator()
846 elif not self.__righthandside_isValid:
847 self.getRightHandSide()
848 return (self.__operator,self.__righthandside)
849
850 def solve(self,**options):
851 """
852 @brief solve the PDE
853
854 @param options
855 """
856 mat,f=self.getSystem()
857 if self.isUsingLumping():
858 out=f/mat
859 self.__copyConstraint(out)
860 else:
861 options[util.TOLERANCE_KEY]=self.getTolerance()
862 options[util.METHOD_KEY]=self.getSolverMethod()
863 options[util.SYMMETRY_KEY]=self.isSymmetric()
864 if self.debug() : print "PDE Debug: solver options: ",options
865 out=mat.solve(f,options)
866 return out
867
868 def getSolution(self,**options):
869 """
870 @brief returns the solution of the PDE
871
872 @param options
873 """
874 if not self.__solution_isValid:
875 if self.debug() : print "PDE Debug: PDE is resolved."
876 self.__solution=self.solve(**options)
877 self.__solution_isValid=True
878 return self.__solution
879 #============ some serivice functions =====================================================
880 def getDomain(self):
881 """
882 @brief returns the domain of the PDE
883 """
884 return self.__domain
885
886 def getNumEquations(self):
887 """
888 @brief returns the number of equations
889 """
890 if self.__numEquations>0:
891 return self.__numEquations
892 else:
893 raise ValueError,"Number of equations is undefined. Please specify argument numEquations."
894
895 def getNumSolutions(self):
896 """
897 @brief returns the number of unknowns
898 """
899 if self.__numSolutions>0:
900 return self.__numSolutions
901 else:
902 raise ValueError,"Number of solution is undefined. Please specify argument numSolutions."
903
904
905 def checkSymmetry(self):
906 """
907 @brief returns if the Operator is symmetric. This is a very expensive operation!!! The symmetry flag is not altered.
908 """
909 raise SystemError,"checkSymmetry is not implemented yet"
910
911 return None
912
913 def getFlux(self,u):
914 """
915 @brief returns the flux J_ij for a given u
916
917 J_ij=A_{ijkl}u_{k,l}+B_{ijk}u_k-X_{ij}
918
919 @param u argument of the operator
920
921 """
922 raise SystemError,"getFlux is not implemented yet"
923 return None
924
925 def applyOperator(self,u):
926 """
927 @brief applies the operator of the PDE to a given solution u in weak from
928
929 @param u argument of the operator
930
931 """
932 return self.getOperator()*escript.Data(u,self.getFunctionSpaceForSolution())
933
934 def getResidual(self,u):
935 """
936 @brief return the residual of u in the weak from
937
938 @param u
939 """
940 return self.applyOperator(u)-self.getRightHandSide()
941
942 class Poisson(LinearPDE):
943 """
944 @brief Class to define a Poisson equstion problem:
945
946 class to define a linear PDE of the form
947
948 -u_{,jj} = f
949
950 with boundary conditons:
951
952 n_j*u_{,j} = 0
953
954 and constraints:
955
956 u=0 where q>0
957
958 """
959
960 def __init__(self,domain=None,f=escript.Data(),q=escript.Data()):
961 LinearPDE.__init__(self,domain=identifyDomain(domain,{"f":f, "q":q}))
962 self._setCoefficient(A=numarray.identity(self.getDomain().getDim()))
963 self.setSymmetryOn()
964 self.setValue(f,q)
965
966 def setValue(self,f=escript.Data(),q=escript.Data()):
967 self._setCoefficient(Y=f,q=q)
968
969
970 # $Log$
971 # Revision 1.2 2004/12/15 07:08:27 jgs
972 # *** empty log message ***
973 #
974 # Revision 1.1.2.2 2004/12/14 03:55:01 jgs
975 # *** empty log message ***
976 #
977 # Revision 1.1.2.1 2004/12/12 22:53:47 gross
978 # linearPDE has been renamed LinearPDE
979 #
980 # Revision 1.1.1.1.2.7 2004/12/07 10:13:08 gross
981 # GMRES added
982 #
983 # Revision 1.1.1.1.2.6 2004/12/07 03:19:50 gross
984 # options for GMRES and PRES20 added
985 #
986 # Revision 1.1.1.1.2.5 2004/12/01 06:25:15 gross
987 # some small changes
988 #
989 # Revision 1.1.1.1.2.4 2004/11/24 01:50:21 gross
990 # Finley solves 4M unknowns now
991 #
992 # Revision 1.1.1.1.2.3 2004/11/15 06:05:26 gross
993 # poisson solver added
994 #
995 # Revision 1.1.1.1.2.2 2004/11/12 06:58:15 gross
996 # 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
997 #
998 # Revision 1.1.1.1.2.1 2004/10/28 22:59:22 gross
999 # finley's RecTest.py is running now: problem in SystemMatrixAdapater fixed
1000 #
1001 # Revision 1.1.1.1 2004/10/26 06:53:56 jgs
1002 # initial import of project esys2
1003 #
1004 # Revision 1.3.2.3 2004/10/26 06:43:48 jgs
1005 # committing Lutz's and Paul's changes to brach jgs
1006 #
1007 # Revision 1.3.4.1 2004/10/20 05:32:51 cochrane
1008 # Added incomplete Doxygen comments to files, or merely put the docstrings that already exist into Doxygen form.
1009 #
1010 # Revision 1.3 2004/09/23 00:53:23 jgs
1011 # minor fixes
1012 #
1013 # Revision 1.1 2004/08/28 12:58:06 gross
1014 # SimpleSolve is not running yet: problem with == of functionsspace
1015 #
1016 #

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26