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

2D Fluid-Grains dam break
=========================
This example simulates a 2-phase dam break, consisting of a column of fluid and
partially immersed grains
We follow the setup described in
X. Sun, M. Sakai, Y. Yamada, Three-dimensional simulation of a solid–liquid flow by the DEM–SPH method, J. Comput. Phys. (2013)

.. code-block:: python
  
  
  # Keywords
  # --------
  # PFEM, fluid-grains, dam break, 2D
  #
  # Description
  # -----------
  # The test demonstrates:
  # - How to perform a PFEM-DEM simulation with moving boundaries.
  # - A validation case for fluid-grains interaction.
  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, scontact, time_integration, pfem
  
  np.random.seed(0)
  

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_f = 0.05  # length of the fluid column
  H_f = 0.1  # height of the fluid column
  l = 0.2  # length of the box
  h = 1.0  # height of the box
  t_gate = 0.05 * L_f  # thickness of the gate
  h_gate = h  # height of the gate
  y_gate = 1e-8  # initial position of the gate
  u_gate = 0.68  # velocity of the gate
  L_g = L_f  # length of the grains column
  H_g = H_f / 3  # height of the grains column
  r_particles = 0.00135  # radius of the grains
  

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

.. code-block:: python
  
  geo_mesh_size = l / 10
  mesh_size = 0.1 * L_f
  smin = 1.4 * r_particles
  smax = 4 * smin
  dmax = 4 * smax
  alpha = 1.3
  

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

.. code-block:: python
  
  gmsh.model.add("ModelInit")
  gmsh.model.occ.addRectangle(l - L_f, 0, 0, L_f, H_f)
  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")
  
  
  def set_geo_mesh(y_gate=1e-4):
      name = gmsh.model.getCurrent()
      gmsh.model.setCurrent("ModelGeo")
      dimtag = gmsh.model.getEntities(2)
      gmsh.model.occ.remove(dimtag, True)
      box = gmsh.model.occ.addRectangle(0, 0, 0, l, h)
      gate = gmsh.model.occ.addRectangle(
          l - L_f - t_gate, y_gate, 0, t_gate, h_gate
      )  # gate
      gmsh.model.occ.cut([(2, box)], [(2, gate)])
      gmsh.model.occ.synchronize()
      gmsh.model.mesh.setSizeCallback(lambda *args: geo_mesh_size)
      gmsh.model.mesh.generate(2)
      geoEntities = gmsh.model.getEntities(1)
      gmsh.write(outputdir + "/lastGeo.msh")
      gmsh.model.setCurrent(name)
      return geoEntities
  
  
  geoEntities = set_geo_mesh()
  gmsh.model.addPhysicalGroup(1, [12], -1, "bottom")
  gmsh.model.addPhysicalGroup(1, [11], -1, "right")
  gmsh.model.addPhysicalGroup(1, [6], -1, "left")
  gmsh.model.addPhysicalGroup(1, [5, 8, 9], -1, "gate")
  

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)
  p0 = gmsh.model.addPhysicalGroup(1, [12], -1, "bottom")
  p1 = gmsh.model.addPhysicalGroup(1, [11], -1, "right")
  p2 = gmsh.model.addPhysicalGroup(1, [6], -1, "left")
  p3 = gmsh.model.addPhysicalGroup(1, [5, 8, 9], -1, "gate")
  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", 10 * mesh_size)
  gmsh.model.mesh.field.setNumber(sizeFieldConstant, "XMin", 0.0)
  gmsh.model.mesh.field.setNumber(sizeFieldConstant, "XMax", l)
  gmsh.model.mesh.field.setNumber(sizeFieldConstant, "YMin", 0.0)
  gmsh.model.mesh.field.setNumber(sizeFieldConstant, "YMax", h)
  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
  mu = 1e-3
  rho_bubble = 1
  sigma = 0.00728
  

DEM parameters
--------------

.. code-block:: python
  
  rhop = 2500
  friction = 0.2  # 0.2
  volume_corr = 0.8  # volume corection on the grains to account for 2D-3D discrepancies
  mass_grains = 0.2
  n_particles = int(mass_grains / (np.pi * r_particles**2 * rhop * L_g))
  

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

.. code-block:: python
  
  cfl = 0.2
  U = 2.5
  U_init = 0.3
  dt = smin / U * cfl
  t = 0
  settle_time = 0.06
  tEnd = 1.0 + settle_time
  

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

.. code-block:: python
  
  f = fluid.FluidProblem(2, g, mu, rho, advection=False)
  f.set_wall_boundary("bottom")
  f.set_wall_boundary("right")
  
  
  def gate_velocity(t):
      if t > settle_time:
          return u_gate
      return 0
  
  
  f.set_wall_boundary("gate", velocity=[0, gate_velocity(t)])
  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("left", velocity_x=0)
  f.set_open_boundary("freeSurface", pressure=0, viscous_flag=False)
  

DEM Problem
-----------

