/[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 1410 by gross, Thu Feb 7 04:24:00 2008 UTC revision 1417 by gross, Mon Feb 25 04:45:48 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 2218  class TransportPDE(object): Line 2219  class TransportPDE(object):
2219    
2220       Mu_{,t} =-(A_{ij}u_{,j})_j-(B_{j}u)_{,j} + C_{j} u_{,j} + Y_i + X_{i,i}       Mu_{,t} =-(A_{ij}u_{,j})_j-(B_{j}u)_{,j} + C_{j} u_{,j} + Y_i + X_{i,i}
2221            
2222         u=r where q>0
2223        
2224       all coefficients are constant over time.       all coefficients are constant over time.
2225    
2226       typical usage:       typical usage:
# Line 2231  class TransportPDE(object): Line 2234  class TransportPDE(object):
2234                u=p.solve(dt)                u=p.solve(dt)
2235    
2236       """       """
2237       def __init__(self,domain,num_equations=1,theta=0.5,trace=True):       def __init__(self,domain,num_equations=1,theta=0.5,useSUPG=False,trace=True):
2238          self.__domain=domain          self.__domain=domain
2239          self.__num_equations=num_equations          self.__num_equations=num_equations
2240          self.__theta=theta          self.__useSUPG=useSUPG
2241          self.__trace=trace          self.__trace=trace
2242            self.__theta=theta
2243          self.__matrix_type=0          self.__matrix_type=0
2244            self.__reduced=False
2245            self.__reassemble=True
2246            if self.__useSUPG:
2247               self.__pde=LinearPDE(domain,numEquations=num_equations,numSolutions=num_equations,debug=trace)
2248               self.__pde.setSymmetryOn()
2249            else:
2250               self.__transport_problem=self.__getNewTransportProblem()
2251          self.setTolerance()          self.setTolerance()
2252          self.__transport_problem=self.__getNewTransportProblem()          self.setReducedOn()
2253            self.__M=escript.Data()
2254            self.__A=escript.Data()
2255            self.__B=escript.Data()
2256            self.__C=escript.Data()
2257            self.__D=escript.Data()
2258            self.__X=escript.Data()
2259            self.__Y=escript.Data()
2260            self.__d=escript.Data()
2261            self.__y=escript.Data()
2262            self.__d_contact=escript.Data()
2263            self.__y_contact=escript.Data()
2264            self.__r=escript.Data()
2265            self.__q=escript.Data()
2266    
2267       def trace(self,text):       def trace(self,text):
2268               if self.__trace: print text               if self.__trace: print text
2269       def getSafeTimeStepSize(self):       def getSafeTimeStepSize(self):
2270          return self.__transport_problem.getSafeTimeStepSize()          if self.__useSUPG:
2271                if self.__reassemble:
2272                   h=self.__domain.getSize()
2273                   dt=None
2274                   if not self.__A.isEmpty():
2275                      dt2=util.inf(h**2*self.__M/util.length(self.__A))
2276                      if dt == None:
2277                         dt = dt2
2278                      else:
2279                         dt=1./(1./dt+1./dt2)
2280                   if not self.__B.isEmpty():
2281                      dt2=util.inf(h*self.__M/util.length(self.__B))
2282                      if dt == None:
2283                         dt = dt2
2284                      else:
2285                         dt=1./(1./dt+1./dt2)
2286                   if not  self.__C.isEmpty():
2287                      dt2=util.inf(h*self.__M/util.length(self.__C))
2288                      if dt == None:
2289                         dt = dt2
2290                      else:
2291                         dt=1./(1./dt+1./dt2)
2292                   if not self.__D.isEmpty():
2293                      dt2=util.inf(self.__M/util.length(self.__D))
2294                      if dt == None:
2295                         dt = dt2
2296                      else:
2297                         dt=1./(1./dt+1./dt2)
2298                   self.__dt = dt/2
2299                return self.__dt
2300            else:
2301                return self.__getTransportProblem().getSafeTimeStepSize()
2302       def getDomain(self):       def getDomain(self):
2303          return self.__domain          return self.__domain
2304       def getTheta(self):       def getTheta(self):
2305          return self.__theta          return self.__theta
2306       def getNumEquations(self):       def getNumEquations(self):
2307          return self.__num_equations          return self.__num_equations
2308         def setReducedOn(self):
2309              if not self.reduced():
2310                  if self.__useSUPG:
2311                     self.__pde.setReducedOrderOn()
2312                  else:
2313                     self.__transport_problem=self.__getNewTransportProblem()
2314              self.__reduced=True
2315         def setReducedOff(self):
2316              if self.reduced():
2317                  if self.__useSUPG:
2318                     self.__pde.setReducedOrderOff()
2319                  else:
2320                     self.__transport_problem=self.__getNewTransportProblem()
2321              self.__reduced=False
2322       def reduced(self):       def reduced(self):
2323               return False           return self.__reduced
2324       def getFunctionSpace(self):       def getFunctionSpace(self):
2325          if self.reduced():          if self.reduced():
2326             return escript.ReducedSolution(self.getDomain())             return escript.ReducedSolution(self.getDomain())
2327          else:          else:
2328             return escript.Solution(self.getDomain())             return escript.Solution(self.getDomain())
2329    
2330       def setTolerance(self,tol=1.e-8):       def setTolerance(self,tol=1.e-8):
2331          self.__tolerance=tol          self.__tolerance=tol
2332            if self.__useSUPG:
2333                  self.__pde.setTolerance(self.__tolerance)
2334    
2335       def __getNewTransportProblem(self):       def __getNewTransportProblem(self):
2336         """         """
# Line 2270  class TransportPDE(object): Line 2342  class TransportPDE(object):
2342                                 self.getNumEquations(), \                                 self.getNumEquations(), \
2343                                 self.getFunctionSpace(), \                                 self.getFunctionSpace(), \
2344                                 self.__matrix_type)                                 self.__matrix_type)
2345              
2346       def setValue(self,M=escript.Data(),A=escript.Data(),B=escript.Data(),C=escript.Data(),D=escript.Data(),X=escript.Data(),Y=escript.Data(),       def __getNewSolutionVector(self):
                   d=escript.Data(),y=escript.Data(),d_contact=escript.Data(),y_contact=escript.Data()):  
