/[escript]/trunk/doc/user/diffusion.tex
ViewVC logotype

Diff of /trunk/doc/user/diffusion.tex

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

revision 573 by lkettle, Thu Mar 2 00:42:53 2006 UTC revision 575 by lkettle, Fri Mar 3 03:33:07 2006 UTC
# Line 26  solver for the temperature diffusion pro Line 26  solver for the temperature diffusion pro
26  The unknown temperature $T$ is a function of its location in the domain and time $t>0$. The governing equation  The unknown temperature $T$ is a function of its location in the domain and time $t>0$. The governing equation
27  in the interior of the domain is given by  in the interior of the domain is given by
28  \begin{equation}  \begin{equation}
29  \rho c\hackscore p T\hackscore{,t} - (\kappa T\hackscore{,i})\hackscore{,i} = q  \rho c\hackscore p T\hackscore{,t} - (\kappa T\hackscore{,i})\hackscore{,i} = q\hackscore H
30  \label{DIFFUSION TEMP EQ 1}  \label{DIFFUSION TEMP EQ 1}
31  \end{equation}  \end{equation}
32  where $\rho c\hackscore p$ and $\kappa$ are given material constants. In case of a composite  where $\rho c\hackscore p$ and $\kappa$ are given material constants. In case of a composite
33  material the parameters depend on their location in the domain. $q$ is  material the parameters depend on their location in the domain. $q\hackscore H$ is
34  a heat source (or sink) within the domain. We are using Einstein summation convention \index{summation convention}  a heat source (or sink) within the domain. We are using Einstein summation convention \index{summation convention}
35  as introduced in \Chap{FirstSteps}. In our case we assume $q$ to be equal to a constant heat production rate  as introduced in \Chap{FirstSteps}. In our case we assume $q\hackscore H$ to be equal to a constant heat production rate
36  $q^{c}$ on a circle or sphere with center $x^c$ and radius $r$ and $0$ elsewhere:  $q^{c}$ on a circle or sphere with center $x^c$ and radius $r$ and $0$ elsewhere:
37  \begin{equation}  \begin{equation}
38  q(x,t)=  q\hackscore H(x,t)=
39  \left\{  \left\{
40  \begin{array}{lcl}  \begin{array}{lcl}
41  q^c  & & \|x-x^c\| \le r \\  q^c  & & \|x-x^c\| \le r \\
# Line 85  T^{(n)}\approx T^{(n-1)}+T\hackscore{,t} Line 85  T^{(n)}\approx T^{(n-1)}+T\hackscore{,t}
85  This is inserted into \eqn{DIFFUSION TEMP EQ 1}. By separating the terms at  This is inserted into \eqn{DIFFUSION TEMP EQ 1}. By separating the terms at
86  $t^{(n)}$ and  $t^{(n-1)}$ one gets for $n=1,2,3\ldots$  $t^{(n)}$ and  $t^{(n-1)}$ one gets for $n=1,2,3\ldots$
87  \begin{equation}  \begin{equation}
88  \frac{\rho c\hackscore p}{h} T^{(n)} - (\kappa T^{(n)}\hackscore{,i})\hackscore{,i} = q +  \frac{\rho c\hackscore p}{h} T^{(n-1)}  \frac{\rho c\hackscore p}{h} T^{(n)} - (\kappa T^{(n)}\hackscore{,i})\hackscore{,i} = q\hackscore H +  \frac{\rho c\hackscore p}{h} T^{(n-1)}
89  \label{DIFFUSION TEMP EQ 7}  \label{DIFFUSION TEMP EQ 7}
90  \end{equation}  \end{equation}
91  where $T^{(0)}=T\hackscore{ref}$ is taken form the initial condition given by \eqn{DIFFUSION TEMP EQ 4}.  where $T^{(0)}=T\hackscore{ref}$ is taken form the initial condition given by \eqn{DIFFUSION TEMP EQ 4}.
# Line 102  first implement a solver for the  bounda Line 102  first implement a solver for the  bounda
102  \subsection{\label{DIFFUSION HELM SEC}Helmholtz Problem}  \subsection{\label{DIFFUSION HELM SEC}Helmholtz Problem}
103  The partial differential equation to be solved for $T^{(n)}$ has the form  The partial differential equation to be solved for $T^{(n)}$ has the form
104  \begin{equation}  \begin{equation}
105  \omega u  - (\kappa u\hackscore{,i})\hackscore{,i} = f  \omega T^{(n)}  - (\kappa T^{(n)}\hackscore{,i})\hackscore{,i} = f
106  \label{DIFFUSION HELM EQ 1}  \label{DIFFUSION HELM EQ 1}
107  \end{equation}  \end{equation}
108  where $u$ plays the role of $T^{(n)}$ and we set  and we set
109  \begin{equation}  \begin{equation}
110  \omega=\frac{\rho c\hackscore p}{h} \mbox{ and } f=q+\frac{\rho c\hackscore p}{h}T^{(n-1)} \;.  \omega=\frac{\rho c\hackscore p}{h} \mbox{ and } f=q\hackscore H +\frac{\rho c\hackscore p}{h}T^{(n-1)} \;.
111  \label{DIFFUSION HELM EQ 1b}  \label{DIFFUSION HELM EQ 1b}
112  \end{equation}  \end{equation}
113  With $g=\eta T\hackscore{ref}$ the radiation condition defined by \eqn{DIFFUSION TEMP EQ 2222}  With $g=\eta T\hackscore{ref}$ the radiation condition defined by \eqn{DIFFUSION TEMP EQ 2222}
114  takes the form  takes the form
115  \begin{equation}  \begin{equation}
116  \kappa u\hackscore{,i} n\hackscore{i} =  g - \eta u\mbox{ on } \Gamma  \kappa T^{(n)}\hackscore{,i} n\hackscore{i} =  g - \eta T^{(n)}\mbox{ on } \Gamma
117  \label{DIFFUSION HELM EQ 2}  \label{DIFFUSION HELM EQ 2}
118  \end{equation}  \end{equation}
119  The partial differential \eqn{DIFFUSION HELM EQ 1} together with boundary conditions of \eqn{DIFFUSION HELM EQ 2}  The partial differential \eqn{DIFFUSION HELM EQ 1} together with boundary conditions of \eqn{DIFFUSION HELM EQ 2}
# Line 148  which we have already used in \Chap{Firs Line 148  which we have already used in \Chap{Firs
148  we will make direct use of the general \LinearPDE class.  we will make direct use of the general \LinearPDE class.
149    
150  By inspecting the Helmholtz equation \index{Helmholtz equation}  By inspecting the Helmholtz equation \index{Helmholtz equation}
151    (\ref{DIFFUSION HELM EQ 1}) and boundary condition (\ref{DIFFUSION HELM EQ 2}) and
152    substituting $u$ for $T^{(n)}$
153  we can easily assign values to the coefficients in the  we can easily assign values to the coefficients in the
154  general PDE of the \LinearPDE class:  general PDE of the \LinearPDE class:
155  \begin{equation}\label{DIFFUSION HELM EQ 3}  \begin{equation}\label{DIFFUSION HELM EQ 3}
# Line 157  d=\eta & y= g &  \\ Line 159  d=\eta & y= g &  \\
159  \end{array}  \end{array}
160  \end{equation}  \end{equation}
161  $\delta\hackscore{ij}$ is the Kronecker symbol \index{Kronecker symbol} defined by $\delta\hackscore{ij}=1$ for  $\delta\hackscore{ij}$ is the Kronecker symbol \index{Kronecker symbol} defined by $\delta\hackscore{ij}=1$ for
162  $i=j$ and $0$ otherwise. Undefined coefficients are assumed to be not present\footnote{There is a difference  $i=j$ and $0$ otherwise. Undefined coefficients are assumed to be not present.\footnote{There is a difference
163  in \escript of being not present and set to zero. As not present coefficients are not processed,  in \escript of being not present and set to zero. As not present coefficients are not processed,
164  it is more efficient to leave a coefficient undefined instead of assigning zero to it.}  it is more efficient to leave a coefficient undefined instead of assigning zero to it.}
165    In this diffusion example we do not need to define a characteristic function $q$ because the
166    boundary conditions we consider in \eqn{DIFFUSION HELM EQ 2} are just the natural boundary
167    conditions which are already defined in the \LinearPDE class (shown in \eqn{LINEARPDE.SINGLE.2}).
168    
169  Defining and solving the Helmholtz equation is very easy now:  Defining and solving the Helmholtz equation is very easy now:
170  \begin{python}  \begin{python}
# Line 215  Now we want to see how we actually solve Line 220  Now we want to see how we actually solve
220  on a rectangular domain  on a rectangular domain
221  of length $l\hackscore{0}=5$ and height $l\hackscore{1}=1$. We choose a simple test solution such that we  of length $l\hackscore{0}=5$ and height $l\hackscore{1}=1$. We choose a simple test solution such that we
222  can verify the returned solution against the exact answer. Actually, we  can verify the returned solution against the exact answer. Actually, we
223  we take $u=x\hackscore{0}$ and then calculate the right hand side terms $f$ and $g$ such that  take $T=x\hackscore{0}$ (here $q\hackscore H = 0$) and then calculate the right hand side terms $f$ and $g$ such that
224  the test solution becomes the solution of the problem. If we assume $\kappa$ as being constant,  the test solution becomes the solution of the problem. If we assume $\kappa$ as being constant,
225  an easy calculation shows that we have to choose $f=\omega \cdot x\hackscore{0}$. On the boundary we get  an easy calculation shows that we have to choose $f=\omega \cdot x\hackscore{0}$. On the boundary we get
226  $\kappa n\hackscore{i} u\hackscore{,i}=\kappa n\hackscore{0}$.    $\kappa n\hackscore{i} u\hackscore{,i}=\kappa n\hackscore{0}$.  
# Line 242  mypde.setValue(A=kappa*kronecker(mydomai Line 247  mypde.setValue(A=kappa*kronecker(mydomai
247  #... calculate error of the PDE solution ...  #... calculate error of the PDE solution ...
248  u=mypde.getSolution()  u=mypde.getSolution()
249  print "error is ",Lsup(u-x[0])  print "error is ",Lsup(u-x[0])
250    saveVTK("x0.xml",sol=u)
251  \end{python}  \end{python}
252    To visualise the solution `x0.~xml' just use the command
253    \begin{python}
254    mayavi -d u.xml -m SurfaceMap &
255    \end{python}
256    and it is easy to see that the solution $T=x\hackscore{0}$ is calculated.
257    
258  The script is similar to the script \file{poisson.py} dicussed in \Chap{FirstSteps}.  The script is similar to the script \file{poisson.py} dicussed in \Chap{FirstSteps}.
259  \code{mydomain.getNormal()} returns the outer normal field on the surface of the domain. The function \function{Lsup}  \code{mydomain.getNormal()} returns the outer normal field on the surface of the domain. The function \function{Lsup}
260  imported by the \code{from esys.escript import *} statement and returns the maximum absolute value of its argument.  imported by the \code{from esys.escript import *} statement and returns the maximum absolute value of its argument.
# Line 279  and coefficients not changing over time Line 291  and coefficients not changing over time
291  $Y$ is reset as it depends on the temperature of the previous time step. This allows the PDE solver to reuse information  $Y$ is reset as it depends on the temperature of the previous time step. This allows the PDE solver to reuse information
292  from previous time steps as much as possible.  from previous time steps as much as possible.
293    
294  The heat source \var{q} which is defined in \eqn{DIFFUSION TEMP EQ 1b} is \var{qc}  The heat source $q\hackscore H$ which is defined in \eqn{DIFFUSION TEMP EQ 1b} is \var{qc}
295  in an area defined as a circle of radius \var{r} and center \var{xc} and zero outside this circle.  in an area defined as a circle of radius \var{r} and center \var{xc} and zero outside this circle.
296  \var{q0} is a fixed constant. The following script defines \var{q} as desired:    \var{q0} is a fixed constant. The following script defines $q\hackscore H$ as desired:  
297  \begin{python}  \begin{python}
298  from esys.escript import length,whereNegative  from esys.escript import length,whereNegative
299  xc=[0.02,0.002]  xc=[0.02,0.002]
300  r=0.001  r=0.001
301  x=mydomain.getX()  x=mydomain.getX()
302  q=q0*whereNegative(length(x-xc)-r)  qH=q0*whereNegative(length(x-xc)-r)
303  \end{python}  \end{python}
304  \var{x} is a \Data class object of  \var{x} is a \Data class object of
305  the \escript module defining locations in the \Domain \var{mydomain}.  the \escript module defining locations in the \Domain \var{mydomain}.
# Line 334  mypde.setSymmetryOn() Line 346  mypde.setSymmetryOn()
346  mypde.setValue(A=kappa*kronecker(mydomain),D=rhocp/h,d=eta,y=eta*Tref)  mypde.setValue(A=kappa*kronecker(mydomain),D=rhocp/h,d=eta,y=eta*Tref)
347  # ... set heat source: ....  # ... set heat source: ....
348  x=mydomain.getX()  x=mydomain.getX()
349  q=qc*whereNegative(length(x-xc)-r)  qH=qc*whereNegative(length(x-xc)-r)
350  # ... set initial temperature ....  # ... set initial temperature ....
351  T=Tref  T=Tref
352  # ... start iteration:  # ... start iteration:
# Line 342  while t<tend: Line 354  while t<tend:
354        i+=1        i+=1
355        t+=h        t+=h
356        print "time step :",t        print "time step :",t
357        mypde.setValue(Y=q+rhocp/h*T)        mypde.setValue(Y=qH+rhocp/h*T)
358        T=mypde.getSolution()        T=mypde.getSolution()
359        saveVTK("T.%d.xml"%i,temp=T)        saveVTK("T.%d.xml"%i,temp=T)
360  \end{python}  \end{python}

Legend:
Removed from v.573  
changed lines
  Added in v.575

  ViewVC Help
Powered by ViewVC 1.1.26