/[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 1400 by gross, Thu Jan 24 06:04:31 2008 UTC revision 1809 by ksteube, Thu Sep 25 06:43:44 2008 UTC
# Line 1  Line 1 
1    
2    ########################################################
3  #  #
4  # $Id$  # Copyright (c) 2003-2008 by University of Queensland
5  #  # Earth Systems Science Computational Center (ESSCC)
6  #######################################################  # http://www.uq.edu.au/esscc
 #  
 #           Copyright 2003-2007 by ACceSS MNRF  
 #       Copyright 2007 by University of Queensland  
 #  
 #                http://esscc.uq.edu.au  
 #        Primary Business: Queensland, Australia  
 #  Licensed under the Open Software License version 3.0  
 #     http://www.opensource.org/licenses/osl-3.0.php  
7  #  #
8  #######################################################  # Primary Business: Queensland, Australia
9    # Licensed under the Open Software License version 3.0
10    # http://www.opensource.org/licenses/osl-3.0.php
11  #  #
12    ########################################################
13    
14    __copyright__="""Copyright (c) 2003-2008 by University of Queensland
15    Earth Systems Science Computational Center (ESSCC)
16    http://www.uq.edu.au/esscc
17    Primary Business: Queensland, Australia"""
18    __license__="""Licensed under the Open Software License version 3.0
19    http://www.opensource.org/licenses/osl-3.0.php"""
20    __url__="http://www.uq.edu.au/esscc/escript-finley"
21    
22  """  """
23  The module provides an interface to define and solve linear partial  The module provides an interface to define and solve linear partial
# Line 34  to define of solve these sepecial PDEs. Line 39  to define of solve these sepecial PDEs.
39  @var __date__: date of the version  @var __date__: date of the version
40  """  """
41    
42    import math
43  import escript  import escript
44  import util  import util
45  import numarray  import numarray
46    
47  __author__="Lutz Gross, l.gross@uq.edu.au"  __author__="Lutz Gross, l.gross@uq.edu.au"
 __copyright__="""  Copyright (c) 2006 by ACcESS MNRF  
                     http://www.access.edu.au  
                 Primary Business: Queensland, Australia"""  
 __license__="""Licensed under the Open Software License version 3.0  
              http://www.opensource.org/licenses/osl-3.0.php"""  
 __url__="http://www.iservo.edu.au/esys"  
 __version__="$Revision$"  
 __date__="$Date$"  
48    
49    
50  class IllegalCoefficient(ValueError):  class IllegalCoefficient(ValueError):
# Line 444  class LinearPDE(object): Line 442  class LinearPDE(object):
442     @cvar CR: The conjugate residual method     @cvar CR: The conjugate residual method
443     @cvar CGS: The conjugate gardient square method     @cvar CGS: The conjugate gardient square method
444     @cvar BICGSTAB: The stabilized BiConjugate Gradient method.     @cvar BICGSTAB: The stabilized BiConjugate Gradient method.
445       @cvar TFQMR: Transport Free Quasi Minimal Residual method.
446       @cvar MINRES: Minimum residual method.
447     @cvar SSOR: The symmetric overrealaxtion method     @cvar SSOR: The symmetric overrealaxtion method
448     @cvar ILU0: The incomplete LU factorization preconditioner  with no fill in     @cvar ILU0: The incomplete LU factorization preconditioner  with no fill in
449     @cvar ILUT: The incomplete LU factorization preconditioner with will in     @cvar ILUT: The incomplete LU factorization preconditioner with will in
# Line 489  class LinearPDE(object): Line 489  class LinearPDE(object):
489     AMG= 22     AMG= 22
490     RILU = 23     RILU = 23
491     TRILINOS = 24     TRILINOS = 24
492       NONLINEAR_GMRES = 25
493       TFQMR = 26
494       MINRES = 27
495    
496     SMALL_TOLERANCE=1.e-13     SMALL_TOLERANCE=1.e-13
497     __PACKAGE_KEY="package"     __PACKAGE_KEY="package"
# Line 947  class LinearPDE(object): Line 950  class LinearPDE(object):
950         sets a new solver         sets a new solver
951    
952         @param solver: sets a new solver method.         @param solver: sets a new solver method.
953         @type solver: one of L{DEFAULT}, L{ITERATIVE} L{DIRECT}, L{CHOLEVSKY}, L{PCG}, L{CR}, L{CGS}, L{BICGSTAB}, L{SSOR}, L{GMRES}, L{PRES20}, L{LUMPING}, L{AMG}         @type solver: one of L{DEFAULT}, L{ITERATIVE} L{DIRECT}, L{CHOLEVSKY}, L{PCG}, L{CR}, L{CGS}, L{BICGSTAB}, L{SSOR}, L{GMRES}, L{TFQMR}, L{MINRES}, L{PRES20}, L{LUMPING}, L{AMG}
954         @param preconditioner: sets a new solver method.         @param preconditioner: sets a new solver method.
955         @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}
956         """         """
# Line 977  class LinearPDE(object): Line 980  class LinearPDE(object):
980         elif m[0]==self.ITERATIVE: method= "ITERATIVE"         elif m[0]==self.ITERATIVE: method= "ITERATIVE"
981         elif m[0]==self.CHOLEVSKY: method= "CHOLEVSKY"         elif m[0]==self.CHOLEVSKY: method= "CHOLEVSKY"
982         elif m[0]==self.PCG: method= "PCG"         elif m[0]==self.PCG: method= "PCG"
983           elif m[0]==self.TFQMR: method= "TFQMR"
984           elif m[0]==self.MINRES: method= "MINRES"
985         elif m[0]==self.CR: method= "CR"         elif m[0]==self.CR: method= "CR"
986         elif m[0]==self.CGS: method= "CGS"         elif m[0]==self.CGS: method= "CGS"
987         elif m[0]==self.BICGSTAB: method= "BICGSTAB"         elif m[0]==self.BICGSTAB: method= "BICGSTAB"
# Line 1746  class LinearPDE(object): Line 1751  class LinearPDE(object):
1751                      d_reduced_times_e=escript.Data()                      d_reduced_times_e=escript.Data()
1752    
1753                   self.__operator=self.__getNewRightHandSide()                   self.__operator=self.__getNewRightHandSide()
1754                   if hasattr(self.getDomain(), "addPDEToLumpedSystem") :                   if False and hasattr(self.getDomain(), "addPDEToLumpedSystem") :
1755                      self.getDomain().addPDEToLumpedSystem(self.__operator, D_times_e, d_times_e)                      self.getDomain().addPDEToLumpedSystem(self.__operator, D_times_e, d_times_e)
1756                      self.getDomain().addPDEToLumpedSystem(self.__operator, D_reduced_times_e, d_reduced_times_e)                      self.getDomain().addPDEToLumpedSystem(self.__operator, D_reduced_times_e, d_reduced_times_e)
1757                   else:                   else:
# Line 2211  def LinearPDESystem(domain,debug=False): Line 2216  def LinearPDESystem(domain,debug=False):
2216     @rtype: L{LinearPDE}     @rtype: L{LinearPDE}
2217     """     """
2218     return LinearPDE(domain,numEquations=domain.getDim(),numSolutions=domain.getDim(),debug=debug)     return LinearPDE(domain,numEquations=domain.getDim(),numSolutions=domain.getDim(),debug=debug)
2219    
2220    class TransportPDE(object):
2221         """
2222         Warning: This is still a very experimental. The class is still changing!
2223    
2224         Mu_{,t} =-(A_{ij}u_{,j})_j-(B_{j}u)_{,j} + C_{j} u_{,j} + Y_i + X_{i,i}
2225        
2226         u=r where q>0
2227        
2228         all coefficients are constant over time.
2229    
2230         typical usage:
2231    
2232             p=TransportPDE(dom)
2233             p.setValue(M=Scalar(1.,Function(dom),C=Scalar(1.,Function(dom)*[-1.,0.])
2234             p.setInitialSolution(u=exp(-length(dom.getX()-[0.1,0.1])**2)
2235             t=0
2236             dt=0.1
2237             while (t<1.):
2238                  u=p.solve(dt)
2239    
2240         """
2241         def __init__(self,domain,num_equations=1,theta=0.5,useSUPG=False,trace=True):
2242            self.__domain=domain
2243            self.__num_equations=num_equations
2244            self.__useSUPG=useSUPG
2245            self.__trace=trace
2246            self.__theta=theta
2247            self.__matrix_type=0
2248            self.__reduced=True
2249            self.__reassemble=True
2250            if self.__useSUPG:
2251               self.__pde=LinearPDE(domain,numEquations=num_equations,numSolutions=num_equations,debug=trace)
2252               self.__pde.setSymmetryOn()
2253               self.__pde.setReducedOrderOn()
2254            else:
2255               self.__transport_problem=self.__getNewTransportProblem()
2256            self.setTolerance()
2257            self.__M=escript.Data()
2258            self.__A=escript.Data()
2259            self.__B=escript.Data()
2260            self.__C=escript.Data()
2261            self.__D=escript.Data()
2262            self.__X=escript.Data()
2263            self.__Y=escript.Data()
2264            self.__d=escript.Data()
2265            self.__y=escript.Data()
2266            self.__d_contact=escript.Data()
2267            self.__y_contact=escript.Data()
2268            self.__r=escript.Data()
2269            self.__q=escript.Data()
2270    
2271         def trace(self,text):
2272                 if self.__trace: print text
2273         def getSafeTimeStepSize(self):
2274            if self.__useSUPG:
2275                if self.__reassemble:
2276                   h=self.__domain.getSize()
2277                   dt=None
2278                   if not self.__A.isEmpty():
2279                      dt2=util.inf(h**2*self.__M/util.length(self.__A))
2280                      if dt == None:
2281                         dt = dt2
2282                      else:
2283                         dt=1./(1./dt+1./dt2)
2284                   if not self.__B.isEmpty():
2285                      dt2=util.inf(h*self.__M/util.length(self.__B))
2286                      if dt == None:
2287                         dt = dt2
2288                      else:
2289                         dt=1./(1./dt+1./dt2)
2290                   if not  self.__C.isEmpty():
2291                      dt2=util.inf(h*self.__M/util.length(self.__C))
2292                      if dt == None:
2293                         dt = dt2
2294                      else:
2295                         dt=1./(1./dt+1./dt2)
2296                   if not self.__D.isEmpty():
2297                      dt2=util.inf(self.__M/util.length(self.__D))
2298                      if dt == None:
2299                         dt = dt2
2300                      else:
2301                         dt=1./(1./dt+1./dt2)
2302                   self.__dt = dt/2
2303                return self.__dt
2304            else:
2305                return self.__getTransportProblem().getSafeTimeStepSize()
2306         def getDomain(self):
2307            return self.__domain
2308         def getTheta(self):
2309            return self.__theta
2310         def getNumEquations(self):
2311            return self.__num_equations
2312         def setReducedOn(self):
2313              if not self.reduced():
2314                  if self.__useSUPG:
2315                     self.__pde.setReducedOrderOn()
2316                  else:
2317                     self.__transport_problem=self.__getNewTransportProblem()
2318              self.__reduced=True
2319         def setReducedOff(self):
2320              if self.reduced():
2321                  if self.__useSUPG:
2322                     self.__pde.setReducedOrderOff()
2323                  else:
2324                     self.__transport_problem=self.__getNewTransportProblem()
2325              self.__reduced=False
2326         def reduced(self):
2327             return self.__reduced
2328         def getFunctionSpace(self):
2329            if self.reduced():
2330               return escript.ReducedSolution(self.getDomain())
2331            else:
2332               return escript.Solution(self.getDomain())
2333    
2334         def setTolerance(self,tol=1.e-8):
2335            self.__tolerance=tol
2336            if self.__useSUPG:
2337                  self.__pde.setTolerance(self.__tolerance)
2338    
2339         def __getNewTransportProblem(self):
2340           """
2341           returns an instance of a new operator
2342           """
2343           self.trace("New Transport problem is allocated.")
2344           return self.getDomain().newTransportProblem( \
2345                                   self.getTheta(),
2346                                   self.getNumEquations(), \
2347                                   self.getFunctionSpace(), \
2348                                   self.__matrix_type)
2349              
2350         def __getNewSolutionVector(self):
2351             if self.getNumEquations() ==1 :
2352                    out=escript.Data(0.0,(),self.getFunctionSpace())
2353             else:
2354                    out=escript.Data(0.0,(self.getNumEquations(),),self.getFunctionSpace())
2355             return out
2356    
2357         def __getTransportProblem(self):
2358           if self.__reassemble:
2359                 self.__source=self.__getNewSolutionVector()
2360                 self.__transport_problem.reset()
2361                 self.getDomain().addPDEToTransportProblem(
2362                             self.__transport_problem,
2363                             self.__source,
2364                             self.__M,
2365                             self.__A,
2366                             self.__B,
2367                             self.__C,
2368                             self.__D,
2369                             self.__X,
2370                             self.__Y,
2371                             self.__d,
2372                             self.__y,
2373                             self.__d_contact,
2374                             self.__y_contact)
2375                 self.__transport_problem.insertConstraint(self.__source,self.__q,self.__r)
2376                 self.__reassemble=False
2377           return self.__transport_problem
2378         def setValue(self,M=None, A=None, B=None, C=None, D=None, X=None, Y=None,
2379                      d=None, y=None, d_contact=None, y_contact=None, q=None, r=None):
2380                 if not M==None:
2381                      self.__reassemble=True
2382                      self.__M=M
2383                 if not A==None:
2384                      self.__reassemble=True
2385                      self.__A=A
2386                 if not B==None:
2387                      self.__reassemble=True
2388                      self.__B=B
2389                 if not C==None:
2390                      self.__reassemble=True
2391                      self.__C=C
2392                 if not D==None:
2393                      self.__reassemble=True
2394                      self.__D=D
2395                 if not X==None:
2396                      self.__reassemble=True
2397                      self.__X=X
2398                 if not Y==None:
2399                      self.__reassemble=True
2400                      self.__Y=Y
2401                 if not d==None:
2402                      self.__reassemble=True
2403                      self.__d=d
2404                 if not y==None:
2405                      self.__reassemble=True
2406                      self.__y=y
2407                 if not d_contact==None:
2408                      self.__reassemble=True
2409                      self.__d_contact=d_contact
2410                 if not y_contact==None:
2411                      self.__reassemble=True
2412                      self.__y_contact=y_contact
2413                 if not q==None:
2414                      self.__reassemble=True
2415                      self.__q=q
2416                 if not r==None:
2417                      self.__reassemble=True
2418                      self.__r=r
2419    
2420         def setInitialSolution(self,u):
2421                 if self.__useSUPG:
2422                     self.__u=util.interpolate(u,self.getFunctionSpace())
2423                 else:
2424                     self.__transport_problem.setInitialValue(util.interpolate(u,self.getFunctionSpace()))
2425    
2426         def solve(self,dt,**kwarg):
2427               if self.__useSUPG:
2428                    if self.__reassemble:
2429                        self.__pde.setValue(D=self.__M,d=self.__d,d_contact=self.__d_contact,q=self.__q) # ,r=self.__r)
2430                        self.__reassemble=False
2431                    dt2=self.getSafeTimeStepSize()
2432                    nn=max(math.ceil(dt/self.getSafeTimeStepSize()),1.)
2433                    dt2=dt/nn
2434                    nnn=0
2435                    u=self.__u
2436                    self.trace("number of substeps is %d."%nn)
2437                    while nnn<nn :
2438                        self.__setSUPG(u,u,dt2/2)
2439                        u_half=self.__pde.getSolution(verbose=True)
2440                        self.__setSUPG(u,u_half,dt2)
2441                        u=self.__pde.getSolution(verbose=True)
2442                        nnn+=1
2443                    self.__u=u
2444                    return self.__u
2445               else:
2446                   kwarg["tolerance"]=self.__tolerance
2447                   tp=self.__getTransportProblem()
2448                   return tp.solve(self.__source,dt,kwarg)
2449         def __setSUPG(self,u0,u,dt):
2450                g=util.grad(u)
2451                X=0
2452                Y=self.__M*u0
2453                X=0
2454                self.__pde.setValue(r=u0)
2455                if not self.__A.isEmpty():
2456                   X=X+dt*util.matrixmult(self.__A,g)
2457                if not self.__B.isEmpty():
2458                   X=X+dt*self.__B*u
2459                if not  self.__C.isEmpty():
2460                   Y=Y+dt*util.inner(self.__C,g)
2461                if not self.__D.isEmpty():
2462                   Y=Y+dt*self.__D*u
2463                if not self.__X.isEmpty():
2464                   X=X+dt*self.__X
2465                if not self.__Y.isEmpty():
2466                   Y=Y+dt*self.__Y
2467                self.__pde.setValue(X=X,Y=Y)
2468                if not self.__y.isEmpty():
2469                   self.__pde.setValue(y=dt*self.__y)
2470                if not self.__y_contact.isEmpty():
2471                   self.__pde.setValue(y=dt*self.__y_contact)
2472                self.__pde.setValue(r=u0)

Legend:
Removed from v.1400  
changed lines
  Added in v.1809

  ViewVC Help
Powered by ViewVC 1.1.26