/[escript]/branches/symbolic_from_3470/escript/py_src/nonlinearPDE.py
ViewVC logotype

Contents of /branches/symbolic_from_3470/escript/py_src/nonlinearPDE.py

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3811 - (show annotations)
Wed Feb 8 01:13:28 2012 UTC (7 years ago) by gross
File MIME type: text/x-python
File size: 18917 byte(s)
some methods added to NonLinearPDE
1 # -*- coding: utf-8 -*-
2
3 ########################################################
4 #
5 # Copyright (c) 2003-2012 by University of Queensland
6 # Earth Systems Science Computational Center (ESSCC)
7 # http://www.uq.edu.au/esscc
8 #
9 # Primary Business: Queensland, Australia
10 # Licensed under the Open Software License version 3.0
11 # http://www.opensource.org/licenses/osl-3.0.php
12 #
13 ########################################################
14
15 __copyright__="""Copyright (c) 2003-2012 by University of Queensland
16 Earth Systems Science Computational Center (ESSCC)
17 http://www.uq.edu.au/esscc
18 Primary Business: Queensland, Australia"""
19 __license__="""Licensed under the Open Software License version 3.0
20 http://www.opensource.org/licenses/osl-3.0.php"""
21 __url__="https://launchpad.net/escript-finley"
22
23 """
24 :var __author__: name of author
25 :var __copyright__: copyrights
26 :var __license__: licence agreement
27 :var __url__: url entry point on documentation
28 :var __version__: version
29 :var __date__: date of the version
30 """
31
32 import numpy
33 from time import time
34 from esys.escript.linearPDEs import LinearPDE, IllegalCoefficient, IllegalCoefficientValue
35 from esys.escript import util, Data, Symbol, Evaluator
36
37 __author__="Cihan Altinay"
38
39 class NonlinearPDE(object):
40 """
41 This class is used to define a general nonlinear, steady, second order PDE
42 for an unknown function *u* on a given domain defined through a `Domain`
43 object.
44
45 For a single PDE having a solution with a single component the nonlinear
46 PDE is defined in the following form:
47
48 *-div(X) + Y = 0*
49
50 where *X*,*Y*=f(*u*,*grad(u)*), *div(F)* denotes the divergence of *F* and
51 *grad(F)* is the spatial derivative of *F*.
52
53 The coefficients *X* (rank 1) and *Y* (scalar) have to be specified through
54 `Symbol` objects.
55
56 The following natural boundary conditions are considered:
57
58 *n[j]*X[j] - y = 0*
59
60 where *n* is the outer normal field. Notice that the coefficient *X*
61 is defined in the PDE. The coefficient *y* is a scalar `Symbol`.
62
63 Constraints for the solution prescribe the value of the solution at certain
64 locations in the domain. They have the form
65
66 *u=u0* where *q>0*
67
68 *u0* and *q* are each scalar where *q* is the characteristic function
69 defining where the constraint is applied. The constraints override any
70 other condition set by the PDE or the boundary condition.
71
72 For a system of PDEs and a solution with several components, *u* is rank
73 one, while the PDE coefficient *X* is rank two and *Y* and *y* is rank one.
74
75 The PDE is solved by linearising the coefficients and iteratively solving
76 the corresponding linear PDE until the error is smaller than a tolerance
77 or a maximum number of iterations is reached.
78
79 Typical usage::
80
81 u = Symbol('u', dim=dom.getDim())
82 p = NonlinearPDE(dom, u)
83 p.setValue(X=grad(u), Y=5*u)
84 v = p.getSolution(u=0.)
85 """
86
87 def __init__(self,domain,u,debug=False):
88 """
89 Initializes a new nonlinear PDE.
90
91 :param domain: domain of the PDE
92 :type domain: `Domain`
93 :param u: The symbol for the unknown PDE function u.
94 :type u: `Symbol`
95 :param debug: if True debug information is printed
96 """
97 self.__COEFFICIENTS = [ "X", "X_reduced", "Y", "Y_reduced", "y", "y_reduced", "y_contact", "y_contact_reduced", "q" ]
98 self.coeffs={}
99 self._u=u
100 self._debug=debug
101 self.dim = domain.getDim()
102 if u.getRank()==0:
103 numEquations=1
104 else:
105 numEquations=u.getShape()[0]
106 numSolutions=numEquations
107 self.lpde=LinearPDE(domain,numEquations,numSolutions,debug)
108
109 def __str__(self):
110 """
111 Returns the string representation of the PDE.
112
113 :return: a simple representation of the PDE
114 :rtype: ``str``
115 """
116 return "<NonlinearPDE %d>"%id(self)
117
118 def getSolutionSymbol(self):
119 """
120 Returns the symbol of the PDE solution
121
122 :return: returns the symbol of the PDE solution
123 :rtype: `Symbol`
124 """
125 return self._u
126
127 def trace(self, text):
128 """
129 Prints the text message if debug mode is switched on.
130
131 :param text: message to be printed
132 :type text: ``string``
133 """
134 if self._debug:
135 print("%s: %s"%(str(self), text))
136
137 def getSolverOptions(self):
138 return self.lpde.getSolverOptions()
139
140 def getSolution(self, **subs):
141 """
142 Returns the solution of the PDE.
143
144 :param subs: Substitutions for all symbols used in the coefficients
145 including the initial value for *u*.
146 :return: the solution
147 :rtype: `Data`
148 """
149
150 u_sym=self._u.atoms().pop().name
151 if not subs.has_key(u_sym):
152 raise KeyError("Initial value for '%s' missing."%u_sym)
153
154 ui=subs[u_sym]
155 if isinstance(ui,float) or isinstance(ui,int):
156 ui=Data(ui, self.lpde.getFunctionSpaceForSolution())
157 subs[u_sym]=ui
158
159 coeffs={}
160 # separate symbolic expressions from other coefficients
161 expressions={}
162 for n, e in self.coeffs.items():
163 if hasattr(e, "atoms"):
164 expressions[n]=e
165 else:
166 coeffs[n]=e
167
168 # evaluate symbolic expressions
169 ev=Evaluator(*expressions.values())
170 ev.subs(**subs)
171 self.trace("Starting expression evaluation.")
172 ttt0=time()
173 res=ev.evaluate()
174 ttt1=time()
175 self.trace("Expressions evaluated. Took %f seconds."%(ttt1-ttt0))
176 names=expressions.keys()
177 for i in range(len(names)):
178 coeffs[names[i]]=res[i]
179
180 # solve linear PDE
181 self.lpde.setValue(**coeffs)
182 u_new=self.lpde.getSolution()
183 self.trace("Initial error: %g"%(util.Lsup(u_new)))
184 self.trace("Initial RHS: %g"%(util.Lsup(self.lpde.getRightHandSide())))
185 n=1
186 # perform Newton iterations until error is small enough or
187 # maximum number of iterations reached
188 while util.Lsup(u_new)>self.lpde.getSolverOptions().getTolerance() and \
189 n<self.lpde.getSolverOptions().getIterMax():
190 delta_u=u_new
191 ev.subs(**{u_sym:ui-u_new})
192 res=ev.evaluate()
193 for i in range(len(names)):
194 coeffs[names[i]]=res[i]
195 self.trace("Lsup(%s)=%s"%(names[i],util.Lsup(res[i])))
196
197 self.lpde.setValue(**coeffs)
198 u_new=self.lpde.getSolution()
199 ui=ui-delta_u
200 n=n+1
201 self.trace("Error after %d iterations: %g"%(n,util.Lsup(u_new)))
202 self.trace("RHS after %d iterations: %g"%(n,util.Lsup(self.lpde.getRightHandSide())))
203
204 self.trace("Final error after %d iterations: %g"%(n,util.Lsup(u_new)))
205 return u_new
206
207 def getNumSolutions(self):
208 """
209 Returns the number of the solution components
210 :rtype: ``int``
211 """
212 s=self._u.getShape()
213 if len(s) > 0:
214 return s[0]
215 else:
216 return 1
217
218 def getShapeOfCoefficient(self,name):
219 """
220 Returns the shape of the coefficient ``name``.
221
222 :param name: name of the coefficient enquired
223 :type name: ``string``
224 :return: the shape of the coefficient ``name``
225 :rtype: ``tuple`` of ``int``
226 :raise IllegalCoefficient: if ``name`` is not a coefficient of the PDE
227 """
228 numSol=self.getNumSolutions()
229 dim = self.dim
230 if name=="X" or name=="X_reduced":
231 if numSol > 1:
232 return (numSol,dim)
233 else:
234 return (dim,)
235 elif name=="q" :
236 if numSol > 1:
237 return (numSol,)
238 else:
239 return ()
240 elif name=="Y" or name=="Y_reduced":
241 if numSol > 1:
242 return (numSol,)
243 else:
244 return ()
245 elif name=="y" or name=="y_reduced":
246 if numSol > 1:
247 return (numSol,)
248 else:
249 return ()
250 elif name=="y_contact" or name=="y_contact_reduced":
251 if numSol > 1:
252 return (numSol,)
253 else:
254 return ()
255 else:
256 raise IllegalCoefficient("Attempt to request unknown coefficient %s"%name)
257
258 def createCoefficient(self, name):
259 """
260 create a new coefficient ``name`` as Symbol
261
262 :param name: name of the coefficient requested
263 :type name: ``string``
264 :return: the value of the coefficient
265 :rtype: `Symbol`
266 :raise IllegalCoefficient: if ``name`` is not a coefficient of the PDE
267 """
268 s=self.getShapeOfCoefficient(name)
269 return Symbol(name, s)
270
271 def getCoefficient(self, name):
272 """
273 Returns the value of the coefficient ``name`` as Symbol
274
275 :param name: name of the coefficient requested
276 :type name: ``string``
277 :return: the value of the coefficient
278 :rtype: `Symbol`
279 :raise IllegalCoefficient: if ``name`` is not a coefficient of the PDE
280 """
281 if self.coeffs.has_key(name):
282 return self.coeffs[name]
283 else:
284 return IllegalCoefficient("Attempt to request undefined coefficient %s"%name)
285
286
287 def setValue(self,**coefficients):
288 """
289 Sets new values to one or more coefficients.
290
291 :param coefficients: new values assigned to coefficients
292 :keyword X: value for coefficient ``X``
293 :type X: `Symbol` or any type that can be cast to a `Data` object
294 :keyword Y: value for coefficient ``Y``
295 :type Y: `Symbol` or any type that can be cast to a `Data` object
296 :keyword y: value for coefficient ``y``
297 :type y: `Symbol` or any type that can be cast to a `Data` object
298 :keyword y_contact: value for coefficient ``y_contact``
299 :type y_contact: `Symbol` or any type that can be cast to a `Data`
300 object
301 :keyword q: mask for location of constraints
302 :type q: `Symbol` or any type that can be cast to a `Data` object
303
304 :raise IllegalCoefficient: if an unknown coefficient keyword is used
305 :raise IllegalCoefficientValue: if a supplied coefficient value has an
306 invalid shape
307 """
308
309 u=self._u
310 for name,val in coefficients.iteritems():
311 shape=util.getShape(val)
312 if not shape == self.getShapeOfCoefficient(name):
313 IllegalCoefficientValue("%s has shape %s but must have shape %d"%(name, self.getShapeOfCoefficient(name), shape ) )
314 rank=len(shape)
315 if name=="X" or name=="X_reduced":
316 # DX/Du = del_X/del_u + del_X/del_grad(u)*del_grad(u)/del_u
317 # \ / \ /
318 # B A
319
320 if rank != u.getRank()+1:
321 raise IllegalCoefficientValue("%s must have rank %d"%(name,u.getRank()+1))
322 if not hasattr(val, 'diff'):
323 A=numpy.zeros(shape+u.grad().getShape())
324 B=numpy.zeros(shape+u.getShape())
325 elif u.getRank()==0:
326 X=self._removeFsFromGrad(val)
327 dXdu=X.diff(u)
328 dgdu=u.grad().diff(u)
329 A=numpy.empty(shape+dgdu.getShape(), dtype=object)
330 B=numpy.empty(shape, dtype=object)
331 ttt0=time()
332 for i in numpy.ndindex(shape):
333 for j in numpy.ndindex(dgdu.getShape()):
334 y=dXdu[i]
335 x=dgdu[j]
336 # expand() and coeff() are very expensive so
337 # we set the unwanted factors to zero to extract
338 # the one we need
339 for jj in numpy.ndindex(dgdu.getShape()):
340 if j==jj: continue
341 y=y.subs(dgdu[jj], 0)
342 B[i]=y.subs(x,0) # terms in u and constants
343 A[i,j]=y.subs(x,1)-B[i]
344 A=Symbol(A)
345 B=Symbol(B)
346 self.trace("Computing A, B took %f seconds."%(time()-ttt0))
347 else: #u.getRank()==1
348 X=self._removeFsFromGrad(val)
349 dXdu=X.diff(u)
350 dgdu=u.grad().diff(u).transpose(2)
351 I,J=shape
352 K,L=u.grad().getShape()
353 A=numpy.empty((I,J,K,L), dtype=object)
354 B=numpy.empty((I,J,K), dtype=object)
355 ttt0=time()
356 for i,j,k,l in numpy.ndindex(I,J,K,L):
357 if dgdu[k,k,l]==0:
358 A[i,j,k,l]=0
359 B[i,j,k]=0
360 else:
361 y=dXdu[i,j,k]
362 x=dgdu[k,k,l]
363 for kk,ll in numpy.ndindex(K,L):
364 if k==kk and l==ll: continue
365 y=y.subs(dgdu[kk,kk,ll], 0)
366 B[i,j,k]=y.subs(x,0) # terms in u and constants
367 A[i,j,k,l]=y.subs(x,1)-B[i,j,k]
368
369 A=Symbol(A)
370 B=Symbol(B)
371 self.trace("Computing A, B took %f seconds."%(time()-ttt0))
372
373 if name=='X_reduced':
374 self.coeffs['A_reduced']=A
375 self.coeffs['B_reduced']=B
376 self.coeffs['X_reduced']=val
377 else:
378 self.coeffs['A']=A
379 self.coeffs['B']=B
380 self.coeffs['X']=val
381 elif name=="Y" or name=="Y_reduced":
382 # DY/Du = del_Y/del_u + del_Y/del_grad(u)*del_grad(u)/del_u
383 # \ / \ /
384 # D C
385 if rank != u.getRank():
386 raise IllegalCoefficientValue("%s must have rank %d"%(name,u.getRank()))
387 if not hasattr(val, 'diff'):
388 C=numpy.zeros(shape+u.grad().getShape())
389 D=numpy.zeros(shape+u.getShape())
390 elif u.getRank()==0:
391 Y=self._removeFsFromGrad(val)
392 dYdu=Y.diff(u)
393 dgdu=u.grad().diff(u)
394 C=numpy.empty(dgdu.getShape(), dtype=object)
395 for j in numpy.ndindex(dgdu.getShape()):
396 y=dYdu
397 x=dgdu[j]
398 # expand() and coeff() are very expensive so
399 # we set the unwanted factors to zero to extract
400 # the one we need
401 for jj in numpy.ndindex(dgdu.getShape()):
402 if j==jj: continue
403 y=y.subs(dgdu[jj], 0)
404 D=y.subs(x,0) # terms in u and constants
405 C[j]=y.subs(x,1)-D
406 C=Symbol(C)
407 else: #u.getRank()==1
408 Y=self._removeFsFromGrad(val)
409 dYdu=Y.diff(u)
410 dgdu=u.grad().diff(u).transpose(2)
411 I,=shape
412 J,K=u.grad().getShape()
413 C=numpy.empty((I,J,K), dtype=object)
414 D=numpy.empty((I,J), dtype=object)
415 ttt0=time()
416 for i,j,k in numpy.ndindex(I,J,K):
417 if dgdu[j,j,k]==0:
418 C[i,j,k]=0
419 D[i,j]=0
420 else:
421 y=dYdu[i,j]
422 x=dgdu[j,j,k]
423 for jj,kk in numpy.ndindex(J,K):
424 if j==jj and k==kk: continue
425 y=y.subs(dgdu[jj,jj,kk], 0)
426 D[i,j]=y.subs(x,0) # terms in u and constants
427 C[i,j,k]=y.subs(x,1)-D[i,j]
428 C=Symbol(C)
429 D=Symbol(D)
430 self.trace("Computing C, D took %f seconds."%(time()-ttt0))
431
432 if name=='Y_reduced':
433 self.coeffs['C_reduced']=C
434 self.coeffs['D_reduced']=D
435 self.coeffs['Y_reduced']=val
436 else:
437 self.coeffs['C']=C
438 self.coeffs['D']=D
439 self.coeffs['Y']=val
440 elif name=="y" or name=="y_reduced":
441 y=val
442 if rank != u.getRank():
443 raise IllegalCoefficientValue("%s must have rank %d"%(name,u.getRank()))
444 if not hasattr(y, 'diff'):
445 d=numpy.zeros(u.getShape())
446 else:
447 d=y.diff(u)
448 if name=='y_reduced':
449 self.coeffs['d_reduced']=d
450 self.coeffs['y_reduced']=y
451 else:
452 self.coeffs['d']=d
453 self.coeffs['y']=y
454 elif name=="y_contact" or name=="y_contact_reduced":
455 y_contact=val
456 if rank != u.getRank():
457 raise IllegalCoefficientValue("%s must have rank %d"%(name,u.getRank()))
458 if not hasattr(y_contact, 'diff'):
459 d_contact=numpy.zeros(u.getShape())
460 else:
461 d_contact=y_contact.diff(u)
462 if name=='y_contact_reduced':
463 self.coeffs['d_contact_reduced']=d_contact
464 self.coeffs['y_contact_reduced']=y_contact
465 else:
466 self.coeffs['d_contact']=d_contact
467 self.coeffs['y_contact']=y_contact
468 elif name=="q":
469 if rank != u.getRank():
470 raise IllegalCoefficientValue("q must have rank %d"%u.getRank())
471 self.coeffs['q']=val
472 else:
473 raise IllegalCoefficient("Attempt to set unknown coefficient %s"%name)
474
475 def _removeFsFromGrad(self, sym):
476 """
477 Returns sym with all occurrences grad_n(a,b,c) replaced by grad_n(a,b).
478 That is, all functionspace parameters are removed.
479 This method is used in setValue() in order to find coefficients of
480 grad(u) with the coeff() function.
481 """
482 from esys.escript import symfn
483 gg=sym.atoms(symfn.grad_n)
484 for g in gg:
485 if len(g.args)==3:
486 r=symfn.grad_n(*g.args[:2])
487 sym=sym.subs(g, r)
488 return sym
489

  ViewVC Help
Powered by ViewVC 1.1.26