heat.tex 17.9 KB
 Thomas Witkowski committed Oct 23, 2008 1 2 3 \section{Time dependent problem} \label{s:heat}  Thomas Witkowski committed Jun 07, 2010 4 5 This is an example for a time dependent scalar problem. The problem is described by the heat equation  Thomas Witkowski committed Oct 23, 2008 6 7 8 9 10 11 12 \begin{eqnarray} \label{eq:heat} \partial_t u - \Delta u &=& f ~~~~~~~ \mbox{in } \Omega \subset \mathbb{R}^{dim}\times\left({t^{begin}},t^{end}\right)\\ u &=& g ~~~~~~~ \mbox{on } \partial\Omega \times\left({t^{begin}},t^{end}\right)\\ u &=& u_0 ~~~~~ \mbox{on } \Omega \times \left({t^{begin}}\right). \end{eqnarray}  Thomas Witkowski committed Jun 07, 2010 13 14 15 16 We solve the problem in the time interval $\left(t^{begin},t^{end}\right)$ with Dirichlet boundary conditions on $\partial\Omega$. The problem is constructed, such that the exact solution is $u(x,t) = \mbox{sin}(\pi t)e^{-10x^2}$. So we set  Thomas Witkowski committed Oct 23, 2008 17 18 19 20 21 22 \begin{eqnarray} f(x,t) &=& \pi \mbox{cos}(\pi t)e^{-10x^2}-\left(400x^2-20dow\right)\mbox{sin}(\pi t)e^{-10x^2}\\ g(x,t) &=& \mbox{sin}(\pi t)e^{-10x^2} \\ u_0(x) &=& \mbox{sin}(\pi t^{begin})e^{-10x^2}. \end{eqnarray}  Thomas Witkowski committed Jun 07, 2010 23 24 We use a variable time discretization scheme. Equation (\ref{eq:heat}) is approximated by  Thomas Witkowski committed Oct 23, 2008 25 26 27  \frac{u^{new} - u^{old}}{\tau} - \left( \theta\Delta u^{new} + (1-\theta)\Delta u^{old} \right) = f(\cdot, t^{old}+\theta\tau).  Thomas Witkowski committed Jun 07, 2010 28 29 30 31 32 33 34 35 36 $\tau = t^{new} - t^{old}$ is the timestep size between the old and the new problem time. $u^{new}$ is the (searched) solution at $t=t^{new}$. $u^{old}$ is the solution at $t=t^{old}$, which is already known from the last timestep. The parameter $\theta$ determines the implicit and explicit treatment of $\Delta u$. For $\theta = 0$ we have the forward explicit Euler scheme, for $\theta = 1$ the backward implicit Euler scheme. $\theta = 0.5$ results in the Crank-Nicholson scheme. If we bring all terms that depend on $u^{old}$ to the right hand side, the equation reads  Thomas Witkowski committed Oct 23, 2008 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100  \frac{u^{new}}{\tau} - \theta\Delta u^{new} = \frac{u^{old}}{\tau} + (1-\theta)\Delta u^{old} + f(\cdot, t^{old}+\theta\tau). %\begin{figure}[h] %\center %\includegraphics[width=0.5\textwidth]{fig/theta.pdf} %\caption{Time discretization with the $\theta$ scheme.} %\label{f:theta} %\end{figure} \begin{table} \center \begin{tabular}{|cc|c|c|c|} \hline & & {\bf 1d} & {\bf 2d} & {\bf 3d} \\ \hline {\bf source code} & \tt src/ & \multicolumn{3}{|c|}{\tt heat.cc} \\ \hline {\bf parameter file} & \tt init/ & \tt heat.dat.1d & \tt heat.dat.2d & \tt heat.dat.3d \\ \hline {\bf macro file} & \tt macro/ & \tt macro.stand.1d & \tt macro.stand.2d & \tt macro.stand.3d \\ \hline {\bf output files} & \tt output/ & \multicolumn{3}{|c|}{\tt heat\_.mesh, heat\_.dat} \\ \hline \end{tabular} \caption{Files of the {\tt heat} example. In the output file names, {\tt } is replaced by the time.} \end{table} \subsection{Source code} Now, we describe the crucial parts of the source code. First, the functions $f$ and $g$ are defined. In contrast to the ellipt example, the functions now are time dependent. This is implemented by deriving the function classes also from class \verb+TimedObject+. This class provides a pointer to the current time, as well as corresponding setting and getting methods. The usage of a pointer to a real value allows to manage the current time in one location. All objects that deal with the same time, point to the same value. In our example, $f$ is evaluated at $t=t^{old}+\theta\tau$, while $g$ (the Dirichlet boundary function for $u^{new}$) is evaluated at $t=t^{new}$. Function $g$ is implemented as follows: \begin{lstlisting}{} class G : public AbstractFunction >, public TimedObject { public: double operator()(const WorldVector& x) const { return sin(M_PI * (*timePtr)) * exp(-10.0 *(x * x)); } }; \end{lstlisting} The variable $timePtr$ is a base class member of \verb+TimedObject+. This pointer has to be set once before $g$ is evaluated the first time. Implementation of function $f$ is done in the same way. \begin{figure} \center \includegraphics[width=0.5\textwidth]{fig/heat_uml.pdf} \caption{UML diagram for class Heat.} \label{f:heat_uml} \end{figure} Now, we begin with the implementation of class \verb+Heat+, that represents the instationary problem. In Figure \ref{f:heat_uml}, its class diagram is shown. \verb+Heat+ is derived from class \verb+ProblemInstatScal+ which leads to following properties: \begin{description} \item \verb+Heat+ implements the \verb+ProblemTimeInterface+, so the adaptation loop can set the current time and schedule timesteps. \item \verb+Heat+ implements \verb+ProblemStatBase+ in the role as initial (stationary) problem. The adaptation loop can compute the initial solution through this interface. The single iteration steps can be overloaded by sub classes of \verb+ProblemInstatScal+. Actually, the initial solution is computed through the method \verb+solveInitialProblem+ of \verb+ProblemTimeInterface+. But this method is implemented by \verb+ProblemInstatScal+ interpreting itself as initial stationary problem. \item \verb+Heat+ knows another implementation of \verb+ProblemStatBase+: This other implementation represents a stationary problem which is solved within each timestep.\end{description} The first lines of class \verb+Heat+ are: \begin{lstlisting}{} class Heat : public ProblemInstatScal { public:  Thomas Witkowski committed Feb 22, 2010 101  Heat(ProblemScal &heatSpace)  Thomas Witkowski committed Oct 23, 2008 102 103 104 105 106  : ProblemInstatScal("heat", heatSpace) { theta = -1.0; GET_PARAMETER(0, name + "->theta", "%f", &theta); TEST_EXIT(theta >= 0)("theta not set!\n");  Thomas Witkowski committed Mar 14, 2011 107  theta1 = theta - 1.0;  Thomas Witkowski committed May 26, 2009 108  }  Thomas Witkowski committed Oct 23, 2008 109 110 111 112 113 114 115 \end{lstlisting} The argument \verb+heatSpace+ is a pointer to the stationary problem which is solved each timestep. It is directly handed to the base class constructor of \verb+ProblemInstatScal+. In the body of the constructor, $\theta$ is read from the parameter file and stored in a member variable. The member variable \verb+theta1+ stores the value of $\theta - 1$. A pointer to this value is used later as factor in the $\theta$-scheme. The next lines show the implementation of the time interface. \begin{lstlisting}{}  Thomas Witkowski committed May 26, 2009 116 117  void setTime(AdaptInfo *adaptInfo) {  Thomas Witkowski committed Feb 22, 2010 118  rhsTime = adaptInfo->getTime() - (1 - theta) * adaptInfo->getTimestep();  Thomas Witkowski committed Oct 23, 2008 119  boundaryTime = adaptInfo->getTime();  Thomas Witkowski committed May 26, 2009 120  }  Thomas Witkowski committed Oct 23, 2008 121   Thomas Witkowski committed May 26, 2009 122 123  void closeTimestep(AdaptInfo *adaptInfo) {  Thomas Witkowski committed Oct 23, 2008 124 125  ProblemInstatScal::closeTimestep(adaptInfo); WAIT;  Thomas Witkowski committed May 26, 2009 126  }  Thomas Witkowski committed Oct 23, 2008 127 128 \end{lstlisting}  Thomas Witkowski committed Jun 07, 2010 129 130 131 The method \verb+setTime+ is called by the adaptation loop to inform the problem about the current time. The right hand side function $f$ will be evaluated at $t^{old}+\theta\tau = t^{new} - (1-\theta)\tau$,  Thomas Witkowski committed Mar 14, 2011 132 133 the Dirichlet boundary function $g$ at $t^{new}$, which is the current time.  Thomas Witkowski committed Jun 07, 2010 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150  The method \verb+closeTimestep+ is called at the end of each timestep by the adaptation loop. In the default implementation of \verb+ProblemInstatScal::closeTimestep+, the solution is written to output files, if specified in the parameter file. Note that the base class implementation of a method must be explicitly called, if the method is overwritten in a sub class. The macro \verb+WAIT+ waits until the \verb+return+ key is pressed by the user, if the corresponding entry in the parameter file is set to 1. The macro \verb+WAIT_REALLY+ would wait, independent of parameter settings. If \verb+closeTimestep+ wouldn't be overloaded here, the default implementation without the \verb+WAIT+ statement would be called after each timestep. Now, the implementation of the \verb+ProblemStatBase+ interface begins. As mentioned above, the instationary problem plays the role of the initial problem by implementing this interface.  Thomas Witkowski committed Oct 23, 2008 151 152 153 154 155  \begin{lstlisting}{} void solve(AdaptInfo *adaptInfo) { problemStat->getSolution()->interpol(exactSolution);  Thomas Witkowski committed May 26, 2009 156  }  Thomas Witkowski committed Oct 23, 2008 157 158 159 160 161 162 163 164 165  void estimate(AdaptInfo *adaptInfo) { double errMax, errSum; errSum = Error::L2Err(*exactSolution, *(problemStat->getSolution()), 0, &errMax, false); adaptInfo->setEstSum(errSum, 0); adaptInfo->setEstMax(errMax, 0);  Thomas Witkowski committed May 26, 2009 166  }  Thomas Witkowski committed Oct 23, 2008 167 168 \end{lstlisting}  Thomas Witkowski committed Jun 07, 2010 169 170 171 172 173 174 175 176 177 178 Here, only the solve and the estimate step are overloaded. For the other steps, there are empty default implementations in \verb+ProblemInstatScal+. Since the mesh is not adapted in the initial problem, the initial adaptation loop will stop after one iteration. In the solve step, the exact solution is interpolated on the macro mesh and stored in the solution vector of the stationary problem. In the estimate step, the L2 error is computed. The maximal element error and the sum over all element errors are stored in \verb+adaptInfo+. To make the exact solution known to the problem, we need a setting function:  Thomas Witkowski committed Oct 23, 2008 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203  \begin{lstlisting}{} void setExactSolution(AbstractFunction > *fct) { exactSolution = fct; } \end{lstlisting} Now, we define some getting functions and the private member variables: \begin{lstlisting}{} double *getThetaPtr() { return θ }; double *getTheta1Ptr() { return &theta1; }; double *getRHSTimePtr() { return &rhsTime; }; double *getBoundaryTimePtr() { return &boundaryTime; }; private: double theta; double theta1; double rhsTime; double boundaryTime; AbstractFunction > *exactSolution; }; \end{lstlisting}  Thomas Witkowski committed Jun 07, 2010 204 205 The definition of class \verb+Heat+ is now finished. In the following, the main program is described.  Thomas Witkowski committed Oct 23, 2008 206 207 208 209 210 211 212 213 214 215 216  \begin{lstlisting}{} int main(int argc, char** argv) { // ===== check for init file ===== TEST_EXIT(argc == 2)("usage: heat initfile\n"); // ===== init parameters ===== Parameters::init(false, argv[1]); // ===== create and init stationary problem =====  Thomas Witkowski committed Feb 22, 2010 217 218  ProblemScal heatSpace("heat->space"); heatSpace.initialize(INIT_ALL);  Thomas Witkowski committed Oct 23, 2008 219 220  // ===== create instationary problem =====  Thomas Witkowski committed Feb 22, 2010 221 222  Heat heat(heatSpace); heat.initialize(INIT_ALL);  Thomas Witkowski committed Oct 23, 2008 223 224 \end{lstlisting}  Thomas Witkowski committed Jun 07, 2010 225 226 227 228 229 230 So far, the stationary space problem \verb+heatSpace+ and the instationary problem \verb+heat+ were created and initialized. \verb+heatSpace+ is an instance of \verb+ProblemScal+. \verb+heat+ is an instance of the class \verb+Heat+ we defined above. \verb+heatSpace+ is given to \verb+heat+ as its stationary problem.  Thomas Witkowski committed Oct 23, 2008 231   Thomas Witkowski committed Jun 07, 2010 232 233 The next step is the creation of the needed \verb+AdaptInfo+ objects and of the instationary adaptation loop:  Thomas Witkowski committed Oct 23, 2008 234 235 236  \begin{lstlisting}{} // create adapt info for heat  Thomas Witkowski committed Feb 22, 2010 237  AdaptInfo adaptInfo("heat->adapt");  Thomas Witkowski committed Oct 23, 2008 238 239  // create initial adapt info  Thomas Witkowski committed Feb 22, 2010 240  AdaptInfo adaptInfoInitial("heat->initial->adapt");  Thomas Witkowski committed Oct 23, 2008 241 242  // create instationary adapt  Thomas Witkowski committed Feb 22, 2010 243 244 245 246 247  AdaptInstationary adaptInstat("heat->adapt", heatSpace, adaptInfo, heat, adaptInfoInitial);  Thomas Witkowski committed Oct 23, 2008 248 249 \end{lstlisting}  Thomas Witkowski committed Jun 07, 2010 250 251 252 253 254 The object \verb+heatSpace+ is handed as \verb+ProblemIterationInterface+ (implemented by class\\ \verb+ProblemScal+) to the adaptation loop. \verb+heat+ is interpreted as \verb+ProblemTimeInterface+ (implemented by class \verb+ProblemInstatScal+).  Thomas Witkowski committed Oct 23, 2008 255   Thomas Witkowski committed Mar 14, 2011 256 The functions $f$ and $g$ are declared in the following way:  Thomas Witkowski committed Oct 23, 2008 257 258 \begin{lstlisting}{} // ===== create boundary functions =====  Thomas Witkowski committed May 26, 2009 259  G *boundaryFct = new G;  Thomas Witkowski committed Feb 22, 2010 260 261 262  boundaryFct->setTimePtr(heat.getBoundaryTimePtr()); heat.setExactSolution(boundaryFct); heatSpace.addDirichletBC(1, boundaryFct);  Thomas Witkowski committed Oct 23, 2008 263 264  // ===== create rhs functions =====  Thomas Witkowski committed Feb 22, 2010 265  int degree = heatSpace.getFESpace()->getBasisFcts()->getDegree();  Thomas Witkowski committed May 26, 2009 266  F *rhsFct = new F(degree);  Thomas Witkowski committed Feb 22, 2010 267  rhsFct->setTimePtr(heat.getRHSTimePtr());  Thomas Witkowski committed Oct 23, 2008 268 \end{lstlisting}  Thomas Witkowski committed Jun 07, 2010 269 270 271 272 The functions interpreted as \verb+TimedObject+s are linked with the corresponding time pointers by \verb+setTimePtr+. The boundary function is handed to \verb+heat+ as exact solution and as Dirichlet boundary function with identifier $1$ to \verb+heatSpace+.  Thomas Witkowski committed Oct 23, 2008 273 274 275 276 277 278 279 280 281  Now, we define the operators: \begin{lstlisting}{} // ===== create operators ===== double one = 1.0; double zero = 0.0; // create laplace  Thomas Witkowski committed Jun 07, 2010 282  Operator A(heatSpace.getFESpace());  Thomas Witkowski committed Feb 22, 2010 283 284 285 286 287 288  A.addSecondOrderTerm(new Laplace_SOT); A.setUhOld(heat.getOldSolution()); if (*(heat.getThetaPtr()) != 0.0) heatSpace.addMatrixOperator(A, heat.getThetaPtr(), &one); if (*(heat.getTheta1Ptr()) != 0.0) heatSpace.addVectorOperator(A, heat.getTheta1Ptr(), &zero);  Thomas Witkowski committed Oct 23, 2008 289 290 \end{lstlisting}  Thomas Witkowski committed Feb 22, 2010 291 292 293 294 295 296 297 298 299 300 Operator \verb+A+ represents $-\Delta u$. It is used as matrix operator on the left hand side with factor $\theta$ and as vector operator on the right hand side with factor $-(1-\theta) = \theta - 1$. These assemble factors are the second arguments of \verb+addMatrixOperator+ and \verb+addVectorOperator+. The third argument is the factor used for estimation. In this example, the estimator will consider the operator only on the left hand side with factor $1$. On the right hand side the operator is applied to the solution of the last timestep. So the old solution is handed to the operator by \verb+setUhOld+.  Thomas Witkowski committed Oct 23, 2008 301 302 303  \begin{lstlisting}{} // create zero order operator  Thomas Witkowski committed Jun 07, 2010 304  Operator C(heatSpace.getFESpace());  Thomas Witkowski committed Feb 22, 2010 305 306  C.addZeroOrderTerm(new Simple_ZOT); C.setUhOld(heat.getOldSolution());  Thomas Witkowski committed Mar 14, 2011 307 308  heatSpace.addMatrixOperator(C, heat.getInvTau(), heat.getInvTau()); heatSpace.addVectorOperator(C, heat.getInvTau(), heat.getInvTau());  Thomas Witkowski committed Oct 23, 2008 309 310 \end{lstlisting}  Thomas Witkowski committed Feb 22, 2010 311 312 The \verb+Simple_ZOT+ of operator \verb+C+ represents the zero order terms for the time discretization. On both sides of the equation $u$  Thomas Witkowski committed Mar 14, 2011 313 314 315 316 are added with the factor $\frac{1}{\tau}$ for both, assembling and error estimation. The inverse of the current timestep is returned by the function \verb+getInvTau()+, which is a member of the class \verb+ProblemInstat+.  Thomas Witkowski committed Oct 23, 2008 317   Thomas Witkowski committed Feb 22, 2010 318 319 Finally, the operator for the right hand side function $f$ is added and the adaptation loop is started:  Thomas Witkowski committed Oct 23, 2008 320 321 322  \begin{lstlisting}{} // create RHS operator  Thomas Witkowski committed Jun 07, 2010 323  Operator F(heatSpace.getFESpace());  Thomas Witkowski committed Feb 22, 2010 324 325  F.addZeroOrderTerm(new CoordsAtQP_ZOT(rhsFct)); heatSpace.addVectorOperator(F);  Thomas Witkowski committed Oct 23, 2008 326 327  // ===== start adaption loop =====  Thomas Witkowski committed Feb 22, 2010 328  int errorCode = adaptInstat.adapt();  Thomas Witkowski committed Oct 23, 2008 329 330 331 } \end{lstlisting}  Thomas Witkowski committed Feb 22, 2010 332 333 334 335 336 337 \verb+CoordsAtQP_ZOT+ is a zero order term that evaluates a given function $fct$ at all needed quadrature points. At the left hand side, it would represent the term $fct(x, t)\cdot u$, on the right hand side, just $fct(x, t)$. Note that the old solution isn't given to the operator here. Otherwise the term would represent $fct(x, t) \cdot u^{old}$ on the right hand side.  Thomas Witkowski committed Oct 23, 2008 338 339 340 341  \subsection{Parameter file} \label{s:heat parameter}  Thomas Witkowski committed Feb 22, 2010 342 343 In this section, we show only the relevant parts of the parameter file \verb+heat.dat.2d+.  Thomas Witkowski committed Oct 23, 2008 344 345 346 347 348 349 350 351 352 353 354 355 356 357 358 359 360 361 362 363 364 365 366 367 368  First the parameter $\theta$ for the time discretization is defined: \begin{lstlisting}{} heat->theta: 1.0 \end{lstlisting} Then we define the initial timestep and the time interval: \begin{lstlisting}{} heat->adapt->timestep: 0.1 heat->adapt->start time: 0.0 heat->adapt->end time: 1.0 \end{lstlisting} Now, tolerances are determined: \begin{lstlisting}{} heat->adapt->tolerance: 0.01 heat->adapt->rel space error: 0.5 heat->adapt->rel time error: 0.5 heat->adapt->time theta 1: 1.0 heat->adapt->time theta 2: 0.3 \end{lstlisting}  Thomas Witkowski committed Feb 22, 2010 369 370 371 372 373 374 375 376 377 The total tolerance is divided in a space tolerance and a time tolerance. The space tolerance is the maximal allowed space error, given by the product of \verb+tolerance+ and \verb+rel space error+. It is reached by adaptive mesh refinements. The time tolerance is the maximal allowed error, due to the timestep size. It is given by the product of \verb+tolerance+ and \verb+rel time error+ and \verb+time theta 1+. It is relevant, only if an implicit time strategy with adaptive timestep size is used. The parameter \verb+time theta 2+ is used to enlarge the timestep, if the estimated time error falls beneath a given threshold.  Thomas Witkowski committed Oct 23, 2008 378 379 380 381 382 383 384  \begin{lstlisting}{} heat->adapt->strategy: 1 heat->adapt->time delta 1: 0.7071 heat->adapt->time delta 2: 1.4142 \end{lstlisting}  Thomas Witkowski committed Feb 22, 2010 385 386 387 388 389 390 391 392 393 394 395 If \verb+strategy+ is $0$, an explicit time strategy with fixed timestep size is used. A value of $1$ stands for the implicit strategy. The time tolerance is reached by successively multiplying the timestep with \verb+time delta 1+. If the estimated timestep error is smaller than the product of \verb+tolerance+ and \verb+rel time error+ and \verb+time theta 2+ at the end of a timestep, the timestep size is multiplied by \verb+time delta 2+. The following lines determine, whether coarsening is allowed in regions with sufficient small errors, and how many refinements or coarsenings are performed for marked elements.  Thomas Witkowski committed Oct 23, 2008 396 397 398 399 400 401 402 403 404 405 406 407 408 409 410 411 412 413 414 415 416 417 \begin{lstlisting}{} heat->adapt->coarsen allowed: 1 heat->adapt->refine bisections: 2 heat->adapt->coarsen bisections: 2 \end{lstlisting} Now, the output behavior is determined: \begin{lstlisting}{} heat->space->output->filename: output/heat heat->space->output->AMDiS format: 1 heat->space->output->AMDiS mesh ext: .mesh heat->space->output->AMDiS data ext: .dat heat->space->output->write every i-th timestep: 10 heat->space->output->append index: 1 heat->space->output->index length: 6 heat->space->output->index decimals: 3 \end{lstlisting}  Thomas Witkowski committed Feb 22, 2010 418 419 420 421 422 423 424 In this example, all output filenames start with prefix \verb+output/heat+ and end with the extensions \verb+.mesh+ and \verb+.dat+. Output is written after every $10$th timestep. The time of the single solution is added after the filename prefix with 6 letters, three of them are decimals. The solution for $t=0$ e.g. would be written in the files \verb+output/heat00.000.mesh+ and \verb+output/heat00.000.dat+.  Thomas Witkowski committed Oct 23, 2008 425   Thomas Witkowski committed Feb 22, 2010 426 427 428 Finally, we set parameter \verb+WAIT+ to $1$. So each call of the macro \verb+WAIT+ in the application will lead to an interruption of the program, until the \verb+return+ key is pressed.  Thomas Witkowski committed Oct 23, 2008 429 430 431 432 433 434  \begin{lstlisting}{} WAIT: 1 \end{lstlisting} \subsection{Macro file}  Thomas Witkowski committed Feb 22, 2010 435 436 We again use the macro file \verb+macro/macro.stand.2d+, which was described in Section \ref{s:ellipt macro}.  Thomas Witkowski committed Oct 23, 2008 437 438  \subsection{Output}  Thomas Witkowski committed Feb 22, 2010 439 440 441 442 443 444 As mentioned above, the output files look like \verb+output/heat00.000.mesh+ and \verb+output/heat00.000.dat+. Depending on the corresponding value in the parameter file only the solution after every $i$-th timestep is written. In Figure \ref{f:heat}, the solution at three timesteps is visualized.  Thomas Witkowski committed Oct 23, 2008 445 446 447 448 449 450 451 452 453 454 455 456  \begin{figure} \center \begin{tabular}{ccc} \includegraphics[width=0.3\textwidth]{fig/heat1.jpg}& \includegraphics[width=0.3\textwidth]{fig/heat2.jpg}& \includegraphics[width=0.3\textwidth]{fig/heat3.jpg}\\ $t \approx 0.2$ & $t \approx 0.47$ & $t \approx 0.81$ \end{tabular} \caption{The solution at three different timesteps.} \label{f:heat} \end{figure}