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

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

Parent Directory Parent Directory | Revision Log Revision Log


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

  ViewVC Help
Powered by ViewVC 1.1.26