2347           if self.getNumEquations() ==1 :           if self.getNumEquations() ==1 :
2348                  self.__source=escript.Data(0.0,(),self.getFunctionSpace())                  out=escript.Data(0.0,(),self.getFunctionSpace())
2349           else:           else:
2350                   self.__source=escript.Data(0.0,(self.getNumEquations(),),self.getFunctionSpace())                  out=escript.Data(0.0,(self.getNumEquations(),),self.getFunctionSpace())
2351           self.getDomain().addPDEToTransportProblem(           return out
                      self.__transport_problem,  
                      self.__source,  
                      M,A,B,C,D,X,Y,d,y,d_contact,y_contact)  
2352    
2353       def setInitialSolution(self,u):       def __getTransportProblem(self):
2354               self.__transport_problem.setInitialValue(util.interpolate(u,self.getFunctionSpace()))         if self.__reassemble:
2355                 self.__source=self.__getNewSolutionVector()
2356                 self.__transport_problem.reset()
2357                 self.getDomain().addPDEToTransportProblem(
2358                             self.__transport_problem,
2359                             self.__source,
2360                             self.__M,
2361                             self.__A,
2362                             self.__B,
2363                             self.__C,
2364                             self.__D,
2365                             self.__X,
2366                             self.__Y,
2367                             self.__d,
2368                             self.__y,
2369                             self.__d_contact,
2370                             self.__y_contact)
2371                 self.__transport_problem.insertConstraint(self.__source,self.__q,self.__r)
2372                 self.__reassemble=False
2373           return self.__transport_problem
2374         def setValue(self,M=None, A=None, B=None, C=None, D=None, X=None, Y=None,
2375                      d=None, y=None, d_contact=None, y_contact=None, q=None, r=None):
2376                 if not M==None:
2377                      self.__reassemble=True
2378                      self.__M=M
2379                 if not A==None:
2380                      self.__reassemble=True
2381                      self.__A=A
2382                 if not B==None:
2383                      self.__reassemble=True
2384                      self.__B=B
2385                 if not C==None:
2386                      self.__reassemble=True
2387                      self.__C=C
2388                 if not D==None:
2389                      self.__reassemble=True
2390                      self.__D=D
2391                 if not X==None:
2392                      self.__reassemble=True
2393                      self.__X=X
2394                 if not Y==None:
2395                      self.__reassemble=True
2396                      self.__Y=Y
2397                 if not d==None:
2398                      self.__reassemble=True
2399                      self.__d=d
2400                 if not y==None:
2401                      self.__reassemble=True
2402                      self.__y=y
2403                 if not d_contact==None:
2404                      self.__reassemble=True
2405                      self.__d_contact=d_contact
2406                 if not y_contact==None:
2407                      self.__reassemble=True
2408                      self.__y_contact=y_contact
2409                 if not q==None:
2410                      self.__reassemble=True
2411                      self.__q=q
2412                 if not r==None:
2413                      self.__reassemble=True
2414                      self.__r=r
2415    
2416       def solve(self,dt):       def setInitialSolution(self,u):
2417             return self.__transport_problem.solve(self.__source,dt,{"verbose" : True , "tolerance" : self.__tolerance})               if self.__useSUPG:
2418                     self.__u=u
2419                 else:
2420                     self.__transport_problem.setInitialValue(util.interpolate(u,self.getFunctionSpace()))
2421    
2422         def solve(self,dt,**kwarg):
2423               if self.__useSUPG:
2424                    if self.__reassemble:
2425                        self.__pde.setValue(D=self.__M,d=self.__d,d_contact=self.__d_contact,q=self.__q,r=self.__r)
2426                        self.__reassemble=False
2427                    dt2=self.getSafeTimeStepSize()
2428                    nn=max(math.ceil(dt/self.getSafeTimeStepSize()),1.)
2429                    dt2=dt/nn
2430                    nnn=0
2431                    u=self.__u
2432                    self.trace("number of substeps is %d."%nn)
2433                    while nnn<nn :
2434                        self.__setSUPG(u,u,dt2/2)
2435                        u_half=self.__pde.getSolution(verbose=True)
2436                        self.__setSUPG(u,u_half,dt2)
2437                        u=self.__pde.getSolution(verbose=True)
2438                        nnn+=1
2439                    self.__u=u
2440                    return self.__u
2441               else:
2442                   kwarg["tolerance"]=self.__tolerance
2443                   tp=self.__getTransportProblem()
2444                   return tp.solve(self.__source,dt,kwarg)
2445         def __setSUPG(self,u0,u,dt):
2446                g=util.grad(u)
2447                X=0
2448                Y=self.__M*u0
2449                X=0
2450                if not self.__A.isEmpty():
2451                   X=X+dt*util.matrixmult(self.__A,g)
2452                if not self.__B.isEmpty():
2453                   X=X+dt*self.__B*u
2454                if not  self.__C.isEmpty():
2455                   Y=Y+dt*util.inner(self.__C,g)
2456                if not self.__D.isEmpty():
2457                   Y=Y+dt*self.__D*u
2458                if not self.__X.isEmpty():
2459                   X=X+dt*self.__X
2460                if not self.__Y.isEmpty():
2461                   Y=Y+dt*self.__Y
2462                self.__pde.setValue(X=X,Y=Y)
2463                if not self.__y.isEmpty():
2464                   self.__pde.setValue(y=dt*self.__y)
2465                if not self.__y_contact.isEmpty():
2466                   self.__pde.setValue(y=dt*self.__y_contact)
2467      
Legend:
Removed from v.1410  
changed lines
  Added in v.1417
  ViewVC Help
Powered by ViewVC 1.1.26