Time-Dependent Problems (01-cathedral-ie)

Git reference: Tutorial example 01-cathedral-ie.

Model problem

This section describes the implementation of a simple time-dependent heat transfer model that describes, in a naive approximation, how the St. Vitus cathedral in Prague responds to changes in the surrounding air temperature during one 24-hour cycle. The geometry is shown below:

We will solve the standard heat transfer equation

(1)c \varrho\frac{\partial T}{\partial t} - \lambda \Delta T = 0

equipped with a Dirichlet condition

T = T_{init}

on the bottom edge \Gamma_{ground} and a Newton condition

\frac{\partial T}{\partial \nu} = \alpha(T_{ext}(t) - T)

on the rest of the boundary \Gamma_{air}. Here, c is the heat capacity of the material, \varrho the material density, \lambda the thermal conductivity, T_{init} the fixed temperature on the ground (same as the initial temperature of the building), and \alpha the heat transfer coefficient between the building and the surrounding air. The surrounding air temperature T_{ext} is time-dependent of the form

T_{ext}(t) = T_{init} + 10\sin(2\pi t/T_{final}),

where T_{final} is 24 hours (translated into seconds).

Equation (1) is equipped with an initial condition of the form

T(x,y,0) = T_{init}(x,y) \ \ \ \mbox{in} \ \Omega.

For simplicity we will use the implicit Euler method with a constant time step \tau, which transforms equation (1) into

c \varrho\frac{T^{n+1} - T^n}{\tau} - \lambda \Delta T^{n+1} = 0.

Weak formulation

The corresponding weak formulation is

\int_{\Omega} c \varrho\frac{T^{n+1}}{\tau}v + \int_{\Omega} \lambda \nabla T^{n+1}\cdot \nabla v + \int_{\Gamma_{air}} \alpha \lambda T^{n+1}v = \int_{\Omega} c \varrho\frac{T^{n}}{\tau}v + \int_{\Gamma_{air}} \alpha \lambda T_{ext}(t^{n+1})v.

In this example we use string boundary markers:

// Boundary markers.
const std::string BDY_GROUND = "Boundary ground";
const std::string BDY_AIR = "Boundary air";

Boundary condition types are defined using the BCTypes class:

// Enter boundary markers.
BCTypes bc_types;
bc_types.add_bc_dirichlet(BDY_GROUND);
bc_types.add_bc_newton(BDY_AIR);

Values for Dirichlet boundary conditions are set via the BCValues class:

// Enter Dirichlet boundary values.
BCValues bc_values;
bc_values.add_const(BDY_GROUND, TEMP_INIT);

Then the space for the temperature T is set up:

// Initialize an H1 space with default shepeset.
H1Space space(&mesh, bc_types, essential_bc_values, P_INIT);
int ndof = Space::get_num_dofs(&space);
info("ndof = %d.", ndof);

Defining weak forms and accessing external functions

Bilinear and linear forms are defined as follows:

template<typename Real, typename Scalar>
Scalar bilinear_form(int n, double *wt, Func<Scalar> *u_ext[], Func<Real> *u, Func<Real> *v, Geom<Real> *e, ExtData<Scalar> *ext)
{
  return HEATCAP * RHO * int_u_v<Real, Scalar>(n, wt, u, v) / TAU +
         LAMBDA * int_grad_u_grad_v<Real, Scalar>(n, wt, u, v);
}

template<typename Real, typename Scalar>
Scalar linear_form(int n, double *wt, Func<Scalar> *u_ext[], Func<Real> *v, Geom<Real> *e, ExtData<Scalar> *ext)
{
  Func<Real> *temp_prev = ext->fn[0];
  return HEATCAP * RHO * int_u_v<Real, Scalar>(n, wt, temp_prev, v) / TAU;
}

template<typename Real, typename Scalar>
Scalar bilinear_form_surf(int n, double *wt, Func<Scalar> *u_ext[], Func<Real> *u, Func<Real> *v, Geom<Real> *e, ExtData<Scalar> *ext)
{
  return LAMBDA * ALPHA * int_u_v<Real, Scalar>(n, wt, u, v);
}

template<typename Real, typename Scalar>
Scalar linear_form_surf(int n, double *wt, Func<Scalar> *u_ext[], Func<Real> *v, Geom<Real> *e, ExtData<Scalar> *ext)
{
  return LAMBDA * ALPHA * temp_ext(TIME) * int_v<Real, Scalar>(n, wt, v);
}

Notice how the previous time level temperature is accessed:

Func<Real> *temp_prev = ext->fn[0];

Setting initial condition

Next we need to initialize the previous time level solution tsln with the initial condition T_{init}. Besides holding the finite element solution, the Solution class can be forced to return zero, to return a constant, or to return an arbitrary function using the methods set_zero(), set_const() and set_exact(), respectively. Here we simply call set_const() and supply the initial temperature:

// Set constant initial condition.
Solution tsln(&mesh, TEMP_INIT);

Registering external functions in weak forms

The weak forms are registered as follows:

// Initialize weak formulation.
WeakForm wf;
wf.add_matrix_form(callback(bilinear_form));
wf.add_matrix_form_surf(callback(bilinear_form_surf), BDY_AIR);
wf.add_vector_form(callback(linear_form), HERMES_ANY, &tsln);
wf.add_vector_form_surf(callback(linear_form_surf), BDY_AIR);

Notice how the previous time level solution ‘tsln’ is registered. A few lines above we saw how it is accessed from inside the weak form.

Initializing the discrete problem

Next, the DiscreteProblem class and the matrix solver structures are initialized:

// Initialize the FE problem.
bool is_linear = true;
DiscreteProblem dp(&wf, &space, is_linear);

// Set up the solver, matrix, and rhs according to the solver selection.
SparseMatrix* matrix = create_matrix(matrix_solver);
Vector* rhs = create_vector(matrix_solver);
Solver* solver = create_linear_solver(matrix_solver, matrix, rhs);

Assembling and the ‘rhs_only’ flag

We are now ready to start the time stepping. Since the stiffness matrix does not depend on the solution, it only needs to be assembled once in the first time step. For all remaining time steps it will be the same, and we just need to re-construct the load vector. This is done via the Boolean variable rhsonly which is set to false before the time stepping begins. For completeness, we show the entire time stepping loop below:

// Time stepping:
int ts = 1;
bool rhs_only = false;
do
{
  info("---- Time step %d, time %3.5f, ext_temp %g", ts, current_time, temp_ext(current_time));

  // First time assemble both the stiffness matrix and right-hand side vector,
  // then just the right-hand side vector.
  if (rhs_only == false) info("Assembling the stiffness matrix and right-hand side vector.");
  else info("Assembling the right-hand side vector (only).");
  dp.assemble(matrix, rhs, rhs_only);
  rhs_only = true;

  // Solve the linear system and if successful, obtain the solution.
  info("Solving the matrix problem.");
  if(solver->solve()) Solution::vector_to_solution(solver->get_solution(), &space, &tsln);
  else error ("Matrix solver failed.\n");

  // Visualize the solution.
  char title[100];
  sprintf(title, "Time %3.2f, exterior temperature %3.5f", current_time, temp_ext(current_time));
  Tview.set_title(title);
  Tview.show(&tsln);

  // Update global time.
  current_time += time_step;

  // Increase time step counter
  ts++;
}
while (current_time < T_FINAL);