Time-Dependent Heat Conduction

Solve transient heat problems with ODE solvers

Run AllReset

This tutorial follows MFEM's Example 16 (ex16.cpp), which solves a time-dependent nonlinear heat equation using the method of lines, implicit ODE solvers, and nonlinear operators.

Problem: Solve $\frac{\partial u}{\partial t} = \nabla \cdot \bigl((\kappa + \alpha u) \nabla u\bigr)$ with natural (Neumann) boundary conditions on a star-shaped domain.

After FEM spatial discretization, this becomes the ODE system:

$\frac{d\mathbf{u}}{dt} = M^{-1}(-K\mathbf{u})$

where $M$ is the mass matrix and $K$ is the diffusion operator with temperature-dependent diffusivity $(\kappa + \alpha u)$ .

Step 1: Load the Mesh

MFEM Example 16 defaults to the star mesh. We load it from MFEM mesh format, just like in the FEM tutorial.

// C++ equivalent (ex16.cpp lines 95, 151-152):
const char *mesh_file = "../data/star.mesh";
Mesh *mesh = new Mesh(mesh_file, 1, 1);
int dim = mesh->Dimension();

Code [1]Run

# MFEM Example 16 Step 1: Load the star mesh
# In C++: Mesh *mesh = new Mesh(mesh_file, 1, 1);

mesh = mfem.mesh.Mesh.fromString(`${STAR_MESH}`)

dim = mesh.dimension
num_elements = mesh.numElements
num_vertices = mesh.numVertices

disp("=== Star Mesh Loaded ===")
disp(strcat("Dimension: ", num2str(dim)))
disp(strcat("Elements: ", num2str(num_elements)))
disp(strcat("Vertices: ", num2str(num_vertices)))

Step 2: Refine the Mesh

Increase resolution with uniform refinement. Each level splits every element into 4 (in 2D).

// C++ equivalent (ex16.cpp lines 162-165):
for (int lev = 0; lev < ref_levels; lev++)
{
   mesh->UniformRefinement();
}

Code [2]Run

# MFEM Example 16 Step 2: Refine mesh
# Default: ref_levels = 2

ref_levels = 2

for l = 1:ref_levels
    mesh.refineUniform()
    disp(strcat("After level ", num2str(l), ": ", num2str(mesh.numElements), " elements"))
end

disp("")
disp(strcat("Final mesh: ", num2str(mesh.numElements), " elements, ", num2str(mesh.numVertices), " vertices"))

Step 3: Visualize the Mesh

Render the refined star mesh with VTK to see the spatial domain.

Code [3]Run

# Visualize the refined star mesh
vtk.colormap("viridis")
vtk.edges("on")
vtk.mesh(mesh)
title("Star Mesh (MFEM ex16)")

Step 4: Define the Finite Element Space

Create an H1 (continuous) finite element space with quadratic elements (order 2), matching ex16 defaults.

// C++ equivalent (ex16.cpp lines 169-173):
H1_FECollection fe_coll(order, dim);
FiniteElementSpace fespace(mesh, &fe_coll);
int fe_size = fespace.GetTrueVSize();

Code [4]Run

# MFEM Example 16 Step 4: Create H1 finite element space
# order = 2 (quadratic), vdim = 1 (scalar temperature field)

order = 2
fespace = mfem.fespace.createH1(mesh, order, 1)
ndofs = fespace.ndofs

disp("=== Finite Element Space ===")
disp(strcat("Polynomial order: ", num2str(order)))
disp(strcat("Number of temperature unknowns: ", num2str(ndofs)))

Step 5: Set Initial Conditions

The initial temperature is $u = 2$ in the hot center ( $|\mathbf{x}| < 0.5$ ) and $u = 1$ elsewhere. All boundaries use natural (Neumann) conditions — no essential BCs.

// C++ equivalent (ex16.cpp lines 179-182, 387-397):
real_t InitialTemperature(const Vector &x) {
   if (x.Norml2() < 0.5) return 2.0;
   else return 1.0;
}
FunctionCoefficient u_0(InitialTemperature);
u_gf.ProjectCoefficient(u_0);

