| 1 |
|
| 2 |
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
| 3 |
% Copyright (c) 2003-2018 by The University of Queensland |
| 4 |
% http://www.uq.edu.au |
| 5 |
% |
| 6 |
% Primary Business: Queensland, Australia |
| 7 |
% Licensed under the Apache License, version 2.0 |
| 8 |
% http://www.apache.org/licenses/LICENSE-2.0 |
| 9 |
% |
| 10 |
% Development until 2012 by Earth Systems Science Computational Center (ESSCC) |
| 11 |
% Development 2012-2013 by School of Earth Sciences |
| 12 |
% Development from 2014 by Centre for Geoscience Computing (GeoComp) |
| 13 |
% |
| 14 |
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
| 15 |
|
| 16 |
\section{The Diffusion Problem} |
| 17 |
\label{DIFFUSION CHAP} |
| 18 |
|
| 19 |
\subsection{\label{DIFFUSION OUT SEC}Outline} |
| 20 |
In this section we discuss how to solve a time-dependent temperature |
| 21 |
diffusion\index{diffusion equation} PDE for a given block of material. |
| 22 |
Within the block there is a heat source which drives temperature diffusion. |
| 23 |
On the surface, energy can radiate into the surrounding environment. |
| 24 |
\fig{DIFFUSION FIG 1} shows the configuration. |
| 25 |
|
| 26 |
\begin{figure} |
| 27 |
\centerline{\includegraphics[width=\figwidth]{DiffusionDomain}} |
| 28 |
\caption{Temperature Diffusion Problem with Circular Heat Source} |
| 29 |
\label{DIFFUSION FIG 1} |
| 30 |
\end{figure} |
| 31 |
|
| 32 |
In \Sec{DIFFUSION TEMP SEC} we present the relevant model. |
| 33 |
A time integration scheme is introduced to calculate the temperature at given time nodes $t^{(n)}$. |
| 34 |
At each time step a Helmholtz equation\index{Helmholtz equation} must be solved. |
| 35 |
The implementation of a Helmholtz equation solver will be discussed in \Sec{DIFFUSION HELM SEC}. |
| 36 |
In Section~\ref{DIFFUSION TRANS SEC} this solver is used to build a solver for |
| 37 |
the temperature diffusion problem. |
| 38 |
|
| 39 |
\subsection{\label{DIFFUSION TEMP SEC}Temperature Diffusion} |
| 40 |
The unknown temperature $T$ is a function of its location in the domain and time $t>0$. |
| 41 |
The governing equation in the interior of the domain is given by |
| 42 |
\begin{equation} |
| 43 |
\rho c_p T_{,t} - (\kappa T_{,i})_{,i} = q_H |
| 44 |
\label{DIFFUSION TEMP EQ 1} |
| 45 |
\end{equation} |
| 46 |
where $\rho c_p$ and $\kappa$ are given material constants. |
| 47 |
In case of a composite material parameters depend on their location in the domain. |
| 48 |
$q_H$ is a heat source (or sink) within the domain. |
| 49 |
We use the Einstein summation convention\index{summation convention} as introduced in \Chap{FirstSteps}. |
| 50 |
In our case we assume $q_H$ to be equal to a constant heat |
| 51 |
production rate $q^{c}$ on a circle or sphere with center $x^c$ and radius $r$, and $0$ elsewhere: |
| 52 |
\begin{equation} |
| 53 |
q_H(x,t)=\left\{ |
| 54 |
\begin{array}{l l} |
| 55 |
q^c & \quad\text{if $\|x-x^c\| \le r$}\\ |
| 56 |
0 & \quad\text{else}\\ |
| 57 |
\end{array} \right. |
| 58 |
\label{DIFFUSION TEMP EQ 1b} |
| 59 |
\end{equation} |
| 60 |
for all $x$ in the domain and time $t>0$. |
| 61 |
|
| 62 |
On the surface of the domain we specify a radiation condition which prescribes |
| 63 |
the normal component of the flux $\kappa T_{,i}$ to be proportional |
| 64 |
to the difference of the current temperature to the surrounding temperature $T_{ref}$: |
| 65 |
\begin{equation} |
| 66 |
\kappa T_{,i} n_i = \eta (T_{ref}-T) |
| 67 |
\label{DIFFUSION TEMP EQ 2} |
| 68 |
\end{equation} |
| 69 |
$\eta$ is a given material coefficient depending on the material of the block and the surrounding medium. |
| 70 |
$n_i$ is the $i$-th component of the outer normal field\index{outer normal field} at the surface of the domain. |
| 71 |
|
| 72 |
To solve the time-dependent \eqn{DIFFUSION TEMP EQ 1} the initial temperature at time $t=0$ has to be given. |
| 73 |
Here we assume that the initial temperature is the surrounding temperature: |
| 74 |
\begin{equation} |
| 75 |
T(x,0)=T_{ref} |
| 76 |
\label{DIFFUSION TEMP EQ 4} |
| 77 |
\end{equation} |
| 78 |
for all $x$ in the domain. Note that the initial conditions satisfy the |
| 79 |
boundary condition defined by \eqn{DIFFUSION TEMP EQ 2}. |
| 80 |
|
| 81 |
The temperature is calculated at discrete time nodes $t^{(n)}$ where |
| 82 |
$t^{(0)}=0$ and $t^{(n)}=t^{(n-1)}+h$, where $h>0$ is the step size which is assumed to be constant. |
| 83 |
In the following, the upper index ${(n)}$ refers to a value at time $t^{(n)}$. |
| 84 |
The simplest and most robust scheme to approximate the time derivative of the |
| 85 |
temperature is the backward Euler\index{backward Euler} scheme. |
| 86 |
The backward Euler scheme is based on the Taylor expansion of $T$ at time $t^{(n)}$: |
| 87 |
\begin{equation} |
| 88 |
T^{(n)}\approx T^{(n-1)}+T_{,t}^{(n)}(t^{(n)}-t^{(n-1)}) |
| 89 |
=T^{(n-1)} + h \cdot T_{,t}^{(n)} |
| 90 |
\label{DIFFUSION TEMP EQ 6} |
| 91 |
\end{equation} |
| 92 |
This is inserted into \eqn{DIFFUSION TEMP EQ 1}. By separating the terms at |
| 93 |
$t^{(n)}$ and $t^{(n-1)}$ one gets for $n=1,2,3\ldots$ |
| 94 |
\begin{equation} |
| 95 |
\frac{\rho c_p}{h} T^{(n)} - (\kappa T^{(n)}_{,i})_{,i} = q_H + \frac{\rho c_p}{h} T^{(n-1)} |
| 96 |
\label{DIFFUSION TEMP EQ 7} |
| 97 |
\end{equation} |
| 98 |
where $T^{(0)}=T_{ref}$ is taken form the initial condition given by \eqn{DIFFUSION TEMP EQ 4}. |
| 99 |
Together with the natural boundary condition |
| 100 |
\begin{equation} |
| 101 |
\kappa T_{,i}^{(n)} n_i = \eta (T_{ref}-T^{(n)}) |
| 102 |
\label{DIFFUSION TEMP EQ 2222} |
| 103 |
\end{equation} |
| 104 |
taken from \eqn{DIFFUSION TEMP EQ 2} this forms a boundary value problem that |
| 105 |
has to be solved for each time step. |
| 106 |
As a first step to implement a solver for the temperature diffusion problem we |
| 107 |
will implement a solver for the boundary value problem that has to be solved |
| 108 |
at each time step. |
| 109 |
|
| 110 |
\subsection{\label{DIFFUSION HELM SEC}Helmholtz Problem} |
| 111 |
The partial differential equation to be solved for $T^{(n)}$ has the form |
| 112 |
\begin{equation} |
| 113 |
\omega T^{(n)} - (\kappa T^{(n)}_{,i})_{,i} = f |
| 114 |
\label{DIFFUSION HELM EQ 1} |
| 115 |
\end{equation} |
| 116 |
and we set |
| 117 |
\begin{equation} |
| 118 |
\omega=\frac{\rho c_p}{h} \mbox{ and } f=q_H +\frac{\rho c_p}{h}T^{(n-1)} \;. |
| 119 |
\label{DIFFUSION HELM EQ 1b} |
| 120 |
\end{equation} |
| 121 |
With $g=\eta T_{ref}$ the radiation condition defined by \eqn{DIFFUSION TEMP EQ 2222} takes the form |
| 122 |
\begin{equation} |
| 123 |
\kappa T^{(n)}_{,i} n_{i} = g - \eta T^{(n)}\text{ on $\Gamma$} |
| 124 |
\label{DIFFUSION HELM EQ 2} |
| 125 |
\end{equation} |
| 126 |
The partial differential \eqn{DIFFUSION HELM EQ 1} together with boundary |
| 127 |
conditions of \eqn{DIFFUSION HELM EQ 2} is the Helmholtz equation\index{Helmholtz equation}. |
| 128 |
|
| 129 |
We want to use the \LinearPDE class provided by \escript to define and solve a |
| 130 |
general linear, steady, second order PDE such as the Helmholtz equation. |
| 131 |
For a single PDE the \LinearPDE class supports the following form: |
| 132 |
\begin{equation}\label{LINEARPDE.SINGLE.1 TUTORIAL} |
| 133 |
-(A_{jl} u_{,l})_{,j}+D u = Y |
| 134 |
\end{equation} |
| 135 |
where we show only the coefficients relevant for the problem discussed here. |
| 136 |
For the general form of a single PDE see \eqn{LINEARPDE.SINGLE.1}. |
| 137 |
The coefficients $A$ and $Y$ have to be specified through \Data objects in |
| 138 |
the \Function on the PDE or objects that can be converted into such \Data objects. |
| 139 |
$A$ is a \RankTwo and $D$ and $Y$ are scalar. |
| 140 |
The following natural boundary conditions\index{boundary condition!natural} |
| 141 |
are considered on $\Gamma$: |
| 142 |
\begin{equation}\label{LINEARPDE.SINGLE.2 TUTORIAL} |
| 143 |
n_{j}A_{jl} u_{,l}+d u= y \;. |
| 144 |
\end{equation} |
| 145 |
Notice that the coefficient $A$ is the same as in the PDE~\eqn{LINEARPDE.SINGLE.1 TUTORIAL}. |
| 146 |
The coefficients $d$ and $y$ are each a \Scalar in the \FunctionOnBoundary. |
| 147 |
Constraints\index{constraint} for the solution prescribe the value of the |
| 148 |
solution at certain locations in the domain. They have the form |
| 149 |
\begin{equation}\label{LINEARPDE.SINGLE.3 TUTORIAL} |
| 150 |
u=r \mbox{ where } q>0 |
| 151 |
\end{equation} |
| 152 |
Both $r$ and $q$ are a \Scalar where $q$ is the characteristic function\index{characteristic function} defining where the constraint is applied. |
| 153 |
The constraints defined by \eqn{LINEARPDE.SINGLE.3 TUTORIAL} override any |
| 154 |
other condition set by \eqn{LINEARPDE.SINGLE.1 TUTORIAL} or \eqn{LINEARPDE.SINGLE.2 TUTORIAL}. |
| 155 |
The \Poisson class of the \linearPDEs module, which we have already used in |
| 156 |
\Chap{FirstSteps}, is in fact a subclass of the more general \LinearPDE class. |
| 157 |
The \linearPDEs module provides a \Helmholtz class but we will make direct use |
| 158 |
of the general \LinearPDE class. |
| 159 |
|
| 160 |
By inspecting the Helmholtz equation\index{Helmholtz equation} |
| 161 |
(\ref{DIFFUSION HELM EQ 1}) and boundary condition (\ref{DIFFUSION HELM EQ 2}), |
| 162 |
and substituting $u$ for $T^{(n)}$, we can easily assign values to the |
| 163 |
coefficients in the general PDE of the \LinearPDE class: |
| 164 |
\begin{equation}\label{DIFFUSION HELM EQ 3} |
| 165 |
\begin{array}{llllll} |
| 166 |
A_{ij}=\kappa \delta_{ij} & D=\omega & Y=f \\ |
| 167 |
d=\eta & y= g & \\ |
| 168 |
\end{array} |
| 169 |
\end{equation} |
| 170 |
$\delta_{ij}$ is the Kronecker symbol\index{Kronecker symbol} defined |
| 171 |
by $\delta_{ij}=1$ for $i=j$ and $0$ otherwise. |
| 172 |
Undefined coefficients are assumed to be not present.\footnote{There is a |
| 173 |
difference in \escript for a coefficient to be not present and set to zero. |
| 174 |
Since in the former case the coefficient is not processed, it is more efficient |
| 175 |
to leave it undefined instead of assigning zero to it.} |
| 176 |
In this diffusion example we do not need to define a characteristic function |
| 177 |
$q$ because the boundary conditions we consider in \eqn{DIFFUSION HELM EQ 2} |
| 178 |
are just the natural boundary conditions which are already defined in the |
| 179 |
\LinearPDE class (shown in \eqn{LINEARPDE.SINGLE.2 TUTORIAL}). |
| 180 |
|
| 181 |
The Helmholtz equation can be set up the following way\footnote{Note that this |
| 182 |
is not a complete code. The full source code can be found in ``helmholtz.py''.}: |
| 183 |
\begin{python} |
| 184 |
mypde=LinearPDE(mydomain) |
| 185 |
mypde.setValue(A=kappa*kronecker(mydomain),D=omega,Y=f,d=eta,y=g) |
| 186 |
u=mypde.getSolution() |
| 187 |
\end{python} |
| 188 |
where we assume that \code{mydomain} is a \Domain object and |
| 189 |
\code{kappa}, \code{omega}, \code{eta}, and \code{g} are given scalar values |
| 190 |
typically \code{float} or \Data objects. The \method{setValue} method |
| 191 |
assigns values to the coefficients of the general PDE. |
| 192 |
The \method{getSolution} method solves the PDE and returns the solution |
| 193 |
\code{u} of the PDE. \function{kronecker} is an \escript function returning |
| 194 |
the Kronecker symbol. |
| 195 |
|
| 196 |
The coefficients can be set by several calls to \method{setValue} in arbitrary |
| 197 |
order. |
| 198 |
If a value is assigned to a coefficient several times, the last assigned value |
| 199 |
is used when the solution is calculated: |
| 200 |
\begin{python} |
| 201 |
mypde = LinearPDE(mydomain) |
| 202 |
mypde.setValue(A=kappa*kronecker(mydomain), d=eta) |
| 203 |
mypde.setValue(D=omega, Y=f, y=g) |
| 204 |
mypde.setValue(d=2*eta) # overwrites d=eta |
| 205 |
u=mypde.getSolution() |
| 206 |
\end{python} |
| 207 |
In some cases the solver of the PDE can make use of the fact that the PDE is |
| 208 |
symmetric\index{symmetric PDE}. A PDE is called symmetric if |
| 209 |
\begin{equation}\label{LINEARPDE.SINGLE.4 TUTORIAL} |
| 210 |
A_{jl}=A_{lj}\;. |
| 211 |
\end{equation} |
| 212 |
Note that $D$ and $d$ may have any value and the right hand sides $Y$, $y$ as |
| 213 |
well as the constraints are not relevant. The Helmholtz problem is symmetric. |
| 214 |
The \LinearPDE class provides the method \method{checkSymmetry} to check if |
| 215 |
the given PDE is symmetric. |
| 216 |
\begin{python} |
| 217 |
mypde = LinearPDE(mydomain) |
| 218 |
mypde.setValue(A=kappa*kronecker(mydomain), d=eta) |
| 219 |
print(mypde.checkSymmetry()) # returns True |
| 220 |
mypde.setValue(B=kronecker(mydomain)[0]) |
| 221 |
print(mypde.checkSymmetry()) # returns False |
| 222 |
mypde.setValue(C=kronecker(mydomain)[0]) |
| 223 |
print(mypde.checkSymmetry()) # returns True |
| 224 |
\end{python} |
| 225 |
Unfortunately, calling \method{checkSymmetry} is very expensive and is only |
| 226 |
recommended for testing and debugging purposes. |
| 227 |
The \method{setSymmetryOn} method is used to declare a PDE symmetric: |
| 228 |
\begin{python} |
| 229 |
mypde = LinearPDE(mydomain) |
| 230 |
mypde.setValue(A=kappa*kronecker(mydomain)) |
| 231 |
mypde.setSymmetryOn() |
| 232 |
\end{python} |
| 233 |
Now we want to see how we solve the Helmholtz equation on a |
| 234 |
rectangular domain of length $l_{0}=5$ and height $l_{1}=1$. |
| 235 |
We choose a simple test solution such that we can verify the returned solution |
| 236 |
against the exact answer. We take $T=x_{0}$ (here |
| 237 |
$q_H = 0$) and then calculate right hand side terms $f$ and $g$ |
| 238 |
such that the test solution becomes the solution of the problem. |
| 239 |
If we assume $\kappa$ is constant, an easy calculation shows that we |
| 240 |
have to choose $f=\omega \cdot x_{0}$. |
| 241 |
On the boundary we get $\kappa n_{i} u_{,i}=\kappa n_{0}$, so we set $g=\kappa n_{0}+\eta x_{0}$. |
| 242 |
The following script \file{helmholtz.py}\index{scripts!\file{helmholtz.py}} |
| 243 |
which is available in the \ExampleDirectory implements this test problem using |
| 244 |
the \finley PDE solver: |
| 245 |
\begin{python} |
| 246 |
from esys.escript import * |
| 247 |
from esys.escript.linearPDEs import LinearPDE |
| 248 |
from esys.finley import Rectangle |
| 249 |
from esys.weipa import saveVTK |
| 250 |
# set some parameters |
| 251 |
kappa=1. |
| 252 |
omega=0.1 |
| 253 |
eta=10. |
| 254 |
# generate domain |
| 255 |
mydomain = Rectangle(l0=5., l1=1., n0=50, n1=10) |
| 256 |
# open PDE and set coefficients |
| 257 |
mypde=LinearPDE(mydomain) |
| 258 |
mypde.setSymmetryOn() |
| 259 |
n=mydomain.getNormal() |
| 260 |
x=mydomain.getX() |
| 261 |
mypde.setValue(A=kappa*kronecker(mydomain), D=omega,Y=omega*x[0], \ |
| 262 |
d=eta, y=kappa*n[0]+eta*x[0]) |
| 263 |
# calculate error of the PDE solution |
| 264 |
u=mypde.getSolution() |
| 265 |
print("error is ",Lsup(u-x[0])) |
| 266 |
saveVTK("x0.vtu", sol=u) |
| 267 |
\end{python} |
| 268 |
To visualize the solution `x0.vtu' you can use the command |
| 269 |
\begin{verbatim} |
| 270 |
mayavi2 -d x0.vtu -m Surface |
| 271 |
\end{verbatim} |
| 272 |
and it is easy to see that the solution $T=x_{0}$ is calculated. |
| 273 |
|
| 274 |
The script is similar to the script \file{poisson.py} discussed in \Chap{FirstSteps}. |
| 275 |
\code{mydomain.getNormal()} returns the outer normal field on the surface of the domain. |
| 276 |
The function \function{Lsup} is imported by the \code{from esys.escript import *} |
| 277 |
statement and returns the maximum absolute value of its argument. |
| 278 |
The error shown by the print statement should be in the order of $10^{-7}$. |
| 279 |
As piecewise bi-linear interpolation is used by \finley to approximate the |
| 280 |
solution, and our solution is a linear function of the spatial coordinates, |
| 281 |
one might expect that the error would be zero or in the order of machine |
| 282 |
precision (typically $\approx 10^{-15}$). |
| 283 |
However most PDE packages use an iterative solver which is terminated when a |
| 284 |
given tolerance has been reached. The default tolerance is $10^{-8}$. |
| 285 |
This value can be altered by using the \method{setTolerance} method of the |
| 286 |
\LinearPDE class. |
| 287 |
|
| 288 |
\subsection{The Transition Problem} |
| 289 |
\label{DIFFUSION TRANS SEC} |
| 290 |
Now we are ready to solve the original time-dependent problem. The main |
| 291 |
part of the script is the loop over time $t$ which takes the following form: |
| 292 |
\begin{python} |
| 293 |
t=0 |
| 294 |
T=Tref |
| 295 |
mypde=LinearPDE(mydomain) |
| 296 |
mypde.setValue(A=kappa*kronecker(mydomain), D=rhocp/h, d=eta, y=eta*Tref) |
| 297 |
while t<t_end: |
| 298 |
mypde.setValue(Y=q+rhocp/h*T) |
| 299 |
T=mypde.getSolution() |
| 300 |
t+=h |
| 301 |
\end{python} |
| 302 |
\var{kappa}, \var{rhocp}, \var{eta} and \var{Tref} are input parameters of the model. |
| 303 |
\var{q} is the heat source in the domain and \var{h} is the time step size. |
| 304 |
The variable \var{T} holds the current temperature. |
| 305 |
It is used to calculate the right hand side coefficient \var{f} in the |
| 306 |
Helmholtz \eqn{DIFFUSION HELM EQ 1}. |
| 307 |
The statement \code{T=mypde.getSolution()} overwrites \var{T} with the |
| 308 |
temperature of the new time step $\var{t}+\var{h}$. |
| 309 |
To get this iterative process going we need to specify the initial temperature |
| 310 |
distribution, which is equal to $T_{ref}$. |
| 311 |
The \LinearPDE object \var{mypde} and the coefficients that do not change over |
| 312 |
time are set up before the loop is entered. |
| 313 |
In each time step only the coefficient $Y$ is reset as it depends on the |
| 314 |
temperature of the previous time step. |
| 315 |
This allows the PDE solver to reuse information from previous time steps as |
| 316 |
much as possible. |
| 317 |
|
| 318 |
The heat source $q_H$ which is defined in \eqn{DIFFUSION TEMP EQ 1b} |
| 319 |
is \var{qc} in an area defined as a circle of radius \var{r} and center |
| 320 |
\var{xc}, and zero outside this circle. \var{q0} is a fixed constant. |
| 321 |
The following script defines $q_H$ as desired: |
| 322 |
\begin{python} |
| 323 |
from esys.escript import length,whereNegative |
| 324 |
xc=[0.02, 0.002] |
| 325 |
r=0.001 |
| 326 |
x=mydomain.getX() |
| 327 |
qH=q0*whereNegative(length(x-xc)-r) |
| 328 |
\end{python} |
| 329 |
\var{x} is an \class{esys.escript.Data} object which contains locations in the \Domain \var{mydomain}. |
| 330 |
The \function{length()} function (also from the \escript module) returns the |
| 331 |
Euclidean norm: |
| 332 |
\begin{equation}\label{DIFFUSION HELM EQ 3aba} |
| 333 |
\|y\|=\sqrt{ |
| 334 |
y_{i} |
| 335 |
y_{i} |
| 336 |
} = \function{esys.escript.length}(y) |
| 337 |
\end{equation} |
| 338 |
So \code{length(x-xc)} calculates the distances of the location \var{x} to the |
| 339 |
center of the circle \var{xc} where the heat source is acting. |
| 340 |
Note that the coordinates of \var{xc} are defined as a list of floating point numbers. |
| 341 |
It is automatically converted into a \Data class object before being |
| 342 |
subtracted from \var{x}. |
| 343 |
The function \function{whereNegative} applied to \code{length(x-xc)-r} returns |
| 344 |
a \Data object which is equal to one where the object is negative (inside the |
| 345 |
circle) and zero elsewhere. |
| 346 |
After multiplication with \var{qc} we get a function with the desired property |
| 347 |
of having value \var{qc} inside the circle and zero elsewhere. |
| 348 |
|
| 349 |
Now we can put the components together to create the script \file{diffusion.py} |
| 350 |
which is available in the \ExampleDirectory\index{scripts!\file{diffusion.py}}: |
| 351 |
\begin{python} |
| 352 |
from esys.escript import * |
| 353 |
from esys.escript.linearPDEs import LinearPDE |
| 354 |
from esys.finley import Rectangle |
| 355 |
from esys.weipa import saveVTK |
| 356 |
#... set some parameters ... |
| 357 |
xc=[0.02, 0.002] |
| 358 |
r=0.001 |
| 359 |
qc=50.e6 |
| 360 |
Tref=0. |
| 361 |
rhocp=2.6e6 |
| 362 |
eta=75. |
| 363 |
kappa=240. |
| 364 |
tend=5. |
| 365 |
# ... time, time step size and counter ... |
| 366 |
t=0 |
| 367 |
h=0.1 |
| 368 |
i=0 |
| 369 |
#... generate domain ... |
| 370 |
mydomain = Rectangle(l0=0.05, l1=0.01, n0=250, n1=50) |
| 371 |
#... open PDE ... |
| 372 |
mypde=LinearPDE(mydomain) |
| 373 |
mypde.setSymmetryOn() |
| 374 |
mypde.setValue(A=kappa*kronecker(mydomain), D=rhocp/h, d=eta, y=eta*Tref) |
| 375 |
# ... set heat source: .... |
| 376 |
x=mydomain.getX() |
| 377 |
qH=qc*whereNegative(length(x-xc)-r) |
| 378 |
# ... set initial temperature .... |
| 379 |
T=Tref |
| 380 |
# ... start iteration: |
| 381 |
while t<tend: |
| 382 |
i+=1 |
| 383 |
t+=h |
| 384 |
print("time step:",t) |
| 385 |
mypde.setValue(Y=qH+rhocp/h*T) |
| 386 |
T=mypde.getSolution() |
| 387 |
saveVTK("T.%d.vtu"%i, temp=T) |
| 388 |
\end{python} |
| 389 |
The script will create the files \file{T.1.vtu}, \file{T.2.vtu}, $\ldots$, |
| 390 |
\file{T.50.vtu} in the directory where the script has been started. |
| 391 |
The files contain the temperature distributions at time steps $1, 2, i, |
| 392 |
\ldots, 50$ in the \VTK file format. |
| 393 |
|
| 394 |
\begin{figure} |
| 395 |
\centerline{\includegraphics[width=\figwidth]{DiffusionRes1}} |
| 396 |
\centerline{\includegraphics[width=\figwidth]{DiffusionRes16}} |
| 397 |
\centerline{\includegraphics[width=\figwidth]{DiffusionRes32}} |
| 398 |
\centerline{\includegraphics[width=\figwidth]{DiffusionRes48}} |
| 399 |
\caption{Results of the Temperature Diffusion Problem for Time Steps 1, 16, 32 |
| 400 |
and 48 (top to bottom)} |
| 401 |
\label{DIFFUSION FIG 2} |
| 402 |
\end{figure} |
| 403 |
\fig{DIFFUSION FIG 2} shows the result for some selected time steps. |
| 404 |
An easy way to visualize the results is the command |
| 405 |
\begin{verbatim} |
| 406 |
mayavi2 -d T.1.vtu -m Surface |
| 407 |
\end{verbatim} |
| 408 |
Use the \emph{Configure Data} window in \mayavi to move forward and backward in time. |
| 409 |
|
Parent Directory
| 
Revision Log