.. code-block:: python
  
  dem = scontact.ParticleProblem(2)
  dem.set_friction_coefficient(friction)
  
  # add outer box
  outer_body = dem.add_body((0, 0), 0, 0)
  x0 = np.array([0, 0])
  x1 = np.array([l, 0])
  x2 = np.array([l, h])
  x3 = np.array([0, h])
  dem.add_segment_to_body(x0, x1, outer_body, "bnd", material="wall")
  dem.add_segment_to_body(x1, x2, outer_body, "bnd", material="wall")
  dem.add_segment_to_body(x2, x3, outer_body, "bnd", material="wall")
  dem.add_segment_to_body(x3, x0, outer_body, "bnd", material="wall")
  dem.add_particle_to_body(x0, 0, outer_body)
  dem.add_particle_to_body(x1, 0, outer_body)
  dem.add_particle_to_body(x2, 0, outer_body)
  dem.add_particle_to_body(x3, 0, outer_body)
  
  # add gate
  gate_body = dem.add_body((0, 0), 0, 0)
  x0 = np.array([l - L_f - t_gate, y_gate])
  x1 = np.array([l - L_f, y_gate])
  x2 = np.array([l - L_f, h_gate])
  x3 = np.array([l - L_f - t_gate, h_gate])
  dem.add_segment_to_body(x0, x1, gate_body, "bnd", material="wall")
  dem.add_segment_to_body(x1, x2, gate_body, "bnd", material="wall")
  dem.add_segment_to_body(x2, x3, gate_body, "bnd", material="wall")
  dem.add_segment_to_body(x3, x0, gate_body, "bnd", material="wall")
  dem.add_particle_to_body(x0, 0, gate_body)
  dem.add_particle_to_body(x1, 0, gate_body)
  dem.add_particle_to_body(x2, 0, gate_body)
  dem.add_particle_to_body(x3, 0, gate_body)
  dem.body_velocity()[gate_body, 1] = u_gate
  
  

initial grains deposition
-------------------------

.. code-block:: python
  
  ratio = 1.2
  diameter_max = 2.0 * r_particles
  diameter_min = diameter_max / ratio
  u = np.random.power(1, n_particles)  # probability density function [0,1]
  d = (
      diameter_max * diameter_min / (diameter_max + u * (diameter_min - diameter_max))
  )  # random particle diameter following the given "u" law
  radii = d / 2
  boundaries = np.array([[l - L_g, 0], [l, 0], [l, 2 * H_g], [l - L_g, 2 * H_g]])
  xp, radii = scontact.fastdeposit_2d(boundaries, radii, outputdir + "/deposit.svg")
  eps = 1e-4 * diameter_max
  for x, r in zip(xp, radii):
      dem.add_particle(x + eps, r, np.pi * r**2 * rhop)
  

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

.. code-block:: python
  
  i = 0
  outf = 10
  gmsh.option.setNumber("General.Verbosity", 0)
  mass = np.pi * dem.r() ** 2 * rhop
  
  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 * smin,
          usePreviousMesh=True,
      )
      oparamcoord = oparamcoord.reshape((-1, 3))
      gmsh.write(outputdir + "/lastMesh.msh")
      ordered_node_tags = pfem.prepareMeshForMigflow(
          i, alphaDomainTag, f, nodetag, oelemtag, oparamcoord
      )
  
      f.set_particles(
          dem.delassus() * volume_corr,
          dem.volume() * volume_corr,
          dem.position(),
          dem.velocity(),
          dem.omega() * 0,
          dem.contact_forces(),
      )
      pfem.applySurfaceTension(f, sigma, False, g, rho_bubble)
  
      if i % outf == 0:  # or t > 0.3:
          f.write_mig(outputdir, t)
          dem.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)
  
      # move gate
      if t < settle_time:
          dem.body_velocity()[gate_body, 1] = 0
      else:
          dem.body_velocity()[gate_body, 1] = u_gate
  
      drag = f.compute_node_force()
      external_forces = drag + g * dem.r() ** 2 * np.pi * rhop
      time_integration._advance_particles(dem, external_forces, dt, 1, 1e-4 * r_particles)
  
      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
      dx = dx.reshape((-1,))
  
      # update gate position and geo mesh
      if y_gate > h / 2:
          dem.body_velocity()[gate_body, 1] = 0
      y_gate += dem.body_velocity()[gate_body, 1] * dt
      geoEntities = set_geo_mesh(y_gate)
  
      # advect nodes and project if needed
      gmsh.model.mesh.advectMeshNodes(
          2,
          alphaDomainTag,
          alphaBoundaryTag,
          "ModelGeo",
          ordered_node_tags,
          dx,
          0.01 * smin,
      )
  
      # set new coordinates
      _, newCoords, _ = gmsh.model.mesh.getNodes(2, alphaDomainTag)
      f.set_coordinates(newCoords)
  
      # update time step
      max_U = np.max([U_init, np.max(np.linalg.norm(f.velocity(), axis=1))])
      dt = smin / max_U * cfl
  
      t += dt
      i += 1
  

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

 python3 -m migflow.plot.migplot output --actors fluid particles --fluid-field velocity --show-edges 1 --bounds 0 0.2 0 0.2