Code [5]Run

# MFEM Example 16 Step 5: Set initial conditions
# u = 2 if |x| < 0.5, u = 1 otherwise

# Create grid function by evaluating at each DOF location
u_gf = mfem.gridfunction.fromFunction(fespace, @(x) 2.0*(norm(x) < 0.5) + 1.0*(norm(x) >= 0.5))

disp("=== Initial Condition ===")
disp("u = 2.0 where |x| < 0.5 (hot center)")
disp("u = 1.0 elsewhere (cool boundary)")
disp(strcat("Max u: ", num2str(max(u_gf.getData())), ", Min u: ", num2str(min(u_gf.getData()))))

# Visualize initial temperature
vtk.colormap("hot")
vtk.gridfunction(mesh, u_gf)
title("Initial Temperature (t=0)")

Step 6: Spatial Discretization — The ConductionOperator

In ex16.cpp, the ConductionOperator class encapsulates the spatial discretization. After applying the FEM weak form, we get:

$M \frac{d\mathbf{u}}{dt} = -K\mathbf{u}$

Matrix	Definition	Role
Mass $M$	$M_{ij} = \int_\Omega \phi_i \phi_j \, dx$	Relates time derivatives to nodal values
Stiffness $K$	$K_{ij} = \int_\Omega (\kappa + \alpha u) \nabla\phi_i \cdot \nabla\phi_j \, dx$	Diffusion with temperature-dependent conductivity

The nonlinearity enters through $K$ : the stiffness matrix depends on the current solution $u$ via the diffusivity $k(u) = \kappa + \alpha u$ .

// C++ equivalent (ex16.cpp lines 284-315):
// ConductionOperator constructor:
M = new BilinearForm(&fespace);
M->AddDomainIntegrator(new MassIntegrator());
M->Assemble();

K = new BilinearForm(&fespace);
GridFunctionCoefficient u_coeff(&u_alpha_gf);
K->AddDomainIntegrator(new DiffusionIntegrator(u_coeff));
K->Assemble();

Code [6]Run

# MFEM Example 16 Step 6a: Assemble the mass matrix M
# M_ij = integral(phi_i * phi_j dx)

M_blf = mfem.bilinearForms.create(fespace)
M_blf.addMassIntegrator()
M_blf.assemble()

disp("=== Mass Matrix M ===")
disp("M_ij = integral(phi_i * phi_j dx)")
disp("Assembled successfully")

Code [7]Run

# MFEM Example 16 Step 6b: Assemble the stiffness matrix K
# K_ij = integral((kappa + alpha*u) * grad(phi_i) . grad(phi_j) dx)
#
# For the initial step, use constant diffusivity kappa
# (will be updated with temperature-dependent diffusivity in time loop)

kappa = 0.5    # base diffusivity (ex16 default)
alpha = 0.01   # nonlinearity parameter (ex16 default)

kappa_coeff = mfem.coefficients.constant(kappa)

K_blf = mfem.bilinearForms.create(fespace)
K_blf.addDiffusionIntegrator(kappa_coeff)
K_blf.assemble()

disp("=== Stiffness Matrix K ===")
disp(strcat("kappa = ", num2str(kappa), " (base diffusivity)"))
disp(strcat("alpha = ", num2str(alpha), " (nonlinearity)"))
disp("K_ij = integral(k(u) * grad(phi_i) . grad(phi_j) dx)")
disp("Assembled with constant kappa for initial step")

Step 7: Time Integration

Heat diffusion is stiff — explicit methods require impractically small time steps ( $\Delta t \leq h^2/(2\kappa d)$ ). We use Backward Euler (implicit), which is unconditionally stable:

$(M + \Delta t\, K) \, \mathbf{u}^{n+1} = M \, \mathbf{u}^n$

This is the same linear solve performed by ConductionOperator::ImplicitSolve in ex16.cpp. With constant diffusivity $\kappa$ , the system matrix $T = M + \Delta t\,K$ is constant — we assemble it once and reuse the CG solver.

