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

Diff of /trunk/esys2/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 147 by jgs, Fri Aug 12 01:45:47 2005 UTC
# Line 1  Line 1 
1  # $Id$  # $Id$
2    
3    from escript.escript import *
4    from escript.modelframe import Model,IterationDivergenceError
5    from escript.linearPDEs import AdvectivePDE,LinearPDE
6    import numarray
7    
 from esys.escript import *  
 from esys.modelframe import Model,IterationDivergenceError  
 from esys.linearPDEs import AdvectivePDE,LinearPDE  
8    
9  class TemperatureDiffusion(Model):  class TemperatureAdvection(Model):
10         """ """         """
11    
12             The conservation of internal heat energy is given by
13    
14             \f[
15                 \rho c_p ( T_{,t}+v_{j}T_{,j} )-(\kappa T_{,i})_{,i}=Q,
16             \f]
17             \f[
18                     n_i\kappa T_{,i}=0
19             \f]
20    
21              it is assummed that \f[ \rho c_p \f] is constant in time.
22    
23              solved by Taylor Galerkin method
24    
25           """
26         def __init__(self,debug=False):         def __init__(self,debug=False):
27             Model.__init__(self,debug=debug)             Model.__init__(self,debug=debug)
28             self.declareParameter(domain=None, \             self.declareParameter(domain=None, \
                                  tend=0., \  
                                  dt=0.1, \  
29                                   temperature=1., \                                   temperature=1., \
30                                     velocity=numarray.zeros([3]),
31                                   density=1., \                                   density=1., \
32                                   c_p=1., \                                   heat_capacity=1., \
33                                   thermal_permabilty=1., \                                   thermal_permabilty=1., \
34                                   reference_temperature=0., \                                   # reference_temperature=0., \
35                                   radiation_coefficient=0., \                                   # radiation_coefficient=0., \
36                                   thermal_source=0., \                                   thermal_source=0., \
37                                   location_fixed_temperature=Data(), \                                   fixed_temperature=0.,
38                                   iterate=True, \                                   location_fixed_temperature=Data(),
39                                   tol=1.e-8, \                                   safety_factor=0.1)
40                                   implicit=True)  
41             self.iter=0         def doInitialization(self):
42               self.__pde=LinearPDE(self.domain)
43         def doInitialization(self,t):             self.__pde.setSymmetryOn()
44             self.tn=t             # self.__pde.setReducedOrderOn()
45             self.pde=LinearPDE(self.domain)             # self.__pde.setLumpingOn()
46             self.pde.setSymmetryOn()             self.__pde.setValue(D=self.heat_capacity*self.density)
   
        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)  
   
        def doIterationFinalization(self,dt):  
           self.tn+=dt  
   
        def getSafeTimeStepSize(self,dt):  
            return self.dt  
   
        def finalize(self):  
             return self.tn>=self.tend  
   
 class TemperatureAdvection(Model):  
        """ """  
   
        def __init__(self,debug=False):  
            Model.__init__(self,debug=debug)  
            self.declareParameter(velocity=numarray.zeros([3]))  
   
        def doInitialization(self,t):  
            self.tn=t  
            self.pde=AdvectivePDE(self.domain)  
47    
48         def getSafeTimeStepSize(self,dt):         def getSafeTimeStepSize(self,dt):
49             v=Lsup(self.velocity)             """returns new step size"""
50             if v>0.:             h=self.domain.getSize()
51                 return min(self.dt,Lsup(self.pde.getDomain().getSize())/v)             return self.safety_factor*inf(h**2/(h*abs(self.heat_capacity*self.density)*length(self.velocity)+self.thermal_permabilty))
52             else:  
53                 return self.dt         def G(self,T,alpha):
54               """tangential operator for taylor galerikin"""
55                     g=grad(T)
56               self.__pde.setValue(X=-self.thermal_permabilty*g, \
57  if __name__=="__main__":                                 Y=self.thermal_source-self.__rhocp*inner(self.velocity,g), \
58     from esys.modelframe import Link,Simulation,ExplicitSimulation                                 r=(self.__fixed_T-self.temperature)*alpha,\
59     from esys.visualization import WriteVTK                                 q=self.location_fixed_temperature)
60     from esys.materials import MaterialTable             return self.__pde.getSolution()
61     from esys.geometry import RectangularDomain,ScalarConstrainer            
62     from esys.input import InterpolatedTimeProfile,GausseanProfile  
63           def doStepPostprocessing(self,dt):
64     dom=RectangularDomain()             """perform taylor galerkin step"""
65     constraints=ScalarConstrainer()             T=self.temperature
66     constraints.domain=Link(dom)         self.__rhocp=self.heat_capacity*self.density
67     constraints.top=1             self.__fixed_T=self.fixed_temperature
68     constraints.bottom=1             self.temperature=dt*self.G(dt/2*self.G(T,1./dt)+T,1./dt)+T
69                 self.trace("Temperature range is %e %e"%(inf(self.temperature),sup(self.temperature)))
    mt=MaterialTable()  
     
    pf=InterpolatedTimeProfile()  
    pf.t=[0.,0.25,0.5,0.75]  
    pf.values=[0.,1.,1.,0.]  
     
    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.147
  ViewVC Help
Powered by ViewVC 1.1.26