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

2D dam break simulation with PFEM
=================================
This example simulates a 2D dam break flow using the Particle Finite Element Method (PFEM).
We use the setup from Cremonesi et al., "A State of the Art Review of the Particle Finite Element Method (PFEM)", 2020.

.. code-block:: python
  
  
  # Keywords
  # --------
  # PFEM, dam break, 2D, fluid
  #
  # Description
  # -----------
  # The test demonstrates:
  # - How to perform a PFEM simulation with free surface dynamics.
  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 = 0.146
  b = 0.175
  
  l_box = 4 * L
  h_box = 8 * L
  
  l_fluid = L
  h_fluid = 2 * L

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

.. code-block:: python
  
  geo_mesh_size = l_fluid / 10
  mesh_size = l_fluid / 10
  smin = 0.3 * mesh_size
  smax = mesh_size
  dmax = 5 * smax
  alpha = 1.3
  

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

.. code-block:: python
  
  gmsh.model.add("ModelInit")
  rect = gmsh.model.occ.addRectangle(0, 0, 0, l_fluid, h_fluid)
  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: geo_mesh_size)
  gmsh.model.mesh.generate(2)
  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)
  
  sizeFieldDistFS = gmsh.model.mesh.field.add("AlphaShapeDistance")
  gmsh.model.mesh.field.setNumber(sizeFieldDistFS, "Tag", alphaBoundaryTag)
  gmsh.model.mesh.field.setNumber(sizeFieldDistFS, "SamplingLength", 0.1 * smin)
  sizeFieldRefine = gmsh.model.mesh.field.add("Threshold")
  gmsh.model.mesh.field.setNumber(sizeFieldRefine, "InField", sizeFieldDistFS)
  gmsh.model.mesh.field.setNumber(sizeFieldRefine, "SizeMin", smin)
  gmsh.model.mesh.field.setNumber(sizeFieldRefine, "SizeMax", smax)
  gmsh.model.mesh.field.setNumber(sizeFieldRefine, "DistMin", 0.0)
  gmsh.model.mesh.field.setNumber(sizeFieldRefine, "DistMax", dmax)
  gmsh.model.mesh.computeAlphaShape(
      2,
      alphaDomainTag,
      alphaBoundaryTag,
      "ModelGeo",
      alpha,
      sizeFieldConstant,
      sizeFieldRefine,
      False,
  )
  

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

.. code-block:: python
  
  g = np.array([0.0, -9.81])
  rho = 1000
  rho_bubble = 1
  mu = 1e-3
  sigma = 0.0
  

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

.. code-block:: python
  
  cfl = 0.3
  U = 10.0
  U_init = U
  dt = mesh_size / U * cfl
  t = 0
  tEnd = 1
  

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")
  f.set_wall_boundary("left")
  f.set_strong_boundary("bottom", velocity_y=0)
  f.set_strong_boundary("right", velocity_x=0)
  f.set_strong_boundary("top", velocity_y=0)
  f.set_strong_boundary("left", velocity_x=0)
  f.set_open_boundary("freeSurface", pressure=0, viscous_flag=False)
  

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

.. code-block:: python
  
  i = 0
  outf = 10
  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,
          sizeFieldRefine,
          sizeFieldRefine,
          boundaryTolerance=0.01 * mesh_size,
          usePreviousMesh=True,
      )
      oparamcoord = oparamcoord.reshape((-1, 3))
      gmsh.write(outputdir + "/lastMesh.msh")
      ordered_node_tags = pfem.prepareMeshForMigflow(
          i, alphaDomainTag, f, nodetag, oelemtag, oparamcoord
      )
  
      pfem.applySurfaceTension_v0(
          f, sigma, bubbleCondition=True, g=g, rho_bubble=rho_bubble
      )
  
      if i % outf == 0:
          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
  
      # 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.0 --fluid-vmax 3.0 --show-edges 1

