/[escript]/trunk/modellib/py_src/temperature.py
ViewVC logotype

Diff of /trunk/modellib/py_src/temperature.py

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

revision 127 by jgs, Fri Jul 22 05:11:29 2005 UTC revision 150 by jgs, Thu Sep 15 03:44:45 2005 UTC
# Line 1  Line 1 
1  # $Id$  # $Id$
2    
   
3  from esys.escript import *  from esys.escript import *
4  from esys.modelframe import Model,IterationDivergenceError  from esys.escript.modelframe import Model,IterationDivergenceError
5  from esys.linearPDEs import AdvectivePDE,LinearPDE  from esys.escript.linearPDEs import AdvectivePDE,LinearPDE
6    import numarray
7    
 class TemperatureDiffusion(Model):  
        """ """  
8    
9         def __init__(self,debug=False):  class TemperatureAdvection(Model):
10             Model.__init__(self,debug=debug)         """
            self.declareParameter(domain=None, \  
                                  tend=0., \  
                                  dt=0.1, \  
                                  temperature=1., \  
                                  density=1., \  
                                  c_p=1., \  
                                  thermal_permabilty=1., \  
                                  reference_temperature=0., \  
                                  radiation_coefficient=0., \  
                                  thermal_source=0., \  
                                  location_fixed_temperature=Data(), \  
                                  iterate=True, \  
                                  tol=1.e-8, \  
                                  implicit=True)  
            self.iter=0  
   
        def doInitialization(self,t):  
            self.tn=t  
            self.pde=LinearPDE(self.domain)  
            self.pde.setSymmetryOn()  
   
        def doIterationInitialization(self,dt):  
             self.iter=0  
             self.T_last=self.temperature  
             self.diff=1.e400  
   
        def doIterationStep(self,dt):  
           T=self.temperature  
           diff=self.diff  
           dim=self.pde.getDim()  
           self.iter+=1  
           rhocp=self.density*self.c_p  
           self.pde.setValue(A=self.thermal_permabilty*kronecker(dim), \  
                               D=rhocp/dt, \  
                               Y=self.thermal_source+rhocp/dt*self.T_last, \  
                               d=self.radiation_coefficient, \  
                               y=self.radiation_coefficient*self.reference_temperature, \  
                               q=self.location_fixed_temperature, \  
                               r=self.T_last)  
           if isinstance(self,TemperatureAdvection): self.pde.setValue(C=self.velocity[:dim]*rhocp)  
           self.pde.setTolerance(self.tol*1.e-2)  
           self.temperature=self.pde.getSolution()  
           self.diff=Lsup(T-self.temperature)  
           if diff<=self.diff:  
               raise IterationDivergenceError,"no improvement in the temperature iteration"  
   
        def terminate(self):  
           if not self.implicit:  
               return True  
           elif self.iter<1:  
               return False  
           else:  
              return self.diff<self.tol*Lsup(self.temperature)  
11    
12         def doIterationFinalization(self,dt):         The conservation of internal heat energy is given by
           self.tn+=dt  
13    
14         def getSafeTimeStepSize(self,dt):         M{S{rho} c_p ( dT/dt+v[j]*grad(T)[j])-grad(\kappa grad(T)_{,i}=Q}
            return self.dt  
15    
16         def finalize(self):         M{n_i\kappa T_{,i}=0}
             return self.tn>=self.tend  
17    
18  class TemperatureAdvection(Model):         it is assummed that M{\rho c_p} is constant in time.
19         """ """  
20           solved by Taylor Galerkin method
21    
22           """
23         def __init__(self,debug=False):         def __init__(self,debug=False):
24             Model.__init__(self,debug=debug)             Model.__init__(self,debug=debug)
25             self.declareParameter(velocity=numarray.zeros([3]))             self.declareParameter(domain=None, \
26                                     temperature=1., \
27         def doInitialization(self,t):                                   velocity=numarray.zeros([3]),
28             self.tn=t                                   density=1., \
29             self.pde=AdvectivePDE(self.domain)                                   heat_capacity=1., \
30                                     thermal_permabilty=1., \
31                                     # reference_temperature=0., \
32                                     # radiation_coefficient=0., \
33                                     thermal_source=0., \
34                                     fixed_temperature=0.,
35                                     location_fixed_temperature=Data(),
36                                     safety_factor=0.1)
37    
38           def doInitialization(self):
39               self.__pde=LinearPDE(self.domain)
40               self.__pde.setSymmetryOn()
41               # self.__pde.setReducedOrderOn()
42               self.__pde.setSolverMethod(self.__pde.LUMPING)
43               self.__pde.setValue(D=self.heat_capacity*self.density)
44    
45         def getSafeTimeStepSize(self,dt):         def getSafeTimeStepSize(self,dt):
46             v=Lsup(self.velocity)             """
47             if v>0.:             returns new step size
48                 return min(self.dt,Lsup(self.pde.getDomain().getSize())/v)             """
49             else:             h=self.domain.getSize()
50                 return self.dt             return self.safety_factor*inf(h**2/(h*abs(self.heat_capacity*self.density)*length(self.velocity)+self.thermal_permabilty))
51    
52                 def G(self,T,alpha):
53               """
54  if __name__=="__main__":             tangential operator for taylor galerikin
55     from esys.modelframe import Link,Simulation,ExplicitSimulation             """
56     from esys.visualization import WriteVTK             g=grad(T)
57     from esys.materials import MaterialTable             self.__pde.setValue(X=-self.thermal_permabilty*g, \
58     from esys.geometry import RectangularDomain,ScalarConstrainer                                 Y=self.thermal_source-self.__rhocp*inner(self.velocity,g), \
59     from esys.input import InterpolatedTimeProfile,GausseanProfile                                 r=(self.__fixed_T-self.temperature)*alpha,\
60                                   q=self.location_fixed_temperature)
61     dom=RectangularDomain()             return self.__pde.getSolution()
62     constraints=ScalarConstrainer()            
63     constraints.domain=Link(dom)  
64     constraints.top=1         def doStepPostprocessing(self,dt):
65     constraints.bottom=1             """
66                 perform taylor galerkin step
67     mt=MaterialTable()             """
68                 T=self.temperature
69     pf=InterpolatedTimeProfile()         self.__rhocp=self.heat_capacity*self.density
70     pf.t=[0.,0.25,0.5,0.75]             self.__fixed_T=self.fixed_temperature
71     pf.values=[0.,1.,1.,0.]             self.temperature=dt*self.G(dt/2*self.G(T,1./dt)+T,1./dt)+T
72                 self.trace("Temperature range is %e %e"%(inf(self.temperature),sup(self.temperature)))
    q=GausseanProfile()  
    q.domain=Link(dom)  
    q.width=0.05  
    q.x_c=numarray.array([0.5,0.5,0.5])  
    q.r=0.01  
    q.A=Link(pf,"out")  
     
    tt=TemperatureDiffusion()  
    tt.domain=Link(dom)  
    tt.tend=1.  
    tt.dt=0.1  
    tt.temperature=0.  
    tt.density=Link(mt)  
    tt.c_p=Link(mt)  
    tt.thermal_permabilty=Link(mt)  
    tt.reference_temperature=0.  
    tt.radiation_coefficient=Link(mt)  
    tt.thermal_source=Link(q,"out")  
    tt.location_fixed_temperature=Link(constraints,"location_of_constraint")  
    tt.implicit=True  
     
    vis=WriteVTK()  
    vis.scalar=Link(tt,"temperature")  
     
    s=ExplicitSimulation([dom,constraints,pf,q,Simulation([mt,tt],debug=True),vis],debug=True)  
    # s=Simulation([dom,constraints,pf,q,Simulation([mt,tt]),vis])  
    s.writeXML()  
    s.run()  

Legend:
Removed from v.127  
changed lines
  Added in v.150

  ViewVC Help
Powered by ViewVC 1.1.26