/[escript]/trunk/escript/py_src/nonlinearPDE.py
ViewVC logotype

Contents of /trunk/escript/py_src/nonlinearPDE.py

Parent Directory Parent Directory | Revision Log Revision Log


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

  ViewVC Help
Powered by ViewVC 1.1.26