/[escript]/trunk/escript/py_src/linearPDEs.py
ViewVC logotype

Diff of /trunk/escript/py_src/linearPDEs.py

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 1137 by gross, Thu May 10 08:11:31 2007 UTC revision 1659 by gross, Fri Jul 18 02:28:13 2008 UTC
# Line 1  Line 1 
1    #
2  # $Id$  # $Id$
3    #
4    #######################################################
5    #
6    #           Copyright 2003-2007 by ACceSS MNRF
7    #       Copyright 2007 by University of Queensland
8    #
9    #                http://esscc.uq.edu.au
10    #        Primary Business: Queensland, Australia
11    #  Licensed under the Open Software License version 3.0
12    #     http://www.opensource.org/licenses/osl-3.0.php
13    #
14    #######################################################
15    #
16    
17  """  """
18  The module provides an interface to define and solve linear partial  The module provides an interface to define and solve linear partial
19  differential equations (PDEs) within L{escript}. L{linearPDEs} does not provide any  differential equations (PDEs) within L{escript}. L{linearPDEs} does not provide any
# Line 19  to define of solve these sepecial PDEs. Line 34  to define of solve these sepecial PDEs.
34  @var __date__: date of the version  @var __date__: date of the version
35  """  """
36    
37    import math
38  import escript  import escript
39  import util  import util
40  import numarray  import numarray
# Line 112  class PDECoefficient(object): Line 128  class PDECoefficient(object):
128         @param reduced: indicates if reduced         @param reduced: indicates if reduced
129         @type reduced: C{bool}         @type reduced: C{bool}
130         """         """
         
