Download this testcase.

2D Laplacian test case

This test case simulates a 2D Laplacian using the Finite Element Method (FEM).

Keywords

FEM, Laplacian

Description

The test demonstrates: - How to perform a simple Laplacian simulation.

import gmsh
import numpy as np
from migflow import advdiff
import shutil, os, sys

Output Directory

Create a clean output directory for simulation results.

outputdir = "output" if len(sys.argv) < 2 else sys.argv[1]
shutil.rmtree(outputdir, ignore_errors=True)
os.makedirs(outputdir)
VISUALIZE = False

Laplacian Problem

def solve_laplacian(lc):
    # Mesh generation
    # The domain is created through the python GMSH api is loaded to generate the mesh and extract the boundaries.
    # The refinement is given by the function set_size_call_back. In this example, a constant mesh size is chosen.
    r = 1.0
    gmsh.model.add("disc")
    gmsh.model.occ.add_disk(0, 0, 0, r, r)
    gmsh.model.occ.synchronize()

    def get_line(x0, x1, eps=1e-6):
        r = gmsh.model.get_entities_in_bounding_box(
            x0[0] - eps, x0[1] - eps, -eps, x1[0] + eps, x1[1] + eps, eps, 1
        )
        return list(tag for _, tag in r)

    gmsh.model.add_physical_group(1, get_line([-r, -r], [r, r]), name="boundary")
    gmsh.model.add_physical_group(2, [1], name="domain")
    gmsh.model.mesh.set_size_callback(lambda dim, tag, x, y, z, mesh_size: lc)
    gmsh.model.mesh.generate(2)
    # The Laplacian problem is described by its dimension, its source term.
    problem = advdiff.AdvDiffProblem(
        2, f=1, k=1, cp=0, temporal=False, advection=False
    )
    problem.load_msh(None)
    problem.set_strong_boundary("boundary", solution=0)
    # The numerical parameters is given and the initial conditions are written in the output directory.
    dt = 1
    print(f"Solving Laplacian with lc={lc:.3e}")
    problem.implicit_euler(dt)
    problem.write_mig(outputdir, lc)
    x = problem.coordinates_fields()
    x, y = x[:, 0], x[:, 1]
    u_ref = u_analytic(x, y)
    nodes_vol = problem.node_volume()
    error_nodes = problem.solution() - u_ref
    error = np.sqrt(np.sum((error_nodes**2) * nodes_vol)) / np.sqrt(
        np.sum((u_ref**2) * nodes_vol)
    )
    return error


lc = np.logspace(-2.0, -1, 6)


def u_analytic(x, y):
    r = np.sqrt(x**2 + y**2)
    return (1 - r**2) / 4


error = np.empty(len(lc))
for i, lci in enumerate(lc):
    error[i] = solve_laplacian(lci)

Plotting the convergence

convergence_rate = np.polyfit(np.log(lc), np.log(error), 1)[0]
print(f"Convergence rate: {convergence_rate:.2f}")
if VISUALIZE:
    import matplotlib.pyplot as plt

    fig, ax = plt.subplots()
    ax.loglog(lc, error, "-o", label="L2 Error")
    ax.loglog(lc, (lc**2) * error[0] / (lc[0] ** 2), "--", label="O(h^2)")
    ax.set_xlabel("Mesh size (lc)")
    ax.set_ylabel("L2 Error")
    ax.legend()
    ax.set_title("2D Laplacian Convergence")
    ax.grid(True)
    plt.show()

Plot

python3 -m migflow.plot.migplot output --actors fluid --problem-type advection-diffusion --fluid-field solution --fluid-vmin 0 --fluid-vmax 0.25 --show-edges 1