A Parabolic Problem

Consider the heat equation with constant heating $f=1$ as an example for a time dependent PDE.

\[\begin{aligned} u_t -\Delta u = f &\quad \text{in}\; \Omega \times (0,T)\\ u = 0 &\quad \text{on}\; \partial\Omega \times (0,T)\\ u = 0 &\quad \text{on}\; \Omega \times \{ 0\}. \end{aligned}\]

We start by discretizing the problem uniformly in time with a finite difference approach. Here, we include $\theta \in [0,1]$ in order to mix implicit and explicit time stepping, where $\theta=1$ means fully implicit and $\theta=0$ fully explicit.

\[\frac{u^{k}-u^{k-1}}{\Delta t} - \Delta(\theta u^{k} + (1-\theta )u^{k-1}) = f\]

Then, for every time step $k \in \mathbb{N}$, we need to solve the variational problem: Find $u^{k} \in H_0^1(\Omega)$ such that

\[\int_\Omega \frac{u^{k}-u^{k-1}}{\Delta t} v\, \mathrm{d}x + \int_\Omega \nabla (\theta u^{k} + (1-\theta )u^{k-1}) \cdot \nabla v\, \mathrm{d}x = \int_\Omega f v\, \mathrm{d}x\]

for all $v \in H_0^1(\Omega)$.

We now multiply by $\Delta t$ and move every known quantity to the right side

\[\int_\Omega u^{k} v\, \mathrm{d}x + \theta\Delta t \int_\Omega \nabla u^{k} \cdot \nabla v\, \mathrm{d}x = \int_\Omega u^{k-1} v\, \mathrm{d}x - (1-\theta)\Delta t \int_\Omega \nabla u^{k-1} \cdot \nabla v\, \mathrm{d}x + \int_\Omega f v\, \mathrm{d}x.\]

leading to the discrete formulation in terms of finite elements reading

\[(M + \theta\Delta t L) u^{k} = (M+(\theta-1)\Delta t L) u^{k-1} + Mf.\]

The following function implements a time loop to solve the problem on the interval $[0,T]$ in tsteps timesteps and outputs a .vtu-file for each step to a folder:

using MinFEM

function parabolic(mesh::Mesh, boundaryIndices::Set{Int64},
                    L::AbstractMatrix, M::AbstractMatrix,
                    f::AbstractVector, u::AbstractVector,
                    T::Float64, tsteps::Int; theta=1.0)

    dt = (T-0)/tsteps

    path = "parabolic_problem/"
    mkpath(path)
    write_to_vtk(u, mesh, "u", path*"parabolic_"*lpad(string(0), 3, '0')*".vtu")

    for i = 1:tsteps
        pde = PDESystem(A=(M + theta*dt*L), b=(M - (1.0-theta)*dt*L)*u + dt*M*f,
                        bc=zeros(mesh.nnodes), DI=boundaryIndices)
        solve!(pde)

        u = copy(pde.state)

        write_to_vtk(u, mesh, "u", path*"parabolic_"*lpad(string(i), 3, '0')*".vtu")
    end

end

The rest of the code is then again similar to the Poisson problem:

mesh = import_mesh("../meshes/Lshaped.msh")

boundary = select_boundaries(mesh)
boundary_nodes = extract_nodes(boundary)

source(x) = 1.0
f = evaluate_mesh_function(mesh, source)

initial_condition(x) = 0.0
u = evaluate_mesh_function(mesh, initial_condition)

L = assemble_laplacian(mesh)
M = assemble_massmatrix(mesh)

parabolic(mesh, boundary_nodes, L, M, f, u, 1.5, 100, theta=1.0)

From Paraview, you can access the folder created and open the file parabolic_*..vtu, which will automatically load all time steps. You can then run it via the Play button, as also explained in the tutorial. The final visualization should look similar to the following: