Download :download:`this testcase <2d_laplacian.py>`.

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.

.. code-block:: python
  
  
  import gmsh
  import numpy as np
  from migflow import advdiff
  import shutil, os, sys
  

Output Directory
----------------
Create a clean output directory for simulation results.

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

Laplacian Problem
-----------------

.. code-block:: python
  
  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
------------------------

.. code-block:: python
  
  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
----
.. code-block:: shell

 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


.. code-block:: python
  
  
