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

temp_trunk_copy/escript/py_src/linearPDEs.py revision 1384 by phornby, Fri Jan 11 02:29:38 2008 UTC trunk/escript/py_src/linearPDEs.py revision 1787 by artak, Mon Sep 15 01:36:34 2008 UTC
# Line 34  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 444  class LinearPDE(object): Line 445  class LinearPDE(object):
445     @cvar CR: The conjugate residual method     @cvar CR: The conjugate residual method
446     @cvar CGS: The conjugate gardient square method     @cvar CGS: The conjugate gardient square method
447     @cvar BICGSTAB: The stabilized BiConjugate Gradient method.     @cvar BICGSTAB: The stabilized BiConjugate Gradient method.
448       @cvar TFQMR: Transport Free Quasi Minimal Residual method.
449       @cvar MINRES: Minimum residual method.
450     @cvar SSOR: The symmetric overrealaxtion method     @cvar SSOR: The symmetric overrealaxtion method
451     @cvar ILU0: The incomplete LU factorization preconditioner  with no fill in     @cvar ILU0: The incomplete LU factorization preconditioner  with no fill in
452     @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 492  class LinearPDE(object):
492     AMG= 22     AMG= 22
493     RILU = 23     RILU = 23
494     TRILINOS = 24     TRILINOS = 24
495       NONLINEAR_GMRES = 25
496       TFQMR = 26
497       MINRES = 27
498
499     SMALL_TOLERANCE=1.e-13     SMALL_TOLERANCE=1.e-13
500     __PACKAGE_KEY="package"     __PACKAGE_KEY="package"
# Line 947  class LinearPDE(object): Line 953  class LinearPDE(object):
953         sets a new solver         sets a new solver
954
955         @param solver: sets a new solver method.         @param solver: sets a new solver method.
956         @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}
957         @param preconditioner: sets a new solver method.         @param preconditioner: sets a new solver method.
958         @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}
959         """         """
# Line 977  class LinearPDE(object): Line 983  class LinearPDE(object):
983         elif m[0]==self.ITERATIVE: method= "ITERATIVE"         elif m[0]==self.ITERATIVE: method= "ITERATIVE"
984         elif m[0]==self.CHOLEVSKY: method= "CHOLEVSKY"         elif m[0]==self.CHOLEVSKY: method= "CHOLEVSKY"
985         elif m[0]==self.PCG: method= "PCG"         elif m[0]==self.PCG: method= "PCG"
986           elif m[0]==self.TFQMR: method= "TFQMR"
987           elif m[0]==self.MINRES: method= "MINRES"
988         elif m[0]==self.CR: method= "CR"         elif m[0]==self.CR: method= "CR"
989         elif m[0]==self.CGS: method= "CGS"         elif m[0]==self.CGS: method= "CGS"
990         elif m[0]==self.BICGSTAB: method= "BICGSTAB"         elif m[0]==self.BICGSTAB: method= "BICGSTAB"
# Line 1746  class LinearPDE(object): Line 1754  class LinearPDE(object):
1754                      d_reduced_times_e=escript.Data()                      d_reduced_times_e=escript.Data()
1755
1756                   self.__operator=self.__getNewRightHandSide()                   self.__operator=self.__getNewRightHandSide()
1760                   else:                   else:
# Line 2190  class LameEquation(LinearPDE): Line 2198  class LameEquation(LinearPDE):
2198       else:       else:
2199          raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name          raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name
2200
2201    def LinearSinglePDE(domain,debug=False):
2202       """
2203       defines a single linear PDEs
2204
2205       @param domain: domain of the PDE
2206       @type domain: L{Domain<escript.Domain>}
2207       @param debug: if True debug informations are printed.
2208       @rtype: L{LinearPDE}
2209       """
2210       return LinearPDE(domain,numEquations=1,numSolutions=1,debug=debug)
2211
2212    def LinearPDESystem(domain,debug=False):
2213       """
2214       defines a system of linear PDEs
2215
2216       @param domain: domain of the PDE
2217       @type domain: L{Domain<escript.Domain>}
2218       @param debug: if True debug informations are printed.
2219       @rtype: L{LinearPDE}
2220       """
2221       return LinearPDE(domain,numEquations=domain.getDim(),numSolutions=domain.getDim(),debug=debug)
2222
2223    class TransportPDE(object):
2224         """
2225         Warning: This is still a very experimental. The class is still changing!
2226
2227         Mu_{,t} =-(A_{ij}u_{,j})_j-(B_{j}u)_{,j} + C_{j} u_{,j} + Y_i + X_{i,i}
2228
2229         u=r where q>0
2230
2231         all coefficients are constant over time.
2232
2233         typical usage:
2234
2235             p=TransportPDE(dom)
2236             p.setValue(M=Scalar(1.,Function(dom),C=Scalar(1.,Function(dom)*[-1.,0.])
2237             p.setInitialSolution(u=exp(-length(dom.getX()-[0.1,0.1])**2)
2238             t=0
2239             dt=0.1
2240             while (t<1.):
2241                  u=p.solve(dt)
2242
2243         """
2244         def __init__(self,domain,num_equations=1,theta=0.5,useSUPG=False,trace=True):
2245            self.__domain=domain
2246            self.__num_equations=num_equations
2247            self.__useSUPG=useSUPG
2248            self.__trace=trace
2249            self.__theta=theta
2250            self.__matrix_type=0
2251            self.__reduced=True
2252            self.__reassemble=True
2253            if self.__useSUPG:
2254               self.__pde=LinearPDE(domain,numEquations=num_equations,numSolutions=num_equations,debug=trace)
2255               self.__pde.setSymmetryOn()
2256               self.__pde.setReducedOrderOn()
2257            else:
2258               self.__transport_problem=self.__getNewTransportProblem()
2259            self.setTolerance()
2260            self.__M=escript.Data()
2261            self.__A=escript.Data()
2262            self.__B=escript.Data()
2263            self.__C=escript.Data()
2264            self.__D=escript.Data()
2265            self.__X=escript.Data()
2266            self.__Y=escript.Data()
2267            self.__d=escript.Data()
2268            self.__y=escript.Data()
2269            self.__d_contact=escript.Data()
2270            self.__y_contact=escript.Data()
2271            self.__r=escript.Data()
2272            self.__q=escript.Data()
2273
2274         def trace(self,text):
2275                 if self.__trace: print text
2276         def getSafeTimeStepSize(self):
2277            if self.__useSUPG:
2278                if self.__reassemble:
2279                   h=self.__domain.getSize()
2280                   dt=None
2281                   if not self.__A.isEmpty():
2282                      dt2=util.inf(h**2*self.__M/util.length(self.__A))
2283                      if dt == None:
2284                         dt = dt2
2285                      else:
2286                         dt=1./(1./dt+1./dt2)
2287                   if not self.__B.isEmpty():
2288                      dt2=util.inf(h*self.__M/util.length(self.__B))
2289                      if dt == None:
2290                         dt = dt2
2291                      else:
2292                         dt=1./(1./dt+1./dt2)
2293                   if not  self.__C.isEmpty():
2294                      dt2=util.inf(h*self.__M/util.length(self.__C))
2295                      if dt == None:
2296                         dt = dt2
2297                      else:
2298                         dt=1./(1./dt+1./dt2)
2299                   if not self.__D.isEmpty():
2300                      dt2=util.inf(self.__M/util.length(self.__D))
2301                      if dt == None:
2302                         dt = dt2
2303                      else:
2304                         dt=1./(1./dt+1./dt2)
2305                   self.__dt = dt/2
2306                return self.__dt
2307            else:
2308                return self.__getTransportProblem().getSafeTimeStepSize()
2309         def getDomain(self):
2310            return self.__domain
2311         def getTheta(self):
2312            return self.__theta
2313         def getNumEquations(self):
2314            return self.__num_equations
2315         def setReducedOn(self):
2316              if not self.reduced():
2317                  if self.__useSUPG:
2318                     self.__pde.setReducedOrderOn()
2319                  else:
2320                     self.__transport_problem=self.__getNewTransportProblem()
2321              self.__reduced=True
2322         def setReducedOff(self):
2323              if self.reduced():
2324                  if self.__useSUPG:
2325                     self.__pde.setReducedOrderOff()
2326                  else:
2327                     self.__transport_problem=self.__getNewTransportProblem()
2328              self.__reduced=False
2329         def reduced(self):
2330             return self.__reduced
2331         def getFunctionSpace(self):
2332            if self.reduced():
2333               return escript.ReducedSolution(self.getDomain())
2334            else:
2335               return escript.Solution(self.getDomain())
2336
2337         def setTolerance(self,tol=1.e-8):
2338            self.__tolerance=tol
2339            if self.__useSUPG:
2340                  self.__pde.setTolerance(self.__tolerance)
2341
2342         def __getNewTransportProblem(self):
2343           """
2344           returns an instance of a new operator
2345           """
2346           self.trace("New Transport problem is allocated.")
2347           return self.getDomain().newTransportProblem( \
2348                                   self.getTheta(),
2349                                   self.getNumEquations(), \
2350                                   self.getFunctionSpace(), \
2351                                   self.__matrix_type)
2352
2353         def __getNewSolutionVector(self):
2354             if self.getNumEquations() ==1 :
2355                    out=escript.Data(0.0,(),self.getFunctionSpace())
2356             else:
2357                    out=escript.Data(0.0,(self.getNumEquations(),),self.getFunctionSpace())
2358             return out
2359
2360         def __getTransportProblem(self):
2361           if self.__reassemble:
2362                 self.__source=self.__getNewSolutionVector()
2363                 self.__transport_problem.reset()
2365                             self.__transport_problem,
2366                             self.__source,
2367                             self.__M,
2368                             self.__A,
2369                             self.__B,
2370                             self.__C,
2371                             self.__D,
2372                             self.__X,
2373                             self.__Y,
2374                             self.__d,
2375                             self.__y,
2376                             self.__d_contact,
2377                             self.__y_contact)
2378                 self.__transport_problem.insertConstraint(self.__source,self.__q,self.__r)
2379                 self.__reassemble=False
2380           return self.__transport_problem
2381         def setValue(self,M=None, A=None, B=None, C=None, D=None, X=None, Y=None,
2382                      d=None, y=None, d_contact=None, y_contact=None, q=None, r=None):
2383                 if not M==None:
2384                      self.__reassemble=True
2385                      self.__M=M
2386                 if not A==None:
2387                      self.__reassemble=True
2388                      self.__A=A
2389                 if not B==None:
2390                      self.__reassemble=True
2391                      self.__B=B
2392                 if not C==None:
2393                      self.__reassemble=True
2394                      self.__C=C
2395                 if not D==None:
2396                      self.__reassemble=True
2397                      self.__D=D
2398                 if not X==None:
2399                      self.__reassemble=True
2400                      self.__X=X
2401                 if not Y==None:
2402                      self.__reassemble=True
2403                      self.__Y=Y
2404                 if not d==None:
2405                      self.__reassemble=True
2406                      self.__d=d
2407                 if not y==None:
2408                      self.__reassemble=True
2409                      self.__y=y
2410                 if not d_contact==None:
2411                      self.__reassemble=True
2412                      self.__d_contact=d_contact
2413                 if not y_contact==None:
2414                      self.__reassemble=True
2415                      self.__y_contact=y_contact
2416                 if not q==None:
2417                      self.__reassemble=True
2418                      self.__q=q
2419                 if not r==None:
2420                      self.__reassemble=True
2421                      self.__r=r
2422
2423         def setInitialSolution(self,u):
2424                 if self.__useSUPG:
2425                     self.__u=util.interpolate(u,self.getFunctionSpace())
2426                 else:
2427                     self.__transport_problem.setInitialValue(util.interpolate(u,self.getFunctionSpace()))
2428
2429         def solve(self,dt,**kwarg):
2430               if self.__useSUPG:
2431                    if self.__reassemble:
2432                        self.__pde.setValue(D=self.__M,d=self.__d,d_contact=self.__d_contact,q=self.__q) # ,r=self.__r)
2433                        self.__reassemble=False
2434                    dt2=self.getSafeTimeStepSize()
2435                    nn=max(math.ceil(dt/self.getSafeTimeStepSize()),1.)
2436                    dt2=dt/nn
2437                    nnn=0
2438                    u=self.__u
2439                    self.trace("number of substeps is %d."%nn)
2440                    while nnn<nn :
2441                        self.__setSUPG(u,u,dt2/2)
2442                        u_half=self.__pde.getSolution(verbose=True)
2443                        self.__setSUPG(u,u_half,dt2)
2444                        u=self.__pde.getSolution(verbose=True)
2445                        nnn+=1
2446                    self.__u=u
2447                    return self.__u
2448               else:
2449                   kwarg["tolerance"]=self.__tolerance
2450                   tp=self.__getTransportProblem()
2451                   return tp.solve(self.__source,dt,kwarg)
2452         def __setSUPG(self,u0,u,dt):
2454                X=0
2455                Y=self.__M*u0
2456                X=0
2457                self.__pde.setValue(r=u0)
2458                if not self.__A.isEmpty():
2459                   X=X+dt*util.matrixmult(self.__A,g)
2460                if not self.__B.isEmpty():
2461                   X=X+dt*self.__B*u
2462                if not  self.__C.isEmpty():
2463                   Y=Y+dt*util.inner(self.__C,g)
2464                if not self.__D.isEmpty():
2465                   Y=Y+dt*self.__D*u
2466                if not self.__X.isEmpty():
2467                   X=X+dt*self.__X
2468                if not self.__Y.isEmpty():
2469                   Y=Y+dt*self.__Y
2470                self.__pde.setValue(X=X,Y=Y)
2471                if not self.__y.isEmpty():
2472                   self.__pde.setValue(y=dt*self.__y)
2473                if not self.__y_contact.isEmpty():
2474                   self.__pde.setValue(y=dt*self.__y_contact)
2475                self.__pde.setValue(r=u0)

Legend:
 Removed from v.1384 changed lines Added in v.1787