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

2D Rayleigh-Benard instability test case
========================================
This test case simulates a 2D Rayleigh-Benard instability using the Finite Element Method (FEM) in a periodic domain.

Keywords
--------
FEM, fluid, convection, 2D, heat transfer


Description
-----------
The test demonstrates:
- How to perform a coupled fluid and heat transfer simulation.
- How to set up boundary conditions for temperature and velocity fields.
- How to set up non constant density.

.. code-block:: python
  
  
  import gmsh
  import numpy as np
  from migflow import fluid, 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)

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.

.. code-block:: python
  
  w = 1
  h = 0.25
  lc = 0.005
  gmsh.model.add("box")
  gmsh.model.occ.add_rectangle(-w / 2, 0, 0, w, h)
  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([-w / 2, 0], [w / 2, 0]), name="Bottom")
  gmsh.model.add_physical_group(1, get_line([-w / 2, h], [w / 2, h]), name="Top")
  gmsh.model.add_physical_group(1, get_line([-w / 2, 0], [-w / 2, h]), name="Left")
  gmsh.model.add_physical_group(1, get_line([w / 2, 0], [w / 2, h]), name="Right")
  gmsh.model.add_physical_group(2, [1], name="domain")
  transform = np.eye(4, 4, dtype=np.float64)
  transform[0, -1] = -w
  gmsh.model.mesh.set_periodic(
      1,
      get_line([-w / 2, 0], [-w / 2, h]),
      get_line([w / 2, 0], [w / 2, h]),
      transform.reshape(-1),
  )
  gmsh.model.mesh.set_size_callback(lambda dim, tag, x, y, z, mesh_size: lc)
  gmsh.model.mesh.generate(2)
  

Fluid Problem
-------------
The fluid is described by its dimension, the external volume force applied, its dynamic viscosity and its density.
Pardiso solver is chosen as a direct solver.
As we target a natural convection, the density depends on the temperature. To do so, the density must be seen as a field and not as a constant.

.. code-block:: python
  
  g = np.array([0, -9.81])
  rho = 1000
  nu = 1e-6
  mu = rho * nu
  beta = 2e-4
  f = fluid.FluidProblem(2, g, mu, density_element=b"triangle_p1")
  f.load_msh(None)
  f.set_wall_boundary("Top", velocity=[0, 0])
  f.set_wall_boundary("Bottom", velocity=[0, 0])
  f.set_mean_pressure(0)
  f.interpolate(
      velocity_y=lambda x: 1e-5
      * np.sin(np.pi * x[:, 0] / w)
      * np.sin(np.pi * x[:, 1] / h)
  )

Temperature Problem
-------------------
The advection-diffusion problem is described by its dimension, its source term, its conductivity and its capacity.

.. code-block:: python
  
  cp = rho * 1.0
  k = 1e-2
  Tt = -10  # Top temperature
  Tb = 10  # Bottom temperature
  c = advdiff.AdvDiffProblem(
      2, 0, k, cp, velocity_element=b"triangle_p1"
  )
  c.load_msh(None)
  c.set_dirichlet_boundary("Top", solution=Tt)
  c.set_dirichlet_boundary("Bottom", solution=Tb)

Time integration
----------------
The numerical parameters is given and the initial conditions are written in the output directory.

.. code-block:: python
  
  dt = 1
  tEnd = 500
  outf = 1
  t = 0
  ii = 0
  
  
  def get_fields(fluid, temp):
      y = fluid.coordinates_fields()[fluid.pressure_index()][:, 1]
      y = y.reshape(-1, 1)
      p1_element = fluid.get_p1_element()
      fields = fluid.get_default_export()
      fields["temperature"] = (temp.solution().reshape(-1, 1), p1_element)
      fields["density"] = (fluid.density().reshape(-1, 1), p1_element)
      return fields
      # return {
      #     "pressure": (fluid.pressure(), p1_element),
      # # "velocity": (fluid.velocity(), p1_element),
      # # "density": (fluid.density(), p1_element),
      # # "temperature": (temp.solution().reshape(-1, 1), p1_element),
      #     # "porosity": (fluid.porosity(), p1_element),
      # }
  
  
  while t < tEnd:
      # Problem coupling
      # ----------------
      f.density()[:] = rho * (1 - beta * (c.solution().reshape(-1, 1) - (Tt + Tb) / 2))
      c.velocity()[:] = f.velocity()[:]
      if ii % outf == 0:
          print(f"ii : {ii} ----- t : {t}/{tEnd}")
          f.write_mig(outputdir, t, fields=get_fields(f, c))
      # Solve
      # -----
      c.implicit_euler(dt)
      f.implicit_euler(dt)
      t += dt
      ii += 1
  

Plot
----
.. code-block:: shell

 python3 -m migflow.plot.migplot output --actors fluid --fluid-field density --fluid-vmin 998 --fluid-vmax 1002 --element-type triangle_p1