131         super(PDECoefficient, self).__init__()         super(PDECoefficient, self).__init__()
132         self.what=where         self.what=where
133         self.pattern=pattern         self.pattern=pattern
# Line 365  class LinearPDE(object): Line 380  class LinearPDE(object):
380    
381     The PDE is symmetrical if     The PDE is symmetrical if
382    
383     M{A[i,j]=A[j,i]}  and M{B[j]=C[j]} and M{A_reduced[i,j]=A_reduced[j,i]}  and M{B_reduced[j]=C_reduced[j]     M{A[i,j]=A[j,i]}  and M{B[j]=C[j]} and M{A_reduced[i,j]=A_reduced[j,i]}  and M{B_reduced[j]=C_reduced[j]}
384    
385     For a system of PDEs and a solution with several components the PDE has the form     For a system of PDEs and a solution with several components the PDE has the form
386    
# Line 416  class LinearPDE(object): Line 431  class LinearPDE(object):
431     of the solution at side 1 and at side 0, denotes the jump of M{u} across discontinuity along the normal calcualted by     of the solution at side 1 and at side 0, denotes the jump of M{u} across discontinuity along the normal calcualted by
432     L{jump<util.jump>}.     L{jump<util.jump>}.
433     The coefficient M{d_contact} is a rank two and M{y_contact} is a rank one both in the L{FunctionOnContactZero<escript.FunctionOnContactZero>} or L{FunctionOnContactOne<escript.FunctionOnContactOne>}.     The coefficient M{d_contact} is a rank two and M{y_contact} is a rank one both in the L{FunctionOnContactZero<escript.FunctionOnContactZero>} or L{FunctionOnContactOne<escript.FunctionOnContactOne>}.
434      The coefficient M{d_contact_reduced} is a rank two and M{y_contact_reduced} is a rank one both in the L{ReducedFunctionOnContactZero<escript.ReducedFunctionOnContactZero>} or L{ReducedFunctionOnContactOne<escript.ReducedFunctionOnContactOne>}.     The coefficient M{d_contact_reduced} is a rank two and M{y_contact_reduced} is a rank one both in the L{ReducedFunctionOnContactZero<escript.ReducedFunctionOnContactZero>} or L{ReducedFunctionOnContactOne<escript.ReducedFunctionOnContactOne>}.
435     In case of a single PDE and a single component solution the contact condition takes the form     In case of a single PDE and a single component solution the contact condition takes the form
436    
437     M{n[j]*J0_{j}=n[j]*J1_{j}=(y_contact+y_contact_reduced)-(d_contact+y_contact_reduced)*jump(u)}     M{n[j]*J0_{j}=n[j]*J1_{j}=(y_contact+y_contact_reduced)-(d_contact+y_contact_reduced)*jump(u)}
# Line 444  class LinearPDE(object): Line 459  class LinearPDE(object):
459     @cvar SCSL: SGI SCSL solver library     @cvar SCSL: SGI SCSL solver library
460     @cvar MKL: Intel's MKL solver library     @cvar MKL: Intel's MKL solver library
461     @cvar UMFPACK: the UMFPACK library     @cvar UMFPACK: the UMFPACK library
462       @cvar TRILINOS: the TRILINOS parallel solver class library from Sandia Natl Labs
463     @cvar ITERATIVE: The default iterative solver     @cvar ITERATIVE: The default iterative solver
464     @cvar AMG: algebraic multi grid     @cvar AMG: algebraic multi grid
465     @cvar RILU: recursive ILU     @cvar RILU: recursive ILU
# Line 473  class LinearPDE(object): Line 489  class LinearPDE(object):
489     PASO= 21     PASO= 21
490     AMG= 22     AMG= 22
491     RILU = 23     RILU = 23
492       TRILINOS = 24
493       NONLINEAR_GMRES = 25
494    
495     SMALL_TOLERANCE=1.e-13     SMALL_TOLERANCE=1.e-13
496     __PACKAGE_KEY="package"     __PACKAGE_KEY="package"
# Line 935  class LinearPDE(object): Line 953  class LinearPDE(object):
953         @param preconditioner: sets a new solver method.         @param preconditioner: sets a new solver method.
954         @type preconditioner: one of L{DEFAULT}, L{JACOBI} L{ILU0}, L{ILUT},L{SSOR}, L{RILU}         @type preconditioner: one of L{DEFAULT}, L{JACOBI} L{ILU0}, L{ILUT},L{SSOR}, L{RILU}
955         """         """
956         if solver==None: solve=self.DEFAULT         if solver==None: solver=self.__solver_method
957           if preconditioner==None: preconditioner=self.__preconditioner
958           if solver==None: solver=self.DEFAULT
959         if preconditioner==None: preconditioner=self.DEFAULT         if preconditioner==None: preconditioner=self.DEFAULT
960         if not (solver,preconditioner)==self.getSolverMethod():         if not (solver,preconditioner)==self.getSolverMethod():
961             self.__solver_method=solver             self.__solver_method=solver
# Line 979  class LinearPDE(object): Line 999  class LinearPDE(object):
999         elif p==self.MKL: package= "MKL"         elif p==self.MKL: package= "MKL"
1000         elif p==self.SCSL: package= "SCSL"         elif p==self.SCSL: package= "SCSL"
1001         elif p==self.UMFPACK: package= "UMFPACK"         elif p==self.UMFPACK: package= "UMFPACK"
1002           elif p==self.TRILINOS: package= "TRILINOS"
1003         else : method="unknown"         else : method="unknown"
1004         return "%s solver of %s package"%(method,package)         return "%s solver of %s package"%(method,package)
1005    
# Line 997  class LinearPDE(object): Line 1018  class LinearPDE(object):
1018         sets a new solver package         sets a new solver package
1019    
1020         @param package: sets a new solver method.         @param package: sets a new solver method.
1021         @type package: one of L{DEFAULT}, L{PASO} L{SCSL}, L{MKL}, L{UMFPACK}         @type package: one of L{DEFAULT}, L{PASO} L{SCSL}, L{MKL}, L{UMFPACK}, L{TRILINOS}
1022         """         """
1023         if package==None: package=self.DEFAULT         if package==None: package=self.DEFAULT
1024         if not package==self.getSolverPackage():         if not package==self.getSolverPackage():
# Line 1682  class LinearPDE(object): Line 1703  class LinearPDE(object):
1703                        raise ValueError,"coefficient B in lumped matrix may not be present."                        raise ValueError,"coefficient B in lumped matrix may not be present."
1704                   if not self.getCoefficientOfGeneralPDE("C").isEmpty():                   if not self.getCoefficientOfGeneralPDE("C").isEmpty():
1705                        raise ValueError,"coefficient C in lumped matrix may not be present."                        raise ValueError,"coefficient C in lumped matrix may not be present."
1706                     if not self.getCoefficientOfGeneralPDE("d_contact").isEmpty():
1707                          raise ValueError,"coefficient d_contact in lumped matrix may not be present."
1708                   if not self.getCoefficientOfGeneralPDE("A_reduced").isEmpty():                   if not self.getCoefficientOfGeneralPDE("A_reduced").isEmpty():
1709                        raise ValueError,"coefficient A_reduced in lumped matrix may not be present."                        raise ValueError,"coefficient A_reduced in lumped matrix may not be present."
1710                   if not self.getCoefficientOfGeneralPDE("B_reduced").isEmpty():                   if not self.getCoefficientOfGeneralPDE("B_reduced").isEmpty():
1711                        raise ValueError,"coefficient B_reduced in lumped matrix may not be present."                        raise ValueError,"coefficient B_reduced in lumped matrix may not be present."
1712                   if not self.getCoefficientOfGeneralPDE("C_reduced").isEmpty():                   if not self.getCoefficientOfGeneralPDE("C_reduced").isEmpty():
1713                        raise ValueError,"coefficient C_reduced in lumped matrix may not be present."                        raise ValueError,"coefficient C_reduced in lumped matrix may not be present."
1714                     if not self.getCoefficientOfGeneralPDE("d_contact_reduced").isEmpty():
1715                          raise ValueError,"coefficient d_contact_reduced in lumped matrix may not be present."
1716                   D=self.getCoefficientOfGeneralPDE("D")                   D=self.getCoefficientOfGeneralPDE("D")
1717                     d=self.getCoefficientOfGeneralPDE("d")
1718                     D_reduced=self.getCoefficientOfGeneralPDE("D_reduced")
1719                     d_reduced=self.getCoefficientOfGeneralPDE("d_reduced")
1720                   if not D.isEmpty():                   if not D.isEmpty():
1721                       if self.getNumSolutions()>1:                       if self.getNumSolutions()>1:
1722                          D_times_e=util.matrix_mult(D,numarray.ones((self.getNumSolutions(),)))                          D_times_e=util.matrix_mult(D,numarray.ones((self.getNumSolutions(),)))
# Line 1696  class LinearPDE(object): Line 1724  class LinearPDE(object):
1724                          D_times_e=D                          D_times_e=D
1725                   else:                   else:
1726                      D_times_e=escript.Data()                      D_times_e=escript.Data()
                  d=self.getCoefficientOfGeneralPDE("d")  
1727                   if not d.isEmpty():                   if not d.isEmpty():
1728                       if self.getNumSolutions()>1:                       if self.getNumSolutions()>1:
1729                          d_times_e=util.matrix_mult(d,numarray.ones((self.getNumSolutions(),)))                          d_times_e=util.matrix_mult(d,numarray.ones((self.getNumSolutions(),)))
# Line 1704  class LinearPDE(object): Line 1731  class LinearPDE(object):
1731                          d_times_e=d                          d_times_e=d
1732                   else:                   else:
1733                      d_times_e=escript.Data()                      d_times_e=escript.Data()
1734                   d_contact=self.getCoefficientOfGeneralPDE("d_contact")        
                  if not d_contact.isEmpty():  
                      if self.getNumSolutions()>1:  
                         d_contact_times_e=util.matrixmult(d_contact,numarray.ones((self.getNumSolutions(),)))  
                      else:  
                         d_contact_times_e=d_contact  
                  else:  
                     d_contact_times_e=escript.Data()  
       
                  self.__operator=self.__getNewRightHandSide()  
                  self.getDomain().addPDEToRHS(self.__operator, \  
                                               escript.Data(), \  
                                               D_times_e, \  
                                               d_times_e,\  
                                               d_contact_times_e)  
                  D_reduced=self.getCoefficientOfGeneralPDE("D_reduced")  
1735                   if not D_reduced.isEmpty():                   if not D_reduced.isEmpty():
1736                       if self.getNumSolutions()>1:                       if self.getNumSolutions()>1:
1737                          D_reduced_times_e=util.matrix_mult(D_reduced,numarray.ones((self.getNumSolutions(),)))                          D_reduced_times_e=util.matrix_mult(D_reduced,numarray.ones((self.getNumSolutions(),)))
# Line 1727  class LinearPDE(object): Line 1739  class LinearPDE(object):
1739                          D_reduced_times_e=D_reduced                          D_reduced_times_e=D_reduced
1740                   else:                   else:
1741                      D_reduced_times_e=escript.Data()                      D_reduced_times_e=escript.Data()
                  d_reduced=self.getCoefficientOfGeneralPDE("d_reduced")  
1742                   if not d_reduced.isEmpty():                   if not d_reduced.isEmpty():
1743                       if self.getNumSolutions()>1:                       if self.getNumSolutions()>1:
1744                          d_reduced_times_e=util.matrix_mult(d_reduced,numarray.ones((self.getNumSolutions(),)))                          d_reduced_times_e=util.matrix_mult(d_reduced,numarray.ones((self.getNumSolutions(),)))
# Line 1735  class LinearPDE(object): Line 1746  class LinearPDE(object):
1746                          d_reduced_times_e=d_reduced                          d_reduced_times_e=d_reduced
1747                   else:                   else:
1748                      d_reduced_times_e=escript.Data()                      d_reduced_times_e=escript.Data()
1749                   d_contact_reduced=self.getCoefficientOfGeneralPDE("d_contact_reduced")  
                  if not d_contact_reduced.isEmpty():  
                      if self.getNumSolutions()>1:  
                         d_contact_reduced_times_e=util.matrixmult(d_contact_reduced,numarray.ones((self.getNumSolutions(),)))  
                      else:  
                         d_contact_reduced_times_e=d_contact_reduced  
                  else:  
                     d_contact_reduced_times_e=escript.Data()  
       
1750                   self.__operator=self.__getNewRightHandSide()                   self.__operator=self.__getNewRightHandSide()
1751                   self.getDomain().addPDEToRHS(self.__operator, \                   if False and hasattr(self.getDomain(), "addPDEToLumpedSystem") :
1752                                                escript.Data(), \                      self.getDomain().addPDEToLumpedSystem(self.__operator, D_times_e, d_times_e)
1753                                                D_times_e, \                      self.getDomain().addPDEToLumpedSystem(self.__operator, D_reduced_times_e, d_reduced_times_e)
1754                                                d_times_e,\                   else:
1755                                                d_contact_times_e)                      self.getDomain().addPDEToRHS(self.__operator, \
1756                   self.getDomain().addPDEToRHS(self.__operator, \                                                   escript.Data(), \
1757                                                escript.Data(), \                                                   D_times_e, \
1758                                                D_reduced_times_e, \                                                   d_times_e,\
1759                                                d_reduced_times_e,\                                                   escript.Data())
1760                                                d_contact_reduced_times_e)                      self.getDomain().addPDEToRHS(self.__operator, \
1761                                                     escript.Data(), \
1762                                                     D_reduced_times_e, \
1763                                                     d_reduced_times_e,\
1764                                                     escript.Data())
1765                        print "RHS:",util.inf(self.__operator),util.sup(self.__operator)
1766                   self.__operator=1./self.__operator                   self.__operator=1./self.__operator
1767                   self.trace("New lumped operator has been built.")                   self.trace("New lumped operator has been built.")
1768                   self.__operator_is_Valid=True                   self.__operator_is_Valid=True
# Line 2185  class LameEquation(LinearPDE): Line 2193  class LameEquation(LinearPDE):
2193       else:       else:
2194          raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name          raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name
2195    
2196    def LinearSinglePDE(domain,debug=False):
2197       """
2198       defines a single linear PDEs
2199    
2200       @param domain: domain of the PDE
2201       @type domain: L{Domain<escript.Domain>}
2202       @param debug: if True debug informations are printed.
2203       @rtype: L{LinearPDE}
2204       """
2205       return LinearPDE(domain,numEquations=1,numSolutions=1,debug=debug)
2206    
2207    def LinearPDESystem(domain,debug=False):
2208       """
2209       defines a system of linear PDEs
2210    
2211       @param domain: domain of the PDE
2212       @type domain: L{Domain<escript.Domain>}
2213       @param debug: if True debug informations are printed.
2214       @rtype: L{LinearPDE}
2215       """
2216       return LinearPDE(domain,numEquations=domain.getDim(),numSolutions=domain.getDim(),debug=debug)
2217    
2218    class TransportPDE(object):
2219         """
2220         Warning: This is still a very experimental. The class is still changing!
2221    
2222         Mu_{,t} =-(A_{ij}u_{,j})_j-(B_{j}u)_{,j} + C_{j} u_{,j} + Y_i + X_{i,i}
2223        
2224         u=r where q>0
2225        
2226         all coefficients are constant over time.
2227    
2228         typical usage:
2229    
2230             p=TransportPDE(dom)
2231             p.setValue(M=Scalar(1.,Function(dom),C=Scalar(1.,Function(dom)*[-1.,0.])
2232             p.setInitialSolution(u=exp(-length(dom.getX()-[0.1,0.1])**2)
2233             t=0
2234             dt=0.1
2235             while (t<1.):
2236                  u=p.solve(dt)
2237    
2238         """
2239         def __init__(self,domain,num_equations=1,theta=0.5,useSUPG=False,trace=True):
2240            self.__domain=domain
2241            self.__num_equations=num_equations
2242            self.__useSUPG=useSUPG
2243            self.__trace=trace
2244            self.__theta=theta
2245            self.__matrix_type=0
2246            self.__reduced=True
2247            self.__reassemble=True
2248            if self.__useSUPG:
2249               self.__pde=LinearPDE(domain,numEquations=num_equations,numSolutions=num_equations,debug=trace)
2250               self.__pde.setSymmetryOn()
2251               self.__pde.setReducedOrderOn()
2252            else:
2253               self.__transport_problem=self.__getNewTransportProblem()
2254            self.setTolerance()
2255            self.__M=escript.Data()
2256            self.__A=escript.Data()
2257            self.__B=escript.Data()
2258            self.__C=escript.Data()
2259            self.__D=escript.Data()
2260            self.__X=escript.Data()
2261            self.__Y=escript.Data()
2262            self.__d=escript.Data()
2263            self.__y=escript.Data()
2264            self.__d_contact=escript.Data()
2265            self.__y_contact=escript.Data()
2266            self.__r=escript.Data()
2267            self.__q=escript.Data()
2268    
2269         def trace(self,text):
2270                 if self.__trace: print text
2271         def getSafeTimeStepSize(self):
2272            if self.__useSUPG:
2273                if self.__reassemble:
2274                   h=self.__domain.getSize()
2275                   dt=None
2276                   if not self.__A.isEmpty():
2277                      dt2=util.inf(h**2*self.__M/util.length(self.__A))
2278                      if dt == None:
2279                         dt = dt2
2280                      else:
2281                         dt=1./(1./dt+1./dt2)
2282                   if not self.__B.isEmpty():
2283                      dt2=util.inf(h*self.__M/util.length(self.__B))
2284                      if dt == None:
2285                         dt = dt2
2286                      else:
2287                         dt=1./(1./dt+1./dt2)
2288                   if not  self.__C.isEmpty():
2289                      dt2=util.inf(h*self.__M/util.length(self.__C))
2290                      if dt == None:
2291                         dt = dt2
2292                      else:
2293                         dt=1./(1./dt+1./dt2)
2294                   if not self.__D.isEmpty():
2295                      dt2=util.inf(self.__M/util.length(self.__D))
2296                      if dt == None:
2297                         dt = dt2
2298                      else:
2299                         dt=1./(1./dt+1./dt2)
2300                   self.__dt = dt/2
2301                return self.__dt
2302            else:
2303                return self.__getTransportProblem().getSafeTimeStepSize()
2304         def getDomain(self):
2305            return self.__domain
2306         def getTheta(self):
2307            return self.__theta
2308         def getNumEquations(self):
2309            return self.__num_equations
2310         def setReducedOn(self):
2311              if not self.reduced():
2312                  if self.__useSUPG:
2313                     self.__pde.setReducedOrderOn()
2314                  else:
2315                     self.__transport_problem=self.__getNewTransportProblem()
2316              self.__reduced=True
2317         def setReducedOff(self):
2318              if self.reduced():
2319                  if self.__useSUPG:
2320                     self.__pde.setReducedOrderOff()
2321                  else:
2322                     self.__transport_problem=self.__getNewTransportProblem()
2323              self.__reduced=False
2324         def reduced(self):
2325             return self.__reduced
2326         def getFunctionSpace(self):
2327            if self.reduced():
2328               return escript.ReducedSolution(self.getDomain())
2329            else:
2330               return escript.Solution(self.getDomain())
2331    
2332         def setTolerance(self,tol=1.e-8):
2333            self.__tolerance=tol
2334            if self.__useSUPG:
2335                  self.__pde.setTolerance(self.__tolerance)
2336    
2337         def __getNewTransportProblem(self):
2338           """
2339           returns an instance of a new operator
2340           """
2341           self.trace("New Transport problem is allocated.")
2342           return self.getDomain().newTransportProblem( \
2343                                   self.getTheta(),
2344                                   self.getNumEquations(), \
2345                                   self.getFunctionSpace(), \
2346                                   self.__matrix_type)
2347              
2348         def __getNewSolutionVector(self):
2349             if self.getNumEquations() ==1 :
2350                    out=escript.Data(0.0,(),self.getFunctionSpace())
2351             else:
2352                    out=escript.Data(0.0,(self.getNumEquations(),),self.getFunctionSpace())
2353             return out
2354    
2355         def __getTransportProblem(self):
2356           if self.__reassemble:
2357                 self.__source=self.__getNewSolutionVector()
2358                 self.__transport_problem.reset()
2359                 self.getDomain().addPDEToTransportProblem(
2360                             self.__transport_problem,
2361                             self.__source,
2362                             self.__M,
2363                             self.__A,
2364                             self.__B,
2365                             self.__C,
2366                             self.__D,
2367                             self.__X,
2368                             self.__Y,
2369                             self.__d,
2370                             self.__y,
2371                             self.__d_contact,
2372                             self.__y_contact)
2373                 self.__transport_problem.insertConstraint(self.__source,self.__q,self.__r)
2374                 self.__reassemble=False
2375           return self.__transport_problem
2376         def setValue(self,M=None, A=None, B=None, C=None, D=None, X=None, Y=None,
2377                      d=None, y=None, d_contact=None, y_contact=None, q=None, r=None):
2378                 if not M==None:
2379                      self.__reassemble=True
2380                      self.__M=M
2381                 if not A==None:
2382                      self.__reassemble=True
2383                      self.__A=A
2384                 if not B==None:
2385                      self.__reassemble=True
2386                      self.__B=B
2387                 if not C==None:
2388                      self.__reassemble=True
2389                      self.__C=C
2390                 if not D==None:
2391                      self.__reassemble=True
2392                      self.__D=D
2393                 if not X==None:
2394                      self.__reassemble=True
2395                      self.__X=X
2396                 if not Y==None:
2397                      self.__reassemble=True
2398                      self.__Y=Y
2399                 if not d==None:
2400                      self.__reassemble=True
2401                      self.__d=d
2402                 if not y==None:
2403                      self.__reassemble=True
2404                      self.__y=y
2405                 if not d_contact==None:
2406                      self.__reassemble=True
2407                      self.__d_contact=d_contact
2408                 if not y_contact==None:
2409                      self.__reassemble=True
2410                      self.__y_contact=y_contact
2411                 if not q==None:
2412                      self.__reassemble=True
2413                      self.__q=q
2414                 if not r==None:
2415                      self.__reassemble=True
2416                      self.__r=r
2417    
2418         def setInitialSolution(self,u):
2419                 if self.__useSUPG:
2420                     self.__u=util.interpolate(u,self.getFunctionSpace())
2421                 else:
2422                     self.__transport_problem.setInitialValue(util.interpolate(u,self.getFunctionSpace()))
2423    
2424         def solve(self,dt,**kwarg):
2425               if self.__useSUPG:
2426                    if self.__reassemble:
2427                        self.__pde.setValue(D=self.__M,d=self.__d,d_contact=self.__d_contact,q=self.__q) # ,r=self.__r)
2428                        self.__reassemble=False
2429                    dt2=self.getSafeTimeStepSize()
2430                    nn=max(math.ceil(dt/self.getSafeTimeStepSize()),1.)
2431                    dt2=dt/nn
2432                    nnn=0
2433                    u=self.__u
2434                    self.trace("number of substeps is %d."%nn)
2435                    while nnn<nn :
2436                        self.__setSUPG(u,u,dt2/2)
2437                        u_half=self.__pde.getSolution(verbose=True)
2438                        self.__setSUPG(u,u_half,dt2)
2439                        u=self.__pde.getSolution(verbose=True)
2440                        nnn+=1
2441                    self.__u=u
2442                    return self.__u
2443               else:
2444                   kwarg["tolerance"]=self.__tolerance
2445                   tp=self.__getTransportProblem()
2446                   return tp.solve(self.__source,dt,kwarg)
2447         def __setSUPG(self,u0,u,dt):
2448                g=util.grad(u)
2449                X=0
2450                Y=self.__M*u0
2451                X=0
2452                self.__pde.setValue(r=u0)
2453                if not self.__A.isEmpty():
2454                   X=X+dt*util.matrixmult(self.__A,g)
2455                if not self.__B.isEmpty():
2456                   X=X+dt*self.__B*u
2457                if not  self.__C.isEmpty():
2458                   Y=Y+dt*util.inner(self.__C,g)
2459                if not self.__D.isEmpty():
2460                   Y=Y+dt*self.__D*u
2461                if not self.__X.isEmpty():
2462                   X=X+dt*self.__X
2463                if not self.__Y.isEmpty():
2464                   Y=Y+dt*self.__Y
2465                self.__pde.setValue(X=X,Y=Y)
2466                if not self.__y.isEmpty():
2467                   self.__pde.setValue(y=dt*self.__y)
2468                if not self.__y_contact.isEmpty():
2469                   self.__pde.setValue(y=dt*self.__y_contact)
2470                self.__pde.setValue(r=u0)

Legend:
Removed from v.1137  
changed lines
  Added in v.1659

  ViewVC Help
Powered by ViewVC 1.1.26