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 |
|