/[escript]/trunk/doc/cookbook/example01.tex
ViewVC logotype

Diff of /trunk/doc/cookbook/example01.tex

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

revision 2957 by artak, Tue Mar 2 01:28:19 2010 UTC revision 2958 by artak, Tue Mar 2 06:00:39 2010 UTC
# Line 24  The first model consists of two blocks o Line 24  The first model consists of two blocks o
24  Initial temperature in \textit{Block 1} is \verb|T1| and in \textit{Block 2} is \verb|T2|.  Initial temperature in \textit{Block 1} is \verb|T1| and in \textit{Block 2} is \verb|T2|.
25  We assume that the system is insulated.  We assume that the system is insulated.
26  What would happen to the temperature distribution in each block over time?  What would happen to the temperature distribution in each block over time?
27  Intuition tells us that heat will be transported from the hotter block to the cooler until both  Intuition tells us that heat will be transported from the hotter block to the cooler one until both
28  blocks have the same temperature.  blocks have the same temperature.
29    
30  \subsection{1D Heat Diffusion Equation}  \subsection{1D Heat Diffusion Equation}
# Line 40  parameters are \textbf{constant}. Line 40  parameters are \textbf{constant}.
40  The heat source is defined by the right hand side of \refEq{eqn:hd} as $q\hackscore{H}$; this can take the form of a constant or a function of time and space. For example $q\hackscore{H} = q\hackscore{0}e^{-\gamma t}$ where we have the output of our heat source decaying with time. There are also two partial derivatives in \refEq{eqn:hd}; $\frac{\partial T}{\partial t}$ describes the change in temperature with time while $\frac{\partial ^2 T}{\partial x^2}$ is the spatial change of temperature. As there is only a single spatial dimension to our problem, our temperature solution $T$ is only dependent on the time $t$ and our position along the iron bar $x$.  The heat source is defined by the right hand side of \refEq{eqn:hd} as $q\hackscore{H}$; this can take the form of a constant or a function of time and space. For example $q\hackscore{H} = q\hackscore{0}e^{-\gamma t}$ where we have the output of our heat source decaying with time. There are also two partial derivatives in \refEq{eqn:hd}; $\frac{\partial T}{\partial t}$ describes the change in temperature with time while $\frac{\partial ^2 T}{\partial x^2}$ is the spatial change of temperature. As there is only a single spatial dimension to our problem, our temperature solution $T$ is only dependent on the time $t$ and our position along the iron bar $x$.
41    
42  \subsection{PDEs and the General Form}  \subsection{PDEs and the General Form}
43  Potentially, it is now possible to solve PDE \refEq{eqn:hd} analytically and this would produce an exact solution to our problem. However, it is not always possible or practical to solve a problem this way. Alternatively, computers can be used to solve these kinds of problems. To do this, a numerical approach is required to discretised  Potentially, it is possible to solve PDE \refEq{eqn:hd} analytically and obtain an exact solution to our problem. However, it is not always possible or practical to solve a problem this way. Alternatively, computers can be used to solve these kinds of problems. To do this, a numerical approach is required to discretise
44  the PDE \refEq{eqn:hd} in time and space so finally we are left with a finite number of equations for a finite number of spatial and time steps in the model. While discretization introduces approximations and a degree of error, we find that a sufficiently sampled model is generally accurate enough for the requirements of the modeller.  the PDE \refEq{eqn:hd} in time and space, then we left with a finite number of equations for a finite number of spatial points and time steps in the model. While discretisation introduces approximations and a degree of error, a sufficiently sampled model is generally accurate enough to satisfy the requirements of the modeller.
45    
46  Firstly, we will discretise the PDE \refEq{eqn:hd} in the time direction which will  Firstly, we discretise the PDE \refEq{eqn:hd} in time. This leaves us with a steady linear PDE which involves spatial derivatives only and needs to be solved in each time step to progress in time. \esc can help us here.
 leave as with a steady linear PDE which is involving spatial derivatives only and needs to be solved in each time  
 step to progress in time - \esc can help us here.  
