The Transition Problem

Now we are ready to solve the original time-dependent problem. The main part of the script is the loop over time $t$ which takes the following form:

kappa, rhocp, eta and Tref are input parameters of the model. q is the heat source in the domain and h is the time step size. The variable T holds the current temperature. It is used to calculate the right hand side coefficient f in the Helmholtz equation in Equation (1.22). The statement T=mypde.getSolution() overwrites T with the temperature of the new time step $\var{t}+\var{h}$. To get this iterative process going we need to specify the initial temperature distribution, which equal to $T\hackscore{ref}$. The LinearPDE class object mypde and coefficients not changing over time are set up before the loop over time is entered. In each time step only the coefficient $Y$ is reset as it depends on the temperature of the previous time step. This allows the PDE solver to reuse information from previous time steps as much as possible.

The heat source $q\hackscore H$ which is defined in Equation (1.16) is qc in an area defined as a circle of radius r and center xc and zero outside this circle. q0 is a fixed constant. The following script defines $q\hackscore H$ as desired:

x is a Data class object of the esys.escript module defining locations in the Domain mydomain. The length() imported from the esys.escript module returns the Euclidean norm:

$$\Vert y\Vert=\sqrt{ y\hackscore{i} y\hackscore{i} } = \function{esys.escript.length}(y)$$ (29)

So length(x-xc) calculates the distances of the location x to the center of the circle xc where the heat source is acting. Note that the coordinates of xc are defined as a list of floating point numbers. It is automatically converted into a Data class object before being subtracted from x. The function whereNegative applied to length(x-xc)-r, returns a Data object which is equal to one where the object is negative (inside the circle) and zero elsewhere. After multiplication with qc we get a function with the desired property of having value qc inside the circle and zero elsewhere.

Now we can put the components together to create the script diffusion.py which is available in the example directory: :

The script will create the files T.1.xml, T.2.xml, $\ldots$, T.50.xml in the directory where the script has been started. The files give the temperature distributions at time steps $1$, $2$, $\ldots$, $50$ in the VTK file format.

Figure 1.6: Results of the Temperature Diffusion Problem for Time Steps $1$ $16$, $32$ and $48$.

Figure 1.6 shows the result for some selected time steps. An easy way to visualize the results is the command

mayavi -d T.1.xml -m SurfaceMap &

Use the Configure Data window in mayavi to move forward and and backwards in time.