Tutorial 1: Create 2D dam break

In this tutorial, we will create the initial state for a two-phase gas-liquid dam break simulation using GeoParticle. The expected result is as follows:

What we will learn in this tutorial:

• Creating 2D geometries: ThickRectangle and FilledRectangle

• Performing operations: subtract and shift

• Accessing Geometry class member: xs, ys, and flatten_coords

Specify the parameters

First, we need to import the necessary libraries and set up the simulation parameters:

import numpy as np
from lammps import lammps
import geoparticle as gp

dl = 1e-4 # particle spacing
l_box = 0.01 # lengt of the box
h_box = 0.006 # height of the box
l_water = 0.003 # length of the water column
h_water = 0.0045 # height of the water column
n_thick = 2 # number of particles in thickness direction

Obtain all the particle coordinates

Create the geometry of the solid wall:

wall = gp.ThickRectangle(
    length=l_box,  # inner length
    width=h_box,  # inner height
    n_thick=n_thick,  # number of particle layers
    dl=dl,  # particle spacing
    plane='XOY',  # plane of the rectangle, 'XOY' is the default
    name='wall'  # name of the geometry for reference
)

One may ask where the wall is located. All the geometries are created by placing their centers (for solids of revolution) or vertices (for other solids like the block) at the anchor. See the documentation of each geometry for its anchor definition. By default, the anchor is set to (0, 0, 0). If we do not specify the anchor, the shift method can be used to move the geometry to the desired location.

Now we should move the water region dl to the right and dl up to avoid overlapping with the wall, as well as the gas. To create the gas, we can create a large rectangle and subtract the water region from it:

Since a 2D geometry is created in the XOY plane by default, we do not need to specify the plane here.

water = gp.FilledRectangle(
    length=l_water, width=h_water, dl=dl, name='water',
    anchor=(dl, dl, 0)
)
gas = gp.FilledRectangle(
    l_box - 2 * dl, h_box - 2 * dl, dl, name='gas'
).shift(x=dl, y=dl)
gas = gas.subtract(
    water,  # which geometry to subtract
    rmax=1e-6  # maximum distance to look for overlapping
)

In addition to using gas = gas.subtract(water, rmax) for subtraction, three other methods are provided:

• gas = gp.Subtract(gas, water, rmax)

• gas = gas - water

• gas -= water

The latter two methods use the default rmax value, i.e., 1e-5 (m), defined in the Geometry class.

Now we have obtained all the particle coordinates.

Create particles in LAMMPS

We should first create a LAMMPS instance:

lmp = lammps(cmdargs=['-screen', 'none', '-log', 'none'])

Now create the simulation box. We can leave some buffer between the geometries and the box boundaries. Access the x-, y-, and z-coordinates of the wall geometry using the xs, ys, and zs properties.

fac_buf = 1
xlo = wall.xs.min() - wall.xs.max() * fac_buf
xhi = wall.xs.max() * (1 + fac_buf)
ylo = wall.ys.min() - wall.ys.max() * fac_buf
yhi = wall.ys.max() * (1 + fac_buf)
zlo = -3e-4
zhi = 3e-4
lmp.commands_string(f"""
dimension    2
atom_style   rheo
units                si
newton               on
boundary         f f p
comm_modify  vel yes
region       simulation_box block {xlo} {xhi} {ylo} {yhi} {zlo} {zhi}
create_box   3 simulation_box
""")

Now we can create particles for each geometry. The required method is create_atoms, which takes the number of atoms, their IDs, types, and coordinates as input arguments. We can use the flatten_coords property of each geometry to get the coordinates in the required format.

n_atoms_wall = wall.size
lmp.create_atoms(
    n_atoms_wall, np.arange(n_atoms_wall) + 1 + lmp.get_natoms(),
    np.full(n_atoms_wall, 1, dtype=int), wall.flatten_coords)
n_atoms_water = water.size
lmp.create_atoms(
    n_atoms_water, np.arange(n_atoms_water) + 1 + lmp.get_natoms(),
    np.full(n_atoms_water, 2, dtype=int), water.flatten_coords)
n_atoms_gas = gas.size
lmp.create_atoms(
    n_atoms_gas, np.arange(n_atoms_gas) + 1 + lmp.get_natoms(),
    np.full(n_atoms_gas, 3, dtype=int), gas.flatten_coords)
n_atoms_all = lmp.get_natoms()
log_n_atom = f'n_atoms_wall: {n_atoms_wall},' \
             f' n_atom_water: {n_atoms_water}' \
             f' n_atoms_all: {n_atoms_all}.'
print(log_n_atom)

Check the overlapping particles

Geoparticle already provides check_overlap method to check the overlapping particles, which is done automatically. However, it only checks the overlapping particles within each geometry. To check the entire system, one way is to compute the union of all geometries and use the check_overlap method. Here, we alternatively use LAMMPS commands for the check. We first delete the overlapping atoms. If the number of atoms remains the same after deleting overlapping atoms, then there is no overlapping atoms

# check atom overlapping
lmp.commands_string(f"""
pair_style      zero {dl * 2}
pair_coeff      * *
neighbor        {dl * 0.1} bin
delete_atoms    overlap {dl * 0.8} all all
""")
n_atoms_all_now = lmp.get_natoms()
if n_atoms_all == n_atoms_all_now:
    print('Congrats, no atoms are overlapped!')
else:
    raise RuntimeError(f'{n_atoms_all - n_atoms_all_now} atoms are overlapped!')

Save the data file for LAMMPS for ensuing visualization and simulation. Note that we must assign mass to each particle before writing the data file.

m0_fluid = dl ** 2 * 993
m0_gas = dl ** 2 * 1.1
lmp.commands_string(f"""
mass 1*2 {m0_fluid}
mass 3 {m0_gas}
run 0
write_data gas_liquid_dam2D.data
""")