// C++ equivalent (ex16.cpp ConductionOperator):
T = Add(1.0, Mmat, dt, Kmat);  // T = M + dt*K
T_solver.SetOperator(*T);
Mmat.Mult(u, z);               // RHS z = M * u^n
T_solver.Mult(z, k);           // solve T * u^{n+1} = z

Code [8]Run

# MFEM Example 16 Step 7: Backward Euler time stepping
# At each step: solve (M + dt*K) * u^{n+1} = M * u^n

dt = 0.01
t_final = 0.1
num_steps = ceil(t_final / dt)

# Finalize M for matrix-vector products
M_blf.finalize()

# Build T = M + dt*K (constant kappa for all steps)
dt_kappa_coeff = mfem.coefficients.constant(dt * kappa)
T_blf = mfem.bilinearForms.create(fespace)
T_blf.addMassIntegrator()
T_blf.addDiffusionIntegrator(dt_kappa_coeff)
T_blf.assemble()
T_blf.finalize()
T_ptr = T_blf.getSystemMatrixPointer()

# CG solver for T (reuse across steps since T is constant)
T_solver = mfem.linear.CGSolver.create()
T_solver.setRelTol(1e-10)
T_solver.setAbsTol(0.0)
T_solver.setMaxIter(200)
T_solver.setPrintLevel(1)
T_solver.setOperator(T_ptr)

# Allocate RHS grid function
rhs_gf = mfem.gridfunction.GridFunction.create(fespace)

disp("=== Time Integration ===")
disp(strcat("Backward Euler, dt = ", num2str(dt), ", ", num2str(num_steps), " steps"))
disp(strcat("kappa = ", num2str(kappa)))

t = 0
for step = 1:num_steps
    # RHS: rhs = M * u^n  (BilinearForm instance method)
    M_blf.mult(u_gf, rhs_gf)

    # Solve: (M + dt*K) * u^{n+1} = rhs
    T_solver.mult(rhs_gf.getPointer(), u_gf.getPointer())

    t = t + dt
end

disp(strcat("Reached t = ", num2str(t)))
disp(strcat("Max temperature: ", num2str(max(u_gf.getData()))))
disp(strcat("Min temperature: ", num2str(min(u_gf.getData()))))

Step 8: Visualize the Evolved Solution

After time stepping with Backward Euler, the initial hot region ( $u=2$ at center) has diffused outward. The peak temperature decreases and the temperature field becomes smoother. We render the scalar temperature field (left) alongside the mesh with edges (right).

// C++ equivalent (ex16.cpp lines 254-258):
u_gf.SetFromTrueDofs(u);
sout << "solution\n" << *mesh << u_gf << flush;

Code [9]Run

# Side-by-side: temperature field (left) and mesh (right)

subplot(1, 2, 1)
vtk.colormap("hot")
vtk.gridfunction(mesh, u_gf)
title(strcat("Temperature at t = ", num2str(t)))

subplot(1, 2, 2)
vtk.colormap("viridis")
vtk.edges("on")
vtk.mesh(mesh)
title("Star Mesh")

Step 9: 3D Visualization with AMR Hex Mesh

MFEM Example 16 supports various mesh types including 3D hexahedral meshes with adaptive mesh refinement (AMR). The amr-hex.mesh is a non-conforming mesh with hanging nodes from local refinement. We place the hot spot ( $u=2$ ) at the origin corner so the temperature gradient is visible on three faces of the cube.

// C++ equivalent:
// ex16 -m ../data/amr-hex.mesh -o 2 -r 0

Code [10]Run

# Load 3D AMR hex mesh from MFEM
mesh3d = mfem.mesh.Mesh.fromString(`${AMR_HEX_MESH}`)

disp("=== 3D AMR Hex Mesh ===")
disp(strcat("Dimension: ", num2str(mesh3d.dimension)))
disp(strcat("Elements: ", num2str(mesh3d.numElements)))
disp(strcat("Vertices: ", num2str(mesh3d.numVertices)))

