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

Contents of /trunk/escriptcore/py_src/heat.py

Parent Directory Parent Directory | Revision Log Revision Log


Revision 5706 - (show annotations)
Mon Jun 29 03:41:36 2015 UTC (4 years, 2 months ago) by sshaw
File MIME type: text/x-python
File size: 6160 byte(s)
all python files now force use of python3 prints and division syntax to stop sneaky errors appearing in py3 environs
1
2 ##############################################################################
3 #
4 # Copyright (c) 2003-2015 by The University of Queensland
5 # http://www.uq.edu.au
6 #
7 # Primary Business: Queensland, Australia
8 # Licensed under the Open Software License version 3.0
9 # http://www.opensource.org/licenses/osl-3.0.php
10 #
11 # Development until 2012 by Earth Systems Science Computational Center (ESSCC)
12 # Development 2012-2013 by School of Earth Sciences
13 # Development from 2014 by Centre for Geoscience Computing (GeoComp)
14 #
15 ##############################################################################
16
17 from __future__ import print_function, division
18
19 __copyright__="""Copyright (c) 2003-2015 by The University of Queensland
20 http://www.uq.edu.au
21 Primary Business: Queensland, Australia"""
22 __license__="""Licensed under the Open Software License version 3.0
23 http://www.opensource.org/licenses/osl-3.0.php"""
24 __url__="https://launchpad.net/escript-finley"
25
26 """
27 Some models for heat advection-diffusion
28
29 :var __author__: name of author
30 :var __copyright__: copyrights
31 :var __license__: licence agreement
32 :var __url__: url entry point on documentation
33 :var __version__: version
34 :var __date__: date of the version
35 """
36
37 __author__="Lutz Gross, l.gross@uq.edu.au"
38
39 from . import escriptcpp as escore
40 from . import linearPDEs as lpe
41 from . import util
42
43 class TemperatureCartesian(lpe.TransportPDE):
44 """
45 Represents and solves the temperature advection-diffusion problem
46
47 *rhocp(T_{,t} + v_i T_{,i} - ( k T_{,i})_i = Q*
48
49 *k T_{,i}*n_i=surface_flux* and *T_{,t} = 0* where ``given_T_mask``>0.
50
51 If surface_flux is not given 0 is assumed.
52
53 Typical usage::
54
55 sp = TemperatureCartesian(domain)
56 sp.setTolerance(1.e-4)
57 t = 0
58 T = ...
59 sp.setValues(rhocp=..., v=..., k=..., given_T_mask=...)
60 sp.setInitialTemperature(T)
61 while t < t_end:
62 sp.setValue(Q=...)
63 T = sp.getTemperature(dt)
64 t += dt
65 """
66 def __init__(self,domain,**kwargs):
67 """
68 Initializes the temperature advection-diffusion problem.
69
70 :param domain: domain of the problem
71 :note: the approximation order is switched to reduced if the approximation order is nnot linear (equal to 1).
72 """
73 lpe.TransportPDE.__init__(self,domain,numEquations=1, **kwargs)
74 order=escore.Solution(domain).getApproximationOrder()
75 if order>1:
76 if escore.ReducedSolution(domain).getApproximationOrder()>1: raise ValueError("Reduced order needs to be equal to 1.")
77 self.setReducedOrderOn()
78 else:
79 self.setReducedOrderOff()
80 self.__rhocp=None
81 self.__v=None
82
83 def setInitialTemperature(self,T):
84 """
85 Same as `setInitialSolution`.
86 """
87 self.setInitialSolution(T)
88
89 def setValue(self,rhocp=None,v=None,k=None,Q=None,surface_flux=None,given_T_mask=None):
90 if rhocp!=None:
91 self.__rhocp=rhocp
92 if v!=None:
93 self.__v=v
94 if rhocp!=None:
95 super(TemperatureCartesian,self).setValue(M=self.__rhocp)
96 if (rhocp!=None or v!=None) and self.__rhocp!=None and self.__v!=None:
97 super(TemperatureCartesian,self).setValue(C=-self.__rhocp*self.__v)
98 if k!=None:
99 super(TemperatureCartesian,self).setValue(A=-k*util.kronecker(self.getDomain()))
100 if Q!=None:
101 super(TemperatureCartesian,self).setValue(Y=Q)
102 if surface_flux!=None:
103 super(TemperatureCartesian,self).setValue(y=surface_flux)
104 if given_T_mask!=None:
105 super(TemperatureCartesian,self).setValue(q=given_T_mask)
106
107 def getTemperature(self,dt,**kwargs):
108 """
109 Same as `getSolution`.
110 """
111 return self.getSolution(dt,**kwargs)
112
113
114 class Tracer(lpe.TransportPDE):
115 """
116 Represents and solves the tracer problem
117
118 *C_{,t} + v_i C_{,i} - ( k T_{,i})_i) = 0*
119
120 *C_{,t} = 0* where ``given_C_mask``>0.
121 *C_{,i}*n_i=0*
122
123 Typical usage::
124
125 sp = Tracer(domain)
126 sp.setTolerance(1.e-4)
127 t = 0
128 T = ...
129 sp.setValues(given_C_mask=...)
130 sp.setInitialTracer(C)
131 while t < t_end:
132 sp.setValue(v=...)
133 dt.getSaveTimeStepSize()
134 C = sp.getTracer(dt)
135 t += dt
136 """
137 def __init__(self,domain,useBackwardEuler=False,**kwargs):
138 """
139 Initializes the Tracer advection problem
140
141 :param domain: domain of the problem
142 :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.
143 :type useBackwardEuler: ``bool``
144 :note: the approximation order is switched to reduced if the approximation order is nnot linear (equal to 1).
145 """
146 lpe.TransportPDE.__init__(self,domain,numEquations=1,useBackwardEuler=useBackwardEuler,**kwargs)
147 order=escore.Solution(domain).getApproximationOrder()
148 if order>1:
149 if escore.ReducedSolution(domain).getApproximationOrder()>1: raise ValueError("Reduced order needs to be equal to 1.")
150 self.setReducedOrderOn()
151 else:
152 self.setReducedOrderOff()
153 super(Tracer,self).setValue(M=1.)
154
155 def setInitialTracer(self,C):
156 """
157 Same as `setInitialSolution`.
158 """
159 self.setInitialSolution(C)
160
161 def setValue(self,v=None,given_C_mask=None, k=None):
162 if v!=None:
163 super(Tracer,self).setValue(C=-v)
164 if k!=None:
165 super(Tracer,self).setValue(A=-k*util.kronecker(self.getDomain()))
166 if given_C_mask!=None:
167 super(Tracer,self).setValue(q=given_C_mask)
168
169 def getTracer(self,dt,**kwargs):
170 """
171 Same as `getSolution`.
172 """
173 return self.getSolution(dt,**kwargs)
174

  ViewVC Help
Powered by ViewVC 1.1.26