xref: /dports/math/alberta/alberta3-920315ae1bbfd1b1fb6672d916619ac37a411e95/doc/heat.tex

\section{Heat equation}
\label{S:heat}%
\idx{implementation of model problems!heat equation}%
\idx{heat equation!implementation}

In this section we describe a model implementation for the (linear)
heat equation
\begin{alignat*}{2}
 \partial_t u - \Delta u &= f &\qquad&
                            \mbox{in } \Omega \subset \R^d \times (0,T),\\
                       u &= g &&\mbox{on } \partial \Omega \times (0,T),\\
                       u &= u_0 && \mbox{on } \Omega\times\{0\}.
\end{alignat*}

We describe here only differences to the implementation of the linear
Poisson problem.  For common (or similar) routines we refer to Section
\ref{S:poisson-impl}.

\subsection{Global variables}
Additionally to the finite element space \code{fe\_space}, the matrix
\code{matrix}, the vectors \code{u\_h} and \code{f\_h} and the
bit-mask \code{dirichlet\_mask} for marking Dirichlet
boundary-segments, we need a vector for storage of the solution
$U_{n}$ from the last time step.  This one is implemented as a global
variable, too. All these global variables are initialized in
\code{main()}.
%%
\bv\begin{lstlisting}
static DOF_REAL_VEC *u_old;
\end{lstlisting}\ev
%%
A global pointer to the \code{ADAPT\_INSTAT} structure is used for
access in the \code{build()} and \code{estimate()} routines, see below.
\bv\begin{lstlisting}
static ADAPT_INSTAT   *adapt_instat;
\end{lstlisting}\ev
Finally, a global variable \code{theta} is used for storing the
parameter $\theta$ and \code{err\_L2} for storing the actual $L^2$
error between true and discrete solution in the actual time step.
\bv\begin{lstlisting}
static REAL theta = 0.5;   /*---  parameter of the time discretization   ---*/
static REAL err_L2  = 0.0; /*---  spatial error in a single time step    ---*/
\end{lstlisting}\ev

\subsection{The main program for the heat equation}%
\label{S:heat_main}

The main function initializes all program parameters from file and
command line (compare \secref{S:parse_parameters}), generates a mesh and a
finite element space, the DOF matrix and vectors, and allocates and
fills the parameter structure \code{ADAPT\_INSTAT} for the adaptive
method for time dependent problems. This structure is accessed by
\code{get\_adapt\_instat()} which already initializes besides the
function pointers all members of this structure from the program
parameters, compare Sections \ref{S:get_adapt} and
\ref{S:ellipt_main}.

The (initial) time step size, read from the parameter file, is reduced
when an initial global mesh refinement is performed. This reduction is
automatically adapted to the order of time discretization (2nd order
when $\theta=0.5$, 1st order otherwise) and space discretization.  For
stability reasons, the time step size is scaled by a factor $10^{-3}$
if $\theta < 0.5$, see also \secref{S:heat_time_discrete}.

Finally, the function pointers for the \code{adapt\_instat()} structure
are adjusted to the problem dependent routines for the heat equation
and the complete numerical simulation is performed by a call to
\code{adapt\_method\_instat()}.