# Create FE space and project initial temperature
fespace3d = mfem.fespace.createH1(mesh3d, 2, 1)
# Hot corner at origin — visible on 3 faces of the cube
u3d = mfem.gridfunction.fromFunction(fespace3d, @(x) 2.0*(norm(x) < 0.5) + 1.0*(norm(x) >= 0.5))

disp(strcat("DOFs: ", num2str(fespace3d.ndofs)))
disp(strcat("Initial max: ", num2str(max(u3d.getData())), ", min: ", num2str(min(u3d.getData()))))

# Time stepping: Backward Euler on the 3D mesh
kappa3d = 0.5
dt3d = 0.01
t3d_final = 0.1

M3d = mfem.bilinearForms.create(fespace3d)
M3d.addMassIntegrator()
M3d.assemble()
M3d.finalize()

dt_kappa3d = mfem.coefficients.constant(dt3d * kappa3d)
T3d = mfem.bilinearForms.create(fespace3d)
T3d.addMassIntegrator()
T3d.addDiffusionIntegrator(dt_kappa3d)
T3d.assemble()
T3d.finalize()

solver3d = mfem.linear.CGSolver.create()
solver3d.setRelTol(1e-10)
solver3d.setAbsTol(0.0)
solver3d.setMaxIter(200)
solver3d.setPrintLevel(0)
solver3d.setOperator(T3d.getSystemMatrixPointer())

rhs3d = mfem.gridfunction.GridFunction.create(fespace3d)

t3d = 0
for step = 1:ceil(t3d_final / dt3d)
    M3d.mult(u3d, rhs3d)
    solver3d.mult(rhs3d.getPointer(), u3d.getPointer())
    t3d = t3d + dt3d
end

disp(strcat("Reached t = ", num2str(t3d)))
disp(strcat("Final max: ", num2str(max(u3d.getData())), ", min: ", num2str(min(u3d.getData()))))

# Side-by-side: solution (left) and mesh (right)
subplot(1, 2, 1)
vtk.colormap("hot")
vtk.gridfunction(mesh3d, u3d)
title(strcat("3D Heat at t = ", num2str(t3d)))

subplot(1, 2, 2)
vtk.colormap("viridis")
vtk.edges("on")
vtk.mesh(mesh3d)
title("AMR Hex Mesh")

Key Concepts

Concept	Description
Method of Lines	Discretize space first (FEM), leaving ODEs in time. Separates spatial accuracy from temporal integration.
Stiff Problems	Problems where explicit methods need impractically small $\Delta t$ . Heat diffusion is inherently stiff.
Backward Euler	Implicit, L-stable, order 1. Solves $(M + \Delta t K) u^{n+1} = M u^n$ at each step. Unconditionally stable.
SDIRK Methods	Singly Diagonally Implicit Runge-Kutta. Higher-order implicit methods (ex16 default). Each stage reuses $T = M + \gamma \Delta t K$ .
ConductionOperator	In ex16.cpp, encapsulates $M$ , $K(u)$ , and solvers. Provides `ImplicitSolve()` for implicit ODE solvers.

Summary

We have walked through MFEM Example 16 step by step:

Load Mesh — Star mesh from MFEM mesh format
Refine — 2 levels of uniform refinement
Visualize — VTK rendering of the mesh
FE Space — H1 continuous, quadratic elements (order 2)
Initial Conditions — $u=2$ (hot center), $u=1$ (cool exterior)
Spatial Discretization — Mass matrix $M$ , stiffness matrix $K$
Time Integration — Backward Euler with CG solver: $(M + \Delta t K) u^{n+1} = M u^n$
Visualization — Evolved temperature field with VTK
3D Example — AMR hex mesh with non-conforming elements

Workbench

Clear

No variables in workbench

Next Steps

Finite Element Method

Learn FEM fundamentals: mesh generation, basis functions, assembly

VTK.js Visualization

Advanced 3D visualization for temperature fields and meshes

Equana