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

Darcy Flow Through a Square Arrangement of Grains
=================================================
This example simulates Darcy flow through a square lattice arrangement
of circular grains in two dimensions. A fluid velocity is imposed across
the periodic domain and the resulting flow field through the porous medium
is computed.

Keywords
--------
FEM, DEM, Darcy, Porous Media

.. code-block:: python
  
  import os
  
  import gmsh
  import shutil
  import numpy as np
  from migflow import fluid, scontact
  
  np.random.seed(42)
  
  outputdir = "output"
  shutil.rmtree(outputdir, ignore_errors=True)
  os.makedirs(outputdir)
  

Physical parameters
-------------------

.. code-block:: python
  
  g = np.array([0, 0])  # gravity
  rho = 1000  # fluid density
  rhop = 1500  # grains density
  mu = 1  # dynamic viscosity
  compacity = 0.5  # Compactness
  porosity = 1 - compacity  # porosity
  V = 0.00001 / porosity  # fluid velocity
  

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

.. code-block:: python
  
  r = 40e-3  # grains radius
  L = 1  # domain width
  h = r / 2  # mesh size
  # h = 1/20                                    # mesh size
  origin = np.array([-L / 2, -L / 2])  # mesh origin
  
  gmsh.initialize()
  gmsh.model.add("mesh")
  gmsh.model.occ.add_rectangle(origin[0], origin[1], 0, 2 * L, 2 * L)
  gmsh.model.occ.synchronize()
  
  
  def get_line(x0, x1, eps=1e-6):
      bbentities = 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 dim, tag in bbentities)
  
  
  gmsh.model.add_physical_group(
      1, get_line([origin[0], origin[1]], [origin[0] + 2 * L, origin[1]]), name="Bottom"
  )
  gmsh.model.add_physical_group(
      1,
      get_line([origin[0], origin[1] + 2 * L], [origin[0] + 2 * L, origin[1] + 2 * L]),
      name="Top",
  )
  gmsh.model.add_physical_group(
      1, get_line([origin[0], origin[1]], [origin[0], origin[1] + 2 * L]), name="Left"
  )
  gmsh.model.add_physical_group(
      1,
      get_line([origin[0] + 2 * L, origin[1]], [origin[0] + 2 * L, origin[1] + 2 * L]),
      name="Right",
  )
  gmsh.model.add_physical_group(2, [1], name="domain")
  gmsh.model.mesh.set_size_callback(lambda dim, tag, x, y, z, lc: h)
  gmsh.model.mesh.generate(2)
  

Fluid Problem
-------------

.. code-block:: python
  
  f = fluid.FluidProblem(2, g, mu, rho)
  f.load_msh(None)
  f.set_open_boundary("Left", velocity=[V, 0])
  f.set_open_boundary("Right", pressure=0.0)
  f.set_wall_boundary("Bottom", velocity=[0, 0])
  f.set_wall_boundary("Top", velocity=[0, 0])
  f.interpolate(velocity_x=V)
  

Particle Problem
----------------

.. code-block:: python
  
  p = scontact.ParticleProblem(2)
  
  x_g = np.arange(0 + 2 * r, 2 * L - 2 * r, 3 * r) + origin[0]
  y_g = np.arange(0 + 2 * r, 2 * L - 2 * r, 3 * r) + origin[1]
  x_g, y_g = np.meshgrid(x_g, y_g)
  
  x_g = x_g.flat
  y_g = y_g.flat
  polygons = np.array([x_g, y_g]).T.reshape(-1, 1, 2)
  grains_size = np.random.uniform(0.6 * r, 1 * r, polygons.shape[0])
  polygons = np.repeat(polygons, 4, axis=1)
  polygons[:, 0, :] += np.array([-grains_size, -grains_size]).T
  polygons[:, 1, :] += np.array([+grains_size, -grains_size]).T
  polygons[:, 2, :] += np.array([+grains_size, +grains_size]).T
  polygons[:, 3, :] += np.array([-grains_size, +grains_size]).T
  
  volume = (2 * grains_size) ** 2
  density = np.full_like(volume, rho)
  
  for poly, ri in zip(polygons, grains_size):
      x_center = np.mean(poly, axis=0)
      p.add_particle(x_center, ri, ri**2 * 4 * rhop)
  
  csr_ptr, csr_element, csr_xi, csr_surface = f.get_mesh_polygons_overlap(polygons)
  
  if False:
      import matplotlib as mpl
      import matplotlib.pyplot as plt
      from matplotlib.patches import Circle
  
      fig, ax = plt.subplots()
      ax.set_aspect("equal")
      triangulation = mpl.tri.Triangulation(
          f.coordinates()[:, 0], f.coordinates()[:, 1], f.elements()
      )
      ax.triplot(triangulation, linewidth=0.2, color="k")
      phi0 = 1 - csr_xi[..., 0] - csr_xi[..., 1]
      phi1 = csr_xi[..., 0]
      phi2 = csr_xi[..., 1]
      phi = np.array([phi0, phi1, phi2]).T  # qp, node
      xelt = f.coordinates()[f.elements()[csr_element], :2]  # qp, node, dim
      xqp = np.einsum("ijk, ij -> ik", xelt, phi)  # qp, dim
      for i, x_i in enumerate(xqp):
          p_area = csr_surface[i]
          r = np.sqrt(p_area / np.pi)
          c = Circle(x_i[:2], r / 2, facecolor="black")
          ax.add_patch(c)
      ax.set_xlim(origin[0], origin[0] + 2 * L)
      ax.set_ylim(origin[1], origin[1] + 2 * L)
      for poly in polygons:
          polygon = plt.Polygon(poly, fill=None, edgecolor="r")
          ax.add_patch(polygon)
      plt.show()
  

Numerical parameters
--------------------

.. code-block:: python
  
  t = 0
  i = 0
  dt = 1e-3
  tEnd = 1.0
  outf = 10
  

Time Integration
----------------

.. code-block:: python
  
  
  
  def get_fields(fluid: fluid.FluidProblem):
      """Return derived output fields for visualization."""
      p1_element = fluid.get_p1_element()
      grad_v = fluid.fields_gradient()[:, :2, :]
      vorticity = (grad_v[:, 1, 0] - grad_v[:, 0, 1]).reshape(-1, 1) * fluid.porosity()
      return {
          "pressure": (fluid.pressure(), p1_element),
          "velocity": (fluid.velocity(), p1_element),
          "porosity": (fluid.porosity(), p1_element),
          "vorticity": (vorticity, p1_element),
      }
  
  
  csr_velocity = np.zeros_like(csr_xi)
  csr_gamma = 1e2 * np.ones(csr_xi.shape[0])
  csr_contact = np.zeros_like(csr_xi)
  f.set_bodies(
      density,
      volume,
      csr_ptr,
      csr_element,
      csr_xi,
      csr_surface,
      csr_velocity,
      csr_gamma,
      csr_contact,
  )
  
  
  while t < tEnd:
      if i % outf == 0:
          print(f"t/tEnd={t:1.1f}/{tEnd:1.1f}, i={i:1d}")
          f.write_mig(outputdir, t, get_fields(f))
          p.write_mig(outputdir, t)
      csr_force = f.get_bodies_csr_force().copy()
      f.set_bodies(
          density,
          volume,
          csr_ptr,
          csr_element,
          csr_xi,
          csr_surface,
          csr_velocity,
          csr_gamma,
          -csr_force,
      )
      f.implicit_euler(dt)
      t += dt
      i += 1
  

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

 python3 -m migflow.plot.migplot output --actors fluid particles --fluid-field porosity