47    
48  For the discretization in time we will use is the Backwards Euler approximation scheme\footnote{see \url{http://en.wikipedia.org/wiki/Euler_method}}. It bases on the  For time discretization we use the Backwards Euler approximation scheme\footnote{see \url{http://en.wikipedia.org/wiki/Euler_method}}. It is based on the
49  approximation  approximation
50  \begin{equation}  \begin{equation}
51  \frac{\partial T(t)}{\partial t} \approx \frac{T(t)-T(t-h)}{h}  \frac{\partial T(t)}{\partial t} \approx \frac{T(t)-T(t-h)}{h}
# Line 72  Substituting \refEq{eqn:Tbeuler} into \r Line 70  Substituting \refEq{eqn:Tbeuler} into \r
70  \frac{\rho c\hackscore p}{h} (T^{(n)} - T^{(n-1)}) - \kappa \frac{\partial^{2} T^{(n)}}{\partial x^{2}} = q\hackscore H  \frac{\rho c\hackscore p}{h} (T^{(n)} - T^{(n-1)}) - \kappa \frac{\partial^{2} T^{(n)}}{\partial x^{2}} = q\hackscore H
71  \label{eqn:hddisc}  \label{eqn:hddisc}
72  \end{equation}  \end{equation}
73  Notice that we evaluate the spatial derivative term at current time $t^{(n)}$ - therefore the name \textbf{backward Euler} scheme. Alternatively, one can use evaluate the spatial derivative term at the previous time $t^{(n-1)}$. This  Notice that we evaluate the spatial derivative term at current time $t^{(n)}$ - therefore the name \textbf{backward Euler} scheme. Alternatively, one can evaluate the spatial derivative term at the previous time $t^{(n-1)}$. This
74  approach is called the \textbf{forward Euler} scheme. This scheme can provide some computational advantages which  approach is called the \textbf{forward Euler} scheme. This scheme can provide some computational advantages, which
75  we are not discussed here but has the major disadvantage that depending on the  are not discussed here. However, this scheme has a major disadvantage, namely depending on the
76  material parameter as well as the discretization of the spatial derivative term the time step size $h$ needs to be chosen sufficiently small to achieve a stable temperature when progressing in time. The term \textit{stable} means  material parameter as well as the discretization of the spatial derivative term, the time step size $h$ needs to be chosen sufficiently small to achieve a stable temperature when progressing in time. The term \textit{stable} means
77  that the approximation of the temperature will not grow beyond its initial bounds and becomes non-physical.  that the approximation of the temperature will not grow beyond its initial bounds and become non-physical.
78  The backward Euler which we use here is unconditionally stable meaning that under the assumption of  The backward Euler scheme, which we use here, is unconditionally stable meaning that under the assumption of
79  physically correct problem set-up the temperature approximation remains physical for all times.  physically correct problem set-up the temperature approximation remains physical for all time steps.
80  The user needs to keep in mind that the discretization error introduced by \refEq{eqn:beuler}  The user needs to keep in mind that the discretization error introduced by \refEq{eqn:beuler}
81  is sufficiently small so a good approximation of the true temperature is calculated. It is  is sufficiently small, thus a good approximation of the true temperature is computed. It is
82  therefore crucial that the user remains critical about his/her results and for instance compares  therefore crucial that the user remains critical about his/her results and for instance compares
83  the results for different time and spatial step sizes.  the results for different time and spatial step sizes.
84    
85  To get the temperature $T^{(n)}$ at time $t^{(n)}$ we need to solve the linear  To get the temperature $T^{(n)}$ at time $t^{(n)}$ we need to solve the linear
86  differential equation \refEq{eqn:hddisc} which is only including spatial derivatives. To solve this problem  differential equation \refEq{eqn:hddisc} which only includes spatial derivatives. To solve this problem
87  we want to to use \esc.  we want to use \esc.
88    
89  \esc interfaces with any given PDE via a general form. For the purpose of this introduction we will illustrate a simpler version of the full linear PDE general form which is available in the \esc user's guide. A simplified form that suits our heat diffusion problem\footnote{In the form of the \esc users guide which using the Einstein convention is written as  In \esc any given PDE described by general form. For the purpose of this introduction we illustrate a simpler version of the general form for full linear PDEs which is available in the \esc user's guide. A simplified form that suits our heat diffusion problem\footnote{The form in the \esc users guide which uses the Einstein convention is written as
90  $-(A\hackscore{jl} u\hackscore{,l})\hackscore{,j}+D u =Y$}  $-(A\hackscore{jl} u\hackscore{,l})\hackscore{,j}+D u =Y$}
91  is described by;  is described by;
92  \begin{equation}\label{eqn:commonform nabla}  \begin{equation}\label{eqn:commonform nabla}
# Line 111  we rearrange \refEq{eqn:hddisc}; Line 109  we rearrange \refEq{eqn:hddisc};
109  \end{equation}  \end{equation}
110  The PDE is now in a form that satisfies \refEq{eqn:commonform nabla} which is required for \esc to solve our PDE. This can be done by generating a solution for successive increments in the time nodes $t^{(n)}$ where  The PDE is now in a form that satisfies \refEq{eqn:commonform nabla} which is required for \esc to solve our PDE. This can be done by generating a solution for successive increments in the time nodes $t^{(n)}$ where
111  $t^{(0)}=0$ and  $t^{(n)}=t^{(n-1)}+h$ where $h>0$ is the step size and assumed to be constant.  $t^{(0)}=0$ and  $t^{(n)}=t^{(n-1)}+h$ where $h>0$ is the step size and assumed to be constant.
112  In the following the upper index ${(n)}$ refers to a value at time $t^{(n)}$. Finally, by comparing \refEq{eqn:hdgenf} with \refEq{eqn:commonform} it can be seen that;  In the following the upper index ${(n)}$ refers to a value at time $t^{(n)}$. Finally, by comparing \refEq{eqn:hdgenf} with \refEq{eqn:commonform} one can see that;
113  \begin{equation}\label{ESCRIPT SET}  \begin{equation}\label{ESCRIPT SET}
114  u=T^{(n)};  u=T^{(n)};
115  A = \kappa; D = \frac{\rho c \hackscore{p}}{h}; f = q \hackscore{H} + \frac{\rho c\hackscore p}{h} T^{(n-1)}  A = \kappa; D = \frac{\rho c \hackscore{p}}{h}; f = q \hackscore{H} + \frac{\rho c\hackscore p}{h} T^{(n-1)}
# Line 120  A = \kappa; D = \frac{\rho c \hackscore{ Line 118  A = \kappa; D = \frac{\rho c \hackscore{
118  \subsection{Boundary Conditions}  \subsection{Boundary Conditions}
119  \label{SEC BOUNDARY COND}  \label{SEC BOUNDARY COND}
120  With the PDE sufficiently modified, consideration must now be given to the boundary conditions of our model. Typically there are two main types of boundary conditions known as \textbf{Neumann} and \textbf{Dirichlet} boundary conditions\footnote{More information on Boundary Conditions is available at Wikipedia \url{http://en.wikipedia.org/wiki/Boundary_conditions}}, respectively.  With the PDE sufficiently modified, consideration must now be given to the boundary conditions of our model. Typically there are two main types of boundary conditions known as \textbf{Neumann} and \textbf{Dirichlet} boundary conditions\footnote{More information on Boundary Conditions is available at Wikipedia \url{http://en.wikipedia.org/wiki/Boundary_conditions}}, respectively.
121  A \textbf{Dirichlet boundary condition} is conceptually simpler and is used to prescribe a known value to the unknown - in our example the temperature - on parts of the boundary or on the entire boundary of the region of interest.  A \textbf{Dirichlet boundary condition} is conceptually simpler and is used to prescribe a known value to the unknown solution (in our example the temperature) on parts of the boundary or on the entire boundary of the region of interest.
122  We discuss Dirichlet boundary condition in our second example presented in Section~\ref{Sec:1DHDv0}.  We discuss Dirichlet boundary condition in our second example presented in Section~\ref{Sec:1DHDv0}.
123    
124  We make the model assumption that the system is insulated so we need  We make the model assumption that the system is insulated so we need
# Line 137  or in a more general case as Line 135  or in a more general case as
135  \end{equation}  \end{equation}
136  where $n$  is the outer normal field \index{outer normal field} at the surface of the domain.  where $n$  is the outer normal field \index{outer normal field} at the surface of the domain.
137  The $\cdot$ (dot) refers to the  dot product of the vectors $\nabla T$ and $n$. In fact, the term $\nabla T \cdot n$ is the normal derivative of  The $\cdot$ (dot) refers to the  dot product of the vectors $\nabla T$ and $n$. In fact, the term $\nabla T \cdot n$ is the normal derivative of
138  the temperature $T$. Other notations which are used are\footnote{The \esc notation for the normal  the temperature $T$. Other notations used here are\footnote{The \esc notation for the normal
139  derivative is $T\hackscore{,i} n\hackscore i$.};  derivative is $T\hackscore{,i} n\hackscore i$.};
140  \begin{equation}  \begin{equation}
141  \nabla T \cdot n  = \frac{\partial T}{\partial n} \; .  \nabla T \cdot n  = \frac{\partial T}{\partial n} \; .
# Line 145  derivative is $T\hackscore{,i} n\hacksco Line 143  derivative is $T\hackscore{,i} n\hacksco
143  A condition of the type \refEq{NEUMAN 1} defines a \textbf{Neuman boundary condition} for the PDE.  A condition of the type \refEq{NEUMAN 1} defines a \textbf{Neuman boundary condition} for the PDE.
144    
145  The PDE \refEq{eqn:hdgenf}  The PDE \refEq{eqn:hdgenf}
146  and the Neuman boundary condition~\ref{eqn:hdgenf} (potentially together with the Dirichlet boundary condition set)  define a \textbf{boundary value problem}.  and the Neuman boundary condition~\ref{eqn:hdgenf} (potentially together with the Dirichlet boundary conditions)  define a \textbf{boundary value problem}.
147  It is a nature of a boundary value problem that it allows to make statements on the solution in the  It is the nature of a boundary value problem to allow making statements about the solution in the
148  interior of the domain from information known on the boundary only. In most cases  interior of the domain from information known on the boundary only. In most cases
149  we use the term partial differential equation but in fact mean a boundary value problem.  we use the term partial differential equation but in fact it is a boundary value problem.
150  It is important to keep in mind that boundary conditions need to be complete and consistent in the sense that  It is important to keep in mind that boundary conditions need to be complete and consistent in the sense that
151  at any point on the boundary either a Dirichlet or a Neuman boundary condition must be set.  at any point on the boundary either a Dirichlet or a Neuman boundary condition must be set.
152    
# Line 165  condition translates into Line 163  condition translates into
163  \begin{equation}\label{NEUMAN 2b}  \begin{equation}\label{NEUMAN 2b}
164  \kappa \frac{\partial T}{\partial x} = 0  \kappa \frac{\partial T}{\partial x} = 0
165  \end{equation}  \end{equation}
166  for the boundary of the domain. This is identical to the Neuman boundary condition we want to set. \esc will take care of this condition for us. We will discuss the Dirichlet boundary condition later.  for the boundary of the domain. This is identical to the Neuman boundary condition we want to set. \esc will take care of this condition for us. We discuss the Dirichlet boundary condition later.
167    
168  \subsection{Outline of the Implementation}  \subsection{Outline of the Implementation}
169  \label{sec:outline}  \label{sec:outline}
170  To solve the heat diffusion equation (equation \refEq{eqn:hd}) we will write a simple \pyt script. At this point we assume that you have some basic understanding of the \pyt programming language. If not there are some pointers and links available in Section \ref{sec:escpybas}. The script we will discuss later in details will have four major steps. Firstly we need to define the domain where we want to  To solve the heat diffusion equation (equation \refEq{eqn:hd}) we write a simple \pyt script. At this point we assume that you have some basic understanding of the \pyt programming language. If not, there are some pointers and links available in Section \ref{sec:escpybas}. The script we discuss later in details have four major steps. Firstly, we need to define the domain where we want to
171  calculate the temperature. For our problem this is the joint blocks of granite which has a rectangular shape. Secondly we need to define the PDE  calculate the temperature. For our problem this is the joint blocks of granite which has a rectangular shape. Secondly, we need to define the PDE to solve in each time step to get the updated temperature. Thirdly, we need to define the coefficients of the PDE and finally we need to solve the PDE. The last two steps need to be repeated until the final time marker has been reached. The work flow is as follows:
 we need to solve in each time step to get the updated temperature. Thirdly we need to define the the coefficients of the PDE and finally we need to solve the PDE. The last two steps need to be repeated until the final time marker has been reached. As a work flow this takes the form;  
172  \begin{enumerate}  \begin{enumerate}
173   \item create domain   \item create domain
174   \item create PDE   \item create PDE
# Line 183  we need to solve in each time step to ge Line 180  we need to solve in each time step to ge
180  \end{enumerate}  \end{enumerate}
181  \item end of calculation  \item end of calculation
182  \end{enumerate}  \end{enumerate}
183  In the terminology of \pyt the domain and PDE are represented by \textbf{objects}. The nice feature of an object is that it defined by it usage and features  In the terminology of \pyt the domain and PDE are represented by \textbf{objects}. The nice feature of an object is that it is defined by its usage and features
184  rather than its actual representation. So we will create a domain object to describe the geometry of the two  rather than its actual representation. So we will create a domain object to describe the geometry of the two
185  granite blocks. The main feature  granite blocks. The main feature
186  of the object we will use is the fact that we can define PDEs and spatially distributed values such as the temperature  of the object we use, is the fact that we can define PDEs and spatially distributed values such as the temperature
187  on a domain. In fact the domain object has many more features - most of them you will  on a domain. In fact the domain object has many more features - most of them you will
188  never use and do not need to understand. Similar a PDE object is defined by the fact that we can define the coefficients of the PDE and solve the PDE. At a  never use and do not need to understand. Similarly, to define a PDE object we use the fact that one needs only to define the coefficients of the PDE and solve the PDE. At a later stage you may use more advanced features of the PDE class, but you need to worry about them only at the point when you use them.
 later stage you may use more advanced features of the PDE class but you need to worry about them only at the point when you use them.  
189    
190    
191  \begin{figure}[t]  \begin{figure}[t]
Legend:
Removed from v.2957  
changed lines
  Added in v.2958
  ViewVC Help
Powered by ViewVC 1.1.26