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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3974 - (hide annotations)
Wed Sep 19 06:27:37 2012 UTC (6 years, 6 months ago) by caltinay
File MIME type: text/x-python
File size: 60489 byte(s)
Update to latest trunk.

1 caltinay 3534 # -*- coding: utf-8 -*-
2    
3     ########################################################
4     #
5 caltinay 3809 # Copyright (c) 2003-2012 by University of Queensland
6 caltinay 3534 # 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 caltinay 3809 __copyright__="""Copyright (c) 2003-2012 by University of Queensland
16 caltinay 3534 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 caltinay 3849 __url__="https://launchpad.net/escript-finley"
22 caltinay 3534
23     """
24     :var __author__: name of author
25     :var __copyright__: copyrights
26     :var __license__: licence agreement
27 caltinay 3849 :var __url__: url entry point on documentation
28 caltinay 3534 :var __version__: version
29     :var __date__: date of the version
30     """
31    
32     import numpy
33 caltinay 3859 import sympy
34 caltinay 3631 from time import time
35 caltinay 3534 from esys.escript.linearPDEs import LinearPDE, IllegalCoefficient, IllegalCoefficientValue
36 caltinay 3857 from esys.escript import util, Data, Symbol, Evaluator, getTotalDifferential, \
37     isSymbol
38 caltinay 3534
39 caltinay 3849 __author__="Cihan Altinay, Lutz Gross"
40 caltinay 3534
41 gross 3839 class iteration_steps_maxReached(Exception):
42 caltinay 3849 """
43     Exception thrown if the maximum number of iteration steps is reached.
44     """
45     pass
46 gross 3839 class DivergenceDetected(Exception):
47 caltinay 3849 """
48     Exception thrown if Newton-Raphson did not converge.
49     """
50     pass
51 gross 3839 class InadmissiblePDEOrdering(Exception):
52 caltinay 3849 """
53     Exception thrown if the ordering of the PDE equations should be revised.
54     """
55     pass
56    
57 caltinay 3855 def concatenateRow(*args):
58     """
59     """
60     if len(args)==1:
61     return args[0]
62    
63     if len(args)>2:
64     return concatenateRow(args[0], concatenateRow(*args[1:]))
65    
66     lshape=util.getShape(args[0])
67     rshape=util.getShape(args[1])
68     if len(lshape)==0:
69     if len(rshape)==0:
70     shape=(2,)
71     res=numpy.zeros(shape, dtype=object)
72     res[0]=args[0]
73     res[1]=args[1]
74     elif len(rshape)==1:
75     shape=(rshape[0]+1,)
76     res=numpy.zeros(shape, dtype=object)
77     res[0]=args[0]
78     res[1:]=args[1]
79     elif len(lshape)==1:
80     if len(rshape)==1:
81     shape=(2,)+lshape
82     res=numpy.zeros(shape, dtype=object)
83     res[0]=args[0]
84     res[1]=args[1]
85     else:
86     shape=(rshape[0]+1,)+lshape
87     res=numpy.zeros(shape, dtype=object)
88     res[0]=args[0]
89     res[1:]=args[1]
90     else:
91     if len(rshape)==1:
92     shape=(lshape[0]+1,)+lshape[1:]
93     res=numpy.zeros(shape, dtype=object)
94     res[:lshape[0]]=args[0]
95     res[lshape[0]]=args[1]
96     else:
97     shape=(lshape[0]+rshape[0],)+lshape[1:]
98     res=numpy.zeros(shape, dtype=object)
99     res[:lshape[0]]=args[0]
100     res[lshape[0]:]=args[1]
101    
102 caltinay 3864 subs=args[0].getDataSubstitutions().copy()
103     subs.update(args[1].getDataSubstitutions())
104 caltinay 3872 dim=args[1].getDim() if args[0].getDim()<0 else args[0].getDim()
105     return Symbol(res, dim=dim, subs=subs)
106 caltinay 3855
107 caltinay 3534 class NonlinearPDE(object):
108     """
109     This class is used to define a general nonlinear, steady, second order PDE
110     for an unknown function *u* on a given domain defined through a `Domain`
111     object.
112    
113     For a single PDE having a solution with a single component the nonlinear
114     PDE is defined in the following form:
115    
116     *-div(X) + Y = 0*
117    
118     where *X*,*Y*=f(*u*,*grad(u)*), *div(F)* denotes the divergence of *F* and
119     *grad(F)* is the spatial derivative of *F*.
120    
121     The coefficients *X* (rank 1) and *Y* (scalar) have to be specified through
122     `Symbol` objects.
123    
124     The following natural boundary conditions are considered:
125    
126     *n[j]*X[j] - y = 0*
127    
128     where *n* is the outer normal field. Notice that the coefficient *X*
129     is defined in the PDE. The coefficient *y* is a scalar `Symbol`.
130    
131     Constraints for the solution prescribe the value of the solution at certain
132     locations in the domain. They have the form
133    
134 gross 3844 *u=r* where *q>0*
135 caltinay 3534
136 gross 3844 *r* and *q* are each scalar where *q* is the characteristic function
137 caltinay 3534 defining where the constraint is applied. The constraints override any
138     other condition set by the PDE or the boundary condition.
139    
140     For a system of PDEs and a solution with several components, *u* is rank
141     one, while the PDE coefficient *X* is rank two and *Y* and *y* is rank one.
142    
143     The PDE is solved by linearising the coefficients and iteratively solving
144     the corresponding linear PDE until the error is smaller than a tolerance
145     or a maximum number of iterations is reached.
146    
147     Typical usage::
148    
149     u = Symbol('u', dim=dom.getDim())
150     p = NonlinearPDE(dom, u)
151 gross 3811 p.setValue(X=grad(u), Y=5*u)
152 caltinay 3534 v = p.getSolution(u=0.)
153     """
154 gross 3839 DEBUG0=0 # no debug info
155     DEBUG1=1 # info on Newton-Raphson solver printed
156 caltinay 3849 DEBUG2=2 # in addition info from linear solver is printed
157     DEBUG3=3 # in addition info on linearization is printed
158     DEBUG4=4 # in addition info on LinearPDE handling is printed
159 caltinay 3859
160 gross 3856 ORDER=0
161 caltinay 3849 def __init__(self, domain, u, debug=DEBUG0):
162 caltinay 3534 """
163     Initializes a new nonlinear PDE.
164    
165     :param domain: domain of the PDE
166     :type domain: `Domain`
167     :param u: The symbol for the unknown PDE function u.
168     :type u: `Symbol`
169 gross 3839 :param debug: level of debug information to be printed
170 caltinay 3534 """
171 gross 3856 self.__COEFFICIENTS = [ "X", "X_reduced", "Y", "Y_reduced", "y", "y_reduced", "y_contact", "y_contact_reduced", "y_dirac"]
172 caltinay 3859
173 gross 3844 self._r=Data()
174 gross 3839 self._set_coeffs={}
175     self._unknown=u
176 caltinay 3534 self._debug=debug
177 gross 3839 self._rtol=1e-6
178     self._atol=0.
179     self._iteration_steps_max=100
180 gross 3877 self._omega_min=0.0001
181 gross 3839 self._quadratic_convergence_limit=0.2
182     self._simplified_newton_limit=0.1
183    
184 gross 3811 self.dim = domain.getDim()
185 caltinay 3534 if u.getRank()==0:
186     numEquations=1
187     else:
188     numEquations=u.getShape()[0]
189     numSolutions=numEquations
190 gross 3839 self._lpde=LinearPDE(domain,numEquations,numSolutions,self._debug > self.DEBUG4 )
191 caltinay 3849
192 caltinay 3534 def __str__(self):
193     """
194 caltinay 3849 Returns the string representation of the PDE
195 caltinay 3534
196     :return: a simple representation of the PDE
197     :rtype: ``str``
198     """
199 caltinay 3849 return "<%s %d>"%(self.__class__.__name__, id(self))
200    
201 caltinay 3831 def getUnknownSymbol(self):
202 gross 3811 """
203 caltinay 3831 Returns the symbol of the PDE unknown
204 caltinay 3534
205 caltinay 3831 :return: the symbol of the PDE unknown
206 gross 3811 :rtype: `Symbol`
207     """
208 gross 3839 return self._unknown
209 caltinay 3849
210 gross 3839 def trace1(self, text):
211 caltinay 3534 """
212 caltinay 3849 Prints the text message if the debug level is greater than DEBUG0
213 caltinay 3534
214     :param text: message to be printed
215     :type text: ``string``
216     """
217 gross 3839 if self._debug > self.DEBUG0:
218 caltinay 3534 print("%s: %s"%(str(self), text))
219 caltinay 3849
220 gross 3839 def trace3(self, text):
221     """
222 caltinay 3849 Prints the text message if the debug level is greater than DEBUG3
223 caltinay 3534
224 gross 3839 :param text: message to be printed
225     :type text: ``string``
226     """
227     if self._debug > self.DEBUG2:
228     print("%s: %s"%(str(self), text))
229    
230 caltinay 3831 def getLinearSolverOptions(self):
231     """
232 caltinay 3849 Returns the options of the linear PDE solver class
233 caltinay 3831 """
234     return self._lpde.getSolverOptions()
235 caltinay 3849
236 gross 3844 def getLinearPDE(self):
237     """
238 caltinay 3849 Returns the linear PDE used to calculate the Newton-Raphson update
239    
240 gross 3844 :rtype: `LinearPDE`
241     """
242     return self._lpde
243 caltinay 3849
244 caltinay 3831 def setOptions(self, **opts):
245     """
246     Allows setting options for the nonlinear PDE.
247     The supported options are:
248     'tolerance' - error tolerance for the Newton method
249 caltinay 3849 'iteration_steps_max' - maximum number of Newton iterations
250     'omega_min' - minimum relaxation factor
251     'atol' - solution norms less than 'atol' are assumed to be 'atol'.
252     This can be useful if one of your solutions is expected to
253     be zero.
254     'quadratic_convergence_limit' - if the norm of the Newton-Raphson
255     correction is reduced by less than 'quadratic_convergence_limit'
256     between two iteration steps quadratic convergence is assumed
257     'simplified_newton_limit' - if the norm of the defect is reduced by
258     less than 'simplified_newton_limit' between two iteration
259     steps and quadratic convergence is detected the iteration
260     swiches to the simplified Newton-Raphson scheme
261 caltinay 3831 """
262     for key in opts:
263     if key=='tolerance':
264 gross 3839 self._rtol=opts[key]
265     elif key=='iteration_steps_max':
266     self._iteration_steps_max=opts[key]
267     elif key=='omega_min':
268     self._omega_min=opts[key]
269     elif key=='atol':
270 caltinay 3849 self._atol=opts[key]
271 gross 3839 elif key=='quadratic_convergence_limit':
272 caltinay 3849 self._quadratic_convergence_limit=opts[key]
273 gross 3839 elif key=='simplified_newton_limit':
274 caltinay 3849 self._simplified_newton_limit=opts[key]
275 caltinay 3831 else:
276     raise KeyError("Invalid option '%s'"%key)
277    
278 caltinay 3534 def getSolution(self, **subs):
279     """
280     Returns the solution of the PDE.
281     :param subs: Substitutions for all symbols used in the coefficients
282 caltinay 3849 including the initial value for the unknown *u*.
283 caltinay 3534 :return: the solution
284     :rtype: `Data`
285     """
286 caltinay 3859 # get the initial value for the iteration process
287    
288     # collect components of unknown in u_syms
289     u_syms=[]
290     simple_u=False
291     for i in numpy.ndindex(self._unknown.getShape()):
292     u_syms.append(Symbol(self._unknown[i]).atoms(sympy.Symbol).pop().name)
293     if len(set(u_syms))==1: simple_u=True
294 gross 3858
295 caltinay 3859 e=Evaluator(self._unknown)
296     for sym in u_syms:
297     if not subs.has_key(sym):
298     raise KeyError("Initial value for '%s' missing."%sym)
299     if not isinstance(subs[sym], Data):
300     subs[sym]=Data(subs[sym], self._lpde.getFunctionSpaceForSolution())
301     e.subs(**{sym:subs[sym]})
302     ui=e.evaluate()
303 caltinay 3534
304 gross 3839 # modify ui so it meets the constraints:
305     q=self._lpde.getCoefficient("q")
306     if not q.isEmpty():
307 caltinay 3849 if hasattr(self, "_r"):
308     r=self._r
309 gross 3856 if isSymbol(r):
310 caltinay 3849 r=Evaluator(r).evaluate(**subs)
311     elif not isinstance(r, Data):
312     r=Data(r, self._lpde.getFunctionSpaceForSolution())
313     elif r.isEmpty():
314     r=0
315     else:
316     r=0
317     ui = q * r + (1-q) * ui
318    
319 caltinay 3534 # separate symbolic expressions from other coefficients
320 caltinay 3831 constants={}
321 caltinay 3534 expressions={}
322 gross 3839 for n, e in self._set_coeffs.items():
323 gross 3856 if isSymbol(e):
324 caltinay 3534 expressions[n]=e
325     else:
326 caltinay 3831 constants[n]=e
327 caltinay 3534
328 caltinay 3831 # set constant PDE values now
329     self._lpde.setValue(**constants)
330 gross 3839 self._lpde.getSolverOptions().setAbsoluteTolerance(0.)
331     self._lpde.getSolverOptions().setVerbosity(self._debug > self.DEBUG1)
332 caltinay 3849 #=====================================================================
333 caltinay 3820 # perform Newton iterations until error is small enough or
334     # maximum number of iterations reached
335     n=0
336 gross 3839 omega=1. # relaxation factor
337     use_simplified_Newton=False
338     defect_norm=None
339     delta_norm=None
340     converged=False
341 caltinay 3849
342 caltinay 3859 #subs[u_sym]=ui
343     if simple_u:
344     subs[u_syms[0]]=ui
345     else:
346     for i in range(len(u_syms)):
347     subs[u_syms[i]]=ui[i]
348    
349 gross 3839 while not converged:
350     if n > self._iteration_steps_max:
351 caltinay 3849 raise iteration_steps_maxReached("maximum number of iteration steps reached, giving up.")
352     self.trace1(5*"="+" iteration step %d "%n + 5*"=")
353 gross 3839 # calculate the correction delta_u
354     if n == 0:
355 caltinay 3849 self._updateLinearPDE(expressions, subs)
356 gross 3839 defect_norm=self._getDefectNorm(self._lpde.getRightHandSide())
357     LINTOL=0.1
358     else:
359 caltinay 3849 if not use_simplified_Newton:
360     self._updateMatrix(expressions, subs)
361 gross 3839 if q_u == None:
362 caltinay 3849 LINTOL = 0.1 * min(qtol/defect_norm)
363 gross 3839 else:
364 caltinay 3849 LINTOL = 0.1 * max( q_u**2, min(qtol/defect_norm) )
365 gross 3839 LINTOL=max(1e-4,min(LINTOL, 0.1))
366 gross 3860 #LINTOL=1.e-5
367 gross 3839 self._lpde.getSolverOptions().setTolerance(LINTOL)
368 caltinay 3849 self.trace1("PDE is solved with rel. tolerance = %e"%LINTOL)
369 gross 3839 delta_u=self._lpde.getSolution()
370 caltinay 3849
371 gross 3839 #check for reduced defect:
372 caltinay 3849 omega=min(2*omega, 1.) # raise omega
373 gross 3839 defect_reduced=False
374 gross 3877 ui_old=ui
375 gross 3839 while not defect_reduced:
376 gross 3877 ui=ui_old-delta_u * omega
377 caltinay 3859 if simple_u:
378     subs[u_syms[0]]=ui
379     else:
380     for i in range(len(u_syms)):
381     subs[u_syms[i]]=ui[i]
382 caltinay 3849 self._updateRHS(expressions, subs)
383 gross 3839 new_defect_norm=self._getDefectNorm(self._lpde.getRightHandSide())
384 gross 3860 defect_reduced=False
385     for i in xrange(len( new_defect_norm)):
386     if new_defect_norm[i] < defect_norm[i]: defect_reduced=True
387    
388     #print new_defect_norm
389     #q_defect=max(self._getSafeRatio(new_defect_norm, defect_norm))
390 caltinay 3849 # if defect_norm==0 and new_defect_norm!=0
391     # this will be util.DBLE_MAX
392 gross 3860 #self.trace1("Defect reduction = %e with relaxation factor %e."%(q_defect, omega))
393     if not defect_reduced:
394 gross 3839 omega*=0.5
395     if omega < self._omega_min:
396 caltinay 3849 raise DivergenceDetected("Underrelaxtion failed to reduce defect, giving up.")
397 gross 3860 self.trace1("Defect reduction with relaxation factor %e."%(omega, ))
398 caltinay 3849 delta_norm, delta_norm_old = self._getSolutionNorm(delta_u) * omega, delta_norm
399 gross 3839 defect_norm, defect_norm_old = new_defect_norm, defect_norm
400     u_norm=self._getSolutionNorm(ui, atol=self._atol)
401 caltinay 3849 # set the tolerance on equation level:
402     qtol=self._getSafeRatio(defect_norm_old * u_norm * self._rtol, delta_norm)
403 gross 3839 # if defect_norm_old==0 and defect_norm_old!=0 this will be util.DBLE_MAX
404     # -> the ordering of the equations is not appropriate.
405 caltinay 3849 # if defect_norm_old==0 and defect_norm_old==0 this is zero so
406     # convergence can happen for defect_norm==0 only.
407 gross 3839 if not max(qtol)<util.DBLE_MAX:
408 caltinay 3849 raise InadmissiblePDEOrdering("Review ordering of PDE equations.")
409     # check stopping criteria
410     if not delta_norm_old == None:
411     q_u=max(self._getSafeRatio(delta_norm, delta_norm_old))
412     # if delta_norm_old==0 and delta_norm!=0
413     # this will be util.DBLE_MAX
414     if q_u <= self._quadratic_convergence_limit and not omega<1. :
415     quadratic_convergence=True
416     self.trace1("Quadratic convergence detected (rate %e <= %e)"%(q_u, self._quadratic_convergence_limit))
417 caltinay 3820
418 caltinay 3849 converged = all( [ defect_norm[i] <= qtol[i] for i in range(len(qtol)) ])
419    
420     else:
421     self.trace1("No quadratic convergence detected (rate %e > %e, omega=%e)"%(q_u, self._quadratic_convergence_limit,omega ))
422     quadratic_convergence=False
423     converged=False
424 gross 3839 else:
425 caltinay 3849 q_u=None
426     converged=False
427     quadratic_convergence=False
428 gross 3839 if self._debug > self.DEBUG0:
429     for i in range(len(u_norm)):
430 caltinay 3849 self.trace1("Component %s: u: %e, du: %e, defect: %e, qtol: %e"%(i, u_norm[i], delta_norm[i], defect_norm[i], qtol[i]))
431     if converged: self.trace1("Iteration has converged.")
432 gross 3839 # Can we switch to simplified Newton?
433     if quadratic_convergence:
434     q_defect=max(self._getSafeRatio(defect_norm, defect_norm_old))
435     if q_defect < self._simplified_newton_limit:
436     use_simplified_Newton=True
437     self.trace1("Simplified Newton-Raphson is applied (rate %e < %e)."%(q_defect, self._simplified_newton_limit))
438     n+=1
439 caltinay 3849 self.trace1(5*"="+" Newton-Raphson iteration completed after %d steps "%n + 5*"=")
440 caltinay 3820 return ui
441 caltinay 3849
442 gross 3839 def _getDefectNorm(self, f):
443     """
444     calculates the norm of the defect ``f``
445 caltinay 3849
446 gross 3839 :param f: defect vector
447     :type f: `Data` of rank 0 or 1.
448     :return: component-by-component norm of ``f``
449     :rtype: `numpy.array` of rank 1
450     :raise ValueError: if shape if ``f`` is incorrect.
451 caltinay 3849
452 gross 3839 """
453     # this still needs some work!!!
454     out=[]
455     s=f.getShape()
456     if len(s) == 0:
457 caltinay 3849 out.append(util.Lsup(f))
458 gross 3839 elif len(s) == 1:
459 caltinay 3849 for i in range(s[0]):
460 gross 3839 out.append(util.Lsup(f[i]))
461     else:
462 caltinay 3849 raise ValueError("Illegal shape of defect vector: %s"%s)
463 gross 3839 return numpy.array(out)
464 caltinay 3849
465 gross 3839 def _getSolutionNorm(self, u, atol=0.):
466     """
467 caltinay 3849 calculates the norm of the solution ``u``
468    
469     :param u: solution vector
470 gross 3839 :type f: `Data` of rank 0 or 1.
471     :return: component-by-component norm of ``u``
472     :rtype: `numpy.array` of rank 1
473     :raise ValueError: if shape if ``u`` is incorrect.
474 caltinay 3849
475 gross 3839 """
476     out=[]
477     s=u.getShape()
478     if len(s) == 0:
479 caltinay 3849 out.append(max(util.Lsup(u),atol) )
480 gross 3839 elif len(s) == 1:
481 caltinay 3849 for i in range(s[0]):
482 gross 3839 out.append(max(util.Lsup(u[i]), atol))
483     else:
484 caltinay 3849 raise ValueError("Illegal shape of solution: %s"%s)
485 gross 3839 return numpy.array(out)
486 caltinay 3849
487 gross 3839 def _getSafeRatio(self, a , b):
488     """
489 caltinay 3849 Returns the component-by-component ratio of ''a'' and ''b''
490    
491     If for a component ``i`` the values ``a[i]`` and ``b[i]`` are both
492     equal to zero their ratio is set to zero.
493     If ``b[i]`` equals zero but ``a[i]`` is positive the ratio is set to
494     `util.DBLE_MAX`.
495    
496 gross 3839 :param a: numerator
497     :param b: denominator
498     :type a: `numpy.array` of rank 1 with non-negative entries.
499     :type b: `numpy.array` of rank 1 with non-negative entries.
500 caltinay 3849 :return: ratio of ``a`` and ``b``
501 gross 3839 :rtype: `numpy.array`
502     """
503     out=0.
504     if a.shape !=b.shape:
505     raise ValueError("shapes must match.")
506     s=a.shape
507     if len(s) != 1:
508 caltinay 3849 raise ValueError("rank one is expected.")
509     out=numpy.zeros(s)
510     for i in range(s[0]):
511     if abs(b[i]) > 0:
512     out[i]=a[i]/b[i]
513     elif abs(a[i]) > 0:
514     out[i] = util.DBLE_MAX
515     else:
516     out[i] = 0
517 gross 3839 return out
518 caltinay 3849
519 gross 3811 def getNumSolutions(self):
520     """
521     Returns the number of the solution components
522     :rtype: ``int``
523     """
524 gross 3839 s=self._unknown.getShape()
525 gross 3811 if len(s) > 0:
526 caltinay 3832 return s[0]
527 gross 3811 else:
528 caltinay 3832 return 1
529 caltinay 3831
530 gross 3811 def getShapeOfCoefficient(self,name):
531     """
532 caltinay 3849 Returns the shape of the coefficient ``name``
533 caltinay 3534
534 gross 3811 :param name: name of the coefficient enquired
535     :type name: ``string``
536     :return: the shape of the coefficient ``name``
537     :rtype: ``tuple`` of ``int``
538     :raise IllegalCoefficient: if ``name`` is not a coefficient of the PDE
539 caltinay 3534 """
540 gross 3811 numSol=self.getNumSolutions()
541     dim = self.dim
542     if name=="X" or name=="X_reduced":
543 caltinay 3849 if numSol > 1:
544     return (numSol,dim)
545 caltinay 3832 else:
546     return (dim,)
547 gross 3839 elif name=="r" or name == "q" :
548 caltinay 3849 if numSol > 1:
549     return (numSol,)
550 caltinay 3832 else:
551     return ()
552 gross 3811 elif name=="Y" or name=="Y_reduced":
553 caltinay 3832 if numSol > 1:
554 caltinay 3849 return (numSol,)
555 caltinay 3832 else:
556     return ()
557 gross 3811 elif name=="y" or name=="y_reduced":
558 caltinay 3849 if numSol > 1:
559     return (numSol,)
560 gross 3811 else:
561 caltinay 3849 return ()
562 gross 3811 elif name=="y_contact" or name=="y_contact_reduced":
563 caltinay 3849 if numSol > 1:
564     return (numSol,)
565 gross 3811 else:
566 caltinay 3849 return ()
567 gross 3839 elif name=="y_dirac":
568 caltinay 3849 if numSol > 1:
569     return (numSol,)
570 gross 3839 else:
571 caltinay 3849 return ()
572 gross 3811 else:
573 caltinay 3849 raise IllegalCoefficient("Attempt to request unknown coefficient %s"%name)
574 caltinay 3534
575 gross 3811 def createCoefficient(self, name):
576     """
577 caltinay 3849 Creates a new coefficient ``name`` as Symbol
578 gross 3811
579 caltinay 3534 :param name: name of the coefficient requested
580     :type name: ``string``
581     :return: the value of the coefficient
582 gross 3845 :rtype: `Symbol` or `Data` (for name = "q")
583 caltinay 3534 :raise IllegalCoefficient: if ``name`` is not a coefficient of the PDE
584 caltinay 3849 """
585 gross 3844 if name == "q":
586 caltinay 3849 return self._lpde.createCoefficient("q")
587 gross 3839 else:
588     s=self.getShapeOfCoefficient(name)
589 gross 3844 return Symbol(name, s, dim=self.dim)
590 caltinay 3831
591 gross 3811 def getCoefficient(self, name):
592 caltinay 3534 """
593 gross 3811 Returns the value of the coefficient ``name`` as Symbol
594 caltinay 3534
595 gross 3811 :param name: name of the coefficient requested
596 caltinay 3534 :type name: ``string``
597 gross 3811 :return: the value of the coefficient
598     :rtype: `Symbol`
599     :raise IllegalCoefficient: if ``name`` is not a coefficient of the PDE
600 caltinay 3534 """
601 gross 3839 if self._set_coeffs.has_key(name):
602     return self._set_coeffs[name]
603 gross 3844 elif name == "r":
604 caltinay 3849 if hasattr(self, "_r"):
605     return self._r
606     else:
607     raise IllegalCoefficient("Attempt to request undefined coefficient %s"%name)
608 gross 3844 elif name == "q":
609 caltinay 3849 return self._lpde.getCoefficient("q")
610 caltinay 3831 else:
611     raise IllegalCoefficient("Attempt to request undefined coefficient %s"%name)
612 caltinay 3534
613     def setValue(self,**coefficients):
614     """
615     Sets new values to one or more coefficients.
616    
617     :param coefficients: new values assigned to coefficients
618 gross 3844
619     :param coefficients: new values assigned to coefficients
620 caltinay 3534 :keyword X: value for coefficient ``X``
621     :type X: `Symbol` or any type that can be cast to a `Data` object
622     :keyword Y: value for coefficient ``Y``
623     :type Y: `Symbol` or any type that can be cast to a `Data` object
624     :keyword y: value for coefficient ``y``
625     :type y: `Symbol` or any type that can be cast to a `Data` object
626     :keyword y_contact: value for coefficient ``y_contact``
627 gross 3839 :type y_contact: `Symbol` or any type that can be cast to a `Data` object
628     :keyword y_dirac: value for coefficient ``y_dirac``
629     :type y_dirac: `Symbol` or any type that can be cast to a `Data` object
630     :keyword q: mask for location of constraint
631 caltinay 3849 :type q: any type that can be cast to a `Data` object
632 gross 3839 :keyword r: value of solution prescribed by constraint
633 caltinay 3849 :type r: `Symbol` or any type that can be cast to a `Data` object
634 caltinay 3534 :raise IllegalCoefficient: if an unknown coefficient keyword is used
635     :raise IllegalCoefficientValue: if a supplied coefficient value has an
636     invalid shape
637 caltinay 3849 """
638    
639 gross 3839 u=self._unknown
640 caltinay 3534 for name,val in coefficients.iteritems():
641     shape=util.getShape(val)
642 gross 3811 if not shape == self.getShapeOfCoefficient(name):
643 caltinay 3855 raise IllegalCoefficientValue("%s has shape %s but must have shape %s"%(name, shape, self.getShapeOfCoefficient(name)))
644 caltinay 3534 rank=len(shape)
645 gross 3839 if name == "q":
646 caltinay 3849 self._lpde.setValue(q=val)
647 gross 3839 elif name == "r":
648 gross 3844 self._r=val
649 gross 3839 elif name=="X" or name=="X_reduced":
650 caltinay 3534 if rank != u.getRank()+1:
651 caltinay 3617 raise IllegalCoefficientValue("%s must have rank %d"%(name,u.getRank()+1))
652 caltinay 3831 T0=time()
653 caltinay 3854 B,A=getTotalDifferential(val, u, 1)
654 caltinay 3617 if name=='X_reduced':
655 gross 3839 self.trace3("Computing A_reduced, B_reduced took %f seconds."%(time()-T0))
656     self._set_coeffs['A_reduced']=A
657     self._set_coeffs['B_reduced']=B
658     self._set_coeffs['X_reduced']=val
659 caltinay 3617 else:
660 gross 3839 self.trace3("Computing A, B took %f seconds."%(time()-T0))
661     self._set_coeffs['A']=A
662     self._set_coeffs['B']=B
663     self._set_coeffs['X']=val
664 caltinay 3617 elif name=="Y" or name=="Y_reduced":
665 caltinay 3534 if rank != u.getRank():
666 caltinay 3617 raise IllegalCoefficientValue("%s must have rank %d"%(name,u.getRank()))
667 caltinay 3831 T0=time()
668 caltinay 3854 D,C=getTotalDifferential(val, u, 1)
669 caltinay 3617 if name=='Y_reduced':
670 gross 3839 self.trace3("Computing C_reduced, D_reduced took %f seconds."%(time()-T0))
671     self._set_coeffs['C_reduced']=C
672     self._set_coeffs['D_reduced']=D
673     self._set_coeffs['Y_reduced']=val
674 caltinay 3617 else:
675 gross 3839 self.trace3("Computing C, D took %f seconds."%(time()-T0))
676     self._set_coeffs['C']=C
677     self._set_coeffs['D']=D
678     self._set_coeffs['Y']=val
679 caltinay 3832 elif name in ("y", "y_reduced", "y_contact", "y_contact_reduced", \
680 gross 3839 "y_dirac"):
681 caltinay 3534 y=val
682     if rank != u.getRank():
683 caltinay 3617 raise IllegalCoefficientValue("%s must have rank %d"%(name,u.getRank()))
684 caltinay 3534 if not hasattr(y, 'diff'):
685     d=numpy.zeros(u.getShape())
686     else:
687     d=y.diff(u)
688 gross 3839 self._set_coeffs[name]=y
689     self._set_coeffs['d'+name[1:]]=d
690 caltinay 3534 else:
691 caltinay 3617 raise IllegalCoefficient("Attempt to set unknown coefficient %s"%name)
692 caltinay 3859
693 gross 3856 def getSensitivity(self, f, g=None, **subs):
694 caltinay 3831 """
695 caltinay 3859 Calculates the sensitivity of the solution of an input factor ``f``
696     in direction ``g``.
697    
698     :param f: the input factor to be investigated. ``f`` may be of rank 0
699     or 1.
700 gross 3856 :type f: `Symbol`
701 caltinay 3859 :param g: the direction(s) of change.
702     If not present, it is *g=eye(n)* where ``n`` is the number of
703     components of ``f``.
704 gross 3856 :type g: ``list`` or single of ``float``, `ndarray` or `Data`.
705     :param subs: Substitutions for all symbols used in the coefficients
706 caltinay 3859 including unknown *u* and the input factor ``f`` to be
707     investigated
708 gross 3856 :return: the sensitivity
709 caltinay 3859 :rtype: `Data` with shape *u.getShape()+(len(g),)* if *len(g)>1* or
710     *u.getShape()* if *len(g)==1*
711 caltinay 3831 """
712 gross 3856 s_f=f.getShape()
713     if len(s_f) == 0:
714 caltinay 3859 len_f=1
715     elif len(s_f) == 1:
716     len_f=s_f[0]
717     else:
718     raise ValueError("rank of input factor must be zero or one.")
719    
720     if not g == None:
721     if len(s_f) == 0:
722     if not isinstance(g, list): g=[g]
723     else:
724     if isinstance(g, list):
725     if len(g) == 0:
726     raise ValueError("no direction given.")
727     if len(getShape(g[0])) == 0: g=[g] # if g[0] is a scalar we assume that the list g is to be interprested a data object
728 gross 3856 else:
729 caltinay 3859 g=[g]
730     # at this point g is a list of directions:
731     len_g=len(g)
732     for g_i in g:
733     if not getShape(g_i) == s_f:
734     raise ValueError("shape of direction (=%s) must match rank of input factor (=%s)"%(getShape(g_i) , s_f) )
735     else:
736     len_g=len_f
737    
738     #*** at this point g is a list of direction or None and len_g is the
739     # number of directions to be investigated.
740    
741     # now we make sure that the operator in the lpde is set (it is the same
742     # as for the Newton-Raphson scheme)
743 gross 3856 # if the solution etc are cached this could be omitted:
744     constants={}
745     expressions={}
746     for n, e in self._set_coeffs.items():
747 caltinay 3859 if n not in self.__COEFFICIENTS:
748     if isSymbol(e):
749     expressions[n]=e
750     else:
751     constants[n]=e
752 gross 3856 self._lpde.setValue(**constants)
753     self._updateMatrix(self, expressions, subs)
754     #=====================================================================
755     self._lpde.getSolverOptions().setAbsoluteTolerance(0.)
756     self._lpde.getSolverOptions().setTolerance(self._rtol)
757     self._lpde.getSolverOptions().setVerbosity(self._debug > self.DEBUG1)
758     #=====================================================================
759     #
760 caltinay 3859 # evaluate the derivatives of X, etc with respect to f:
761     #
762 gross 3856 ev=Evaluator()
763     names=[]
764     if hasattr(self, "_r"):
765 caltinay 3859 if isSymbol(self._r):
766     names.append('r')
767     ev.addExpression(self._r.diff(f))
768     for n in self._set_coeffs.keys():
769     if n in self.__COEFFICIENTS and isSymbol(self._set_coeffs[n]):
770 gross 3856 if n=="X" or n=="X_reduced":
771     T0=time()
772     B,A=getTotalDifferential(self._set_coeffs[n], f, 1)
773     if n=='X_reduced':
774     self.trace3("Computing A_reduced, B_reduced took %f seconds."%(time()-T0))
775     names.append('A_reduced'); ev.addExpression(A)
776     names.append('B_reduced'); ev.addExpression(B)
777 caltinay 3859 else:
778     self.trace3("Computing A, B took %f seconds."%(time()-T0))
779 gross 3856 names.append('A'); ev.addExpression(A)
780     names.append('B'); ev.addExpression(B)
781     elif n=="Y" or n=="Y_reduced":
782 caltinay 3859 T0=time()
783     D,C=getTotalDifferential(self._set_coeffs[n], f, 1)
784 gross 3856 if n=='Y_reduced':
785     self.trace3("Computing C_reduced, D_reduced took %f seconds."%(time()-T0))
786     names.append('C_reduced'); ev.addExpression(C)
787     names.append('D_reduced'); ev.addExpression(D)
788     else:
789 caltinay 3859 self.trace3("Computing C, D took %f seconds."%(time()-T0))
790 gross 3856 names.append('C'); ev.addExpression(C)
791     names.append('D'); ev.addExpression(D)
792     elif n in ("y", "y_reduced", "y_contact", "y_contact_reduced", "y_dirac"):
793 caltinay 3859 names.append('d'+name[1:]); ev.addExpression(self._set_coeffs[name].diff(f))
794 gross 3856 relevant_symbols['d'+name[1:]]=self._set_coeffs[name].diff(f)
795     res=ev.evaluate()
796     if len(names)==1: res=[res]
797 caltinay 3864 self.trace3("RHS expressions evaluated in %f seconds."%(time()-T0))
798 gross 3856 for i in range(len(names)):
799     self.trace3("util.Lsup(%s)=%s"%(names[i],util.Lsup(res[i])))
800     coeffs_f=dict(zip(names,res))
801 caltinay 3859 #
802 caltinay 3831
803 caltinay 3859 # now we are ready to calculate the right hand side coefficients into
804     # args by multiplication with g and grad(g).
805     if len_g >1:
806     if self.getNumSolutions() == 1:
807     u_g=Data(0., (len_g,), self._lpde.getFunctionSpaceForSolution())
808     else:
809     u_g=Data(0., (self.getNumSolutions(), len_g), self._lpde.getFunctionSpaceForSolution())
810    
811     for i in xrange(len_g):
812     # reset coefficients may be set at previous calls:
813     args={}
814     for n in self.__COEFFICIENTS: args[n]=Data()
815     args['r']=Data()
816     if g == None: # g_l=delta_{il} and len_f=len_g
817     for n,v in coeffs_f:
818     name=None
819     if len_f > 1:
820     val=v[:,i]
821     else:
822     val=v
823     if n.startswith("d"):
824     name='y'+n[1:]
825     elif n.startswith("D"):
826     name='Y'+n[1:]
827     elif n.startswith("r"):
828     name='r'+n[1:]
829     if name: args[name]=val
830 gross 3856 else:
831 caltinay 3859 g_i=g[i]
832     for n,v in coeffs_f:
833     name=None
834     if n.startswith("d"):
835     name = 'y'+n[1:]
836     val = self.__mm(v, g_i)
837     elif n.startswith("r"):
838     name= 'r'
839     val = self.__mm(v, g_i)
840     elif n.startswith("D"):
841     name = 'Y'+n[1:]
842     val = self.__mm(v, g_i)
843     elif n.startswith("B") and isinstance(g_i, Data):
844     name = 'Y'+n[1:]
845     val = self.__mm(v, grad(g_i))
846     elif n.startswith("C"):
847     name = 'X'+n[1:]
848     val = matrix_multiply(v, g_i)
849     elif n.startswith("A") and isinstance(g_i, Data):
850     name = 'X'+n[1:]
851     val = self.__mm(v, grad(g_i))
852     if name:
853     if args.has_key(name):
854     args[name]+=val
855     else:
856     args[name]=val
857     self._lpde.setValue(**args)
858     u_g_i=self._lpde.getSolution()
859    
860     if len_g >1:
861     if self.getNumSolutions() == 1:
862     u_g[i]=-u_g_i
863     else:
864     u_g[:,i]=-u_g_i
865     else:
866     u_g=-u_g_i
867    
868 gross 3856 return u_g
869    
870     def __mm(self, m, v):
871 caltinay 3859 """
872     a sligtly crude matrix*matrix multiplication
873     m is A-coefficient, u is vector, v is grad vector: A_ijkl*v_kl
874     m is B-coefficient, u is vector, v is vector: B_ijk*v_k
875     m is C-coefficient, u is vector, v is grad vector: C_ikl*v_kl
876     m is D-coefficient, u is vector, v is vector: D_ij*v_j
877     m is A-coefficient, u is scalar, v is grad vector: A_jkl*v_kl
878     m is B-coefficient, u is scalar, v is vector: B_jk*v_k
879     m is C-coefficient, u is scalar, v is grad vector: C_kl*v_kl
880     m is D-coefficient, u is scalar, v is vector: D_j*v_j
881     m is A-coefficient, u is vector, v is grad scalar: A_ijl*v_l
882     m is B-coefficient, u is vector, v is scalar: B_ij*v
883     m is C-coefficient, u is vector, v is grad scalar: C_il*v_l
884     m is D-coefficient, u is vector, v is scalar: D_i*v
885     m is A-coefficient, u is scalar, v is grad scalar: A_jl*v_l
886     m is B-coefficient, u is scalar, v is scalar: B_j*v
887     m is C-coefficient, u is scalar, v is grad scalar: C_l*v_l
888     m is D-coefficient, u is scalar, v is scalar: D*v
889     """
890     s_m=getShape(m)
891     s_v=getShape(v)
892 gross 3856
893 caltinay 3859 # m is B-coefficient, u is vector, v is scalar: B_ij*v
894     # m is D-coefficient, u is vector, v is scalar: D_i*v
895     # m is B-coefficient, u is scalar, v is scalar: B_j*v
896     # m is D-coefficient, u is scalar, v is scalar: D*v
897     if s_m == () or s_v == ():
898     return m*v
899 gross 3856
900 caltinay 3859 # m is D-coefficient, u is scalar, v is vector: D_j*v_j
901     # m is C-coefficient, u is scalar, v is grad scalar: C_l*v_l
902     if len(s_m) == 1:
903     return inner(m,v)
904 gross 3856
905 caltinay 3859 # m is D-coefficient, u is vector, v is vector: D_ij*v_j
906     # m is B-coefficient, u is scalar, v is vector: B_jk*v_k
907     # m is C-coefficient, u is vector, v is grad scalar: C_il*v_l
908     # m is A-coefficient, u is scalar, v is grad scalar: A_jl*v_l
909     if len(s_m) == 2 and len(s_v) == 1:
910     return matrix_mult(m,v)
911    
912     # m is C-coefficient, u is scalar, v is grad vector: C_kl*v_kl
913     if len(s_m) == 2 and len(s_v) == 2:
914     return inner(m,v)
915    
916     # m is B-coefficient, u is vector, v is vector: B_ijk*v_k
917     # m is A-coefficient, u is vector, v is grad scalar: A_ijl*v_l
918     if len(s_m) == 3 and len(s_v) == 1:
919     return matrix_mult(m,v)
920    
921     # m is A-coefficient, u is scalar, v is grad vector: A_jkl*v_kl
922     # m is C-coefficient, u is vector, v is grad vector: C_ikl*v_kl
923     if len(s_m) == 3 and len(s_v) == 2:
924     return tensor_mult(m,v)
925    
926     # m is A-coefficient, u is vector, v is grad vector: A_ijkl*v_kl
927     if len(s_m) == 4 and len(s_v) == 2:
928     return tensor_mult(m,v)
929    
930 caltinay 3849 def _updateRHS(self, expressions, subs):
931     """
932     """
933 caltinay 3831 ev=Evaluator()
934     names=[]
935     for name in expressions:
936     if name in self.__COEFFICIENTS:
937     ev.addExpression(expressions[name])
938     names.append(name)
939     if len(names)==0:
940     return
941 gross 3839 self.trace3("Starting expression evaluation.")
942 caltinay 3864 T0=time()
943 caltinay 3831 ev.subs(**subs)
944     res=ev.evaluate()
945     if len(names)==1: res=[res]
946 caltinay 3864 self.trace3("RHS expressions evaluated in %f seconds."%(time()-T0))
947 caltinay 3831 for i in range(len(names)):
948 gross 3839 self.trace3("util.Lsup(%s)=%s"%(names[i],util.Lsup(res[i])))
949 gross 3856 args=dict(zip(names,res))
950     # reset coefficients may be set at previous calls:
951 caltinay 3859 for n in self.__COEFFICIENTS:
952     if not args.has_key(n): args[n]=Data()
953     if not args.has_key('r'): args['r']=Data()
954 gross 3856 self._lpde.setValue(**args)
955 caltinay 3831
956 caltinay 3849 def _updateMatrix(self, expressions, subs):
957     """
958     """
959 caltinay 3831 ev=Evaluator()
960     names=[]
961     for name in expressions:
962     if not name in self.__COEFFICIENTS:
963     ev.addExpression(expressions[name])
964     names.append(name)
965     if len(names)==0:
966     return
967 gross 3839 self.trace3("Starting expression evaluation.")
968 caltinay 3864 T0=time()
969 caltinay 3831 ev.subs(**subs)
970     res=ev.evaluate()
971     if len(names)==1: res=[res]
972 caltinay 3864 self.trace3("Matrix expressions evaluated in %f seconds."%(time()-T0))
973 caltinay 3831 for i in range(len(names)):
974 gross 3839 self.trace3("util.Lsup(%s)=%s"%(names[i],util.Lsup(res[i])))
975 caltinay 3831 self._lpde.setValue(**dict(zip(names,res)))
976    
977 caltinay 3849 def _updateLinearPDE(self, expressions, subs):
978     self._updateMatrix(expressions, subs)
979     self._updateRHS(expressions, subs)
980 caltinay 3831
981 caltinay 3617
982 gross 3844 class VariationalProblem(object):
983     """
984     This class is used to define a general constraint vartional problem
985 caltinay 3849 for an unknown function *u* and (spatially variable) parameter *p* on a
986     given domain defined through a `Domain` object. *u* may be a scalar or a
987     vector. *p* which may be a scalar or a vector may not be present.
988    
989     The solution *u* and the paremeter *p* are given as the solution of the
990     minimization problem:
991    
992 gross 3844 *min_{u,p} int(H) + int(h)*
993    
994 caltinay 3849 where int{H} refers to integration over the domain and
995     *H*=f(*x*,*u*,*grad(u)*,*p*, *grad(p)*) is a function which may depend on
996     the location *x* within the domain and is a function of the solution *u*
997     and the parameter *p* and their gradients *grad(u)* and *grad(p)*,
998     respectively.
999     Similarly, int{H} refers to integration over the boundary of the domain and
1000     *h=f(*x*,*u*, *p*) is a function which may depend on the location *x*
1001     within the domain boundary and is a function of the solution *u* and the
1002     parameter *p*.
1003 gross 3844
1004 caltinay 3849 If *p* is present, *u* is the solution of a PDE with coefficients depending
1005     on the parameter *p*. The PDE defines a constraint for the variational
1006     problem. It is assumed that, if *p* is present, for any given parameter *p*
1007     a unique solution *u* of the PDE exists.
1008    
1009 gross 3844 For a single PDE having a solution with a single component the nonlinear
1010     PDE is defined in the following form:
1011    
1012     *-div(X) + Y = 0*
1013    
1014 caltinay 3849 where *X*,*Y*=f(*x*,*u*,*grad(u)*, *p*, *grad(p)*), *div(F)* denotes the
1015     divergence of *F* and *grad(F)* is the spatial derivative of *F*.
1016 gross 3844
1017     The coefficients *X* (rank 1) and *Y* (scalar) have to be specified through
1018     `Symbol` objects.
1019    
1020     The following natural boundary conditions are considered:
1021    
1022     *n[j]*X[j] - y = 0*
1023    
1024     where *n* is the outer normal field. Notice that the coefficient *X*
1025     is defined in the PDE. The coefficient *y* is a scalar `Symbol`.
1026    
1027     Constraints for the solution prescribe the value of the solution at certain
1028     locations in the domain. They have the form
1029    
1030     *u=r* where *q>0*
1031    
1032     *r* and *q* are each scalar where *q* is the characteristic function
1033     defining where the constraint is applied. The constraints override any
1034     other condition set by the PDE or the boundary condition.
1035    
1036     For a system of PDEs and a solution with several components, *u* is rank
1037     one, while the PDE coefficient *X* is rank two and *Y* and *y* is rank one.
1038    
1039     The PDE is solved by linearising the coefficients and iteratively solving
1040     the corresponding linear PDE until the error is smaller than a tolerance
1041     or a maximum number of iterations is reached.
1042    
1043     Typical usage:
1044    
1045     s=Symbol('s', dim=dom.getDim())
1046     T = Symbol('T', dim=dom.getDim())
1047     log_k = Symbol('log_k', dim=dom.getDim())
1048 caltinay 3849 v = VariationalProblem(dom, u=T, p=log_k)
1049 caltinay 3855 v.setValue(X=exp(-log_k)*grad(T), Y=s, h=0.3*(T-0.3)**2)
1050     T, log_k, l = v.getSolution(T=0., log_k=1, s=2.45)
1051     sT,S_log_k=v.getSensitivity(s, direction=1, T=T, log_k=log_k, s=2.45)
1052 caltinay 3849
1053 gross 3844 """
1054     DEBUG0=0 # no debug info
1055     DEBUG1=1 # info on Newton-Raphson solver printed
1056 caltinay 3849 DEBUG2=2 # in addition info from linear solver is printed
1057     DEBUG3=3 # in addition info on linearization is printed
1058     DEBUG4=4 # in addition info on LinearPDE handling is printed
1059     def __init__(self, domain, u, p=None, debug=DEBUG0):
1060 gross 3844 """
1061 caltinay 3849 Initializes a new variational problem.
1062 gross 3844
1063     :param domain: domain of the PDE
1064     :type domain: `Domain`
1065     :param u: The symbol for the unknown PDE function u.
1066     :type u: `Symbol`
1067 caltinay 3849 :param p: The symbol for the parameter p.
1068     :type p: `Symbol`
1069 gross 3844 :param debug: level of debug information to be printed
1070     """
1071     self.__COEFFICIENTS = [ "X", "X_reduced", "Y", "Y_reduced", "y", "y_reduced", "y_contact", "y_contact_reduced", "y_dirac", \
1072     "H", "h_reduced", "h", "h_reduced", "h_contact", "h_contact_reduced", "h_dirac" ]
1073    
1074     self._set_coeffs={}
1075     self._unknown=u
1076     self._parameter=p
1077     self._debug=debug
1078     self.dim = domain.getDim()
1079 caltinay 3849
1080 gross 3844 if self._parameter == None:
1081 caltinay 3849 U=u
1082     else:
1083 gross 3869 self._lagrangean=Symbol("lambda%s"%id(self), self._unknown.getShape(), dim=self.dim)
1084 caltinay 3855 U=concatenateRow(self._parameter, self._unknown, self._lagrangean)
1085 gross 3844 self.__PDE=NonlinearPDE(domain, u=U, debug=debug)
1086 caltinay 3849
1087 gross 3844 def __str__(self):
1088     """
1089 caltinay 3849 Returns the string representation of the problem.
1090 gross 3844
1091 caltinay 3849 :return: a simple representation of the variational problem
1092 gross 3844 :rtype: ``str``
1093     """
1094 caltinay 3849 return "<%s %d>"%(self.__class__.__name__, id(self))
1095    
1096 gross 3844 def trace1(self, text):
1097     """
1098 caltinay 3849 Prints the text message if the debug level is greater than DEBUG0
1099 gross 3844
1100     :param text: message to be printed
1101     :type text: ``string``
1102     """
1103     if self._debug > self.DEBUG0:
1104     print("%s: %s"%(str(self), text))
1105 caltinay 3849
1106 gross 3844 def trace3(self, text):
1107     """
1108 caltinay 3849 Prints the text message if the debug level is greater than DEBUG3
1109 gross 3844
1110     :param text: message to be printed
1111     :type text: ``string``
1112     """
1113     if self._debug > self.DEBUG2:
1114     print("%s: %s"%(str(self), text))
1115 caltinay 3849
1116     def getNonlinearPDE(self):
1117     """
1118     Returns the `NonlinearPDE` used to solve the variational problem
1119    
1120     :return: underlying nonlinear PDE
1121     :rtype: `NonlinearPDE`
1122     """
1123     return self.__PDE
1124    
1125 gross 3844 def getNumSolutions(self):
1126     """
1127 caltinay 3849 Returns the number of solution components
1128 gross 3844 :rtype: ``int``
1129     """
1130     s=self._unknown.getShape()
1131     if len(s) > 0:
1132     return s[0]
1133     else:
1134     return 1
1135 caltinay 3849
1136 gross 3844 def getNumParameters(self):
1137     """
1138 caltinay 3849 Returns the number of parameter components. If no parameter is present
1139 gross 3844 zero is returned.
1140 caltinay 3849
1141 gross 3844 :rtype: ``int``
1142     """
1143 caltinay 3855 if self._parameter == None:
1144 caltinay 3849 return 0
1145     else:
1146 caltinay 3855 s=self._parameter.getShape()
1147 caltinay 3849 if len(s) > 0:
1148     return s[0]
1149     else:
1150     return 1
1151    
1152     def getShapeOfCoefficient(self, name):
1153 gross 3844 """
1154 caltinay 3849 Returns the shape of the coefficient ``name``
1155 gross 3844
1156     :param name: name of the coefficient enquired
1157     :type name: ``string``
1158     :return: the shape of the coefficient ``name``
1159     :rtype: ``tuple`` of ``int``
1160     :raise IllegalCoefficient: if ``name`` is not a coefficient of the PDE
1161     """
1162     numSol=self.getNumSolutions()
1163     numParams=self.getNumParameters()
1164     dim = self.dim
1165     if (name=="X" or name=="X_reduced") and numParams>0:
1166 caltinay 3849 if numSol > 1:
1167     return (numSol,dim)
1168 gross 3844 else:
1169     return (dim,)
1170     elif name=="r" or name == "q" :
1171 caltinay 3849 if numSol > 1:
1172     return (numSol,)
1173 gross 3844 else:
1174     return ()
1175     elif ( name=="rp" or name == "qp") and numParams>0:
1176 caltinay 3849 if numParams > 1:
1177     return (numParams,)
1178 gross 3844 else:
1179     return ()
1180     elif ( name=="Y" or name=="Y_reduced") and numParams>0:
1181     if numSol > 1:
1182 caltinay 3849 return (numSol,)
1183 gross 3844 else:
1184     return ()
1185     elif ( name=="y" or name=="y_reduced") and numParams>0:
1186 caltinay 3849 if numSol > 1:
1187     return (numSol,)
1188 gross 3844 else:
1189 caltinay 3849 return ()
1190 gross 3844 elif ( name=="y_contact" or name=="y_contact_reduced") and numParams>0:
1191 caltinay 3849 if numSol > 1:
1192     return (numSol,)
1193 gross 3844 else:
1194 caltinay 3849 return ()
1195 gross 3844 elif name=="y_dirac" and numParams>0:
1196 caltinay 3849 if numSol > 1:
1197     return (numSol,)
1198 gross 3844 else:
1199 caltinay 3849 return ()
1200 gross 3844 elif name=="H" or name=="H_reduced":
1201     return ()
1202     elif name=="h" or name=="h_reduced":
1203 caltinay 3849 return ()
1204 gross 3844 elif name=="h_contact" or name=="h_contact_reduced":
1205 caltinay 3849 return ()
1206 gross 3844 elif name=="h_dirac":
1207 caltinay 3849 return ()
1208 gross 3844 else:
1209 caltinay 3849 raise IllegalCoefficient("Attempt to request unknown coefficient %s"%name)
1210 gross 3844
1211     def createCoefficient(self, name):
1212     """
1213 caltinay 3849 Creates a new coefficient ``name`` as Symbol
1214 gross 3844
1215     :param name: name of the coefficient requested
1216     :type name: ``string``
1217     :return: the value of the coefficient
1218 caltinay 3849 :rtype: `Symbol` or `Data`
1219 gross 3844 :raise IllegalCoefficient: if ``name`` is not a coefficient of the PDE
1220 caltinay 3849 """
1221 gross 3844 numSol=self.getNumSolutions()
1222 caltinay 3849 numParams=self.getNumParameters()
1223 gross 3844 if name == "q":
1224 caltinay 3849 if numParams > 0:
1225 gross 3844 return self.getNonlinearPDE().createCoefficient("q")[numParams:numParams+numSol]
1226     else:
1227 caltinay 3849 return self.getNonlinearPDE().createCoefficient("q")
1228 gross 3844 elif name == "qp":
1229 caltinay 3849 if numParams > 0:
1230 gross 3844 return self.getNonlinearPDE().createCoefficient("q")[:numParams]
1231     else:
1232 caltinay 3849 raise IllegalCoefficient("Attempt to request coefficient %s"%name)
1233 gross 3844 else:
1234     s=self.getShapeOfCoefficient(name)
1235     return Symbol(name, s, dim=self.dim)
1236    
1237     def getCoefficient(self, name):
1238     """
1239     Returns the value of the coefficient ``name`` as Symbol
1240    
1241     :param name: name of the coefficient requested
1242     :type name: ``string``
1243     :return: the value of the coefficient
1244     :rtype: `Symbol`
1245     :raise IllegalCoefficient: if ``name`` is not a coefficient of the PDE
1246     """
1247     if self._set_coeffs.has_key(name):
1248     return self._set_coeffs[name]
1249     else:
1250     raise IllegalCoefficient("Attempt to request undefined coefficient %s"%name)
1251 caltinay 3849
1252 gross 3845 def __getNonlinearPDECoefficient(self, extension, capson=False, order=0):
1253 caltinay 3849 """
1254     """
1255     if capson:
1256     H_key="H"+extension
1257     X_key="X"+extension
1258     Y_key="Y"+extension
1259     Z_key="Z"+extension
1260     else:
1261     H_key="h"+extension
1262     X_key="x"+extension
1263     Y_key="y"+extension
1264     Z_key="z"+extension
1265 gross 3844
1266 gross 3858 Z=0
1267 caltinay 3849 if self._set_coeffs.has_key(H_key): Z+=self._set_coeffs[H_key]
1268 caltinay 3855 if self._set_coeffs.has_key(X_key): Z+=util.inner(self._set_coeffs[X_key], util.grad(self._lagrangean))
1269     if self._set_coeffs.has_key(Y_key): Z+=util.inner(self._set_coeffs[Y_key], self._lagrangean)
1270 caltinay 3849
1271 caltinay 3855 if self.getNumParameters() > 0:
1272 caltinay 3849 if order == 0:
1273     Yp=getTotalDifferential(Z, self._parameter, order=0)
1274     Yu=getTotalDifferential(Z, self._unknown, order=0)
1275     Yl=getTotalDifferential(Z, self._lagrangean, order=0)
1276     Y=concatenateRow(Yp, Yl, Yu) # order different from solution!
1277     else:
1278     Yp,Xp=getTotalDifferential(Z, self._parameter, order=1)
1279     Yu,Xu=getTotalDifferential(Z, self._unknown, order=1)
1280     Yl,Xl=getTotalDifferential(Z, self._lagrangean, order=1)
1281     Y=concatenateRow(Yp, Yl, Yu) # order different from solution!
1282     X=concatenateRow(Xp, Xl, Xu) # order different from solution!
1283     else:
1284     if order == 0:
1285     Y=getTotalDifferential(Z, self._unknown, order=0)
1286     else:
1287     Y,X=getTotalDifferential(Z, self._unknown, order=1)
1288 gross 3858 if order == 0:
1289 caltinay 3859 return Y
1290     else:
1291     return Y,X
1292 caltinay 3849
1293 gross 3844 def setValue(self,**coefficients):
1294     """
1295     Sets new values to one or more coefficients.
1296    
1297     :keyword H: value for coefficient ``H``
1298     :type H: `Symbol`
1299     :keyword h: value for coefficient ``h``
1300     :type h: `Symbol`
1301     :keyword h_contact: value for coefficient `h_contact``
1302     :type h_contact: `Symbol`
1303     :keyword h_dirac: value for coefficient ``y_dirac``
1304     :type h_dirac: `Symbol`
1305     :keyword X: value for coefficient ``X``
1306     :type X: `Symbol` or any type that can be cast to a `Data` object
1307 gross 3845 :keyword Y: value for coefficient ``Y``=r
1308 gross 3844 :type Y: `Symbol` or any type that can be cast to a `Data` object
1309     :keyword y: value for coefficient ``y``
1310     :type y: `Symbol` or any type that can be cast to a `Data` object
1311     :keyword y_contact: value for coefficient ``y_contact``
1312     :type y_contact: `Symbol` or any type that can be cast to a `Data` object
1313     :keyword y_dirac: value for coefficient ``y_dirac``
1314     :type y_dirac: `Symbol` or any type that can be cast to a `Data` object
1315     :keyword q: mask for location of constraint
1316 caltinay 3849 :type q: any type that can be casted to a `Data` object
1317 gross 3844 :keyword r: value of solution prescribed by constraint
1318 caltinay 3849 :type r: `Symbol` or any type that can be cast to a `Data` object
1319 gross 3844 :keyword qp: mask for location of constraint fro parameter
1320 caltinay 3849 :type qp: any type that can be casted to a `Data` object
1321 gross 3844 :keyword rp: value of the parameter prescribed by parameter constraint
1322 caltinay 3849 :type rp: `Symbol` or any type that can be cast to a `Data` object
1323    
1324 gross 3844 :raise IllegalCoefficient: if an unknown coefficient keyword is used
1325     :raise IllegalCoefficientValue: if a supplied coefficient value has an
1326     invalid shape
1327     """
1328     numSol=self.getNumSolutions()
1329     numParams=self.getNumParameters()
1330     update=[]
1331     for name,val in coefficients.iteritems():
1332     shape=util.getShape(val)
1333     if not shape == self.getShapeOfCoefficient(name):
1334 caltinay 3862 raise IllegalCoefficientValue("%s has shape %s but must have shape %s"%(name, shape, self.getShapeOfCoefficient(name)))
1335 gross 3844 if name == "q":
1336 gross 3858 self._q = val
1337 caltinay 3849 update.append("q")
1338    
1339 gross 3844 elif name == "qp":
1340 caltinay 3849 if numParams<1:
1341     raise IllegalCoefficientValue("Illegal coefficient %s - no parameter present."%name)
1342 gross 3858 self._qp = val
1343 caltinay 3849 update.append("q")
1344    
1345     elif name == "r":
1346 gross 3858 self._r = val
1347 caltinay 3849 update.append("r")
1348    
1349 gross 3844 elif name == "rp":
1350 caltinay 3849 if numParams<1:
1351     raise IllegalCoefficientValue("Illegal coefficient %s - no parameter present."%name)
1352 gross 3858 self._rp = val
1353 caltinay 3849 update.append("r")
1354    
1355 gross 3844 elif name=="X":
1356 caltinay 3849 if numParams<1:
1357     raise IllegalCoefficientValue("Illegal coefficient %s - no parameter present."%name)
1358 gross 3844 self._set_coeffs['X']=val
1359     update.append("Y")
1360 caltinay 3849
1361 gross 3844 elif name=="X_reduced":
1362 caltinay 3849 if numParams<1:
1363     raise IllegalCoefficientValue("Illegal coefficient %s - no parameter present."%name)
1364 gross 3844 self._set_coeffs['X_reduced']=val
1365     update.append("Y_reduced")
1366 caltinay 3849
1367 gross 3844 elif name=="Y":
1368 caltinay 3849 if numParams<1:
1369     raise IllegalCoefficientValue("Illegal coefficient %s - no parameter present."%name)
1370 gross 3844 self._set_coeffs['Y']=val
1371     update.append("Y")
1372 caltinay 3849
1373 gross 3844 elif name=="Y_reduced":
1374 caltinay 3849 if numParams<1:
1375     raise IllegalCoefficientValue("Illegal coefficient %s - no parameter present."%name)
1376 gross 3844 self._set_coeffs['Y_reduced']=val
1377     update.append("Y_reduced")
1378 caltinay 3849
1379 gross 3844 elif name=="H":
1380     self._set_coeffs['H']=val
1381     update.append("Y")
1382 caltinay 3849
1383 gross 3845 elif name=="H_reduced":
1384 gross 3844 self._set_coeffs['H_reduced']=val
1385     update.append("Y_reduced")
1386 caltinay 3849
1387 gross 3844 elif name=="y":
1388 caltinay 3849 if numParams<1:
1389     raise IllegalCoefficientValue("Illegal coefficient %s - no parameter present."%name)
1390 gross 3844 self._set_coeffs['y']=val
1391     update.append("y")
1392 caltinay 3849
1393 gross 3844 elif name=="y_reduced":
1394 caltinay 3849 if numParams<1:
1395     raise IllegalCoefficientValue("Illegal coefficient %s - no parameter present."%name)
1396 gross 3844 self._set_coeffs['y_reduced']=val
1397     update.append("y_reduced")
1398 caltinay 3849
1399 gross 3844 elif name=="h":
1400     self._set_coeffs['h']=val
1401     update.append("y")
1402 caltinay 3849
1403 gross 3845 elif name=="h_reduced":
1404 gross 3844 self._set_coeffs['h_reduced']=val
1405     update.append("y_reduced")
1406 caltinay 3849
1407 gross 3844 elif name=="y_contact":
1408 caltinay 3849 if numParams<1:
1409     raise IllegalCoefficientValue("Illegal coefficient %s - no parameter present."%name)
1410 gross 3844 self._set_coeffs['y_contact']=val
1411     update.append("y_contact")
1412 caltinay 3849
1413 gross 3844 elif name=="y_contact_reduced":
1414 caltinay 3849 if numParams<1:
1415     raise IllegalCoefficientValue("Illegal coefficient %s - no parameter present."%name)
1416 gross 3844 self._set_coeffs['y_contact_reduced']=val
1417     update.append("y_contact_reduced")
1418 caltinay 3849
1419 gross 3844 elif name=="h_contact":
1420     self._set_coeffs['h_contact']=val
1421     update.append("y_contact")
1422 caltinay 3849
1423 gross 3845 elif name=="h_contact_reduced":
1424 gross 3844 self._set_coeffs['h_contact_reduced']=val
1425     update.append("y_contact_reduced")
1426 caltinay 3849
1427 gross 3844 elif name=="y_dirac":
1428 caltinay 3849 if numParams<1:
1429     raise IllegalCoefficientValue("Illegal coefficient %s - no parameter present."%name)
1430 gross 3844 self._set_coeffs['y_dirac']=val
1431     update.append("y_dirac")
1432 caltinay 3849
1433 gross 3844 elif name=="h_dirac":
1434     self._set_coeffs['h_diract']=val
1435     update.append("y_dirac")
1436 caltinay 3849 else:
1437 gross 3844 raise IllegalCoefficient("Attempt to set unknown coefficient %s"%name)
1438    
1439 gross 3858 # now we can update the coefficients of the nonlinear PDE:
1440     coeff2={}
1441     if "q" in update:
1442 caltinay 3859 if numParams>0:
1443     q=self.getNonlinearPDE().createCoefficient("q")
1444     if hasattr(self, "_qp"): q[:numParams]=self._qp
1445     if hasattr(self, "_q"):
1446     q[numParams:numParams+numSol]=self._q
1447     q[numParams+numSol:]=self._q
1448     else:
1449     q=self._q
1450     coeff2["q"]=q
1451 gross 3858 if "r" in update:
1452 caltinay 3859 if numParams>0:
1453     if hasattr(self, "_rp"):
1454     rp=self._rp
1455 caltinay 3849 else:
1456 caltinay 3859 rp=numpy.zeros(self.getShapeOfCoefficient('rp'))
1457     if hasattr(self, "_r"):
1458     r=self._r
1459     else:
1460     r=numpy.zeros(self.getShapeOfCoefficient('r'))
1461     coeff2["r"]=concatenateRow(rp, r, numpy.zeros(self.getShapeOfCoefficient('r')) )
1462     else:
1463     coeff2["r"]=self._r
1464 caltinay 3849
1465 gross 3858 if "Y" in update:
1466 caltinay 3859 Y,X = self.__getNonlinearPDECoefficient("", capson=True, order=1)
1467     coeff2["Y"]=Y
1468     coeff2["X"]=X
1469 gross 3858 if "y" in update:
1470 caltinay 3859 coeff2["y"] = self.__getNonlinearPDECoefficient("", capson=False)
1471 gross 3858 if "y_contact" in update:
1472 caltinay 3859 coeff2["y_contact"] = self.__getNonlinearPDECoefficient("_contact", capson=False)
1473 gross 3858 if "y_dirac" in update:
1474 caltinay 3859 coeff2["y_dirac"] = self.__getNonlinearPDECoefficient("_dirac", capson=False)
1475 gross 3858 if "Y_reduced" in update:
1476 caltinay 3859 Y,X = self.__getNonlinearPDECoefficient("_reduced",capson=True, order=1)
1477     coeff2["Y_reduced"]=Y
1478     coeff2["X_reduced"]=X
1479 gross 3858 if "y_reduced" in update:
1480 caltinay 3859 coeff2["y_reduced"]= self.__getNonlinearPDECoefficient("_reduced",capson=False)
1481 gross 3858 if "y_contact_reduced" in update:
1482 caltinay 3859 coeff2["y_contact_reduced"] = self.__getNonlinearPDECoefficient("_contact_reduced",capson=False)
1483 caltinay 3849
1484 gross 3858 # and now we can set the PDE coefficients:
1485     self.getNonlinearPDE().setValue(**coeff2)
1486 gross 3845
1487     def getSolution(self, **subs):
1488     """
1489 caltinay 3849 Returns the solution of the variational problem.
1490 gross 3845 :param subs: Substitutions for all symbols used in the coefficients
1491 caltinay 3849 including the initial value for solution *u* and for the
1492     parameter *p* (if present)
1493 gross 3845 :return: parameter, corresponding solution and lagrangean multiplier
1494     :rtype: tuple of `Data` or single `Data` (if no parameter present)
1495     """
1496     numSol=self.getNumSolutions()
1497     numParams=self.getNumParameters()
1498 caltinay 3859
1499 caltinay 3849 if numParams>0:
1500 caltinay 3859 fs=self.getNonlinearPDE().getLinearPDE().getFunctionSpaceForSolution()
1501 caltinay 3855 subs[self._lagrangean.name]=Data(0., self._lagrangean.getShape(), fs)
1502 caltinay 3849 p_sym=self._parameter.atoms().pop().name
1503     if not subs.has_key(p_sym):
1504     raise KeyError("Initial value for '%s' missing."%p_sym)
1505 caltinay 3855 pi=subs[p_sym]
1506 caltinay 3849 if not isinstance(pi,Data):
1507 caltinay 3855 fs=self.getNonlinearPDE().getLinearPDE().getFunctionSpaceForSolution()
1508     pi=Data(pi, fs)
1509     subs[p_sym]=pi
1510 caltinay 3849
1511     Ui=self.getNonlinearPDE().getSolution(**subs)
1512    
1513 caltinay 3859 if numParams > 0:
1514 caltinay 3849 # return parameter, solution, lagrangean multiplier
1515 gross 3869 if numParams == 1:
1516     p=Ui[0]
1517     else:
1518     p=Ui[:numParams]
1519     if numSol == 1:
1520     u=Ui[numParams]
1521     l=Ui[numParams+1]
1522     else:
1523     u=Ui[numParams:numParams+numSol]
1524     l= Ui[numParams+numSol:]
1525     return p,u,l
1526 gross 3845 else:
1527 caltinay 3849 return Ui
1528    

  ViewVC Help
Powered by ViewVC 1.1.26