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

2D lid-driven cavity flow with PFEM
===================================
This example simulates a 2D lid-driven cavity flow using the Particle Finite Element Method (PFEM).


Keywords
--------
PFEM, lid driven cavity, 2D, fluid


Description
-----------
The test demonstrates:
- How to perform a PFEM simulation in a fixed domain (without free surfaces)

.. code-block:: python
  
  import os, sys, shutil
  import numpy as np
  

Gmsh pfem branch
----------------
This test requires the alphashapes branch of GMSH to use the Particle Finite Element (PFEM) module.
PFEM is a lagrangian method where the mesh is regenerated at each time step, making it suitable for problems with large deformations and free surfaces.

.. code-block:: python
  
  sys.path.insert(0, os.environ["GMSH_PFEM_DIR"])
  
  import gmsh
  from migflow import fluid, pfem
  

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)
  
  

Geometrical parameters and mesh generation
------------------------------------------

.. code-block:: python
  
  l_box = 1.0
  h_box = 1.0
  

Mesh parameters
---------------

.. code-block:: python
  
  geo_mesh_size = l_box / 10
  mesh_size = l_box / 40
  alpha = 10.0
  

Initial fluid mesh
------------------

.. code-block:: python
  
  gmsh.model.add("ModelInit")
  gmsh.model.occ.addRectangle(0, 0, 0, l_box, h_box)
  gmsh.model.occ.synchronize()
  gmsh.model.mesh.setSizeCallback(lambda *args: mesh_size)
  gmsh.model.mesh.generate(2)
  nodeTags, coords, _ = gmsh.model.mesh.getNodes()
  

Solid domain
------------

.. code-block:: python
  
  gmsh.model.add("ModelGeo")
  gmsh.model.occ.addRectangle(0, 0, 0, l_box, h_box)
  gmsh.model.occ.synchronize()
  gmsh.model.mesh.setSizeCallback(lambda *args: mesh_size)
  gmsh.model.mesh.generate(2)
  nodeTags, coords, _ = gmsh.model.mesh.getNodes()
  gmsh.model.addPhysicalGroup(1, [1], -1, "bottom")
  gmsh.model.addPhysicalGroup(1, [2], -1, "right")
  gmsh.model.addPhysicalGroup(1, [3], -1, "top")
  gmsh.model.addPhysicalGroup(1, [4], -1, "left")
  geoEntities = gmsh.model.getEntities(1)

PFEM mesh and size fields
-------------------------

.. code-block:: python
  
  gmsh.model.add("ModelFluid")
  alphaDomainTag = gmsh.model.addDiscreteEntity(2, -1, [])
  for dim, tag in geoEntities:
      gmsh.model.addDiscreteEntity(dim, tag, [])
  alphaBoundaryTag = gmsh.model.addDiscreteEntity(1, -1, [])
  gmsh.model.mesh.addNodes(2, alphaDomainTag, nodeTags, coords)
  gmsh.model.addPhysicalGroup(1, [1], -1, "bottom")
  gmsh.model.addPhysicalGroup(1, [2], -1, "right")
  gmsh.model.addPhysicalGroup(1, [3], -1, "top")
  gmsh.model.addPhysicalGroup(1, [4], -1, "left")
  gmsh.model.addPhysicalGroup(1, [alphaBoundaryTag], -1, "freeSurface")
  gmsh.model.addPhysicalGroup(2, [alphaDomainTag], -1, "domain")
  

Size fields
-----------

.. code-block:: python
  
  sizeFieldConstant = gmsh.model.mesh.field.add("Box")
  gmsh.model.mesh.field.setNumber(sizeFieldConstant, "VIn", mesh_size)
  gmsh.model.mesh.field.setNumber(sizeFieldConstant, "VOut", mesh_size)
  gmsh.model.mesh.field.setNumber(sizeFieldConstant, "XMin", 0.0)
  gmsh.model.mesh.field.setNumber(sizeFieldConstant, "XMax", l_box)
  gmsh.model.mesh.field.setNumber(sizeFieldConstant, "YMin", 0.0)
  gmsh.model.mesh.field.setNumber(sizeFieldConstant, "YMax", h_box)
  gmsh.model.mesh.field.setNumber(sizeFieldConstant, "Thickness", 0.001)
  

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

.. code-block:: python
  
  g = np.array([0.0, 0.0])
  rho = 1000
  mu = 1e-3
  Re = 1000
  v_top = Re * mu / (rho * l_box)

Time parameters
---------------

.. code-block:: python
  
  cfl = 0.5
  U = v_top
  U_init = v_top
  dt = mesh_size / U * cfl
  t = 0
  tEnd = 25000.0
  

Fluid problem
-------------

.. code-block:: python
  
  f = fluid.FluidProblem(2, g, mu, rho, advection=False)
  f.set_wall_boundary("bottom")
  f.set_wall_boundary("right")
  f.set_wall_boundary("top", velocity=[v_top, 0])
  f.set_wall_boundary("left")
  f.set_strong_boundary("bottom", velocity=[0, 0])
  f.set_strong_boundary("right", velocity=[0, 0])
  f.set_strong_boundary("top", velocity=[v_top, 0])
  f.set_strong_boundary("left", velocity=[0, 0])
  

Simulation Loop
---------------
Time integration of coupled fluid–particle motion.

.. code-block:: python
  
  i = 0
  outf = 25
  gmsh.option.setNumber("General.Verbosity", 0)
  
  while t < tEnd:
      print(f"{i:4d}, {t:.6g}/{tEnd:.6g}, {dt:.6g}")
      # Update PFEM mesh
      nodetag, oelemtag, oparamcoord, _ = gmsh.model.mesh.computeAlphaShape(
          2,
          alphaDomainTag,
          alphaBoundaryTag,
          "ModelGeo",
          alpha,
          sizeFieldConstant,
          sizeFieldConstant,
          boundaryTolerance=0.01 * mesh_size,
          usePreviousMesh=False,
      )
      oparamcoord = oparamcoord.reshape((-1, 3))
      gmsh.write(outputdir + "/lastMesh.msh")
      ordered_node_tags = pfem.prepareMeshForMigflow(
          i, alphaDomainTag, f, nodetag, oelemtag, oparamcoord
      )
  
      if i % outf == 0:  # or t > 0.3:
          f.write_mig(outputdir, t)
  
      # nodes velocity, be aware that its dimension is (n_nodes, 3) not (n_nodes, 2)
      u_old = np.zeros_like(f.coordinates())
      u_old[:, :2] = f.velocity()
  
      f.implicit_euler(dt)
  
      dx = np.zeros((f.coordinates().shape[0], 3))
      u = np.zeros_like(dx)
      u[:, :2] = f.velocity()
      if i == 0:
          u_old = u
      dx = u * dt + 0.5 * (u - u_old) / dt * dt**2
  
      # Do not move the top boundary nodes
      top_nodes = np.unique(f.mesh_boundaries()[b"top"].flatten())
      dx[top_nodes, :] = 0.0
  
      # advect nodes and project if needed
      gmsh.model.mesh.advectMeshNodes(
          2,
          alphaDomainTag,
          alphaBoundaryTag,
          "ModelGeo",
          ordered_node_tags,
          dx.flatten(),
          0.01 * mesh_size,
      )
  
      # set new coordinates
      _, newCoords, _ = gmsh.model.mesh.getNodes(2, alphaDomainTag)
      f.set_coordinates(newCoords)
  
      t += dt
      i += 1
  

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

 python3 -m migflow.plot.migplot output --actors fluid --fluid-field velocity --fluid-vmin 0 --fluid-vmax 0.001 --show-edges 1

