/[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 3631 - (show annotations)
Thu Oct 20 03:12:46 2011 UTC (7 years, 6 months ago) by caltinay
File MIME type: text/x-python
File size: 16079 byte(s)
Fix compatibility with latest sympy version.

1 # -*- coding: utf-8 -*-
2
3 ########################################################
4 #
5 # Copyright (c) 2003-2010 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-2010 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=r* where *q>0*
67
68 *r* 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=1.)
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.coeffs={}
98 self._u=u
99 self._debug=debug
100 if u.getRank()==0:
101 numEquations=1
102 else:
103 numEquations=u.getShape()[0]
104 numSolutions=numEquations
105 self.lpde=LinearPDE(domain,numEquations,numSolutions,debug)
106
107 def __str__(self):
108 """
109 Returns the string representation of the PDE.
110
111 :return: a simple representation of the PDE
112 :rtype: ``str``
113 """
114 return "<NonlinearPDE %d>"%id(self)
115
116 def trace(self, text):
117 """
118 Prints the text message if debug mode is switched on.
119
120 :param text: message to be printed
121 :type text: ``string``
122 """
123 if self._debug:
124 print("%s: %s"%(str(self), text))
125
126 def getSolverOptions(self):
127 return self.lpde.getSolverOptions()
128
129 def getSolution(self, **subs):
130 """
131 Returns the solution of the PDE.
132
133 :param subs: Substitutions for all symbols used in the coefficients
134 including the initial value for *u*.
135 :return: the solution
136 :rtype: `Data`
137 """
138
139 u_sym=self._u.atoms().pop().name
140 if not subs.has_key(u_sym):
141 raise KeyError("Initial value for '%s' missing."%u_sym)
142
143 ui=subs[u_sym]
144 if isinstance(ui,float) or isinstance(ui,int):
145 ui=Data(ui, self.lpde.getFunctionSpaceForSolution())
146 subs[u_sym]=ui
147
148 coeffs={}
149 # separate symbolic expressions from other coefficients
150 expressions={}
151 for n, e in self.coeffs.items():
152 if hasattr(e, "atoms"):
153 expressions[n]=e
154 else:
155 coeffs[n]=e
156
157 # evaluate symbolic expressions
158 ev=Evaluator(*expressions.values())
159 ev.subs(**subs)
160 self.trace("Starting expression evaluation.")
161 ttt0=time()
162 res=ev.evaluate()
163 ttt1=time()
164 self.trace("Expressions evaluated. Took %f seconds."%(ttt1-ttt0))
165 names=expressions.keys()
166 for i in range(len(names)):
167 coeffs[names[i]]=res[i]
168
169 # solve linear PDE
170 self.lpde.setValue(**coeffs)
171 u_new=self.lpde.getSolution()
172 self.trace("Initial error: %g"%(util.Lsup(u_new)))
173 self.trace("Initial RHS: %g"%(util.Lsup(self.lpde.getRightHandSide())))
174 n=1
175 # perform Newton iterations until error is small enough or
176 # maximum number of iterations reached
177 while util.Lsup(u_new)>self.lpde.getSolverOptions().getTolerance() and \
178 n<self.lpde.getSolverOptions().getIterMax():
179 delta_u=u_new
180 ev.subs(**{u_sym:ui-u_new})
181 res=ev.evaluate()
182 for i in range(len(names)):
183 coeffs[names[i]]=res[i]
184 self.trace("Lsup(%s)=%s"%(names[i],util.Lsup(res[i])))
185
186 self.lpde.setValue(**coeffs)
187 u_new=self.lpde.getSolution()
188 ui=ui-delta_u
189 n=n+1
190 self.trace("Error after %d iterations: %g"%(n,util.Lsup(u_new)))
191 self.trace("RHS after %d iterations: %g"%(n,util.Lsup(self.lpde.getRightHandSide())))
192
193 self.trace("Final error after %d iterations: %g"%(n,util.Lsup(u_new)))
194 return u_new
195
196 def getCoefficient(self, name):
197 """
198 Returns the value of the coefficient ``name``.
199
200 :param name: name of the coefficient requested
201 :type name: ``string``
202 :return: the value of the coefficient
203 :rtype: `Symbol`
204 :raise IllegalCoefficient: if ``name`` is not a coefficient of the PDE
205 """
206 if self.hasCoefficient(name):
207 return self.coeffs[name]
208 raise IllegalCoefficient("Illegal coefficient %s requested."%name)
209
210 def hasCoefficient(self, name):
211 """
212 Returns True if ``name`` is the name of a coefficient.
213
214 :param name: name of the coefficient enquired
215 :type name: ``string``
216 :return: True if ``name`` is the name of a coefficient of the PDE,
217 False otherwise
218 :rtype: ``bool``
219
220 """
221 return self.coeffs.has_key(name)
222
223 def setValue(self,**coefficients):
224 """
225 Sets new values to one or more coefficients.
226
227 :param coefficients: new values assigned to coefficients
228 :keyword X: value for coefficient ``X``
229 :type X: `Symbol` or any type that can be cast to a `Data` object
230 :keyword Y: value for coefficient ``Y``
231 :type Y: `Symbol` or any type that can be cast to a `Data` object
232 :keyword y: value for coefficient ``y``
233 :type y: `Symbol` or any type that can be cast to a `Data` object
234 :keyword y_contact: value for coefficient ``y_contact``
235 :type y_contact: `Symbol` or any type that can be cast to a `Data`
236 object
237 :keyword r: values prescribed to the solution at the locations of
238 constraints
239 :type r: `Symbol` or any type that can be cast to a `Data` object
240 :keyword q: mask for location of constraints
241 :type q: `Symbol` or any type that can be cast to a `Data` object
242
243 :raise IllegalCoefficient: if an unknown coefficient keyword is used
244 :raise IllegalCoefficientValue: if a supplied coefficient value has an
245 invalid shape
246 """
247 u=self._u
248 for name,val in coefficients.iteritems():
249 shape=util.getShape(val)
250 rank=len(shape)
251 if name=="X" or name=="X_reduced":
252 # DX/Du = del_X/del_u + del_X/del_grad(u)*del_grad(u)/del_u
253 # \ / \ /
254 # B A
255
256 if rank != u.getRank()+1:
257 raise IllegalCoefficientValue("%s must have rank %d"%(name,u.getRank()+1))
258 if not hasattr(val, 'diff'):
259 A=numpy.zeros(shape+u.grad().getShape())
260 B=numpy.zeros(shape+u.getShape())
261 elif u.getRank()==0:
262 X=self._removeFsFromGrad(val)
263 dXdu=X.diff(u)
264 dgdu=u.grad().diff(u)
265 A=numpy.empty(shape+dgdu.getShape(), dtype=object)
266 for i in numpy.ndindex(shape):
267 for j in numpy.ndindex(dgdu.getShape()):
268 tmp=dXdu[i].expand().coeff(dgdu[j])
269 A[i,j]=0 if tmp is None else tmp
270 A=Symbol(A).simplify()
271 B=(dXdu-util.matrix_mult(A,dgdu)).simplify()
272 else: #u.getRank()==1
273 X=self._removeFsFromGrad(val)
274 dXdu=X.diff(u)
275 dgdu=u.grad().diff(u).transpose(2)
276 I,J=shape
277 K,L=u.grad().getShape()
278 A=numpy.empty((I,J,K,L), dtype=object)
279 for i in range(I):
280 for j in range(J):
281 for k in range(K):
282 for l in range(L):
283 if dgdu[k,k,l]==0:
284 A[i,j,k,l]=0
285 else:
286 tmp=dXdu[i,j,k].expand().coeff(dgdu[k,k,l])
287 A[i,j,k,l]=0 if tmp is None else tmp
288 #tmp=sympy.collect(dXdu[i,j,k].expand(), dgdu[k,k,l], evaluate=False, exact=True)
289 #try:
290 # A[i,j,k,l]=tmp[dgdu[k,k,l]]
291 #except KeyError:
292 # A[i,j,k,l]=0
293
294 A=Symbol(A).expand() #.simplify()
295 B=(dXdu-A.tensorProduct(dgdu.transpose(1),2)).expand() #simplify()
296
297 if name=='X_reduced':
298 self.coeffs['A_reduced']=A
299 self.coeffs['B_reduced']=B
300 self.coeffs['X_reduced']=val
301 else:
302 self.coeffs['A']=A
303 self.coeffs['B']=B
304 self.coeffs['X']=val
305 elif name=="Y" or name=="Y_reduced":
306 # DY/Du = del_Y/del_u + del_Y/del_grad(u)*del_grad(u)/del_u
307 # \ / \ /
308 # D C
309 if rank != u.getRank():
310 raise IllegalCoefficientValue("%s must have rank %d"%(name,u.getRank()))
311 if not hasattr(val, 'diff'):
312 C=numpy.zeros(shape+u.grad().getShape())
313 D=numpy.zeros(shape+u.getShape())
314 elif u.getRank()==0:
315 Y=self._removeFsFromGrad(val)
316 dYdu=Y.diff(u)
317 dgdu=u.grad().diff(u)
318 C=numpy.empty(dgdu.getShape(), dtype=object)
319 for j in numpy.ndindex(dgdu.getShape()):
320 tmp=dYdu.expand().coeff(dgdu[j])
321 #tmp=dYdu.collect(dgdu[j])
322 C[j]=0 if tmp is None else tmp
323 C=Symbol(C).simplify()
324 D=(dYdu-util.inner(C,dgdu)).simplify()
325 else: #u.getRank()==1
326 Y=self._removeFsFromGrad(val)
327 dYdu=Y.diff(u)
328 dgdu=u.grad().diff(u).transpose(2)
329 I,=shape
330 J,K=u.grad().getShape()
331 C=numpy.empty((I,J,K), dtype=object)
332 for i in range(I):
333 for j in range(J):
334 for k in range(K):
335 if dgdu[j,j,k]==0:
336 C[i,j,k]=0
337 else:
338 tmp=dYdu[i,j].expand().coeff(dgdu[j,j,k])
339 C[i,j,k]=0 if tmp is None else tmp
340 C=Symbol(C).simplify()
341 D=(dYdu-C.tensorProduct(dgdu.transpose(1),2)).simplify()
342
343 if name=='Y_reduced':
344 self.coeffs['C_reduced']=C
345 self.coeffs['D_reduced']=D
346 self.coeffs['Y_reduced']=val
347 else:
348 self.coeffs['C']=C
349 self.coeffs['D']=D
350 self.coeffs['Y']=val
351 elif name=="y" or name=="y_reduced":
352 y=val
353 if rank != u.getRank():
354 raise IllegalCoefficientValue("%s must have rank %d"%(name,u.getRank()))
355 if not hasattr(y, 'diff'):
356 d=numpy.zeros(u.getShape())
357 else:
358 d=y.diff(u)
359 if name=='y_reduced':
360 self.coeffs['d_reduced']=d
361 self.coeffs['y_reduced']=y
362 else:
363 self.coeffs['d']=d
364 self.coeffs['y']=y
365 elif name=="y_contact" or name=="y_contact_reduced":
366 y_contact=val
367 if rank != u.getRank():
368 raise IllegalCoefficientValue("%s must have rank %d"%(name,u.getRank()))
369 if not hasattr(y_contact, 'diff'):
370 d_contact=numpy.zeros(u.getShape())
371 else:
372 d_contact=y_contact.diff(u)
373 if name=='y_contact_reduced':
374 self.coeffs['d_contact_reduced']=d_contact
375 self.coeffs['y_contact_reduced']=y_contact
376 else:
377 self.coeffs['d_contact']=d_contact
378 self.coeffs['y_contact']=y_contact
379 elif name=="q":
380 if rank != u.getRank():
381 raise IllegalCoefficientValue("q must have rank %d"%u.getRank())
382 self.coeffs['q']=val
383 elif name=="r":
384 if rank != u.getRank():
385 raise IllegalCoefficientValue("r must have rank %d"%u.getRank())
386 self.coeffs['r']=val
387 else:
388 raise IllegalCoefficient("Attempt to set unknown coefficient %s"%name)
389
390 def _removeFsFromGrad(self, sym):
391 """
392 Returns sym with all occurrences grad_n(a,b,c) replaced by grad_n(a,b).
393 That is, all functionspace parameters are removed.
394 This method is used in setValue() in order to find coefficients of
395 grad(u) with the coeff() function.
396 """
397 from esys.escript import symfn
398 gg=sym.atoms(symfn.grad_n)
399 for g in gg:
400 if len(g.args)==3:
401 r=symfn.grad_n(*g.args[:2])
402 sym=sym.subs(g, r)
403 return sym
404

  ViewVC Help
Powered by ViewVC 1.1.26