/[escript]/trunk/escript/py_src/heat.py
ViewVC logotype

Diff of /trunk/escript/py_src/heat.py

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

revision 1417 by gross, Mon Feb 25 04:45:48 2008 UTC revision 2548 by jfenwick, Mon Jul 20 06:20:06 2009 UTC
# Line 1  Line 1 
1  # $Id:$  
2  #  ########################################################
 #######################################################  
 #  
 #       Copyright 2008 by University of Queensland  
3  #  #
4  #                http://esscc.uq.edu.au  # Copyright (c) 2003-2009 by University of Queensland
5  #        Primary Business: Queensland, Australia  # Earth Systems Science Computational Center (ESSCC)
6  #  Licensed under the Open Software License version 3.0  # http://www.uq.edu.au/esscc
 #     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__="https://launchpad.net/escript-finley"
21    
22  """  """
23  Some models for heat advection-diffusion  Some models for heat advection-diffusion
# Line 24  Some models for heat advection-diffusion Line 31  Some models for heat advection-diffusion
31  """  """
32    
33  __author__="Lutz Gross, l.gross@uq.edu.au"  __author__="Lutz Gross, l.gross@uq.edu.au"
 __copyright__="""  Copyright (c) 2008 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:$"  
34    
35  # from escript import *  # from escript import *
36  import util  import util
37  from linearPDEs import TransportPDE  from linearPDEs import TransportPDE
38    
39  class TemperatureCartesian(TransportPDE):  class TemperatureCartesian(TransportPDE):
40        """      """
41        solves      Represents and solves the temperature advection-diffusion problem
42    
43        M{rhocp(T_{,t} + v_i T_{,i} - ( k T_{,i})_i = Q}
44    
45        M{k T_{,i}*n_i=surface_flux} and M{T_{,t} = 0} where C{given_T_mask}>0.
46    
47            rhocp(T_{,t} + v_i T_{,i} - ( k T_{,i})_i = Q      If surface_flux is not given 0 is assumed.
48            
49                   k T_{,i}*n_i=surface_flux      Typical usage::
50                
51                   T_{,t} = 0  where given_T_mask>0          sp = TemperatureCartesian(domain)
52            sp.setTolerance(1.e-4)
53        if surface_flux is not give 0 is assumed.          t = 0
54            T = ...
55        typical usage:          sp.setValues(rhocp=...,  v=..., k=..., given_T_mask=...)
56            sp.setInitialTemperature(T)
57              sp=TemperatureCartesian(domain)          while t < t_end:
58              sp.setTolerance(1.e-4)              sp.setValue(Q=...)
59              t=0              T = sp.getTemperature(dt)
60              T=...              t += dt
61              sp.setValues(rhocp = ..,  v=.., k=.., given_T_mask=..)      """
62              sp.setInitialTemperature(T)      def __init__(self,domain,useBackwardEuler=False,**kwargs):
             while t < t_end:  
                 sp.setValue(Q= ...)  
                 T=sp.getTemperature(dt)  
                 t+=dt  
       """  
       def __init__(self,domain,theta=0.5,**kwargs):  
63          """          """
64          initializes tmperature advection-diffuision problem          Initializes the temperature advection-diffusion problem.
65    
66          @param domain: domain of the problem          @param domain: domain of the problem
67          @param theta: method control: theta=1. backward Euler          @param useBackwardEuler: if set the backward Euler scheme is used. Otherwise the Crank-Nicholson scheme is applied. Not that backward Euler scheme will return a safe time step size which is practically infinity as the scheme is unconditional unstable. So other measures need to be applied to control the time step size. The Crank-Nicholson scheme provides a higher accuracy but requires to limit the time step size to be stable.
68                                        theta=0.5: Crank-Nicholson scheme          @type useBackwardEuler: C{bool}
                                       thete=0.: forward Euler (not recommendable)  
69          """          """
70          TransportPDE.__init__(self,domain,num_equations=1,theta=theta,**kwargs)          TransportPDE.__init__(self,domain,numEquations=1,useBackwardEuler=useBackwardEuler,**kwargs)
71          self.setReducedOn()          self.setReducedOrderOn()
72          self.__rhocp=None          self.__rhocp=None
73          self.__v=None          self.__v=None
74        def setInitialTemperature(self,T):  
75        def setInitialTemperature(self,T):
76            """
77            Same as L{setInitialSolution}.
78            """
79          self.setInitialSolution(T)          self.setInitialSolution(T)
80        def setValue(self,rhocp=None,v=None,k=None,Q=None,surface_flux=None,given_T_mask=None):  
81             if rhocp!=None:      def setValue(self,rhocp=None,v=None,k=None,Q=None,surface_flux=None,given_T_mask=None):
82                 self.__rhocp=rhocp          if rhocp!=None:
83             if v!=None:              self.__rhocp=rhocp
84                 self.__v=v          if v!=None:
85             if rhocp!=None:              self.__v=v
86                 super(TemperatureCartesian,self).setValue(M=self.__rhocp)          if rhocp!=None:
87             if (rhocp!=None or v!=None) and self.__rhocp!=None and self.__v!=None:              super(TemperatureCartesian,self).setValue(M=self.__rhocp)
88                 super(TemperatureCartesian,self).setValue(C=-self.__rhocp*self.__v)          if (rhocp!=None or v!=None) and self.__rhocp!=None and self.__v!=None:
89             if k!=None:              super(TemperatureCartesian,self).setValue(C=-self.__rhocp*self.__v)
90                 super(TemperatureCartesian,self).setValue(A=-k*util.kronecker(self.getDomain()))          if k!=None:
91             if Q!=None:              super(TemperatureCartesian,self).setValue(A=-k*util.kronecker(self.getDomain()))
92                 super(TemperatureCartesian,self).setValue(Y=Q)          if Q!=None:
93             if surface_flux!=None:              super(TemperatureCartesian,self).setValue(Y=Q)
94                 super(TemperatureCartesian,self).setValue(y=surface_flux)          if surface_flux!=None:
95             if given_T_mask!=None:              super(TemperatureCartesian,self).setValue(y=surface_flux)
96                 super(TemperatureCartesian,self).setValue(q=given_T_mask)          if given_T_mask!=None:
97                          super(TemperatureCartesian,self).setValue(q=given_T_mask)
98        def getTemperature(self,dt,**kwargs):  
99            return self.solve(dt,**kwargs)      def getTemperature(self,dt,**kwargs):
100            """
101            Same as L{getSolution}.
102            """
103            return self.getSolution(dt,**kwargs)
104    
105    
106    class Tracer(TransportPDE):
107        """
108        Represents and solves the tracer problem
109    
110        M{C_{,t} + v_i C_{,i} - ( k T_{,i})_i) = 0}
111    
112        M{C_{,t} = 0} where C{given_C_mask}>0.
113        M{C_{,i}*n_i=0}
114    
115        Typical usage::
116    
117            sp = Tracer(domain)
118            sp.setTolerance(1.e-4)
119            t = 0
120            T = ...
121            sp.setValues(given_C_mask=...)
122            sp.setInitialTracer(C)
123            while t < t_end:
124                sp.setValue(v=...)
125                dt.getSaveTimeStepSize()
126                C = sp.getTracer(dt)
127                t += dt
128        """
129        def __init__(self,domain,useBackwardEuler=False,**kwargs):
130            """
131            Initializes the Tracer advection problem
132    
133            @param domain: domain of the problem
134            @param useBackwardEuler: if set the backward Euler scheme is used. Otherwise the Crank-Nicholson scheme is applied. Not that backward Euler scheme will return a safe time step size which is practically infinity as the scheme is unconditional unstable. So other measures need to be applied to control the time step size. The Crank-Nicholson scheme provides a higher accuracy but requires to limit the time step size to be stable.
135            @type useBackwardEuler: C{bool}
136            """
137            TransportPDE.__init__(self,domain,numEquations=1,useBackwardEuler=useBackwardEuler,**kwargs)
138            self.setReducedOrderOn()
139            super(Tracer,self).setValue(M=1.)
140    
141        def setInitialTracer(self,C):
142            """
143            Same as L{setInitialSolution}.
144            """
145            self.setInitialSolution(C)
146    
147        def setValue(self,v=None,given_C_mask=None, k=None):
148            if v!=None:
149                super(Tracer,self).setValue(C=-v)
150            if k!=None:
151                super(Tracer,self).setValue(A=-k*util.kronecker(self.getDomain()))
152            if given_C_mask!=None:
153                super(Tracer,self).setValue(q=given_C_mask)
154    
155        def getTracer(self,dt,**kwargs):
156            """
157            Same as L{getSolution}.
158            """
159            return self.getSolution(dt,**kwargs)
160    

Legend:
Removed from v.1417  
changed lines
  Added in v.2548

  ViewVC Help
Powered by ViewVC 1.1.26