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

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

Parent Directory Parent Directory | Revision Log Revision Log


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

  ViewVC Help
Powered by ViewVC 1.1.26