The \code{heat.c} demo-program only implements Dirichlet boundary
conditions by setting all bits of \code{dirichlet\_mask} to $1$. The
implementation of more complicated boundary conditions is exemplified
in the explanation for the \code{ellipt.c} program, see
\secref{S:poisson-impl}.
%%
\bv\begin{lstlisting}
int main(int argc, char **argv)
{
  FUNCNAME("main");
  MACRO_DATA     *data;
  MESH           *mesh;
  const BAS_FCTS *lagrange;
  int             n_refine = 0, p = 1, dim;
  char            filename[PATH_MAX];
  REAL            fac = 1.0;

  /*****************************************************************************
   * first of all, initialize the access to parameters of the init file
   ****************************************************************************/
  parse_parameters(argc, argv, "INIT/heat.dat");

  GET_PARAMETER(1, "mesh dimension", "%d", &dim);
  GET_PARAMETER(1, "macro file name", "%s", filename);
  GET_PARAMETER(1, "global refinements", "%d", &n_refine);
  GET_PARAMETER(1, "parameter theta", "%e", &theta);
  GET_PARAMETER(1, "polynomial degree", "%d", &p);
  GET_PARAMETER(1, "online graphics", "%d", &do_graphics);

  /*****************************************************************************
  *  get a mesh, and read the macro triangulation from file
  ****************************************************************************/
  data = read_macro(filename);
  mesh = GET_MESH(dim, "ALBERTA mesh", data,
                  NULL /* init_node_projection() */,
                  NULL /* init_wall_trafos() */);
  free_macro_data(data);

  init_leaf_data(mesh, sizeof(struct heat_leaf_data),
                 NULL /* refine_leaf_data() */,
                 NULL /* coarsen_leaf_data() */);

  /* Finite element spaces can be added at any time, but it is more
   * efficient to do so before refining the mesh a lot.
   */
  lagrange = get_lagrange(mesh->dim, p);
  TEST_EXIT(lagrange, "no lagrange BAS_FCTS\n");

  fe_space = get_fe_space(mesh, lagrange->name, lagrange, 1, ADM_FLAGS_DFLT);

  global_refine(mesh, n_refine * dim, FILL_NOTHING);

  /*****************************************************************************
   *  initialize the global variables shared across build(), solve()
   *  and estimate().
   ****************************************************************************/
  matrix = get_dof_matrix("A", fe_space, NULL);
  f_h    = get_dof_real_vec("f_h", fe_space);
  u_h    = get_dof_real_vec("u_h", fe_space);
  u_h->refine_interpol = fe_space->bas_fcts->real_refine_inter;
  u_h->coarse_restrict = fe_space->bas_fcts->real_coarse_inter;
  u_old  = get_dof_real_vec("u_old", fe_space);
  u_old->refine_interpol = fe_space->bas_fcts->real_refine_inter;
  u_old->coarse_restrict = fe_space->bas_fcts->real_coarse_inter;
  dof_set(0.0, u_h);      /*  initialize u_h  !                          */

  BNDRY_FLAGS_ALL(dirichlet_mask); /* Only Dirichlet b.c. supported here */

  /*****************************************************************************
   *  init adapt_instat structure
   ****************************************************************************/
  adapt_instat = get_adapt_instat(dim, "heat", "adapt", 2, adapt_instat);

  /* Some animation in between ... */
  if (do_graphics) {
    graphics(mesh, NULL, NULL, NULL, adapt_instat->start_time);
  }

  /*****************************************************************************
   *  adapt time step size to refinement level and polynomial degree,
   *  based on the known L2-error error estimates.
   ****************************************************************************/
  if (theta < 0.5) {
    WARNING("You are using the explicit Euler scheme\n");
    WARNING("Use a sufficiently small time step size!!!\n");
    fac = 1.0e-3;
  }

  if (theta == 0.5) {
    adapt_instat->timestep *= fac*pow(2, -(REAL)(p*(n_refine))/2.0);
  } else {
    adapt_instat->timestep *= fac*pow(2, -(REAL)(p*(n_refine)));
  }
  MSG("using initial timestep size = %.4le\n", adapt_instat->timestep);

  eval_time_u0 = adapt_instat->start_time;

  adapt_instat->adapt_initial->get_el_est = get_el_est;
  adapt_instat->adapt_initial->estimate = est_initial;
  adapt_instat->adapt_initial->solve = interpol_u0;

  adapt_instat->adapt_space->get_el_est   = get_el_est;
  adapt_instat->adapt_space->get_el_estc  = get_el_estc;
  adapt_instat->adapt_space->estimate = estimate;
  adapt_instat->adapt_space->build_after_coarsen = build;
  adapt_instat->adapt_space->solve = solve;

  adapt_instat->init_timestep  = init_timestep;
  adapt_instat->set_time       = set_time;
  adapt_instat->get_time_est   = get_time_est;
  adapt_instat->close_timestep = close_timestep;

  /*****************************************************************************
   * ... off we go ...
   ****************************************************************************/
  adapt_method_instat(mesh, adapt_instat);

  WAIT_REALLY;

  return 0;
}
\end{lstlisting}\ev

\subsection{The parameter file for the heat equation}%
\label{S:heat_param}

The parameter file for the heat equation \code{INIT/heat.dat} (here
for the 2d simulations) is similar to the parameter file for the
Poisson problem. The main differences are additional parameters for
the adaptive procedure, see \secref{S:adapt_instat_in_ALBERTA}. These
additional parameters may also be optimized for 1d, 2d, and 3d.

Via the parameter \code{write~finite~element~data} storage of
meshes and finite element solution for post-processing purposes
can be done. The parameter \code{write~statistical~data} selects
storage of files containing number of DOFs, estimate, error, etc. versus
time. Finally, \code{data~path} can prescribe an existing path
for storing such data.

\bv\begin{verbatim}
mesh dimension:       2
macro file name:      Macro/macro.amc
global refinements:   4
polynomial degree:    1

online graphics:      1

% graphic windows: solution, estimate, mesh, and error if size > 0
graphic windows:       400 400 400 400
% for gltools graphics you can specify the range for the values of
% discrete solution for displaying:  min max
% automatical scaling by display routine if min >= max
gltools range:  -1.0 1.0

solver:                2 % 1: BICGSTAB 2: CG 3: GMRES 4: ODIR 5: ORES
solver max iteration:  1000
solver restart:        10  %  only used for GMRES
solver tolerance:      1.e-12
solver info:           2
solver precon:         1   % 0: no precon 1: diag precon
                           % 2: HB precon 3: BPX precon
                           % 4: SSOR, omega = 1.0, #iter = 3
                           % 5: SSOR, with control over omega and #iter
                           % 6: ILU(k)
precon ssor omega:    1.0  % for precon == 5
precon ssor iter:     1    % for precon == 5
precon ilu(k):        8    % for precon == 6

parameter theta:               1.0
adapt->start_time:             0.0
adapt->end_time:               2.0

