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

Advection–Diffusion or Fluid Transport in a 2D Periodic Channel
===============================================================

This example illustrates how to solve a 2D periodic channel problem
using either the advection–diffusion equation (ADV_DIFF = 1)
or the incompressible Navier–Stokes equations (ADV_DIFF = 0).

The domain is a rectangular channel with periodic lateral boundaries
and Dirichlet conditions at the top and bottom.

Keywords
--------
Advection-Diffusion, Periodic Boundaries, Mesh Adaptation

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

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)
  

Problem Configuration
---------------------

.. code-block:: python
  
  ADV_DIFF = 1  # 1 → advection–diffusion, 0 → fluid simulation
  
  

Helper Class for Mesh Sizing
----------------------------

.. code-block:: python
  
  class MeshSize:
      """Callable mesh-size function for gmsh."""
  
      def __init__(self, lc):
          self._lc = lc
  
      def __call__(self, dim, tag, x, y, z, lc):
          return self._lc
  
  

Geometrical parameters
----------------------

.. code-block:: python
  
  width, height = 4.0, 1.0
  origin = [-width / 2, -height / 2]
  lc = 0.1
  mesh_size = MeshSize(lc)
  

Physical Parameters
-------------------

.. code-block:: python
  
  g = np.array([0, 0])  # gravity (ignored for adv-diff)
  rho = 1.204  # nitrogen density [kg/m³]
  nu = 1.8e-5  # kinematic viscosity [m²/s]
  mu = rho * nu  # dynamic viscosity [kg/(m·s)]
  alpha = 3e-3  # thermal dilatation coefficient
  k = 2e-2  # thermal conductivity [W/mK]
  cp = 1010e-3  # heat capacity [kJ/kg·K]
  
  

Generate Initial Mesh
---------------------

.. code-block:: python
  
  def gen_mesh(width, height, mesh_size, origin=np.array([0, 0])):
      """Generate a rectangular periodic mesh in gmsh.
  
      Parameters
      ----------
      width : float
          Domain width.
      height : float
          Domain height.
      mesh_size : MeshSize
          Callable object controlling element size.
      origin : np.ndarray, optional
          Origin of the rectangle.
      """
      origin = np.asarray(origin)
      gmsh.model.add("box")
      gmsh.model.occ.add_rectangle(origin[0], origin[1], 0, width, height)
      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 [tag for _, tag in r]
  
      h, w = height, width
      l_tag = get_line(origin + [0, 0], origin + [0, h])
      r_tag = get_line(origin + [w, 0], origin + [w, h])
      # Define periodicity (right ↔ left)
      affine_matrix = np.eye(4, 4, dtype=np.float64)
      affine_matrix[0, 3] = width
      gmsh.model.mesh.set_periodic(1, r_tag, l_tag, affine_matrix.reshape(-1))
      gmsh.model.occ.synchronize()
      # Define physical groups
      gmsh.model.add_physical_group(
          1, get_line(origin + [0, 0], origin + [w, 0]), name="Bottom"
      )
      gmsh.model.add_physical_group(
          1, get_line(origin + [0, h], origin + [w, h]), name="Top"
      )
      gmsh.model.add_physical_group(2, [1], name="domain")
      gmsh.model.mesh.set_size_callback(mesh_size)
  
  
  gmsh.initialize()
  gen_mesh(width, height, mesh_size, origin)
  gmsh.model.mesh.generate(2)
  
  

Problem Definition
------------------

.. code-block:: python
  
  def define_problem(adv_diff):
      """Create either an advection–diffusion or fluid problem."""
      if adv_diff:
          # --- Advection–Diffusion problem ---
          c = advdiff.AdvDiffProblem(2, 0, 0, 1.0, explicit=True, cg=False)
          c.load_msh(None)
          c.set_weak_boundary("Bottom")
          c.set_weak_boundary("Top")
          c.velocity()[:, 0] = 1.0  # uniform flow to the right
          x = c.coordinates_fields()
          c.solution()[x[:, 0] > -1e-8] = 1  # initialize concentration at inflow
          return c
      else:
          # --- Fluid problem (reference) ---
          f = fluid.FluidProblem(2, g, [mu, mu], [rho, rho])
          f.load_msh(None)
          f.set_mean_pressure(0)
          f.set_wall_boundary("Bottom")
          f.set_wall_boundary("Top")
          s0 = np.zeros_like(f.solution())
          s0[f.velocity_index()[:, 0]] = 1.0  # initial velocity in x-direction
          f.solution()[:] = s0
          c0 = np.zeros_like(f.coordinates()[:, 0])
          xc0 = f.coordinates()
          c0[xc0[:, 0] > -1e-8] = 1
          f.set_concentration_cg(c0)
          return f
  
  
  p = define_problem(ADV_DIFF)
  

Mesh Adaptation
---------------

.. code-block:: python
  
  mesh_size = MeshSize(0.05)
  gen_mesh(width, height, mesh_size, origin)
  gmsh.model.mesh.clear()
  gmsh.model.mesh.generate(2)
  if ADV_DIFF:
      p.adapt_mesh()
      p.velocity()[:, 0] = 1.0
  else:
      p.adapt_mesh()
  

Simulation Parameters
---------------------

.. code-block:: python
  
  dt = 0.01  # time step [s]
  tEnd = 1.0  # final time [s]
  outf = 1  # output frequency
  

Time Integration Loop
---------------------

.. code-block:: python
  
  t, i = 0, 0
  while t < tEnd:
      print(f"{i:4d}, {t:.6g}/{tEnd:.6g}")
      if i % outf == 0:
          p.write_mig(outputdir, t)
      if ADV_DIFF:
          p.explicit_euler(dt, limiter=True)
      else:
          p.advance_concentration(dt)
      t += dt
      i += 1
  

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 --element-type triangle_p1dg


.. code-block:: python
  
  
