# -*- coding: utf-8 -*- ############################################################################## # # Copyright (c) 2003-2020 by The University of Queensland # http://www.uq.edu.au # # Primary Business: Queensland, Australia # Licensed under the Apache License, version 2.0 # http://www.apache.org/licenses/LICENSE-2.0 # # Development until 2012 by Earth Systems Science Computational Center (ESSCC) # Development 2012-2013 by School of Earth Sciences # Development from 2014 by Centre for Geoscience Computing (GeoComp) # Development from 2019 by School of Earth and Environmental Sciences # ############################################################################## from __future__ import print_function, division __copyright__="""Copyright (c) 2003-2020 by The University of Queensland http://www.uq.edu.au Primary Business: Queensland, Australia""" __license__="""Licensed under the Apache License, version 2.0 http://www.apache.org/licenses/LICENSE-2.0""" __url__="https://launchpad.net/escript-finley" """ :var __author__: name of author :var __copyright__: copyrights :var __license__: licence agreement :var __url__: url entry point on documentation :var __version__: version :var __date__: date of the version """ from .start import HAVE_SYMBOLS import numpy from time import time from . import linearPDEs as lpe from . import util from .escriptcpp import Data if HAVE_SYMBOLS: import sympy import esys.escriptcore.symbolic as symb __author__="Cihan Altinay, Lutz Gross" class iteration_steps_maxReached(Exception): """ Exception thrown if the maximum number of iteration steps is reached. """ pass class DivergenceDetected(Exception): """ Exception thrown if Newton-Raphson did not converge. """ pass class InadmissiblePDEOrdering(Exception): """ Exception thrown if the ordering of the PDE equations should be revised. """ pass def concatenateRow(*args): """ """ if len(args)==1: return args[0] if len(args)>2: return concatenateRow(args[0], concatenateRow(*args[1:])) lshape=util.getShape(args[0]) rshape=util.getShape(args[1]) if len(lshape)==0: if len(rshape)==0: shape=(2,) res=numpy.zeros(shape, dtype=object) res[0]=args[0] res[1]=args[1] elif len(rshape)==1: shape=(rshape[0]+1,) res=numpy.zeros(shape, dtype=object) res[0]=args[0] res[1:]=args[1] elif len(lshape)==1: if len(rshape)==1: shape=(2,)+lshape res=numpy.zeros(shape, dtype=object) res[0]=args[0] res[1]=args[1] else: shape=(rshape[0]+1,)+lshape res=numpy.zeros(shape, dtype=object) res[0]=args[0] res[1:]=args[1] else: if len(rshape)==1: shape=(lshape[0]+1,)+lshape[1:] res=numpy.zeros(shape, dtype=object) res[:lshape[0]]=args[0] res[lshape[0]]=args[1] else: shape=(lshape[0]+rshape[0],)+lshape[1:] res=numpy.zeros(shape, dtype=object) res[:lshape[0]]=args[0] res[lshape[0]:]=args[1] subs=args[0].getDataSubstitutions().copy() subs.update(args[1].getDataSubstitutions()) dim=args[1].getDim() if args[0].getDim()<0 else args[0].getDim() return symb.Symbol(res, dim=dim, subs=subs) class NonlinearPDE(object): """ This class is used to define a general nonlinear, steady, second order PDE for an unknown function *u* on a given domain defined through a `Domain` object. For a single PDE having a solution with a single component the nonlinear PDE is defined in the following form: *-div(X) + Y = 0* where *X*,*Y*=f(*u*,*grad(u)*). *div(F)* denotes the divergence of *F* and *grad(F)* is the spatial derivative of *F*. The coefficients *X* (rank 1) and *Y* (scalar) have to be specified through `Symbol` objects. The following natural boundary conditions are considered: *n[j]*X[j] + y = 0* where *n* is the outer normal field. Notice that the coefficient *X* is defined in the PDE. The coefficient *y* is a scalar `Symbol`. Constraints for the solution prescribe the value of the solution at certain locations in the domain. They have the form *u=r* where *q>0* *r* and *q* are each scalar where *q* is the characteristic function defining where the constraint is applied. The constraints override any other condition set by the PDE or the boundary condition. For a system of PDEs and a solution with several components, *u* is rank one, while the PDE coefficient *X* is rank two and *y* is rank one. The PDE is solved by linearising the coefficients and iteratively solving the corresponding linear PDE until the error is smaller than a tolerance or a maximum number of iterations is reached. Typical usage:: u = Symbol('u', dim=dom.getDim()) p = NonlinearPDE(dom, u) p.setValue(X=grad(u), Y=1+5*u) v = p.getSolution(u=0.) """ DEBUG0=0 # no debug info DEBUG1=1 # info on Newton-Raphson solver printed DEBUG2=2 # in addition info from linear solver is printed DEBUG3=3 # in addition info on linearization is printed DEBUG4=4 # in addition info on LinearPDE handling is printed ORDER=0 def __init__(self, domain, u, debug=DEBUG0): """ Initializes a new nonlinear PDE. :param domain: domain of the PDE :type domain: `Domain` :param u: The symbol for the unknown PDE function u. :type u: `Symbol` :param debug: level of debug information to be printed """ if not HAVE_SYMBOLS: raise RuntimeError("Trying to instantiate a NonlinearPDE but sympy not available") self.__COEFFICIENTS = [ "X", "X_reduced", "Y", "Y_reduced", "y", "y_reduced", "y_contact", "y_contact_reduced", "y_dirac"] self._r=Data() self._set_coeffs={} self._unknown=u self._debug=debug self._rtol=1e-6 self._atol=0. self._iteration_steps_max=100 self._omega_min=0.0001 self._quadratic_convergence_limit=0.2 self._simplified_newton_limit=0.1 self.dim = domain.getDim() if u.getRank()==0: numEquations=1 else: numEquations=u.getShape()[0] numSolutions=numEquations self._lpde=lpe.LinearPDE(domain,numEquations,numSolutions,self._debug > self.DEBUG4 ) def __str__(self): """ Returns the string representation of the PDE :return: a simple representation of the PDE :rtype: ``str`` """ return "<%s %d>"%(self.__class__.__name__, id(self)) def getUnknownSymbol(self): """ Returns the symbol of the PDE unknown :return: the symbol of the PDE unknown :rtype: `Symbol` """ return self._unknown def trace1(self, text): """ Prints the text message if the debug level is greater than DEBUG0 :param text: message to be printed :type text: ``string`` """ if self._debug > self.DEBUG0: print("%s: %s"%(str(self), text)) def trace3(self, text): """ Prints the text message if the debug level is greater than DEBUG3 :param text: message to be printed :type text: ``string`` """ if self._debug > self.DEBUG2: print("%s: %s"%(str(self), text)) def getLinearSolverOptions(self): """ Returns the options of the linear PDE solver class """ return self._lpde.getSolverOptions() def getLinearPDE(self): """ Returns the linear PDE used to calculate the Newton-Raphson update :rtype: `LinearPDE` """ return self._lpde def setOptions(self, **opts): """ Allows setting options for the nonlinear PDE. The supported options are: ``tolerance`` error tolerance for the Newton method ``iteration_steps_max`` maximum number of Newton iterations ``omega_min`` minimum relaxation factor ``atol`` solution norms less than ``atol`` are assumed to be ``atol``. This can be useful if one of your solutions is expected to be zero. ``quadratic_convergence_limit`` if the norm of the Newton-Raphson correction is reduced by less than ``quadratic_convergence_limit`` between two iteration steps quadratic convergence is assumed. ``simplified_newton_limit`` if the norm of the defect is reduced by less than ``simplified_newton_limit`` between two iteration steps and quadratic convergence is detected the iteration switches to the simplified Newton-Raphson scheme. """ for key in opts: if key=='tolerance': self._rtol=opts[key] elif key=='iteration_steps_max': self._iteration_steps_max=opts[key] elif key=='omega_min': self._omega_min=opts[key] elif key=='atol': self._atol=opts[key] elif key=='quadratic_convergence_limit': self._quadratic_convergence_limit=opts[key] elif key=='simplified_newton_limit': self._simplified_newton_limit=opts[key] else: raise KeyError("Invalid option '%s'"%key) def getSolution(self, **subs): """ Returns the solution of the PDE. :param subs: Substitutions for all symbols used in the coefficients including the initial value for the unknown *u*. :return: the solution :rtype: `Data` """ # get the initial value for the iteration process # collect components of unknown in u_syms u_syms=[] simple_u=False for i in numpy.ndindex(self._unknown.getShape()): u_syms.append(symb.Symbol(self._unknown[i]).atoms(sympy.Symbol).pop().name) if len(set(u_syms))==1: simple_u=True e=symb.Evaluator(self._unknown) for sym in u_syms: if not sym in subs: raise KeyError("Initial value for '%s' missing."%sym) if not isinstance(subs[sym], Data): subs[sym]=Data(subs[sym], self._lpde.getFunctionSpaceForSolution()) e.subs(**{sym:subs[sym]}) ui=e.evaluate() # modify ui so it meets the constraints: q=self._lpde.getCoefficient("q") if not q.isEmpty(): if hasattr(self, "_r"): r=self._r if symb.isSymbol(r): r=symb.Evaluator(r).evaluate(**subs) elif not isinstance(r, Data): r=Data(r, self._lpde.getFunctionSpaceForSolution()) elif r.isEmpty(): r=0 else: r=0 ui = q * r + (1-q) * ui # separate symbolic expressions from other coefficients constants={} expressions={} for n, e in sorted(self._set_coeffs.items(), key=lambda x: x[0]): if symb.isSymbol(e): expressions[n]=e else: constants[n]=e # set constant PDE values now self._lpde.setValue(**constants) self._lpde.getSolverOptions().setAbsoluteTolerance(0.) self._lpde.getSolverOptions().setVerbosity(self._debug > self.DEBUG1) #===================================================================== # perform Newton iterations until error is small enough or # maximum number of iterations reached n=0 omega=1. # relaxation factor use_simplified_Newton=False defect_norm=None delta_norm=None converged=False #subs[u_sym]=ui if simple_u: subs[u_syms[0]]=ui else: for i in range(len(u_syms)): subs[u_syms[i]]=ui[i] while not converged: if n > self._iteration_steps_max: raise iteration_steps_maxReached("maximum number of iteration steps reached, giving up.") self.trace1(5*"="+" iteration step %d "%n + 5*"=") # calculate the correction delta_u if n == 0: self._updateLinearPDE(expressions, subs, **constants) defect_norm=self._getDefectNorm(self._lpde.getRightHandSide()) LINTOL=0.1 else: if not use_simplified_Newton: self._updateMatrix(expressions, subs) if q_u is None: LINTOL = 0.1 * min(qtol/defect_norm) else: LINTOL = 0.1 * max( q_u**2, min(qtol/defect_norm) ) LINTOL=max(1e-4,min(LINTOL, 0.1)) #LINTOL=1.e-5 self._lpde.getSolverOptions().setTolerance(LINTOL) self.trace1("PDE is solved with rel. tolerance = %e"%LINTOL) delta_u=self._lpde.getSolution() #check for reduced defect: omega=min(2*omega, 1.) # raise omega defect_reduced=False ui_old=ui while not defect_reduced: ui=ui_old-delta_u * omega if simple_u: subs[u_syms[0]]=ui else: for i in range(len(u_syms)): subs[u_syms[i]]=ui[i] self._updateRHS(expressions, subs, **constants) new_defect_norm=self._getDefectNorm(self._lpde.getRightHandSide()) defect_reduced=False for i in range(len( new_defect_norm)): if new_defect_norm[i] < defect_norm[i]: defect_reduced=True #print new_defect_norm #q_defect=max(self._getSafeRatio(new_defect_norm, defect_norm)) # if defect_norm==0 and new_defect_norm!=0 # this will be util.DBLE_MAX #self.trace1("Defect reduction = %e with relaxation factor %e."%(q_defect, omega)) if not defect_reduced: omega*=0.5 if omega < self._omega_min: raise DivergenceDetected("Underrelaxtion failed to reduce defect, giving up.") self.trace1("Defect reduction with relaxation factor %e."%(omega, )) delta_norm, delta_norm_old = self._getSolutionNorm(delta_u) * omega, delta_norm defect_norm, defect_norm_old = new_defect_norm, defect_norm u_norm=self._getSolutionNorm(ui, atol=self._atol) # set the tolerance on equation level: qtol=self._getSafeRatio(defect_norm_old * u_norm * self._rtol, delta_norm) # if defect_norm_old==0 and defect_norm_old!=0 this will be util.DBLE_MAX # -> the ordering of the equations is not appropriate. # if defect_norm_old==0 and defect_norm_old==0 this is zero so # convergence can happen for defect_norm==0 only. if not max(qtol) %e, omega=%e)"%(q_u, self._quadratic_convergence_limit,omega )) quadratic_convergence=False converged=False else: q_u=None converged=False quadratic_convergence=False if self._debug > self.DEBUG0: for i in range(len(u_norm)): self.trace1("Component %s: u: %e, du: %e, defect: %e, qtol: %e"%(i, u_norm[i], delta_norm[i], defect_norm[i], qtol[i])) if converged: self.trace1("Iteration has converged.") # Can we switch to simplified Newton? if quadratic_convergence: q_defect=max(self._getSafeRatio(defect_norm, defect_norm_old)) if q_defect < self._simplified_newton_limit: use_simplified_Newton=True self.trace1("Simplified Newton-Raphson is applied (rate %e < %e)."%(q_defect, self._simplified_newton_limit)) n+=1 self.trace1(5*"="+" Newton-Raphson iteration completed after %d steps "%n + 5*"=") return ui def _getDefectNorm(self, f): """ calculates the norm of the defect ``f`` :param f: defect vector :type f: `Data` of rank 0 or 1. :return: component-by-component norm of ``f`` :rtype: ``numpy.array`` of rank 1 :raise ValueError: if shape if ``f`` is incorrect. """ # this still needs some work!!! out=[] s=f.getShape() if len(s) == 0: out.append(util.Lsup(f)) elif len(s) == 1: for i in range(s[0]): out.append(util.Lsup(f[i])) else: raise ValueError("Illegal shape of defect vector: %s"%s) return numpy.array(out) def _getSolutionNorm(self, u, atol=0.): """ calculates the norm of the solution ``u`` :param u: solution vector :type u: `Data` of rank 0 or 1. :return: component-by-component norm of ``u`` :rtype: ``numpy.array`` of rank 1 :raise ValueError: if shape of ``u`` is incorrect. """ out=[] s=u.getShape() if len(s) == 0: out.append(max(util.Lsup(u),atol) ) elif len(s) == 1: for i in range(s[0]): out.append(max(util.Lsup(u[i]), atol)) else: raise ValueError("Illegal shape of solution: %s"%s) return numpy.array(out) def _getSafeRatio(self, a , b): """ Returns the component-by-component ratio of ''a'' and ''b'' If for a component ``i`` the values ``a[i]`` and ``b[i]`` are both equal to zero their ratio is set to zero. If ``b[i]`` equals zero but ``a[i]`` is positive the ratio is set to `util.DBLE_MAX`. :param a: numerator :param b: denominator :type a: ``numpy.array`` of rank 1 with non-negative entries. :type b: ``numpy.array`` of rank 1 with non-negative entries. :return: ratio of ``a`` and ``b`` :rtype: ``numpy.array`` """ out=0. if a.shape !=b.shape: raise ValueError("shapes must match.") s=a.shape if len(s) != 1: raise ValueError("rank one is expected.") out=numpy.zeros(s) for i in range(s[0]): if abs(b[i]) > 0: out[i]=a[i]/b[i] elif abs(a[i]) > 0: out[i] = util.DBLE_MAX else: out[i] = 0 return out def getNumSolutions(self): """ Returns the number of the solution components :rtype: ``int`` """ s=self._unknown.getShape() if len(s) > 0: return s[0] else: return 1 def getShapeOfCoefficient(self,name): """ Returns the shape of the coefficient ``name`` :param name: name of the coefficient enquired :type name: ``string`` :return: the shape of the coefficient ``name`` :rtype: ``tuple`` of ``int`` :raise IllegalCoefficient: if ``name`` is not a coefficient of the PDE """ numSol=self.getNumSolutions() dim = self.dim if name=="X" or name=="X_reduced": if numSol > 1: return (numSol,dim) else: return (dim,) elif name=="r" or name == "q" : if numSol > 1: return (numSol,) else: return () elif name=="Y" or name=="Y_reduced": if numSol > 1: return (numSol,) else: return () elif name=="y" or name=="y_reduced": if numSol > 1: return (numSol,) else: return () elif name=="y_contact" or name=="y_contact_reduced": if numSol > 1: return (numSol,) else: return () elif name=="y_dirac": if numSol > 1: return (numSol,) else: return () else: raise lpe.IllegalCoefficient("Attempt to request unknown coefficient %s"%name) def createCoefficient(self, name): """ Creates a new coefficient ``name`` as Symbol :param name: name of the coefficient requested :type name: ``string`` :return: the value of the coefficient :rtype: `Symbol` or `Data` (for name = "q") :raise IllegalCoefficient: if ``name`` is not a coefficient of the PDE """ if name == "q": return self._lpde.createCoefficient("q") else: s=self.getShapeOfCoefficient(name) return symb.Symbol(name, s, dim=self.dim) def getCoefficient(self, name): """ Returns the value of the coefficient ``name`` as Symbol :param name: name of the coefficient requested :type name: ``string`` :return: the value of the coefficient :rtype: `Symbol` :raise IllegalCoefficient: if ``name`` is not a coefficient of the PDE """ if name in self._set_coeffs: return self._set_coeffs[name] elif name == "r": if hasattr(self, "_r"): return self._r else: raise lpe.IllegalCoefficient("Attempt to request undefined coefficient %s"%name) elif name == "q": return self._lpde.getCoefficient("q") else: raise lpe.IllegalCoefficient("Attempt to request undefined coefficient %s"%name) def setValue(self,**coefficients): """ Sets new values to one or more coefficients. :param coefficients: new values assigned to coefficients :param coefficients: new values assigned to coefficients :keyword X: value for coefficient ``X`` :type X: `Symbol` or any type that can be cast to a `Data` object :keyword Y: value for coefficient ``Y`` :type Y: `Symbol` or any type that can be cast to a `Data` object :keyword y: value for coefficient ``y`` :type y: `Symbol` or any type that can be cast to a `Data` object :keyword y_contact: value for coefficient ``y_contact`` :type y_contact: `Symbol` or any type that can be cast to a `Data` object :keyword y_dirac: value for coefficient ``y_dirac`` :type y_dirac: `Symbol` or any type that can be cast to a `Data` object :keyword q: mask for location of constraint :type q: any type that can be cast to a `Data` object :keyword r: value of solution prescribed by constraint :type r: `Symbol` or any type that can be cast to a `Data` object :raise IllegalCoefficient: if an unknown coefficient keyword is used :raise IllegalCoefficientValue: if a supplied coefficient value has an invalid shape """ u=self._unknown for name,val in sorted(coefficients.items(), key=lambda x: x[0]): shape=util.getShape(val) if not shape == self.getShapeOfCoefficient(name): raise lpe.IllegalCoefficientValue("%s has shape %s but must have shape %s"%(name, shape, self.getShapeOfCoefficient(name))) rank=len(shape) if name == "q": self._lpde.setValue(q=val) elif name == "r": self._r=val elif name=="X" or name=="X_reduced": if rank != u.getRank()+1: raise lpe.IllegalCoefficientValue("%s must have rank %d"%(name,u.getRank()+1)) T0=time() B,A=symb.getTotalDifferential(val, u, 1) if name=='X_reduced': self.trace3("Computing A_reduced, B_reduced took %f seconds."%(time()-T0)) self._set_coeffs['A_reduced']=A self._set_coeffs['B_reduced']=B self._set_coeffs['X_reduced']=val else: self.trace3("Computing A, B took %f seconds."%(time()-T0)) self._set_coeffs['A']=A self._set_coeffs['B']=B self._set_coeffs['X']=val elif name=="Y" or name=="Y_reduced": if rank != u.getRank(): raise lpe.IllegalCoefficientValue("%s must have rank %d"%(name,u.getRank())) T0=time() D,C=symb.getTotalDifferential(val, u, 1) if name=='Y_reduced': self.trace3("Computing C_reduced, D_reduced took %f seconds."%(time()-T0)) self._set_coeffs['C_reduced']=C self._set_coeffs['D_reduced']=D self._set_coeffs['Y_reduced']=val else: self.trace3("Computing C, D took %f seconds."%(time()-T0)) self._set_coeffs['C']=C self._set_coeffs['D']=D self._set_coeffs['Y']=val elif name in ("y", "y_reduced", "y_contact", "y_contact_reduced", \ "y_dirac"): y=val if rank != u.getRank(): raise lpe.IllegalCoefficientValue("%s must have rank %d"%(name,u.getRank())) if not hasattr(y, 'diff'): d=numpy.zeros(u.getShape()+u.getShape()) else: d=y.diff(u) self._set_coeffs[name]=y self._set_coeffs['d'+name[1:]]=d else: raise lpe.IllegalCoefficient("Attempt to set unknown coefficient %s"%name) def getSensitivity(self, f, g=None, **subs): """ Calculates the sensitivity of the solution of an input factor ``f`` in direction ``g``. :param f: the input factor to be investigated. ``f`` may be of rank 0 or 1. :type f: `Symbol` :param g: the direction(s) of change. If not present, it is *g=eye(n)* where ``n`` is the number of components of ``f``. :type g: ``list`` or single of ``float``, ``numpy.array`` or `Data`. :param subs: Substitutions for all symbols used in the coefficients including unknown *u* and the input factor ``f`` to be investigated :return: the sensitivity :rtype: `Data` with shape *u.getShape()+(len(g),)* if *len(g)>1* or *u.getShape()* if *len(g)==1* """ s_f=f.getShape() if len(s_f) == 0: len_f=1 elif len(s_f) == 1: len_f=s_f[0] else: raise ValueError("rank of input factor must be zero or one.") if not g is None: if len(s_f) == 0: if not isinstance(g, list): g=[g] else: if isinstance(g, list): if len(g) == 0: raise ValueError("no direction given.") 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 else: g=[g] # at this point g is a list of directions: len_g=len(g) for g_i in g: if not getShape(g_i) == s_f: raise ValueError("shape of direction (=%s) must match rank of input factor (=%s)"%(getShape(g_i) , s_f) ) else: len_g=len_f #*** at this point g is a list of direction or None and len_g is the # number of directions to be investigated. # now we make sure that the operator in the lpde is set (it is the same # as for the Newton-Raphson scheme) # if the solution etc are cached this could be omitted: constants={} expressions={} for n, e in sorted(self._set_coeffs.items(), key=lambda x: x[0]): if n not in self.__COEFFICIENTS: if symb.isSymbol(e): expressions[n]=e else: constants[n]=e self._lpde.setValue(**constants) self._updateMatrix(self, expressions, subs) #===================================================================== self._lpde.getSolverOptions().setAbsoluteTolerance(0.) self._lpde.getSolverOptions().setTolerance(self._rtol) self._lpde.getSolverOptions().setVerbosity(self._debug > self.DEBUG1) #===================================================================== # # evaluate the derivatives of X, etc with respect to f: # ev=symb.Evaluator() names=[] if hasattr(self, "_r"): if symb.isSymbol(self._r): names.append('r') ev.addExpression(self._r.diff(f)) for n in sorted(self._set_coeffs.keys()): if n in self.__COEFFICIENTS and symb.isSymbol(self._set_coeffs[n]): if n=="X" or n=="X_reduced": T0=time() B,A=symb.getTotalDifferential(self._set_coeffs[n], f, 1) if n=='X_reduced': self.trace3("Computing A_reduced, B_reduced took %f seconds."%(time()-T0)) names.append('A_reduced'); ev.addExpression(A) names.append('B_reduced'); ev.addExpression(B) else: self.trace3("Computing A, B took %f seconds."%(time()-T0)) names.append('A'); ev.addExpression(A) names.append('B'); ev.addExpression(B) elif n=="Y" or n=="Y_reduced": T0=time() D,C=symb.getTotalDifferential(self._set_coeffs[n], f, 1) if n=='Y_reduced': self.trace3("Computing C_reduced, D_reduced took %f seconds."%(time()-T0)) names.append('C_reduced'); ev.addExpression(C) names.append('D_reduced'); ev.addExpression(D) else: self.trace3("Computing C, D took %f seconds."%(time()-T0)) names.append('C'); ev.addExpression(C) names.append('D'); ev.addExpression(D) elif n in ("y", "y_reduced", "y_contact", "y_contact_reduced", "y_dirac"): names.append('d'+name[1:]); ev.addExpression(self._set_coeffs[name].diff(f)) relevant_symbols['d'+name[1:]]=self._set_coeffs[name].diff(f) res=ev.evaluate() if len(names)==1: res=[res] self.trace3("RHS expressions evaluated in %f seconds."%(time()-T0)) if self._debug > self.DEBUG2: for i in range(len(names)): self.trace3("util.Lsup(%s)=%s"%(names[i],util.Lsup(res[i]))) coeffs_f=dict(zip(names,res)) # # now we are ready to calculate the right hand side coefficients into # args by multiplication with g and grad(g). if len_g >1: if self.getNumSolutions() == 1: u_g=Data(0., (len_g,), self._lpde.getFunctionSpaceForSolution()) else: u_g=Data(0., (self.getNumSolutions(), len_g), self._lpde.getFunctionSpaceForSolution()) for i in range(len_g): # reset coefficients may be set at previous calls: args={} for n in self.__COEFFICIENTS: args[n]=Data() args['r']=Data() if g is None: # g_l=delta_{il} and len_f=len_g for n,v in coeffs_f: name=None if len_f > 1: val=v[:,i] else: val=v if n.startswith("d"): name='y'+n[1:] elif n.startswith("D"): name='Y'+n[1:] elif n.startswith("r"): name='r'+n[1:] if name: args[name]=val else: g_i=g[i] for n,v in coeffs_f: name=None if n.startswith("d"): name = 'y'+n[1:] val = self.__mm(v, g_i) elif n.startswith("r"): name= 'r' val = self.__mm(v, g_i) elif n.startswith("D"): name = 'Y'+n[1:] val = self.__mm(v, g_i) elif n.startswith("B") and isinstance(g_i, Data): name = 'Y'+n[1:] val = self.__mm(v, grad(g_i)) elif n.startswith("C"): name = 'X'+n[1:] val = matrix_multiply(v, g_i) elif n.startswith("A") and isinstance(g_i, Data): name = 'X'+n[1:] val = self.__mm(v, grad(g_i)) if name: if name in args: args[name]+=val else: args[name]=val self._lpde.setValue(**args) u_g_i=self._lpde.getSolution() if len_g >1: if self.getNumSolutions() == 1: u_g[i]=-u_g_i else: u_g[:,i]=-u_g_i else: u_g=-u_g_i return u_g def __mm(self, m, v): """ a sligtly crude matrix*matrix multiplication m is A-coefficient, u is vector, v is grad vector: A_ijkl*v_kl m is B-coefficient, u is vector, v is vector: B_ijk*v_k m is C-coefficient, u is vector, v is grad vector: C_ikl*v_kl m is D-coefficient, u is vector, v is vector: D_ij*v_j m is A-coefficient, u is scalar, v is grad vector: A_jkl*v_kl m is B-coefficient, u is scalar, v is vector: B_jk*v_k m is C-coefficient, u is scalar, v is grad vector: C_kl*v_kl m is D-coefficient, u is scalar, v is vector: D_j*v_j m is A-coefficient, u is vector, v is grad scalar: A_ijl*v_l m is B-coefficient, u is vector, v is scalar: B_ij*v m is C-coefficient, u is vector, v is grad scalar: C_il*v_l m is D-coefficient, u is vector, v is scalar: D_i*v m is A-coefficient, u is scalar, v is grad scalar: A_jl*v_l m is B-coefficient, u is scalar, v is scalar: B_j*v m is C-coefficient, u is scalar, v is grad scalar: C_l*v_l m is D-coefficient, u is scalar, v is scalar: D*v """ s_m=getShape(m) s_v=getShape(v) # m is B-coefficient, u is vector, v is scalar: B_ij*v # m is D-coefficient, u is vector, v is scalar: D_i*v # m is B-coefficient, u is scalar, v is scalar: B_j*v # m is D-coefficient, u is scalar, v is scalar: D*v if s_m == () or s_v == (): return m*v # m is D-coefficient, u is scalar, v is vector: D_j*v_j # m is C-coefficient, u is scalar, v is grad scalar: C_l*v_l if len(s_m) == 1: return inner(m,v) # m is D-coefficient, u is vector, v is vector: D_ij*v_j # m is B-coefficient, u is scalar, v is vector: B_jk*v_k # m is C-coefficient, u is vector, v is grad scalar: C_il*v_l # m is A-coefficient, u is scalar, v is grad scalar: A_jl*v_l if len(s_m) == 2 and len(s_v) == 1: return matrix_mult(m,v) # m is C-coefficient, u is scalar, v is grad vector: C_kl*v_kl if len(s_m) == 2 and len(s_v) == 2: return inner(m,v) # m is B-coefficient, u is vector, v is vector: B_ijk*v_k # m is A-coefficient, u is vector, v is grad scalar: A_ijl*v_l if len(s_m) == 3 and len(s_v) == 1: return matrix_mult(m,v) # m is A-coefficient, u is scalar, v is grad vector: A_jkl*v_kl # m is C-coefficient, u is vector, v is grad vector: C_ikl*v_kl if len(s_m) == 3 and len(s_v) == 2: return tensor_mult(m,v) # m is A-coefficient, u is vector, v is grad vector: A_ijkl*v_kl if len(s_m) == 4 and len(s_v) == 2: return tensor_mult(m,v) def _updateRHS(self, expressions, subs, **constants): """ """ ev=symb.Evaluator() names=[] for name in sorted(expressions): if name in self.__COEFFICIENTS: ev.addExpression(expressions[name]) names.append(name) if len(names)==0: return self.trace3("Starting expression evaluation.") T0=time() ev.subs(**subs) res=ev.evaluate() if len(names)==1: res=[res] self.trace3("RHS expressions evaluated in %f seconds."%(time()-T0)) if self._debug > self.DEBUG2: for i in range(len(names)): self.trace3("util.Lsup(%s)=%s"%(names[i],util.Lsup(res[i]))) args=dict(zip(names,res)) # reset coefficients may be set at previous calls: for n in self.__COEFFICIENTS: if not n in args and not n in constants: args[n]=Data() args=dict(list(args.items())+list(constants.items())) if not 'r' in args: args['r']=Data() self._lpde.setValue(**args) def _updateMatrix(self, expressions, subs): """ """ ev=symb.Evaluator() names=[] for name in sorted(expressions): if not name in self.__COEFFICIENTS: ev.addExpression(expressions[name]) names.append(name) if len(names)==0: return self.trace3("Starting expression evaluation.") T0=time() ev.subs(**subs) res=ev.evaluate() if len(names)==1: res=[res] self.trace3("Matrix expressions evaluated in %f seconds."%(time()-T0)) if self._debug > self.DEBUG2: for i in range(len(names)): self.trace3("util.Lsup(%s)=%s"%(names[i],util.Lsup(res[i]))) self._lpde.setValue(**dict(zip(names,res))) def _updateLinearPDE(self, expressions, subs, **constants): self._updateMatrix(expressions, subs) self._updateRHS(expressions, subs, **constants)