adapt->tolerance:              1.0e-3
adapt->timestep:               1.0e-1
adapt->rel_initial_error:      0.5
adapt->rel_space_error:        0.5
adapt->rel_time_error:         0.5
adapt->strategy:               0   % 0=explicit, 1=implicit
adapt->max_iteration:          10
adapt->info:                   2

adapt->initial->strategy:      3    % 0=none, 1=GR, 2=MS, 3=ES, 4=GERS
adapt->initial->MS_gamma:      0.5
adapt->initial->max_iteration: 10
adapt->initial->info:          2

adapt->space->strategy:        3    % 0=none, 1=GR, 2=MS, 3=ES, 4=GERS
adapt->space->ES_theta:        0.9
adapt->space->ES_theta_c:      0.2
adapt->space->max_iteration:   10
adapt->space->coarsen_allowed: 1   % 0|1
adapt->space->info:            2

estimator C0:                  0.1
estimator C1:                  0.1
estimator C2:                  0.1
estimator C3:                  0.1

write finite element data:     1  % write data for post-processing or not
write statistical data:        0  % write statistical data or not
data path:                     ./data  % path for data to be written

WAIT:                          0
\end{verbatim}\ev

\begin{figure}[htbp]
%\vspace*{-4mm}
\begin{center}
\includegraphics[width=0.45\hsize]{EPS/tau-2d}
\hfil
\includegraphics[width=0.45\hsize]{EPS/n_dof-2d}
\end{center}
\vspace*{-4mm}
\caption[Adaptivity: heat-equation, time-step size and number of DOFs, 2d]{Time step
  size (left) and number of DOFs for different polynomial degree
  (right) over time in 2d.}
\label{F:heat2d}
\end{figure}
\begin{figure}[htbp]
\vspace*{-4mm}
\begin{center}
\includegraphics[width=0.45\hsize]{EPS/tau-3d}
\quad
\includegraphics[width=0.45\hsize]{EPS/n_dof-3d}
\end{center}
\vspace*{-4mm}
\caption[Adaptivity: heat-equation, time-step size and number of DOFs, 3d]{Time step size (left) and number of DOFs for different
  polynomial degree (right) over time in 3d.}
\label{F:heat3d}
\end{figure}

Figures \ref{F:heat2d} and \ref{F:heat3d} show the variation of time
step sizes and number of DOFs over time, automatically generated by
the adaptive method in two and three space dimensions for a problem
with time-periodic data. The number of DOFs is depicted for different
spatial discretization order and shows the strong benefit from using a
higher order method. The size of time steps was nearly the same for
all spatial discretizations. Parameters for the adaptive procedure can
be taken from the corresponding parameter files in 2d and 3d in the
distribution.

\subsection{Functions for leaf data}%
\label{S:heat_leaf_data}

For time dependent problems, mesh adaption usually also includes
coarsening of previously (for smaller $t$) refined parts of the mesh.
For storage of local coarsening error estimates, the leaf data structure
is enlarged by a second \code{REAL}. Functions \code{rw\_el\_estc()} and
\code{get\_el\_estc()} are provided for access to that storage location,
in addition to the functions \code{rw\_el\_est()} and \code{get\_el\_est()}
which were already defined in \code{ellipt.c}.
\bv\begin{lstlisting}
struct heat_leaf_data
{
  REAL estimate;            /*  one real for the element indicator          */
  REAL est_c;               /*  one real for the coarsening indicator       */
};

static REAL *rw_el_est(EL *el)
{
  if (IS_LEAF_EL(el))
    return &((struct heat_leaf_data *)LEAF_DATA(el))->estimate;
  else
    return NULL;
}

static REAL get_el_est(EL *el)
{
  if (IS_LEAF_EL(el))
    return ((struct heat_leaf_data *)LEAF_DATA(el))->estimate;
  else
    return 0.0;
}

static REAL *rw_el_estc(EL *el)
{
  if (IS_LEAF_EL(el))
    return &((struct heat_leaf_data *)LEAF_DATA(el))->est_c;
  else
    return NULL;
}

static REAL get_el_estc(EL *el)
{
  if (IS_LEAF_EL(el))
    return ((struct heat_leaf_data *)LEAF_DATA(el))->est_c;
  else
    return 0.0;
}
\end{lstlisting}\ev

\subsection{Data of the differential equation}%
\label{S:heat_data}

Data for the heat equation are the initial values $u_0$, right hand
side $f$, and boundary values $g$. When the true solution $u$ is
known, it can be used for computing the true error between discrete
and exact solution.

The sample problem is defined such that the exact solution is
\[
   u(x,t) = \sin(\pi t) e^{-10 |x|^2} \quad\text{on } (0,1)^d \times [0,1].
\]

