Heat on a Torus

The cells of a power diagram make a mesh for finite volumes: one unknown per cell, and one flux per facet. On a torus, the facets across the period are fluxes like the others. The cell complex gives everything the scheme needs:

• the volumes of the cells;
• for each facet, its area and its two seeds;
• the vector from one seed to the other. With a periodicity, this points to the image of the second seed that the facet faces.

We solve the heat equation on the unit torus,

$${\partial }_{t}u=\mathrm{\Delta }u,$$

with the classical two-point flux: the flux through the facet $f$ between cells $i$ and $j$ is its area over the distance between the two seeds, times the difference of the values,

$$|C_i|\frac{d{u}_i}{dt}=\sum _{f=(i,j)}\frac{|f|}{|p_j^f-p_i|}(u_j-u_i),$$

where $p_j^f$ is the image of the seed $j$ across $f$. On a Voronoi diagram the facet is orthogonal to $p_j^f-p_i$, which makes this flux consistent.

The mesh

The seeds are spread on the torus by a Lloyd quantization with the same periodicity. The positions it returns stay next to where they started, so % 1.0 brings them back into the square.

import numpy as np
import scipy.sparse as sp
import sdot

positions = np.asarray( sdot.quantization( sdot.Box( nb_dims = 2 ), 1500, seed = 0, periodicity = True ) ) % 1.0
pd = sdot.PowerDiagram( positions, boundaries = sdot.box_half_spaces( [ 0, 0 ], [ 1, 1 ] ), periodicity = True )
cx = pd.cell_complex()

The scheme

Each facet adds its transmissibility to the matrix $L$, at the rows and columns of its two cells. A facet between a cell and an image of itself exchanges nothing, so it is skipped. Time is discretized by implicit Euler, $(M+\mathrm{\Delta }tL)u^{k+1}=M{u}^k$, with $M$ the diagonal of the cell volumes.

n = len( positions )
a, b = cx.facets.sides.T
keep = a != b
a, b = a[ keep ], b[ keep ]
t = cx.facets.measures[ keep ] / np.linalg.norm( cx.facets.seed_vectors[ keep ], axis = 1 )
L = sp.coo_matrix( ( np.concatenate( [ t, t, -t, -t ] ),
                     ( np.concatenate( [ a, b, a, b ] ), np.concatenate( [ a, b, b, a ] ) ) ), shape = ( n, n ) ).tocsr()
M = sp.diags( cx.cells.measures )

Below, a hot disc and a cold band diffuse. The disc sits across the corner of the square, so it is a single disc on the torus. Each flux leaves one cell and enters another, so the total heat $\sum _i|C_i|u_i$ is conserved exactly: $0.029209$ at every step.

Open in a separate page ↗

Use the slider at the top of the viewer (or the play button) to go through time.

Convergence

The mode $\cos (2\pi x)$ decays exactly as $e^{-4{\pi }^{2}t}$. Up to $t=0.01$, we compare the decay computed by the scheme with the exact one, using the same implicit Euler steps in both so that only the error in space remains:

cells	decay of the scheme	exact	relative error
250	0.67706	0.67409	$4.4\cdot {10}^{-3}$
1000	0.67483	0.67409	$1.1\cdot {10}^{-3}$
4000	0.67427	0.67409	$2.8\cdot {10}^{-4}$

Each time the number of cells is multiplied by 4, the cell size $h$ is halved and the error is divided by 4: the scheme is of second order in $h$ on these centroidal cells.

The figures and the table are generated by docs/examples/finite_volumes.py.