/[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 3975 - (show annotations)
Thu Sep 20 01:54:06 2012 UTC (7 years ago) by caltinay
File MIME type: text/x-python
File size: 60489 byte(s)
Merged symbolic branch into trunk. Curious what daniel and spartacus have to
say...

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

  ViewVC Help
Powered by ViewVC 1.1.26