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

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

Parent Directory Parent Directory | Revision Log Revision Log


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

  ViewVC Help
Powered by ViewVC 1.1.26