/[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 3789 - (show annotations)
Tue Jan 31 04:55:05 2012 UTC (7 years, 2 months ago) by caltinay
File MIME type: text/x-python
File size: 16737 byte(s)
Fast forward to latest trunk revision 3788.

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 ttt0=time()
280 for i in range(I):
281 for j in range(J):
282 for k in range(K):
283 for l in range(L):
284 if dgdu[k,k,l]==0:
285 A[i,j,k,l]=0
286 else:
287 tmp=dXdu[i,j,k].expand().coeff(dgdu[k,k,l])
288 A[i,j,k,l]=0 if tmp is None else tmp
289 #tmp=sympy.collect(dXdu[i,j,k].expand(), dgdu[k,k,l], evaluate=False, exact=True)
290 #try:
291 # A[i,j,k,l]=tmp[dgdu[k,k,l]]
292 #except KeyError:
293 # A[i,j,k,l]=0
294
295 self.trace("Computing A took %f seconds."%(time()-ttt0))
296 ttt0=time()
297 A=Symbol(A).expand() #.simplify()
298 self.trace("Expanding A took %f seconds."%(time()-ttt0))
299 ttt0=time()
300 B=(dXdu-A.tensorProduct(dgdu.transpose(1),2)).expand() #simplify()
301 self.trace("Expanding B took %f seconds."%(time()-ttt0))
302
303 if name=='X_reduced':
304 self.coeffs['A_reduced']=A
305 self.coeffs['B_reduced']=B
306 self.coeffs['X_reduced']=val
307 else:
308 self.coeffs['A']=A
309 self.coeffs['B']=B
310 self.coeffs['X']=val
311 elif name=="Y" or name=="Y_reduced":
312 # DY/Du = del_Y/del_u + del_Y/del_grad(u)*del_grad(u)/del_u
313 # \ / \ /
314 # D C
315 if rank != u.getRank():
316 raise IllegalCoefficientValue("%s must have rank %d"%(name,u.getRank()))
317 if not hasattr(val, 'diff'):
318 C=numpy.zeros(shape+u.grad().getShape())
319 D=numpy.zeros(shape+u.getShape())
320 elif u.getRank()==0:
321 Y=self._removeFsFromGrad(val)
322 dYdu=Y.diff(u)
323 dgdu=u.grad().diff(u)
324 C=numpy.empty(dgdu.getShape(), dtype=object)
325 for j in numpy.ndindex(dgdu.getShape()):
326 tmp=dYdu.expand().coeff(dgdu[j])
327 #tmp=dYdu.collect(dgdu[j])
328 C[j]=0 if tmp is None else tmp
329 C=Symbol(C).simplify()
330 D=(dYdu-util.inner(C,dgdu)).simplify()
331 else: #u.getRank()==1
332 Y=self._removeFsFromGrad(val)
333 dYdu=Y.diff(u)
334 dgdu=u.grad().diff(u).transpose(2)
335 I,=shape
336 J,K=u.grad().getShape()
337 C=numpy.empty((I,J,K), dtype=object)
338 ttt0=time()
339 for i in range(I):
340 for j in range(J):
341 for k in range(K):
342 if dgdu[j,j,k]==0:
343 C[i,j,k]=0
344 else:
345 tmp=dYdu[i,j].expand().coeff(dgdu[j,j,k])
346 C[i,j,k]=0 if tmp is None else tmp
347 self.trace("Computing C took %f seconds."%(time()-ttt0))
348 ttt0=time()
349 C=Symbol(C).simplify()
350 self.trace("Simplifying C took %f seconds."%(time()-ttt0))
351 ttt0=time()
352 D=(dYdu-C.tensorProduct(dgdu.transpose(1),2)).simplify()
353 self.trace("Simplifying D took %f seconds."%(time()-ttt0))
354
355 if name=='Y_reduced':
356 self.coeffs['C_reduced']=C
357 self.coeffs['D_reduced']=D
358 self.coeffs['Y_reduced']=val
359 else:
360 self.coeffs['C']=C
361 self.coeffs['D']=D
362 self.coeffs['Y']=val
363 elif name=="y" or name=="y_reduced":
364 y=val
365 if rank != u.getRank():
366 raise IllegalCoefficientValue("%s must have rank %d"%(name,u.getRank()))
367 if not hasattr(y, 'diff'):
368 d=numpy.zeros(u.getShape())
369 else:
370 d=y.diff(u)
371 if name=='y_reduced':
372 self.coeffs['d_reduced']=d
373 self.coeffs['y_reduced']=y
374 else:
375 self.coeffs['d']=d
376 self.coeffs['y']=y
377 elif name=="y_contact" or name=="y_contact_reduced":
378 y_contact=val
379 if rank != u.getRank():
380 raise IllegalCoefficientValue("%s must have rank %d"%(name,u.getRank()))
381 if not hasattr(y_contact, 'diff'):
382 d_contact=numpy.zeros(u.getShape())
383 else:
384 d_contact=y_contact.diff(u)
385 if name=='y_contact_reduced':
386 self.coeffs['d_contact_reduced']=d_contact
387 self.coeffs['y_contact_reduced']=y_contact
388 else:
389 self.coeffs['d_contact']=d_contact
390 self.coeffs['y_contact']=y_contact
391 elif name=="q":
392 if rank != u.getRank():
393 raise IllegalCoefficientValue("q must have rank %d"%u.getRank())
394 self.coeffs['q']=val
395 elif name=="r":
396 if rank != u.getRank():
397 raise IllegalCoefficientValue("r must have rank %d"%u.getRank())
398 self.coeffs['r']=val
399 else:
400 raise IllegalCoefficient("Attempt to set unknown coefficient %s"%name)
401
402 def _removeFsFromGrad(self, sym):
403 """
404 Returns sym with all occurrences grad_n(a,b,c) replaced by grad_n(a,b).
405 That is, all functionspace parameters are removed.
406 This method is used in setValue() in order to find coefficients of
407 grad(u) with the coeff() function.
408 """
409 from esys.escript import symfn
410 gg=sym.atoms(symfn.grad_n)
411 for g in gg:
412 if len(g.args)==3:
413 r=symfn.grad_n(*g.args[:2])
414 sym=sym.subs(g, r)
415 return sym
416

  ViewVC Help
Powered by ViewVC 1.1.26