All library subroutines which evaluate a given data function (for
integration, e.g.) are defined for space dependent functions only and
do not know about a time variable. Thus, such a `simple' space dependent
function $f_{\textrm{space}}(x)$ has to be derived from a space--time
dependent function $f(x,t)$.  We do this by keeping the time in a
global variable, and setting
\[
  f_{\textrm{space}}(x) := f(x,t).
\]
\bv\begin{lstlisting}
static REAL eval_time_u = 0.0;
static REAL u(const REAL_D x)
{
  return sin(M_PI*eval_time_u)*exp(-10.0*SCP_DOW(x,x));
}

static REAL eval_time_u0 = 0.0;
static REAL u0(const REAL_D x)
{
  eval_time_u = eval_time_u0;
  return u(x);
}

static REAL eval_time_g = 0.0;
static REAL g(const REAL_D x)              /* boundary values, not optional */
{
  eval_time_u = eval_time_g;
  return u(x);
}

static REAL eval_time_f = 0.0;
static REAL f(const REAL_D x)              /* -Delta u, not optional        */
{
  REAL  r2 = SCP_DOW(x,x), ux  = sin(M_PI*eval_time_f)*exp(-10.0*r2);
  REAL  ut = M_PI*cos(M_PI*eval_time_f)*exp(-10.0*r2);
  return ut - (400.0*r2 - 20.0*DIM)*ux;
}
\end{lstlisting}\ev

  As indicated, the times for evaluation of boundary data and right
  hand side may be chosen independent of each other depending on the
  kind of time discretization. The value of \code{eval\_time\_f} and
  \code{eval\_time\_g} are set by the function \code{set\_time()}.
  Similarly, the evaluation time for the exact solution is set by
  \code{estimate()} where also the evaluation time of $f$ is set for
  the evaluation of the element residual. In order to start the
  simulation not only at $t=0$, we have introduced a variable
  \code{eval\_time\_u0}, which is set in \code{main()} at the
  beginning of the program to the value of
  \code{adapt\_instat->start\_time}.

\subsection{Time discretization}%
\label{S:heat_time_discrete}

The model implementation uses a variable time discretization scheme.
Initial data is interpolated on the initial mesh,
\[
    U_0 = I_0 u_0.
\]
For $\theta\in[0,1]$, the solution $U_{n+1} \approx u(\cdot, t_{n+1})$
is given by $U_{n+1} \in I_{n+1}g(\cdot, t_{n+1}) + \Xc_{n+1}$ such that
\begin{alignat}{1}
\label{E:heat_time_discrete}
  \frac1{\tau_{n+1}} (U_{n+1}, \Phi) + \theta (\nabla U_{n+1}, \nabla\Phi)
    =& \frac1{\tau_{n+1}} (I_{n+1}U_{n}, \Phi)
      - (1-\theta)(\nabla I_{n+1}U_{n}, \nabla\Phi) \\
      &+ (f(\cdot, t_{n}+\theta\tau_{n+1}), \Phi)
    \qquad \text{for all }\Phi\in\Xc_{n+1}. \notag
\end{alignat}
For $\theta=0$, this is the forward (explicit) Euler scheme, for
$\theta=1$ the backward (implicit) Euler scheme. For $\theta=0.5$,
we obtain the Cranck--Nicholson scheme, which is of second order in
time. For $\theta\in[0.5,1.0]$, the scheme is unconditionally stable,
while for $\theta<0.5$ stability is only guaranteed if the time step
size is small enough. For that reason, the time step size is scaled
by an additional factor of $10^{-3}$ in the main program if $\theta<0.5$.
But this might not be enough for guaranteeing stability of the scheme!
We do recommend to use $\theta = 0.5,1$ only.

\subsection{Initial data interpolation}%
\label{S:heat_initial}

Initial data $u_0$ is just interpolated on the initial mesh, thus the
\code{solve()} entry in \code{adapt\_instat->adapt\_initial} will
point to a routine \code{interpol\_u0()} which implements this by the
library interpolation routine. No \code{build()} routine is needed by
the initial mesh adaption procedure.
\bv\begin{lstlisting}
static void interpol_u0(MESH *mesh)
{
  dof_compress(mesh);
  interpol(u0, u_h);

  return;
}
\end{lstlisting}\ev


\subsection{The assemblage of the discrete system}%
\label{S:heat_build}

Using a matrix notation, the discrete problem (\ref{E:heat_time_discrete})
can be written as
\[
   \Big(\frac1{\tau_{n+1}} \boldsymbol{M}
          + \theta \boldsymbol{A} \Big) \boldsymbol{U}_{n+1}
    = \Big(\frac1{\tau_{n+1}} \boldsymbol{M}
   - (1-\theta) \boldsymbol{A} \Big) \boldsymbol{U}_{n} + \boldsymbol{F}_{n+1}.
\]
Here, $\boldsymbol{M}=(\Phi_i, \Phi_j)$ denotes the mass matrix and
$\boldsymbol{A}=(\nabla \Phi_i, \nabla \Phi_j)$ the stiffness matrix
(up to Dirichlet boundary DOFs). The system matrix on the left hand
side is not the same as the one applied to the
old solution on the right hand side. But we want to compute the
contribution of the solution form the old time step $U_n$ to the right
hand side vector efficiently by a simple matrix--vector multiplication
and thus avoiding additional element-wise integration.  For doing this
without storing both matrices $\boldsymbol{M}$ and $\boldsymbol{A}$ we
are using the element-wise strategy explained and used in
Section~\ref{S:nonlin_solve} when assembling the linearized equation
in the Newton iteration for solving the nonlinear reaction--diffusion
equation.

The subroutine \code{assemble()} generates both the system matrix and
the right hand side at the same time. The mesh elements are visited
via the non-recursive mesh traversal routines. On every leaf element,
both the element mass matrix \code{c\_mat} and the element stiffness
matrix \code{a\_mat} are calculated using the \code{el\_matrix\_fct()}
provided by \code{fill\_matrix\_info()}.  For this purpose, two
different operators (the mass and stiffness operators) are defined and
applied on each element.  The stiffness operator uses the same
\code{LALt()} function for the second order term as described in
Section \ref{S:ellipt_build}; the mass operator implements only the
constant zero order coefficient $c=1/\tau_{n+1}$, which is passed in
\code{struct op\_data} and evaluated in the function \code{c()}. The
initialization and access to these operators is done in the same way
as in Section~\ref{S:nonlin_solve} where this is described in detail.
During the non-recursive mesh traversal, the element stiffness matrix
and the mass matrix are computed and added to the global system
matrix.  Then, the contribution to the right hand side vector of the
solution from the old time step is computed by a matrix--vector
product of these element matrices with the local coefficient vector on
the element of $U_n$ and added to the global load vector (see
\tableref{tab:elementblas2} for \code{bi\_mat\_el\_vec()}).

After this step, the the right hand side $f$ and Dirichlet boundary
values $g$ are treated by the standard routines.

\begin{compatibility}
  In contrast to previous \ALBERTA versions, the element-vectors and
  -matrices are no longer flat \code{C}-arrays, but ``cooked''
  data-structures, with some support routines doing basis linear
  algebra. See \secref{S:elvecmat}.
\end{compatibility}

\bv\begin{lstlisting}
struct op_data         /* application data (resp. "user_data") */
{
  REAL_BD Lambda;      /*  the gradient of the barycentric coordinates */
  REAL    det;         /*  |det D F_S|                                 */

  REAL    tau_1;
};

static REAL c(const EL_INFO *el_info, const QUAD *quad, int iq, void *ud)
{
  struct op_data *info = (struct op_data *)ud;

  return info->tau_1*info->det;
}

static void assemble(DOF_MATRIX *matrix, DOF_REAL_VEC *fh, DOF_REAL_VEC *uh,
                     const DOF_REAL_VEC *u_old,
                     REAL theta, REAL tau,
                     REAL (*f)(const REAL_D), REAL (*g)(const REAL_D),
                     const BNDRY_FLAGS dirichlet_mask)
{
  /* Some quantities remembered across calls. Think of this routine
   * like being a "library function" ... The stuff is re-initialized
   * whenever the finite element space changes. We use fe_space->admin
   * to check for changes in the finite element space because
   * DOF_ADMIN's are persisitent within ALBERTA, while fe-space are
   * not.
   */
  static EL_MATRIX_INFO stiff_elmi, mass_elmi;
  static const DOF_ADMIN *admin = NULL;
  static const QUAD *quad = NULL;
  static struct op_data op_data[1]; /* storage for det and Lambda */

  /* Remaining (non-static) variables. */
  const BAS_FCTS *bas_fcts;
  FLAGS          fill_flag;
  REAL           *f_vec;
  int            nbf;
  EL_SCHAR_VEC   *bound;

  /* Initialize persistent variables. */
  if (admin != uh->fe_space->admin) {
    OPERATOR_INFO stiff_opi = { NULL, }, mass_opi = { NULL, };

    admin    = uh->fe_space->admin;

    stiff_opi.row_fe_space   = uh->fe_space;
    stiff_opi.quad[2]        = quad;
    stiff_opi.LALt.real      = LALt;
    stiff_opi.LALt_pw_const  = true;
    stiff_opi.LALt_symmetric = true;
    stiff_opi.user_data      = op_data;

    fill_matrix_info(&stiff_opi, &stiff_elmi);

    mass_opi.row_fe_space = uh->fe_space;
    mass_opi.quad[0]      = quad;
    mass_opi.c.real       = c;
    mass_opi.c_pw_const   = true;
    mass_opi.user_data    = op_data;

    fill_matrix_info(&mass_opi, &mass_elmi);

    quad = get_quadrature(uh->fe_space->bas_fcts->dim,
                          2*uh->fe_space->bas_fcts->degree);
  }

  op_data->tau_1 = 1.0/tau;

  /* Assemble the matrix and the right hand side. The idea is to
   * assemble the local mass and stiffness matrices only once, and to
   * use it to update both, the system matrix and the contribution of
   * the time discretisation to the RHS.
   */
  clear_dof_matrix(matrix);
  dof_set(0.0, fh);
  f_vec = fh->vec;

  bas_fcts = uh->fe_space->bas_fcts;
  nbf      = bas_fcts->n_bas_fcts;

  bound = get_el_schar_vec(bas_fcts);

  fill_flag = CALL_LEAF_EL|FILL_COORDS|FILL_BOUND;
  TRAVERSE_FIRST(uh->fe_space->mesh, -1, fill_flag) {
    const EL_REAL_VEC  *u_old_loc;
    const EL_DOF_VEC   *dof;
    const EL_BNDRY_VEC *bndry_bits;
    const EL_MATRIX *stiff_loc, *mass_loc;

    /* Get the local coefficients of u_old, boundary info, dof-mapping */
    u_old_loc  = fill_el_real_vec(NULL, el_info->el, u_old);
    dof        = get_dof_indices(NULL, uh->fe_space, el_info->el);
    bndry_bits = get_bound(NULL, bas_fcts, el_info);

    /* Initialization of values used by LALt and c. It is not
     * necessary to introduce an extra "init_element()" hook for our
     * OPERATOR_INFO structures; the line below is just what would be
     * contained in that function (compare with ellipt.c).
     *
     * Beware to replace the "..._0cd()" for co-dimension 0 by its
     * proper ..._dim() variant if ever "porting" this stuff to
     * parametric meshes on manifolds.
     */
    op_data->det = el_grd_lambda_0cd(el_info, op_data->Lambda);

    /* Obtain the local (i.e. per-element) matrices. */
    stiff_loc = stiff_elmi.el_matrix_fct(el_info, stiff_elmi.fill_info);
    mass_loc  = mass_elmi.el_matrix_fct(el_info, mass_elmi.fill_info);

    /* Translate the geometric boundary classification into
     * Dirichlet/Neumann/Interior boundary condition
     * interpretation. Inside the loop over the mesh-elements we need
     * only to care about Dirichlet boundary conditions.
     */
    dirichlet_map(bound, bndry_bits, dirichlet_mask);

    /* add theta*a(psi_i,psi_j) + 1/tau*m(4*u^3*psi_i,psi_j) */
    if (theta) {
      add_element_matrix(matrix,
                         theta, stiff_loc, NoTranspose, dof, dof, bound);
    }
    add_element_matrix(matrix, 1.0, mass_loc, NoTranspose, dof, dof, bound);

      /* compute the contributions from the old time-step:
       *
       * f += -(1-theta)*a(u_old,psi_i) + 1/tau*m(u_old,psi_i)
       */
    bi_mat_el_vec(-(1.0 - theta), stiff_loc,
                  1.0, mass_loc, u_old_loc,
                  1.0, fh, dof, bound);

  } TRAVERSE_NEXT();

  free_el_schar_vec(bound);

  /* Indicate that the boundary conditions are built into the matrix,
   * needed e.g. by the hierarchical preconditioners.
   */
  BNDRY_FLAGS_CPY(matrix->dirichlet_bndry, dirichlet_mask);

  /* Add the "force-term" to the right hand side (L2scp_...() is additive) */
  L2scp_fct_bas(f, quad, fh);

  /* Close the system by imposing suitable boundary conditions. Have a
   * look at ellipt.c for how to impose more complicated stuff; here
   * we only use Dirichlet b.c.
   */
  dirichlet_bound(fh, uh, NULL, dirichlet_mask, g);
}
\end{lstlisting}\ev

The \code{build()} routine for one time step of the heat equation is
nearly a dummy routine and just calls the \code{assemble()} routine
described above. In order to avoid holes in vectors and matrices, as a
first step, the mesh is compressed. This guarantees optimal
performance of the BLAS1 routines used in the iterative solvers.

\bv\begin{lstlisting}
static void build(MESH *mesh, U_CHAR flag)
{
  FUNCNAME("build");

  dof_compress(mesh);

  INFO(adapt_instat->adapt_space->info, 2,
    "%d DOFs for %s\n", fe_space->admin->size_used, fe_space->name);

  assemble(matrix, f_h, u_h, u_old, theta, adapt_instat->timestep,
           f, g, dirichlet_mask);
}
\end{lstlisting}\ev

The resulting linear system is solved by calling the \code{oem\_solve\_s()}
library routine. This is done via the \code{solve()} subroutine described
in \secref{S:ellipt_solve}.


\subsection{Error estimation}%
\label{S:heat_estimate}

The initial error $\|U_0-u_0\|_{L^2(\Omega)}$ is calculated exactly
(up to quadrature error) by a call to \code{L2\_err()}. Local error
contributions are written via \code{rw\_el\_est()} to the \code{estimate}
value in \code{struct heat\_leaf\_data}. The \code{err\_max} and
\code{err\_sum} of the \code{ADAPT\_STAT} structure (which will be
\code{adapt\_instat->adapt\_initial}, see below) are set accordingly.
\bv\begin{lstlisting}
static REAL est_initial(MESH *mesh, ADAPT_STAT *adapt)
{
  err_L2 = adapt->err_sum =
    L2_err(u0, u_h, NULL, false, false, rw_el_est, &adapt->err_max);
  return adapt->err_sum;
}
\end{lstlisting}\ev

  In each time step, error estimation is done by the library routine
  \code{heat\_est()}, which generates both time and space
  discretization indicators, compare Section \ref{S:heat_estimate}.
  Similar to the estimator for elliptic problems, a function
  \code{r()} is needed for computing contributions of lower order
  terms and the right hand side. The flag for passing information
  about the discrete solution $U_{n+1}$ or its gradient to \code{r()}
  is set to zero in \code{estimate()} since no lower order term is
  involved.

  Local element indicators are stored to the \code{estimate} or
  \code{est\_c} entries inside the data structure
  \code{struct heat\_leaf\_data} via
  \code{rw\_el\_est()} and \code{rw\_el\_estc()}. The \code{err\_max}
  and \code{err\_sum} entries of \code{adapt->adapt\_space} are set
  accordingly. The temporal error indicator is the return value by
  \code{heat\_est()} and is stored in a global variable for later
  access by \code{get\_time\_est()}. In this example, the true
  solution is known and thus the true error $\|u(\cdot,t_{n+1})
  -U_{n+1}\|_{L^2(\Omega)}$ is calculated additionally for comparison.
\bv\begin{lstlisting}
static REAL r(const EL_INFO *el_info,
              const QUAD *quad, int iq,
              REAL uh_at_qp, const REAL_D grd_uh_at_qp,
              REAL t)
{
  REAL_D      x;

  coord_to_world(el_info, quad->lambda[iq], x);
  eval_time_f = t;

  return -f(x);
}


static REAL estimate(MESH *mesh, ADAPT_STAT *adapt)
{
  FUNCNAME("estimate");
  static REAL    C[4] = {-1.0, 1.0, 1.0, 1.0};
  static REAL_DD A = {{0.0}};
  FLAGS          r_flag = 0;  /* = (INIT_UH|INIT_GRD_UH), if needed by r()  */
  int            n;
  REAL           space_est;

  eval_time_u = adapt_instat->time;

  if (C[0] < 0) {
    C[0] = 1.0;
    GET_PARAMETER(1, "estimator C0", "%f", &C[0]);
    GET_PARAMETER(1, "estimator C1", "%f", &C[1]);
    GET_PARAMETER(1, "estimator C2", "%f", &C[2]);
    GET_PARAMETER(1, "estimator C3", "%f", &C[3]);

    for (n = 0; n < DIM_OF_WORLD; n++) {
      A[n][n] = 1.0;   /* set diogonal of A; all other elements are zero */
    }
  }

  time_est = heat_est(u_h, u_old, adapt_instat, rw_el_est, rw_el_estc,
                      -1 /* quad_degree */,
                      C, (const REAL_D *)A, dirichlet_mask,
                      r, r_flag, NULL /* gn() */, 0 /* gn_flag */);

  space_est = adapt_instat->adapt_space->err_sum;
  err_L2 = L2_err(u, u_h, NULL, false, false, NULL, NULL);

  INFO(adapt_instat->info,2,
    "---8<---------------------------------------------------\n");
  INFO(adapt_instat->info, 2,"time = %.4le with timestep = %.4le\n",
                              adapt_instat->time, adapt_instat->timestep);
  INFO(adapt_instat->info, 2,"estimate   = %.4le, max = %.4le\n", space_est,
                              sqrt(adapt_instat->adapt_space->err_max));
  INFO(adapt_instat->info, 2,"||u-uh||L2 = %.4le, ratio = %.2lf\n", err_L2,
                              err_L2/MAX(space_est,1.e-20));

  return adapt_instat->adapt_space->err_sum;
}
\end{lstlisting}\ev


\subsection{Time steps}%
\label{S:heat_timestep}

Time dependent problems are calculated step by step in single time steps. In a
fully implicit time-adaptive strategy, each time step includes an adaptation
of the time step size as well as an adaptation of the corresponding spatial
discretization. First, the time step size is adapted and then the mesh
adaptation procedure is performed. This second part may again push the
estimate for the time discretization error over the corresponding tolerance.
In this case, the time step size is again reduced and the whole procedure is
iterated until both, time and space discretization error estimates meet the
prescribed tolerances (or until a maximal number of iterations is performed).
For details and other time-adaptive strategies see
\secref{S:adapt_instat_in_ALBERTA}.

Besides the \code{build()}, \code{solve()}, and \code{estimate()}
routines for the adaptation of the initial grid and the grids in each
time steps, additional routines for initializing time steps, setting
time step size of the actual time step, and finalizing a time
step are needed. For adaptive control of the time step size, the
function \code{get\_time\_est()} gives information
about the size of the time discretization error.  The actual time
discretization error is stored in the global variable \code{time\_est}
and its value is set in the function \code{estimate()}.

During the initialization of a new time step in
\code{init\_timestep()}, the discrete solution \code{u\_h} from the
old time step (or from interpolation of initial data) is copied to
\code{u\_old}. In the function \code{set\_time()} evaluation times for
the right hand side $f$ and Dirichlet boundary data $g$ are set
accordingly to the chosen time discretization scheme. Since a time
step can be rejected by the adaptive method by a too large time
discretization error, this function can be called several times during
the computation of a single time step.  On each call, information
about the actual time and time step size is accessed via the entries
\code{time} and \code{timestep} of the \code{adapt\_instat} structure.

After accepting the current time step size and current grid by the
adaptive method, the time step is completed by
\code{close\_time\_step()}.  The variables \code{space\_est},
\code{time\_est}, and \code{err\_L2} now hold the final estimates resp.
error, and \code{u\_h} the accepted finite element solution for this
time step. The final mesh and discrete solution can now be written to
file for post-processing purposes, depending on the parameter value of
\code{write finite element data}. The file name for the mesh/solution
consists of the data path, the base name \code{mesh}/\code{u\_h}, and
the iteration counter of the actual time step. Such a composition can
be easily constructed by the function \code{generate\_filename()},
described in \secref{S:generate_filename}.
Mesh and finite element solution are then exported to file
by the \code{write\_*\_xdr()} routines in a portable binary format.
Using this procedure, the sequence of discrete solutions can easily
be visualized by the program \code{alberta\_movi} which
is an interface to GRAPE and comes with the distribution of \ALBERTA,
compare \secref{S:graph_GRAPE}.

Depending on the parameter value of \code{write statistical data}, the
evolution of estimates, error, number of DOFs, size of time step size,
etc.  are written to files by the function
\code{write\_statistical\_data()}, which is included in \code{heat.c}
but not described here. It produces for each quantity a two-column
data file where the first column contains time and the second column
estimate, error, etc. Such data can easily be evaluated by standard
(graphic) tools.

Finally, a graphical output of the solution and the mesh is generated
via the \code{graphics()} routine already used in the previous examples.

\bv\begin{lstlisting}
static REAL time_est = 0.0;

static REAL get_time_est(MESH *mesh, ADAPT_INSTAT *adapt)
{
  return(time_est);
}

static void init_timestep(MESH *mesh, ADAPT_INSTAT *adapt)
{
  FUNCNAME("init_timestep");

  INFO(adapt_instat->info,1,
    "---8<---------------------------------------------------\n");
  INFO(adapt_instat->info, 1,"starting new timestep\n");

  dof_copy(u_h, u_old);
  return;
}

static void set_time(MESH *mesh, ADAPT_INSTAT *adapt)
{
  FUNCNAME("set_time");

  INFO(adapt->info,1,
    "---8<---------------------------------------------------\n");
  if (adapt->time == adapt->start_time)
  {
    INFO(adapt->info, 1,"start time: %.4le\n", adapt->time);
  }
  else
  {
    INFO(adapt->info, 1,"timestep for (%.4le %.4le), tau = %.4le\n",
                         adapt->time-adapt->timestep, adapt->time,
                         adapt->timestep);
  }

  eval_time_f = adapt->time - (1 - theta)*adapt->timestep;
  eval_time_g = adapt->time;

  return;
}

static void close_timestep(MESH *mesh, ADAPT_INSTAT *adapt)
{
  FUNCNAME("close_timestep");
  static REAL err_max = 0.0;                     /* max space-time error    */
  static REAL est_max = 0.0;                     /* max space-time estimate */
  static int  write_fe_data = 0, write_stat_data = 0;
  static int  step = 0;
  static char path[256] = "./";

  REAL        space_est = adapt->adapt_space->err_sum;
  REAL        tolerance = adapt->rel_time_error*adapt->tolerance;

  err_max = MAX(err_max, err_L2);
  est_max = MAX(est_max, space_est + time_est);

  INFO(adapt->info,1,
    "---8<---------------------------------------------------\n");

  if (adapt->time == adapt->start_time)
  {
    tolerance = adapt->adapt_initial->tolerance;
    INFO(adapt->info,1,"start time: %.4le\n", adapt->time);
  }
  else
  {
    tolerance += adapt->adapt_space->tolerance;
    INFO(adapt->info,1,"timestep for (%.4le %.4le), tau = %.4le\n",
                        adapt->time-adapt->timestep, adapt->time,
                        adapt->timestep);
  }
  INFO(adapt->info,2,"max. est.  = %.4le, tolerance = %.4le\n",
                      est_max, tolerance);
  INFO(adapt->info,2,"max. error = %.4le, ratio = %.2lf\n",
                      err_max, err_max/MAX(est_max,1.0e-20));

  if (!step) {
    GET_PARAMETER(1, "write finite element data", "%d", &write_fe_data);
    GET_PARAMETER(1, "write statistical data", "%d", &write_stat_data);
    GET_PARAMETER(1, "data path", "%s", path);
  }

  /*****************************************************************************
   * write mesh and discrete solution to file for post-processing
   ****************************************************************************/

  if (write_fe_data) {
    const char *fn;

    fn= generate_filename(path, "mesh", step);
    write_mesh_xdr(mesh, fn, adapt->time);
    fn= generate_filename(path, "u_h", step);
    write_dof_real_vec_xdr(u_h, fn);
  }

  step++;

  /*****************************************************************************
   * write data about estimate, error, time step size, etc.
   ****************************************************************************/

  if (write_stat_data) {
    int n_dof = fe_space->admin->size_used;
    write_statistics(path, adapt, n_dof, space_est, time_est, err_L2);
  }

  if (do_graphics) {
    graphics(mesh, u_h, get_el_est, u, adapt->time);
  }

  return;
}
\end{lstlisting}\ev


%%% Local Variables:
%%% mode: latex
%%% TeX-master: "alberta-man"
%%% End: