Skip to content
Snippets Groups Projects

Compare revisions

Changes are shown as if the source revision was being merged into the target revision. Learn more about comparing revisions.

Source

Select target project
No results found
Select Git revision
  • Sparse
  • WallLaw
  • fhennig/pystencils2.0-compat
  • improved_comm
  • master
  • suffa/psm_optimization
  • release/0.2.1
  • release/0.2.10
  • release/0.2.11
  • release/0.2.12
  • release/0.2.13
  • release/0.2.14
  • release/0.2.15
  • release/0.2.2
  • release/0.2.3
  • release/0.2.4
  • release/0.2.5
  • release/0.2.6
  • release/0.2.7
  • release/0.2.8
  • release/0.2.9
  • release/0.3.0
  • release/0.3.1
  • release/0.3.2
  • release/0.3.3
  • release/0.3.4
  • release/0.4.0
  • release/0.4.1
  • release/0.4.2
  • release/0.4.3
  • release/0.4.4
  • release/1.0
  • release/1.0.1
  • release/1.1
  • release/1.1.1
  • release/1.2
  • release/1.3
  • release/1.3.1
  • release/1.3.2
  • release/1.3.3
  • release/1.3.4
  • release/1.3.5
  • release/1.3.6
  • release/1.3.7
44 results

Target

Select target project
  • ravi.k.ayyala/lbmpy
  • brendan-waters/lbmpy
  • anirudh.jonnalagadda/lbmpy
  • jbadwaik/lbmpy
  • alexander.reinauer/lbmpy
  • itischler/lbmpy
  • he66coqe/lbmpy
  • ev81oxyl/lbmpy
  • Bindgen/lbmpy
  • da15siwa/lbmpy
  • holzer/lbmpy
  • RudolfWeeber/lbmpy
  • pycodegen/lbmpy
13 results
Select Git revision
  • GetterSetterAPI
  • HRR
  • HydroPressure
  • InplaceConfig
  • Outflow
  • PhaseField
  • Sparse
  • UBBVelocity
  • UpdateAPISparse
  • WallLaw
  • WetNodeBoundaries
  • csebug
  • feature/sparse
  • feature/try
  • improved_comm
  • install_requires
  • master
  • phaseField
  • relaxationrates
  • test_martin
  • release/0.2.1
  • release/0.2.10
  • release/0.2.11
  • release/0.2.12
  • release/0.2.13
  • release/0.2.14
  • release/0.2.15
  • release/0.2.2
  • release/0.2.3
  • release/0.2.4
  • release/0.2.5
  • release/0.2.6
  • release/0.2.7
  • release/0.2.8
  • release/0.2.9
  • release/0.3.0
  • release/0.3.1
  • release/0.3.2
  • release/0.3.3
  • release/0.3.4
  • release/0.4.0
  • release/0.4.1
  • release/0.4.2
  • release/0.4.3
  • release/0.4.4
  • release/1.0
  • release/1.0.1
  • release/1.1
  • release/1.1.1
  • release/1.2
  • release/1.3
  • release/1.3.1
  • release/1.3.2
  • release/1.3.3
  • release/1.3.4
  • release/1.3.5
  • release/1.3.6
57 results
Show changes
Showing
with 299 additions and 43 deletions
%% Cell type:code id: tags:
``` python
import os
import time
from tqdm.notebook import tqdm
import numpy as np
import math
from scipy.special import erfc
from pystencils.session import *
from lbmpy.session import *
from pystencils.simp import sympy_cse
from pystencils.boundaries import BoundaryHandling
from lbmpy.phasefield_allen_cahn.parameter_calculation import AllenCahnParameters
from lbmpy.phasefield_allen_cahn.contact_angle import ContactAngle
from lbmpy.phasefield_allen_cahn.kernel_equations import *
from lbmpy.advanced_streaming import LBMPeriodicityHandling
from lbmpy.boundaries import NoSlip, LatticeBoltzmannBoundaryHandling
```
%% Cell type:markdown id: tags:
If `pycuda` is installed the simulation automatically runs on GPU
If `cupy` is installed the simulation automatically runs on GPU
%% Cell type:code id: tags:
``` python
try:
import pycuda
import cupy
except ImportError:
pycuda = None
cupy = None
gpu = False
target = ps.Target.CPU
print('No pycuda installed')
print('No cupy installed')
if pycuda:
if cupy:
gpu = True
target = ps.Target.GPU
```
%% Output
No pycuda installed
%% Cell type:markdown id: tags:
# Capillary wave simulated with a phase-field model for immiscible fluids
%% Cell type:markdown id: tags:
## Geometry Setup
First of all the stencils for the phase-field LB step as well as the stencil for the hydrodynamic LB step are defined. According to the stencils the simulation runs either in 2D or 3D
%% Cell type:code id: tags:
``` python
stencil_phase = LBStencil(Stencil.D2Q9)
stencil_hydro = LBStencil(Stencil.D2Q9)
assert(stencil_hydro.D == stencil_phase.D)
dimensions = stencil_phase.D
```
%% Cell type:code id: tags:
``` python
# user defined input
Re = 10 # Reynolds number
domain_width = 50
fluid_depth = 0.5 * domain_width
amplitude = 0.01 * domain_width
relaxation_rate_heavy = 1.99
mobility = 0.02 # phase field mobility
interface_width = 5 # phase field interface width
density_heavy = 1.0 # density of heavy phase
density_ratio = 1000
density_light = density_heavy / density_ratio # density of light phase
kinematic_viscosity_ratio = 1
kinematic_viscosity_heavy = 1 / 3 * (1 / relaxation_rate_heavy - 0.5)
kinematic_viscosity_light = kinematic_viscosity_heavy / kinematic_viscosity_ratio
wavelength = domain_width
wavenumber = 2.0 * np.pi / domain_width
wave_frequency = Re * kinematic_viscosity_heavy / domain_width / amplitude # angular wave frequency
surface_tension = wave_frequency**2 * (density_heavy + density_light) / wavenumber**3
gravitational_acceleration = 0
Pe = domain_width * amplitude * wave_frequency / mobility # Peclet number
Cn = interface_width / domain_width # Cahn number
dynamic_viscosity_heavy = kinematic_viscosity_heavy * density_heavy
relaxation_time_heavy = 3.0 * kinematic_viscosity_heavy
kinematic_viscosity_light = kinematic_viscosity_heavy / kinematic_viscosity_ratio
dynamic_viscosity_light = kinematic_viscosity_light * density_light
relaxation_time_light = 3.0 * kinematic_viscosity_light
timesteps = int(80 / wave_frequency)
data_extract_frequency = int(0.1 / wave_frequency)
vtk_output_frequency = int(1 / wave_frequency)
vtk_output_path = "vtk_out/capillary-wave"
vtk_base_directory = vtk_output_path.split("/")[0] # create directory for vtk-output if it does not yet exist
if not os.path.exists(vtk_base_directory):
os.mkdir(os.getcwd() + "/" + vtk_base_directory)
domain_size = (domain_width, domain_width)
filename = "pf-re-" + str(Re) + "-resolution-" + str(domain_width) + ".txt"
print("timesteps =", timesteps)
print("Re =", Re)
print("Pe =", Pe)
print("Cn =", Cn)
print("domain_width =", domain_width)
print("fluid_depth =", fluid_depth)
print("amplitude =", amplitude)
print("relaxation_rate_heavy =", relaxation_rate_heavy)
print("mobility =", mobility)
print("interface_width =", interface_width)
print("density_heavy =", density_heavy)
print("density_light =", density_light)
print("density_ratio =", density_ratio)
print("kinematic_viscosity_heavy =", kinematic_viscosity_heavy)
print("kinematic_viscosity_light =", kinematic_viscosity_light)
print("kinematic_viscosity_ratio =", kinematic_viscosity_ratio)
print("dynamic_viscosity_heavy =", dynamic_viscosity_heavy)
print("dynamic_viscosity_light =", dynamic_viscosity_light)
print("dynamic_viscosity_ratio =", dynamic_viscosity_heavy/dynamic_viscosity_light)
print("wavelength =", wavelength)
print("wavenumber =", wavenumber)
print("wave_frequency =", wave_frequency)
print("surface_tension =", surface_tension)
print("gravitational_acceleration =", gravitational_acceleration)
```
%% Output
timesteps = 238800
Re = 10
Pe = 0.41876046901171776
Cn = 0.1
domain_width = 50
fluid_depth = 25.0
amplitude = 0.5
relaxation_rate_heavy = 1.99
mobility = 0.02
interface_width = 5
density_heavy = 1.0
density_light = 0.001
density_ratio = 1000
kinematic_viscosity_heavy = 0.0008375209380234356
kinematic_viscosity_light = 0.0008375209380234356
kinematic_viscosity_ratio = 1
dynamic_viscosity_heavy = 0.0008375209380234356
dynamic_viscosity_light = 8.375209380234357e-07
dynamic_viscosity_ratio = 1000.0
wavelength = 50
wavenumber = 0.12566370614359174
wave_frequency = 0.00033500837520937423
surface_tension = 5.6612953740701554e-05
gravitational_acceleration = 0
%% Cell type:code id: tags:
``` python
parameters = AllenCahnParameters(density_heavy=density_heavy, density_light=density_light,
dynamic_viscosity_heavy=dynamic_viscosity_heavy,
dynamic_viscosity_light=dynamic_viscosity_light,
surface_tension=surface_tension, mobility=mobility,
interface_thickness=interface_width)
```
%% Cell type:code id: tags:
``` python
parameters
```
%% Output
<lbmpy.phasefield_allen_cahn.parameter_calculation.AllenCahnParameters at 0x146131580>
%% Cell type:markdown id: tags:
## Fields
As a next step all fields which are needed get defined. To do so we create a `datahandling` object. More details about it can be found in the third tutorial of the [pystencils framework]( http://pycodegen.pages.walberla.net/pystencils/). Basically it holds all fields and manages the kernel runs.
%% Cell type:code id: tags:
``` python
# create a datahandling object
dh = ps.create_data_handling((domain_size), periodicity=(True, False), parallel=False, default_target=target)
# fields
g = dh.add_array("g", values_per_cell=len(stencil_hydro))
dh.fill("g", 0.0, ghost_layers=True)
h = dh.add_array("h",values_per_cell=len(stencil_phase))
dh.fill("h", 0.0, ghost_layers=True)
g_tmp = dh.add_array("g_tmp", values_per_cell=len(stencil_hydro))
dh.fill("g_tmp", 0.0, ghost_layers=True)
h_tmp = dh.add_array("h_tmp",values_per_cell=len(stencil_phase))
dh.fill("h_tmp", 0.0, ghost_layers=True)
u = dh.add_array("u", values_per_cell=dh.dim)
dh.fill("u", 0.0, ghost_layers=True)
rho = dh.add_array("rho", values_per_cell=1)
dh.fill("rho", 1.0, ghost_layers=True)
C = dh.add_array("C")
dh.fill("C", 0.0, ghost_layers=True)
C_tmp = dh.add_array("C_tmp")
dh.fill("C_tmp", 0.0, ghost_layers=True)
```
%% Cell type:code id: tags:
``` python
# set the frequency of the file outputs
vtk_writer = dh.create_vtk_writer(vtk_output_path, ["C", "u", "rho"], ghost_layers=True)
```
%% Cell type:code id: tags:
``` python
rho_L = parameters.symbolic_density_light
rho_H = parameters.symbolic_density_heavy
density = rho_L + C.center * (rho_H - rho_L)
body_force = [0, 0, 0]
body_force[1] = parameters.symbolic_gravitational_acceleration * density
```
%% Cell type:markdown id: tags:
## Definition of the lattice Boltzmann methods
%% Cell type:code id: tags:
``` python
w_c = parameters.symbolic_omega_phi
config_phase = LBMConfig(stencil=stencil_phase, method=Method.MRT, compressible=True,
delta_equilibrium=False,
force=sp.symbols("F_:2"), velocity_input=u,
weighted=True, relaxation_rates=[0, w_c, w_c, 1, 1, 1, 1, 1, 1],
output={'density': C_tmp}, kernel_type='stream_pull_collide')
method_phase = create_lb_method(lbm_config=config_phase)
```
%% Cell type:code id: tags:
``` python
omega = parameters.omega(C)
config_hydro = LBMConfig(stencil=stencil_hydro, method=Method.MRT, compressible=False,
weighted=True, relaxation_rates=[omega, 1, 1, 1],
force=sp.symbols("F_:2"),
output={'velocity': u}, kernel_type='collide_stream_push')
method_hydro = create_lb_method(lbm_config=config_hydro)
```
%% Cell type:markdown id: tags:
## Initialization
%% Cell type:code id: tags:
``` python
h_updates = initializer_kernel_phase_field_lb(method_phase, C, u, h, parameters)
g_updates = initializer_kernel_hydro_lb(method_hydro, 1, u, g)
h_init = ps.create_kernel(h_updates, target=dh.default_target, cpu_openmp=True).compile()
g_init = ps.create_kernel(g_updates, target=dh.default_target, cpu_openmp=True).compile()
```
%% Output
[['spatialInner0'], ['spatialInner1']]
[['spatialInner0'], ['spatialInner1']]
%% Cell type:code id: tags:
``` python
# initialize the phase field
def initialize_phasefield():
Nx = domain_size[0]
Ny = domain_size[1]
for block in dh.iterate(ghost_layers=True, inner_ghost_layers=False):
# get x as cell center coordinate, i.e., including shift with 0.5
x = np.zeros_like(block.midpoint_arrays[0])
x[:, :] = block.midpoint_arrays[0]
# get y as cell center coordinate, i.e., including shift with 0.5
y = np.zeros_like(block.midpoint_arrays[1])
y[:, :] = block.midpoint_arrays[1]
tmp = fluid_depth + amplitude * np.cos(x / (Nx - 1) * np.pi * 2 + np.pi)
# initialize diffuse interface with tanh profile
init_values = 0.5 - 0.5 * np.tanh((y - tmp) / (interface_width / 2))
block["C"][:, :] = init_values
block["C_tmp"][:, :] = init_values
if gpu:
dh.all_to_gpu()
dh.run_kernel(h_init, **parameters.symbolic_to_numeric_map)
dh.run_kernel(g_init)
```
%% Cell type:code id: tags:
``` python
# plot initial profile
x = np.arange(0, domain_size[0], 1) # start,stop,step
y = fluid_depth + amplitude * np.cos(x / (domain_size[0] - 1) * np.pi * 2 + np.pi)
plt.plot(x,y, marker='o')
plt.show()
```
%% Output
%% Cell type:code id: tags:
``` python
# check symmetry of sine-curve
for i in range(0, domain_size[0] - 1):
if (i >= domain_size[0] * 0.5):
continue
if (y[i] != y[domain_size[0] - 1 - i]):
print("asymmetric:", y[i], y[domain_size[0] - 1 - i])
```
%% Output
asymmetric: 25.111260466978155 25.11126046697816
%% Cell type:code id: tags:
``` python
initialize_phasefield()
```
%% Cell type:code id: tags:
``` python
plt.scalar_field(dh.cpu_arrays["C"])
plt.show()
```
%% Output
%% Cell type:code id: tags:
``` python
# plot initial interface profile
x = np.arange(0, domain_size[1]+2, 1) # start,stop,step
y = dh.cpu_arrays["C"][int(domain_size[0] * 0.5), :]
plt.plot(x,y, marker='o')
plt.grid()
plt.show()
```
%% Output
%% Cell type:code id: tags:
``` python
force_h = interface_tracking_force(C, stencil_phase, parameters)
hydro_force = hydrodynamic_force(g, C, method_hydro, parameters, body_force)
```
%% Cell type:markdown id: tags:
## Definition of the LB update rules
%% Cell type:code id: tags:
``` python
lbm_optimisation = LBMOptimisation(symbolic_field=h, symbolic_temporary_field=h_tmp)
allen_cahn_update_rule = create_lb_update_rule(lbm_config=config_phase,
lbm_optimisation=lbm_optimisation)
allen_cahn_update_rule = add_interface_tracking_force(allen_cahn_update_rule, force_h)
ast_kernel = ps.create_kernel(allen_cahn_update_rule, target=dh.default_target, cpu_openmp=True)
kernel_allen_cahn_lb = ast_kernel.compile()
```
%% Output
[['spatialInner0'], ['spatialInner1']]
%% Cell type:code id: tags:
``` python
force_Assignments = hydrodynamic_force_assignments(g, u, C, method_hydro, parameters, body_force)
lbm_optimisation = LBMOptimisation(symbolic_field=g, symbolic_temporary_field=g_tmp)
hydro_lb_update_rule = create_lb_update_rule(lbm_config=config_hydro,
lbm_optimisation=lbm_optimisation)
hydro_lb_update_rule = add_hydrodynamic_force(hydro_lb_update_rule, force_Assignments, C, g, parameters)
ast_kernel = ps.create_kernel(hydro_lb_update_rule, target=dh.default_target, cpu_openmp=True)
kernel_hydro_lb = ast_kernel.compile()
```
%% Output
[['spatialInner0'], ['spatialInner1']]
%% Cell type:markdown id: tags:
## Boundary Conditions
%% Cell type:code id: tags:
``` python
# periodic Boundarys for g, h and C
periodic_BC_C = dh.synchronization_function(C.name, target=dh.default_target, optimization = {"openmp": True})
periodic_BC_g = LBMPeriodicityHandling(stencil=stencil_hydro, data_handling=dh, pdf_field_name=g.name,
streaming_pattern='push')
periodic_BC_h = LBMPeriodicityHandling(stencil=stencil_phase, data_handling=dh, pdf_field_name=h.name,
streaming_pattern='pull')
# No slip boundary for the phasefield lbm
bh_allen_cahn = LatticeBoltzmannBoundaryHandling(method_phase, dh, 'h',
target=dh.default_target, name='boundary_handling_h',
streaming_pattern='pull')
# No slip boundary for the velocityfield lbm
bh_hydro = LatticeBoltzmannBoundaryHandling(method_hydro, dh, 'g' ,
target=dh.default_target, name='boundary_handling_g',
streaming_pattern='push')
contact_angle = BoundaryHandling(dh, C.name, stencil_hydro, target=dh.default_target)
contact = ContactAngle(90, interface_width)
wall = NoSlip()
if dimensions == 2:
bh_allen_cahn.set_boundary(wall, make_slice[:, 0])
bh_allen_cahn.set_boundary(wall, make_slice[:, -1])
bh_hydro.set_boundary(wall, make_slice[:, 0])
bh_hydro.set_boundary(wall, make_slice[:, -1])
contact_angle.set_boundary(contact, make_slice[:, 0])
contact_angle.set_boundary(contact, make_slice[:, -1])
else:
bh_allen_cahn.set_boundary(wall, make_slice[:, 0, :])
bh_allen_cahn.set_boundary(wall, make_slice[:, -1, :])
bh_hydro.set_boundary(wall, make_slice[:, 0, :])
bh_hydro.set_boundary(wall, make_slice[:, -1, :])
contact_angle.set_boundary(contact, make_slice[:, 0, :])
contact_angle.set_boundary(contact, make_slice[:, -1, :])
bh_allen_cahn.prepare()
bh_hydro.prepare()
contact_angle.prepare()
```
%% Cell type:markdown id: tags:
## Full timestep
%% Cell type:code id: tags:
``` python
# definition of the timestep for the immiscible fluids model
def timeloop():
# Solve the interface tracking LB step with boundary conditions
periodic_BC_h()
bh_allen_cahn()
dh.run_kernel(kernel_allen_cahn_lb, **parameters.symbolic_to_numeric_map)
dh.swap("C", "C_tmp")
# apply the three phase-phase contact angle
contact_angle()
# periodic BC of the phase-field
periodic_BC_C()
# solve the hydro LB step with boundary conditions
dh.run_kernel(kernel_hydro_lb, **parameters.symbolic_to_numeric_map)
periodic_BC_g()
bh_hydro()
# compute density (only for vtk output)
# must be done BEFORE swapping fields to avoid having outdated values
compute_density()
# field swaps
dh.swap("h", "h_tmp")
dh.swap("g", "g_tmp")
```
%% Cell type:code id: tags:
``` python
def initialize_hydrostatic_pressure():
for block in dh.iterate(ghost_layers=True, inner_ghost_layers=False):
# get y as cell center coordinate, i.e., including shift with 0.5
y = np.zeros_like(block.midpoint_arrays[1])
y[:, :] = block.midpoint_arrays[1]
# compute hydrostatic density
rho_hydrostatic = 1 + 3 * gravitational_acceleration * (y - fluid_depth)
# subtract 1 because PDFs in incompressible LBM are centered around zero in lbmpy
rho_hydrostatic -= 1
# set equilibrium PDFs with velocity=0 and rho;
for i in range(0, stencil_hydro.Q, 1):
block["g"][:,:,i] = method_hydro.weights[i] * rho_hydrostatic[:,:]
block["g_tmp"][:,:,i] = method_hydro.weights[i] * rho_hydrostatic[:,:]
```
%% Cell type:code id: tags:
``` python
def compute_density():
for block in dh.iterate(ghost_layers=True, inner_ghost_layers=False):
# PDFs in incompressible LBM are centered around zero in lbmpy
# => add 1 to whole sum, i.e., initialize with 1
block["rho"].fill(1);
# compute density
for i in range(block["g"].shape[-1]):
block["rho"][:,:] += block["g"][:,:,i]
```
%% Cell type:code id: tags:
``` python
if (gravitational_acceleration != 0):
initialize_hydrostatic_pressure()
compute_density()
plt.scalar_field(dh.cpu_arrays["rho"])
```
%% Cell type:code id: tags:
``` python
print("================================= start of the simulation ===================================")
start = time.time()
pbar = tqdm(total=timesteps)
timestep = []
surface_position = []
symmetry_norm = []
for i in range(0, timesteps):
sum_c_2 = 0.0
sum_delta_c_2 = 0.0
# write vtk output
if(i % vtk_output_frequency == 0):
if gpu:
dh.to_cpu("C")
dh.to_cpu("u")
dh.to_cpu("rho")
vtk_writer(i)
# extract data (to be written to file)
if(i % data_extract_frequency == 0):
if gpu:
dh.to_cpu("C")
dh.to_cpu("u")
dh.to_cpu("rho")
pbar.update(data_extract_frequency)
timestep.append(i)
ny = domain_size[1]
# get index containing phase field value < 0.5
i1 = np.argmax(dh.cpu_arrays["C"][int(domain_size[0] * 0.5), :] < 0.5)
i0 = i1 - 1 # index containing phase field value >= 0.5
f0 = dh.cpu_arrays["C"][int(domain_size[0] * 0.5), i0] # phase field value >= 0.5
f1 = dh.cpu_arrays["C"][int(domain_size[0] * 0.5), i1] # phase field value < 0.5
# coordinate of cell center is index+0.5-1 (0.5 to get to cell center; -1 to remove ghost layer from index)
y0 = i0 + 0.5 - 1
y1 = i1 + 0.5 - 1
#interpolate
surface_position.append( y0 + (y1 - y0) / (f1 - f0) * (0.5 - f0) )
# evaluate symmetry in x-direction
for y in range(0, domain_size[1] - 1):
for x in range(0, domain_size[0] - 1):
if (x >= domain_size[0] * 0.5):
continue
x_mirrored = domain_size[0] - 1 - x;
sum_c_2 += dh.cpu_arrays["C"][x, y]**2
sum_delta_c_2 += (dh.cpu_arrays["C"][x, y] - dh.cpu_arrays["C"][x_mirrored, y])**2
symmetry_norm.append( (sum_delta_c_2 / sum_c_2)**2 )
timeloop()
pbar.close()
end = time.time()
sim_time = end - start
print("\n")
print("time needed for the calculation: %4.4f" % sim_time , "seconds")
nrOfCells = np.prod(domain_size)
mlups = nrOfCells * timesteps / sim_time * 1e-6
print("MLUPS: %4.2f" % mlups)
```
%% Output
================================= start of the simulation ===================================
time needed for the calculation: 92.7764 seconds
MLUPS: 6.43
%% Cell type:code id: tags:
``` python
if gpu:
dh.to_cpu("C")
dh.to_cpu("u")
dh.to_cpu("rho")
plt.scalar_field(dh.cpu_arrays["C"])
```
%% Output
<matplotlib.image.AxesImage at 0x146e329a0>
%% Cell type:code id: tags:
``` python
# non-dimensionalize time and surface position
t_nd = [value * wave_frequency for value in timestep]
a_nd = [(value - fluid_depth) / amplitude for value in surface_position]
```
%% Cell type:code id: tags:
``` python
# plot surface position over time
plt.xlabel('$t/-$', fontsize=18)
plt.ylabel('$a/-$', fontsize=18)
plt.grid(color='black', linestyle='-', linewidth=0.1)
plt.plot(t_nd, a_nd, color=(0.121, 0.231, 0.4))
plt.show()
```
%% Output
%% Cell type:code id: tags:
``` python
# plot analytical solutions
# analytical model from Prosperetti, Motion of two superposed viscous fluids, 1981, doi: 10.1063/1.863522
nu = kinematic_viscosity_heavy
k = wavenumber
omega_0 = wave_frequency
rho_u = density_heavy
rho_l = density_light
a0 = amplitude
beta = rho_l * rho_u / (rho_l + rho_u)**2
k2nu = k**2 * nu # helper variable
# polynomials from equation (21)
p0 = (1 - 4 * beta) * k2nu**2 + omega_0**2
p1 = 4 * (1 - 3 * beta) * k2nu**1.5
p2 = 2 * (1 - 6 * beta) * k2nu
p3 = -4 * beta * k2nu**(0.5)
p4 = 1
p = [p0, p1, p2, p3, p4]
# get roots of equation (21)
z = np.polynomial.polynomial.polyroots(p)
Z = np.array([(z[1] - z[0]) * (z[2] - z[0]) * (z[3] - z[0]),
(z[0] - z[1]) * (z[2] - z[1]) * (z[3] - z[1]),
(z[0] - z[2]) * (z[1] - z[2]) * (z[3] - z[2]),
(z[0] - z[3]) * (z[1] - z[3]) * (z[2] - z[3])])
t = np.arange(0, timesteps, 1)
# equation (20)
tmp0 = (1 - 4 * beta) * k2nu**2 # helper variable
term0 = 4 * tmp0 / (8 * tmp0 + omega_0**2) * a0 * erfc((k2nu * t)**0.5) # first term in equation (20)
term1 = 0.0
for i in range(0, 4, 1):
tmp1 = z[i]**2 - k2nu # helper variable
term1 += z[i] / Z[i] * omega_0**2 * a0 / tmp1 * np.exp(tmp1 * t) * erfc(z[i] * t**0.5)
a = np.real(term0 + term1)
# non-dimesionalize time and surface position
pos_anal_lin_org = a / a0
t_anal_org = t * omega_0
# non-dimesionalize time and surface position
a_anal_nd = a / amplitude
t_anal_nd = t * wave_frequency
plt.plot(t_anal_nd, a_anal_nd, label="analytical solution")
plt.plot(t_nd, a_nd, label="simulation")
# correction factor from Denner et al., Dispersion and viscous attenuation of capillary waves with finite amplitude, 2016, doi: 10.1140/epjst/e2016-60199-2
# model is now valid for initial amplitudes of about 0.1*wavelength
a0_hat = a0 / wavelength
C = -4.5 * a0_hat**3 + 5.3 * a0_hat**2 + 0.18 * a0_hat + 1
a_anal_nd_corr = a_anal_nd * C
t_anal_nd_corr = t_anal_nd * C
plt.plot(t_anal_nd_corr, a_anal_nd_corr, label="linear analytical solution with correction factor")
plt.xlabel('$t/-$', fontsize=18)
plt.ylabel('$a/-$', fontsize=18)
plt.legend()
plt.grid()
plt.show()
```
%% Output
%% Cell type:code id: tags:
``` python
# plot symmetry norm over time
plt.xlabel('$t/-$', fontsize=18)
plt.ylabel('$symmetry-norm/-$', fontsize=18)
plt.grid(color='black', linestyle='-', linewidth=0.1)
plt.plot(t_nd, symmetry_norm, color=(0.121, 0.231, 0.4))
plt.show()
```
%% Output
%% Cell type:code id: tags:
``` python
# store results to file
import csv
with open(filename, 'w') as f:
writer = csv.writer(f, delimiter='\t')
writer.writerows(zip(t_nd, a_nd, symmetry_norm))
```
......
%% Cell type:code id: tags:
``` python
import os
import time
from tqdm.notebook import tqdm
import numpy as np
import math
from pystencils.session import *
from lbmpy.session import *
from pystencils.simp import sympy_cse
from pystencils.boundaries import BoundaryHandling
from lbmpy.phasefield_allen_cahn.parameter_calculation import AllenCahnParameters
from lbmpy.phasefield_allen_cahn.contact_angle import ContactAngle
from lbmpy.phasefield_allen_cahn.kernel_equations import *
from lbmpy.advanced_streaming import LBMPeriodicityHandling
from lbmpy.boundaries import NoSlip, LatticeBoltzmannBoundaryHandling
```
%% Cell type:markdown id: tags:
If `pycuda` is installed the simulation automatically runs on GPU
If `cupy` is installed the simulation automatically runs on GPU
%% Cell type:code id: tags:
``` python
try:
import pycuda
import cupy
except ImportError:
pycuda = None
cupy = None
gpu = False
target = ps.Target.CPU
print('No pycuda installed')
print('No cupy installed')
if pycuda:
if cupy:
gpu = True
target = ps.Target.GPU
```
%% Output
No pycuda installed
%% Cell type:markdown id: tags:
# Gravity wave simulated with a phase-field model for immiscible fluids
%% Cell type:markdown id: tags:
## Geometry Setup
First of all the stencils for the phase-field LB step as well as the stencil for the hydrodynamic LB step are defined. According to the stencils the simulation runs either in 2D or 3D
%% Cell type:code id: tags:
``` python
stencil_phase = LBStencil(Stencil.D2Q9)
stencil_hydro = LBStencil(Stencil.D2Q9)
assert(stencil_hydro.D == stencil_phase.D)
dimensions = stencil_phase.D
```
%% Cell type:code id: tags:
``` python
# user defined input
Re = 10 # Reynolds number
domain_width = 50
fluid_depth = 0.5 * domain_width
amplitude = 0.01 * domain_width
relaxation_rate_heavy = 1.99
mobility = 0.02 # phase field mobility
interface_width = 5 # phase field interface width
density_heavy = 1.0 # density of heavy phase
density_ratio = 1000
density_light = density_heavy / density_ratio # density of light phase
kinematic_viscosity_ratio = 1
kinematic_viscosity_heavy = 1 / 3 * (1 / relaxation_rate_heavy - 0.5)
kinematic_viscosity_light = kinematic_viscosity_heavy / kinematic_viscosity_ratio
wavelength = domain_width
wavenumber = 2.0 * np.pi / domain_width
wave_frequency = Re * kinematic_viscosity_heavy / domain_width / amplitude # angular wave frequency
surface_tension = 0
gravitational_acceleration = - wave_frequency**2 / wavenumber / np.tanh(wavenumber * fluid_depth)
Pe = domain_width * amplitude * wave_frequency / mobility # Peclet number
Cn = interface_width / domain_width # Cahn number
dynamic_viscosity_heavy = kinematic_viscosity_heavy * density_heavy
relaxation_time_heavy = 3.0 * kinematic_viscosity_heavy
kinematic_viscosity_light = kinematic_viscosity_heavy / kinematic_viscosity_ratio
dynamic_viscosity_light = kinematic_viscosity_light * density_light
relaxation_time_light = 3.0 * kinematic_viscosity_light
timesteps = int(80 / wave_frequency)
data_extract_frequency = int(0.1 / wave_frequency)
vtk_output_frequency = int(1 / wave_frequency)
vtk_output_path = "vtk_out/gravity-wave"
vtk_base_directory = vtk_output_path.split("/")[0] # create directory for vtk-output if it does not yet exist
if not os.path.exists(vtk_base_directory):
os.mkdir(os.getcwd() + "/" + vtk_base_directory)
domain_size = (domain_width, domain_width)
filename = "pf-re-" + str(Re) + "-resolution-" + str(domain_width) + ".txt"
print("timesteps =", timesteps)
print("Re =", Re)
print("Pe =", Pe)
print("Cn =", Cn)
print("domain_width =", domain_width)
print("fluid_depth =", fluid_depth)
print("amplitude =", amplitude)
print("relaxation_rate_heavy =", relaxation_rate_heavy)
print("mobility =", mobility)
print("interface_width =", interface_width)
print("density_heavy =", density_heavy)
print("density_light =", density_light)
print("density_ratio =", density_ratio)
print("kinematic_viscosity_heavy =", kinematic_viscosity_heavy)
print("kinematic_viscosity_light =", kinematic_viscosity_light)
print("kinematic_viscosity_ratio =", kinematic_viscosity_ratio)
print("dynamic_viscosity_heavy =", dynamic_viscosity_heavy)
print("dynamic_viscosity_light =", dynamic_viscosity_light)
print("dynamic_viscosity_ratio =", dynamic_viscosity_heavy/dynamic_viscosity_light)
print("wavelength =", wavelength)
print("wavenumber =", wavenumber)
print("wave_frequency =", wave_frequency)
print("surface_tension =", surface_tension)
print("gravitational_acceleration =", gravitational_acceleration)
print("data_extract_frequency =", data_extract_frequency)
print("vtk_output_frequency =", vtk_output_frequency)
```
%% Output
timesteps = 238800
Re = 10
Pe = 0.41876046901171776
Cn = 0.1
domain_width = 50
fluid_depth = 25.0
amplitude = 0.5
relaxation_rate_heavy = 1.99
mobility = 0.02
interface_width = 5
density_heavy = 1.0
density_light = 0.001
density_ratio = 1000
kinematic_viscosity_heavy = 0.0008375209380234356
kinematic_viscosity_light = 0.0008375209380234356
kinematic_viscosity_ratio = 1
dynamic_viscosity_heavy = 0.0008375209380234356
dynamic_viscosity_light = 8.375209380234357e-07
dynamic_viscosity_ratio = 1000.0
wavelength = 50
wavenumber = 0.12566370614359174
wave_frequency = 0.00033500837520937423
surface_tension = 0
gravitational_acceleration = -8.964447065459421e-07
data_extract_frequency = 298
vtk_output_frequency = 2985
%% Cell type:code id: tags:
``` python
parameters = AllenCahnParameters(density_heavy=density_heavy, density_light=density_light,
dynamic_viscosity_heavy=dynamic_viscosity_heavy,
dynamic_viscosity_light=dynamic_viscosity_light,
gravitational_acceleration=gravitational_acceleration,
surface_tension=surface_tension, mobility=mobility,
interface_thickness=interface_width)
```
%% Cell type:code id: tags:
``` python
parameters
```
%% Output
<lbmpy.phasefield_allen_cahn.parameter_calculation.AllenCahnParameters at 0x12d1bd550>
%% Cell type:markdown id: tags:
## Fields
As a next step all fields which are needed get defined. To do so we create a `datahandling` object. More details about it can be found in the third tutorial of the [pystencils framework]( http://pycodegen.pages.walberla.net/pystencils/). Basically it holds all fields and manages the kernel runs.
%% Cell type:code id: tags:
``` python
# create a datahandling object
dh = ps.create_data_handling((domain_size), periodicity=(True, False), parallel=False, default_target=target)
# fields
g = dh.add_array("g", values_per_cell=len(stencil_hydro))
dh.fill("g", 0.0, ghost_layers=True)
h = dh.add_array("h",values_per_cell=len(stencil_phase))
dh.fill("h", 0.0, ghost_layers=True)
g_tmp = dh.add_array("g_tmp", values_per_cell=len(stencil_hydro))
dh.fill("g_tmp", 0.0, ghost_layers=True)
h_tmp = dh.add_array("h_tmp",values_per_cell=len(stencil_phase))
dh.fill("h_tmp", 0.0, ghost_layers=True)
u = dh.add_array("u", values_per_cell=dh.dim)
dh.fill("u", 0.0, ghost_layers=True)
rho = dh.add_array("rho", values_per_cell=1)
dh.fill("rho", 1.0, ghost_layers=True)
C = dh.add_array("C")
dh.fill("C", 0.0, ghost_layers=True)
C_tmp = dh.add_array("C_tmp")
dh.fill("C_tmp", 0.0, ghost_layers=True)
```
%% Cell type:code id: tags:
``` python
# set the frequency of the file outputs
vtk_writer = dh.create_vtk_writer(vtk_output_path, ["C", "u", "rho"], ghost_layers=True)
```
%% Cell type:code id: tags:
``` python
rho_L = parameters.symbolic_density_light
rho_H = parameters.symbolic_density_heavy
density = rho_L + C.center * (rho_H - rho_L)
body_force = [0, 0, 0]
body_force[1] = parameters.symbolic_gravitational_acceleration * density
```
%% Cell type:markdown id: tags:
## Definition of the lattice Boltzmann methods
%% Cell type:code id: tags:
``` python
w_c = parameters.symbolic_omega_phi
config_phase = LBMConfig(stencil=stencil_phase, method=Method.MRT, compressible=True,
delta_equilibrium=False,
force=sp.symbols("F_:2"), velocity_input=u,
weighted=True, relaxation_rates=[0, w_c, w_c, 1, 1, 1, 1, 1, 1],
output={'density': C_tmp}, kernel_type='stream_pull_collide')
method_phase = create_lb_method(lbm_config=config_phase)
```
%% Cell type:code id: tags:
``` python
omega = parameters.omega(C)
config_hydro = LBMConfig(stencil=stencil_hydro, method=Method.MRT, compressible=False,
weighted=True, relaxation_rates=[omega, 1, 1, 1],
force=sp.symbols("F_:2"),
output={'velocity': u}, kernel_type='collide_stream_push')
method_hydro = create_lb_method(lbm_config=config_hydro)
```
%% Cell type:code id: tags:
``` python
method_hydro
```
%% Output
<lbmpy.methods.momentbased.momentbasedmethod.MomentBasedLbMethod at 0x12d4a1e50>
%% Cell type:markdown id: tags:
## Initialization
%% Cell type:code id: tags:
``` python
h_updates = initializer_kernel_phase_field_lb(method_phase, C, u, h, parameters)
g_updates = initializer_kernel_hydro_lb(method_hydro, 1, u, g)
h_init = ps.create_kernel(h_updates, target=dh.default_target, cpu_openmp=True).compile()
g_init = ps.create_kernel(g_updates, target=dh.default_target, cpu_openmp=True).compile()
```
%% Output
[['spatialInner0'], ['spatialInner1']]
[['spatialInner0'], ['spatialInner1']]
%% Cell type:code id: tags:
``` python
# initialize the phase field
def initialize_phasefield():
Nx = domain_size[0]
Ny = domain_size[1]
for block in dh.iterate(ghost_layers=True, inner_ghost_layers=False):
# get x as cell center coordinate, i.e., including shift with 0.5
x = np.zeros_like(block.midpoint_arrays[0])
x[:, :] = block.midpoint_arrays[0]
# get y as cell center coordinate, i.e., including shift with 0.5
y = np.zeros_like(block.midpoint_arrays[1])
y[:, :] = block.midpoint_arrays[1]
tmp = fluid_depth + amplitude * np.cos(x / (Nx - 1) * np.pi * 2 + np.pi)
# initialize diffuse interface with tanh profile
init_values = 0.5 - 0.5 * np.tanh((y - tmp) / (interface_width / 2))
block["C"][:, :] = init_values
block["C_tmp"][:, :] = init_values
if gpu:
dh.all_to_gpu()
dh.run_kernel(h_init, **parameters.symbolic_to_numeric_map)
dh.run_kernel(g_init)
```
%% Cell type:code id: tags:
``` python
# plot initial profile
x = np.arange(0, domain_size[0], 1) # start,stop,step
y = fluid_depth + amplitude * np.cos(x / (domain_size[0] - 1) * np.pi * 2 + np.pi)
plt.plot(x,y, marker='o')
plt.show()
```
%% Output
%% Cell type:code id: tags:
``` python
initialize_phasefield()
```
%% Cell type:code id: tags:
``` python
plt.scalar_field(dh.cpu_arrays["C"])
plt.show()
```
%% Output
%% Cell type:code id: tags:
``` python
# plot initial interface profile
x = np.arange(0, domain_size[1]+2, 1) # start,stop,step
y = dh.cpu_arrays["C"][int(domain_size[0] * 0.5), :]
plt.plot(x,y, marker='o')
plt.grid()
plt.show()
```
%% Output
%% Cell type:code id: tags:
``` python
force_h = interface_tracking_force(C, stencil_phase, parameters)
hydro_force = hydrodynamic_force(g, C, method_hydro, parameters, body_force)
```
%% Cell type:markdown id: tags:
## Definition of the LB update rules
%% Cell type:code id: tags:
``` python
lbm_optimisation = LBMOptimisation(symbolic_field=h, symbolic_temporary_field=h_tmp)
allen_cahn_update_rule = create_lb_update_rule(lbm_config=config_phase,
lbm_optimisation=lbm_optimisation)
allen_cahn_update_rule = add_interface_tracking_force(allen_cahn_update_rule, force_h)
ast_kernel = ps.create_kernel(allen_cahn_update_rule, target=dh.default_target, cpu_openmp=True)
kernel_allen_cahn_lb = ast_kernel.compile()
```
%% Output
[['spatialInner0'], ['spatialInner1']]
%% Cell type:code id: tags:
``` python
force_Assignments = hydrodynamic_force_assignments(g, u, C, method_hydro, parameters, body_force)
lbm_optimisation = LBMOptimisation(symbolic_field=g, symbolic_temporary_field=g_tmp)
hydro_lb_update_rule = create_lb_update_rule(lbm_config=config_hydro,
lbm_optimisation=lbm_optimisation)
hydro_lb_update_rule = add_hydrodynamic_force(hydro_lb_update_rule, force_Assignments, C, g, parameters)
ast_kernel = ps.create_kernel(hydro_lb_update_rule, target=dh.default_target, cpu_openmp=True)
kernel_hydro_lb = ast_kernel.compile()
```
%% Output
[['spatialInner0'], ['spatialInner1']]
%% Cell type:markdown id: tags:
## Boundary Conditions
%% Cell type:code id: tags:
``` python
# periodic Boundarys for g, h and C
periodic_BC_C = dh.synchronization_function(C.name, target=dh.default_target, optimization = {"openmp": True})
periodic_BC_g = LBMPeriodicityHandling(stencil=stencil_hydro, data_handling=dh, pdf_field_name=g.name,
streaming_pattern='push')
periodic_BC_h = LBMPeriodicityHandling(stencil=stencil_phase, data_handling=dh, pdf_field_name=h.name,
streaming_pattern='pull')
# No slip boundary for the phasefield lbm
bh_allen_cahn = LatticeBoltzmannBoundaryHandling(method_phase, dh, 'h',
target=dh.default_target, name='boundary_handling_h',
streaming_pattern='pull')
# No slip boundary for the velocityfield lbm
bh_hydro = LatticeBoltzmannBoundaryHandling(method_hydro, dh, 'g' ,
target=dh.default_target, name='boundary_handling_g',
streaming_pattern='push')
contact_angle = BoundaryHandling(dh, C.name, stencil_hydro, target=dh.default_target)
contact = ContactAngle(90, interface_width)
wall = NoSlip()
if dimensions == 2:
bh_allen_cahn.set_boundary(wall, make_slice[:, 0])
bh_allen_cahn.set_boundary(wall, make_slice[:, -1])
bh_hydro.set_boundary(wall, make_slice[:, 0])
bh_hydro.set_boundary(wall, make_slice[:, -1])
contact_angle.set_boundary(contact, make_slice[:, 0])
contact_angle.set_boundary(contact, make_slice[:, -1])
else:
bh_allen_cahn.set_boundary(wall, make_slice[:, 0, :])
bh_allen_cahn.set_boundary(wall, make_slice[:, -1, :])
bh_hydro.set_boundary(wall, make_slice[:, 0, :])
bh_hydro.set_boundary(wall, make_slice[:, -1, :])
contact_angle.set_boundary(contact, make_slice[:, 0, :])
contact_angle.set_boundary(contact, make_slice[:, -1, :])
bh_allen_cahn.prepare()
bh_hydro.prepare()
contact_angle.prepare()
```
%% Cell type:markdown id: tags:
## Full timestep
%% Cell type:code id: tags:
``` python
# definition of the timestep for the immiscible fluids model
def timeloop():
# Solve the interface tracking LB step with boundary conditions
periodic_BC_h()
bh_allen_cahn()
dh.run_kernel(kernel_allen_cahn_lb, **parameters.symbolic_to_numeric_map)
dh.swap("C", "C_tmp")
# apply the three phase-phase contact angle
contact_angle()
# periodic BC of the phase-field
periodic_BC_C()
# solve the hydro LB step with boundary conditions
dh.run_kernel(kernel_hydro_lb, **parameters.symbolic_to_numeric_map)
periodic_BC_g()
bh_hydro()
# compute density (only for vtk output)
# must be done BEFORE swapping fields to avoid having outdated values
compute_density()
# field swaps
dh.swap("h", "h_tmp")
dh.swap("g", "g_tmp")
```
%% Cell type:code id: tags:
``` python
def initialize_hydrostatic_pressure():
for block in dh.iterate(ghost_layers=True, inner_ghost_layers=False):
# get y as cell center coordinate, i.e., including shift with 0.5
y = np.zeros_like(block.midpoint_arrays[1])
y[:, :] = block.midpoint_arrays[1]
# compute hydrostatic density
rho_hydrostatic = 1 + 3 * gravitational_acceleration * (y - fluid_depth)
# subtract 1 because PDFs in incompressible LBM are centered around zero in lbmpy
rho_hydrostatic -= 1
# set equilibrium PDFs with velocity=0 and rho;
for i in range(0, stencil_hydro.Q, 1):
block["g"][:,:,i] = method_hydro.weights[i] * rho_hydrostatic[:,:]
block["g_tmp"][:,:,i] = method_hydro.weights[i] * rho_hydrostatic[:,:]
```
%% Cell type:code id: tags:
``` python
def compute_density():
for block in dh.iterate(ghost_layers=True, inner_ghost_layers=False):
# PDFs in incompressible LBM are centered around zero in lbmpy
# => add 1 to whole sum, i.e., initialize with 1
block["rho"].fill(1);
# compute density
for i in range(block["g"].shape[-1]):
block["rho"][:,:] += block["g"][:,:,i]
```
%% Cell type:code id: tags:
``` python
if (gravitational_acceleration != 0):
initialize_hydrostatic_pressure()
compute_density()
plt.scalar_field(dh.cpu_arrays["rho"])
```
%% Output
%% Cell type:code id: tags:
``` python
print("================================= start of the simulation ===================================")
start = time.time()
pbar = tqdm(total=timesteps)
timestep = []
surface_position = []
symmetry_norm = []
mass = []
for i in range(0, timesteps):
sum_c_2 = 0.0
sum_delta_c_2 = 0.0
# write vtk output
if(i % vtk_output_frequency == 0):
if gpu:
dh.to_cpu("C")
dh.to_cpu("u")
dh.to_cpu("rho")
vtk_writer(i)
# extract data (to be written to file)
if(i % data_extract_frequency == 0):
pbar.update(data_extract_frequency)
timestep.append(i)
ny = domain_size[1]
# get index containing phase field value < 0.5
i1 = np.argmax(dh.cpu_arrays["C"][int(domain_size[0] * 0.5), :] < 0.5)
i0 = i1 - 1 # index containing phase field value >= 0.5
f0 = dh.cpu_arrays["C"][int(domain_size[0] * 0.5), i0] # phase field value >= 0.5
f1 = dh.cpu_arrays["C"][int(domain_size[0] * 0.5), i1] # phase field value < 0.5
# coordinate of cell center is index+0.5-1 (0.5 to get to cell center; -1 to remove ghost layer from index)
y0 = i0 + 0.5 - 1
y1 = i1 + 0.5 - 1
#interpolate
surface_position.append( y0 + (y1 - y0) / (f1 - f0) * (0.5 - f0) )
# evaluate symmetry in x-direction
for y in range(0, domain_size[1] - 1):
for x in range(0, domain_size[0] - 1):
if (x >= domain_size[0] * 0.5):
continue
x_mirrored = domain_size[0] - 1 - x;
sum_c_2 += dh.cpu_arrays["C"][x, y]**2
sum_delta_c_2 += (dh.cpu_arrays["C"][x, y] - dh.cpu_arrays["C"][x_mirrored, y])**2
symmetry_norm.append( (sum_delta_c_2 / sum_c_2)**2 )
mass.append(np.sum(dh.cpu_arrays["C"]))
timeloop()
pbar.close()
end = time.time()
sim_time = end - start
print("\n")
print("time needed for the calculation: %4.4f" % sim_time , "seconds")
nrOfCells = np.prod(domain_size)
mlups = nrOfCells * timesteps / sim_time * 1e-6
print("MLUPS: %4.2f" % mlups)
```
%% Output
================================= start of the simulation ===================================
time needed for the calculation: 99.6110 seconds
MLUPS: 5.99
%% Cell type:code id: tags:
``` python
if gpu:
dh.to_cpu("C")
dh.to_cpu("u")
dh.to_cpu("rho")
plt.scalar_field(dh.cpu_arrays["C"])
```
%% Output
<matplotlib.image.AxesImage at 0x12db854f0>
%% Cell type:code id: tags:
``` python
# non-dimensionalize time and surface position
t_nd = [value * wave_frequency for value in timestep]
a_nd = [(value - fluid_depth) / amplitude for value in surface_position]
```
%% Cell type:code id: tags:
``` python
plt.plot(t_nd, mass, color=(0.121, 0.231, 0.4))
plt.show()
```
%% Output
%% Cell type:code id: tags:
``` python
# plot surface position over time
plt.xlabel('$t/-$', fontsize=18)
plt.ylabel('$a/-$', fontsize=18)
plt.grid(color='black', linestyle='-', linewidth=0.1)
plt.plot(t_nd, a_nd, color=(0.121, 0.231, 0.4))
plt.show()
```
%% Output
%% Cell type:code id: tags:
``` python
# plot analytical solutions
x = 0.5 # position where surface elevation is evaluated;
# 0.5 means center of first cell;
# amplitude of analytical solution is at 0 (not in domain center)
t = np.arange(0, timesteps, 1)
a = amplitude * np.exp(-2 * kinematic_viscosity_heavy * wavenumber**2 * t) # damping of wave amplitude
# linear analytical solution
a_anal_lin = a * np.cos(wavenumber * x - wave_frequency * t) + fluid_depth
# non-linear analytical solution
a_anal_nonlin = (a * np.cos(wavenumber * x - wave_frequency * t) +
wavenumber * a**2 * 0.5 *
(1 + 3 / (2 * np.sinh(wavenumber * fluid_depth)**2)) *
1 / np.tanh(wavenumber * fluid_depth) *
np.cos(2 * wavenumber * x - 2 * wave_frequency * t) + fluid_depth)
# non-dimesionalize time and surface position
a_anal_lin_nd = (a_anal_lin - fluid_depth) / amplitude
a_anal_nonlin_nd = (a_anal_nonlin - fluid_depth) / amplitude
t_anal_nd = t * wave_frequency
plt.plot(t_anal_nd, a_anal_lin_nd, label="linear analytical solution")
#plt.plot(t_anal_nd, a_anal_nonlin_nd, label="nonlinear analytical solution")
plt.plot(t_nd, a_nd, label="simulation")
plt.xlabel('$t/-$', fontsize=18)
plt.ylabel('$a/-$', fontsize=18)
plt.legend()
plt.grid()
plt.show()
```
%% Output
%% Cell type:code id: tags:
``` python
# plot symmetry norm over time
plt.xlabel('$t/-$', fontsize=18)
plt.ylabel('$symmetry-norm/-$', fontsize=18)
plt.grid(color='black', linestyle='-', linewidth=0.1)
plt.plot(t_nd, symmetry_norm, color=(0.121, 0.231, 0.4))
plt.show()
```
%% Output
%% Cell type:code id: tags:
``` python
# store results to file
import csv
with open(filename, 'w') as f:
writer = csv.writer(f, delimiter='\t')
writer.writerows(zip(t_nd, a_nd, symmetry_norm))
```
......
......@@ -2,18 +2,23 @@
The cumulant lattice Boltzmann equation in three dimensions: Theory and validation
by Geier, Martin; Schönherr, Martin; Pasquali, Andrea; Krafczyk, Manfred (2015)
Chapter 5.1
:cite:`geier2015` Chapter 5.1
NOTE: for integration tests, the parameter study is greatly shortened, i.e., the runs are shortened in time and
not all resolutions and viscosities are considered. Nevertheless, all values used by Geier et al. are still in
the setup, only commented, and remain ready to be used (check for comments that start with `NOTE`).
"""
import numpy as np
import sympy as sp
import pytest
from pystencils import Target
import sympy as sp
from lbmpy.creationfunctions import update_with_default_parameters
from lbmpy import LatticeBoltzmannStep, LBStencil
from lbmpy.creationfunctions import LBMConfig, LBMOptimisation
from lbmpy.db import LbmpyJsonSerializer
from lbmpy.enums import Method, Stencil
from lbmpy.relaxationrates import (
relaxation_rate_from_lattice_viscosity, relaxation_rate_from_magic_number)
from lbmpy.scenarios import create_fully_periodic_flow
from pystencils import Target, create_data_handling, CreateKernelConfig
def get_exponent_term(l, **kwargs):
......@@ -53,14 +58,13 @@ def get_analytical_solution(l, l_0, u_0, v_0, nu, t):
def plot_y_velocity(vel, **kwargs):
import matplotlib.pyplot as plt
from matplotlib import cm
vel = vel[:, vel.shape[1]//2, :, 1]
vel = vel[:, vel.shape[1] // 2, :, 1]
ranges = [np.arange(n, dtype=float) for n in vel.shape]
x, y = np.meshgrid(*ranges, indexing='ij')
fig = plt.gcf()
ax = fig.gca(projection='3d')
ax.plot_surface(x, y, vel, cmap=cm.coolwarm, rstride=2, cstride=2,
ax.plot_surface(x, y, vel, cmap='coolwarm', rstride=2, cstride=2,
linewidth=0, antialiased=True, **kwargs)
......@@ -71,9 +75,9 @@ def fit_and_get_slope(x_values, y_values):
def get_phase_error(phases, evaluation_interval):
xValues = np.arange(len(phases) * evaluation_interval, step=evaluation_interval)
phaseError = fit_and_get_slope(xValues, phases)
return phaseError
x_values = np.arange(len(phases) * evaluation_interval, step=evaluation_interval)
phase_error = fit_and_get_slope(x_values, phases)
return phase_error
def get_viscosity(abs_values, evaluation_interval, l):
......@@ -92,33 +96,50 @@ def get_amplitude_and_phase(vel_slice):
return fft_abs[max_index], fft_phase[max_index]
def run(l, l_0, u_0, v_0, nu, y_size, periodicity_in_kernel, **kwargs):
def run(l, l_0, u_0, v_0, nu, y_size, lbm_config, lbm_optimisation, config):
inv_result = {
'viscosity_error': np.nan,
'phase_error': np.nan,
'mlups': np.nan,
}
if 'initial_velocity' in kwargs:
del kwargs['initial_velocity']
if lbm_config.initial_velocity:
# no manually specified initial velocity allowed in shear wave setup
lbm_config.initial_velocity = None
print("Running size l=%d nu=%f " % (l, nu), kwargs)
print(f"Running size l={l} nu={nu}")
initial_vel_arr = get_initial_velocity_field(l, l_0, u_0, v_0, y_size)
omega = relaxation_rate_from_lattice_viscosity(nu)
kwargs['relaxation_rates'] = [sp.sympify(r) for r in kwargs['relaxation_rates']]
if 'entropic' not in kwargs or not kwargs['entropic']:
kwargs['relaxation_rates'] = [r.subs(sp.Symbol("omega"), omega) for r in kwargs['relaxation_rates']]
if lbm_config.fourth_order_correction and omega < 1.75:
pytest.skip("Fourth-order correction only allowed for omega >= 1.75.")
lbm_config.relaxation_rates = [sp.sympify(r) for r in lbm_config.relaxation_rates]
lbm_config.relaxation_rates = [r.subs(sp.Symbol("omega"), omega) for r in lbm_config.relaxation_rates]
periodicity_in_kernel = (lbm_optimisation.builtin_periodicity != (False, False, False))
domain_size = initial_vel_arr.shape[:-1]
scenario = create_fully_periodic_flow(initial_vel_arr, periodicity_in_kernel=periodicity_in_kernel, **kwargs)
data_handling = create_data_handling(domain_size, periodicity=not periodicity_in_kernel,
default_ghost_layers=1, parallel=False)
if 'entropic' in kwargs and kwargs['entropic']:
scenario.kernelParams['omega'] = kwargs['relaxation_rates'][0].subs(sp.Symbol("omega"), omega)
scenario = LatticeBoltzmannStep(data_handling=data_handling, name="periodic_scenario", lbm_config=lbm_config,
lbm_optimisation=lbm_optimisation, config=config)
for b in scenario.data_handling.iterate(ghost_layers=False):
np.copyto(b[scenario.velocity_data_name], initial_vel_arr[b.global_slice])
scenario.set_pdf_fields_from_macroscopic_values()
# NOTE: use those values to limit the runtime in integration tests
total_time_steps = 2000 * (l // l_0) ** 2
initial_time_steps = 1100 * (l // l_0) ** 2
eval_interval = 100 * (l // l_0) ** 2
# NOTE: for simulating the real shear-wave scenario from Geier et al. use the following values
# total_time_steps = 20000 * (l // l_0) ** 2
# initial_time_steps = 11000 * (l // l_0) ** 2
# eval_interval = 1000 * (l // l_0) ** 2
total_time_steps = 20000 * (l // l_0) ** 2
initial_time_steps = 11000 * (l // l_0) ** 2
eval_interval = 1000 * (l // l_0) ** 2
scenario.run(initial_time_steps)
if np.isnan(scenario.velocity_slice()).any():
print(" Result", inv_result)
return inv_result
magnitudes = []
......@@ -149,65 +170,81 @@ def run(l, l_0, u_0, v_0, nu, y_size, periodicity_in_kernel, **kwargs):
def create_full_parameter_study():
from pystencils.runhelper import ParameterStudy
params = {
setup_params = {
'l_0': 32,
'u_0': 0.096,
'v_0': 0.1,
'ySize': 1,
'periodicityInKernel': True,
'stencil': 'D3Q27',
'compressible': True,
'y_size': 1
}
ls = [32 * 2 ** i for i in range(0, 5)]
nus = [1e-2, 1e-3, 1e-4, 1e-5]
srt_and_trt_methods = [{'method': method,
'stencil': stencil,
'compressible': comp,
'relaxation_rates': ["omega", str(relaxation_rate_from_magic_number(sp.Symbol("omega")))],
'equilibrium_order': eqOrder,
'continuous_equilibrium': mbEq,
'optimization': {'target': Target.CPU, 'split': True, 'cse_pdfs': True}}
for method in ('srt', 'trt')
for stencil in ('D3Q19', 'D3Q27')
omega, omega_f = sp.symbols("omega, omega_f")
# NOTE: use those values to limit the runtime in integration tests
ls = [32]
nus = [1e-5]
# NOTE: for simulating the real shear-wave scenario from Geier et al. use the following values
# ls = [32 * 2 ** i for i in range(0, 5)]
# nus = [1e-2, 1e-3, 1e-4, 1e-5]
srt_and_trt_methods = [LBMConfig(method=method,
stencil=LBStencil(stencil),
compressible=comp,
relaxation_rates=[omega, str(relaxation_rate_from_magic_number(omega))],
equilibrium_order=eqOrder,
continuous_equilibrium=mbEq)
for method in (Method.SRT, Method.TRT)
for stencil in (Stencil.D3Q19, Stencil.D3Q27)
for comp in (True, False)
for eqOrder in (1, 2, 3)
for mbEq in (True, False)]
methods = [{'method': 'trt', 'relaxation_rates': ["omega", 1]},
{'method': 'mrt3', 'relaxation_rates': ["omega", 1, 1], 'equilibrium_order': 2},
{'method': 'mrt3', 'relaxation_rates': ["omega", 1, 1], 'equilibrium_order': 3},
{'method': 'mrt3', 'cumulant': True, 'relaxation_rates': ["omega", 1, 1], # cumulant
'continuous_equilibrium': True, 'equilibrium_order': 3,
'optimization': {'cse_global': True}},
{'method': 'trt-kbc-n4', 'relaxation_rates': ["omega", "omega_f"], 'entropic': True, # entropic order 2
'continuous_equilibrium': True, 'equilibrium_order': 2},
{'method': 'trt-kbc-n4', 'relaxation_rates': ["omega", "omega_f"], 'entropic': True, # entropic order 3
'continuous_equilibrium': True, 'equilibrium_order': 3},
# entropic cumulant:
{'method': 'trt-kbc-n4', 'relaxation_rates': ["omega", "omega_f"], 'entropic': True,
'cumulant': True, 'continuous_equilibrium': True, 'equilibrium_order': 3},
{'method': 'trt-kbc-n2', 'relaxation_rates': ["omega", "omega_f"], 'entropic': True,
'cumulant': True, 'continuous_equilibrium': True, 'equilibrium_order': 3},
{'method': 'mrt3', 'relaxation_rates': ["omega", "omega_f", "omega_f"], 'entropic': True,
'cumulant': True, 'continuous_equilibrium': True, 'equilibrium_order': 3},
optimization_srt_trt = LBMOptimisation(split=True, cse_pdfs=True)
config = CreateKernelConfig(target=Target.CPU)
methods = [LBMConfig(stencil=LBStencil(Stencil.D3Q27), method=Method.TRT, relaxation_rates=[omega, 1]),
LBMConfig(stencil=LBStencil(Stencil.D3Q27), method=Method.MRT, relaxation_rates=[omega],
equilibrium_order=2),
LBMConfig(stencil=LBStencil(Stencil.D3Q27), method=Method.MRT, relaxation_rates=[omega],
equilibrium_order=3),
LBMConfig(stencil=LBStencil(Stencil.D3Q27), method=Method.CUMULANT, relaxation_rates=[omega], # cumulant
compressible=True, continuous_equilibrium=True, equilibrium_order=3),
LBMConfig(stencil=LBStencil(Stencil.D3Q27), method=Method.CUMULANT, relaxation_rates=[omega], # cumulant with fourth-order correction
compressible=True, continuous_equilibrium=True, fourth_order_correction=0.1),
LBMConfig(stencil=LBStencil(Stencil.D3Q27), method=Method.TRT_KBC_N4, relaxation_rates=[omega, omega_f],
entropic=True, zero_centered=False, continuous_equilibrium=True, equilibrium_order=2),
LBMConfig(stencil=LBStencil(Stencil.D3Q27), method=Method.TRT_KBC_N4, relaxation_rates=[omega, omega_f],
entropic=True, zero_centered=False, continuous_equilibrium=True, equilibrium_order=3),
# entropic cumulant: not supported for the moment
# LBMConfig(method=Method.CUMULANT, relaxation_rates=["omega", "omega_f"], entropic=True, zero_centered=False,
# compressible=True, continuous_equilibrium=True, equilibrium_order=3)
]
parameter_study = ParameterStudy(run, database_connector="shear_wave_db")
parameter_study = ParameterStudy(run, database_connector="shear_wave_db",
serializer_info=('lbmpy_serializer', LbmpyJsonSerializer))
for L in ls:
for nu in nus:
for methodParams in methods + srt_and_trt_methods:
simulation_params = params.copy()
simulation_params.update(methodParams)
if 'optimization' not in simulation_params:
simulation_params['optimization'] = {}
new_params, new_opt_params = update_with_default_parameters(simulation_params,
simulation_params['optimization'], False)
simulation_params = new_params
simulation_params['optimization'] = new_opt_params
simulation_params['L'] = L
for methodParams in srt_and_trt_methods:
simulation_params = setup_params.copy()
simulation_params['lbm_config'] = methodParams
simulation_params['lbm_optimisation'] = optimization_srt_trt
simulation_params['config'] = config
simulation_params['l'] = L
simulation_params['nu'] = nu
l_0 = simulation_params['l_0']
parameter_study.add_run(simulation_params.copy(), weight=(L / l_0) ** 4)
for methodParams in methods:
simulation_params = setup_params.copy()
simulation_params['lbm_config'] = methodParams
simulation_params['lbm_optimisation'] = LBMOptimisation()
simulation_params['config'] = config
simulation_params['l'] = L
simulation_params['nu'] = nu
l_0 = simulation_params['l_0']
parameter_study.add_run(simulation_params.copy(), weight=(L / l_0) ** 4)
......@@ -215,15 +252,25 @@ def create_full_parameter_study():
def test_shear_wave():
pytest.importorskip('pycuda')
pytest.importorskip('cupy')
params = {
'y_size': 1,
'l_0': 32,
'u_0': 0.096,
'v_0': 0.1,
'stencil': 'D3Q19',
'compressible': True,
"optimization": {"target": Target.GPU}
'l': 32,
'nu': 1e-2,
}
run(32, nu=1e-2, equilibrium_order=2, method='srt', y_size=1, periodicity_in_kernel=True,
relaxation_rates=[sp.Symbol("omega"), 5, 5], continuous_equilibrium=True, **params)
kernel_config = CreateKernelConfig(target=Target.GPU)
lbm_config = LBMConfig(method=Method.SRT, relaxation_rates=[sp.Symbol("omega")], stencil=LBStencil(Stencil.D3Q27),
compressible=True, continuous_equilibrium=True, equilibrium_order=2)
run(**params, lbm_config=lbm_config, config=kernel_config, lbm_optimisation=LBMOptimisation())
@pytest.mark.longrun
def test_all_scenarios():
parameter_study = create_full_parameter_study()
parameter_study.run()
Source diff could not be displayed: it is too large. Options to address this: view the blob.
......@@ -415,7 +415,7 @@
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"display_name": "Python 3 (ipykernel)",
"language": "python",
"name": "python3"
},
......@@ -429,7 +429,7 @@
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.8.2"
"version": "3.11.0rc1"
}
},
"nbformat": 4,
......
%% Cell type:code id: tags:
``` python
import pytest
pytest.importorskip('pycuda')
pytest.importorskip('cupy')
```
%% Output
---------------------------------------------------------------------------
Skipped Traceback (most recent call last)
/var/folders/07/0d7kq8fd0sx24cs53zz90_qc0000gp/T/ipykernel_16968/622163826.py in <module>
1 import pytest
----> 2 pytest.importorskip('pycuda')
/opt/local/Library/Frameworks/Python.framework/Versions/3.9/lib/python3.9/site-packages/_pytest/outcomes.py in importorskip(modname, minversion, reason)
210 if reason is None:
211 reason = f"could not import {modname!r}: {exc}"
--> 212 raise Skipped(reason, allow_module_level=True) from None
213 mod = sys.modules[modname]
214 if minversion is None:
Skipped: could not import 'pycuda': No module named 'pycuda'
<module 'cupy' from '/home/markus/.local/lib/python3.11/site-packages/cupy/__init__.py'>
%% Cell type:code id: tags:
``` python
from lbmpy.session import *
from lbmpy.phasefield.n_phase_boyer import *
from lbmpy.phasefield.kerneleqs import *
from lbmpy.phasefield.contact_angle_circle_fitting import *
from scipy.ndimage.filters import gaussian_filter
from scipy.ndimage import gaussian_filter
from pystencils.simp import sympy_cse_on_assignment_list
one = sp.sympify(1)
import pyximport
pyximport.install(language_level=3)
from lbmpy.phasefield.simplex_projection import simplex_projection_2d # NOQA
```
%% Cell type:markdown id: tags:
# Simulation arbitrary surface tension case
%% Cell type:code id: tags:
``` python
n = 4
dx, dt = 1, 1
mobility = 2e-3
domain_size = (150, 150)
ε = one * 4
penalty_factor = 0
stabilization_factor = 10
κ = (one, one/2, one/3, one/4)
sigma_factor = one / 15
σ = sp.ImmutableDenseMatrix(n, n, lambda i,j: sigma_factor* (κ[i] + κ[j]) if i != j else 0 )
```
%% Cell type:code id: tags:
``` python
dh = create_data_handling(domain_size, periodicity=True, default_target=ps.Target.GPU)
c = dh.add_array('c', values_per_cell=n)
c_tmp = dh.add_array_like('c_tmp', 'c')
μ = dh.add_array('mu', values_per_cell=n)
cvec = c.center_vector
μvec = μ.center_vector
```
%% Cell type:code id: tags:
``` python
α, _ = diffusion_coefficients(σ)
f = lambda c: c**2 * ( 1 - c ) **2
a, b = compute_ab(f)
capital_f = capital_f0(cvec, σ) + correction_g(cvec, σ) + stabilization_factor * stabilization_term(cvec, α)
f_bulk = free_energy_bulk(capital_f, b, ε) + penalty_factor * (one - sum(cvec))
f_if = free_energy_interfacial(cvec, σ, a, ε)
f = f_bulk + f_if
```
%% Cell type:code id: tags:
``` python
#f_bulk
```
%% Cell type:code id: tags:
``` python
μ_assignments = mu_kernel(f, cvec, c, μ)
μ_assignments = [Assignment(a.lhs, a.rhs.doit()) for a in μ_assignments]
μ_assignments = sympy_cse_on_assignment_list(μ_assignments)
```
%% Cell type:code id: tags:
``` python
discretize = fd.Discretization2ndOrder(dx=dx, dt=dt)
```
%% Cell type:code id: tags:
``` python
def lapl(e):
return sum(ps.fd.diff(e, d, d) for d in range(dh.dim))
```
%% Cell type:code id: tags:
``` python
rhs = α * μvec
discretized_rhs = [discretize(fd.expand_diff_full( lapl(mobility * rhs_i) + fd.transient(cvec[i], idx=i), functions=μvec))
for i, rhs_i in enumerate(rhs)]
c_assignments = [Assignment(lhs, rhs)
for lhs, rhs in zip(c_tmp.center_vector, discretized_rhs)]
```
%% Cell type:code id: tags:
``` python
#c_assignments
```
%% Cell type:code id: tags:
``` python
μ_sync = dh.synchronization_function(μ.name)
c_sync = dh.synchronization_function(c.name)
optimization = {'cpu_openmp': False, 'cpu_vectorize_info': None}
config = ps.CreateKernelConfig(cpu_openmp=False, target=dh.default_target)
μ_kernel = create_kernel(μ_assignments, config=config).compile()
c_kernel = create_kernel(c_assignments, config=config).compile()
def set_c(slice_obj, values):
for block in dh.iterate(slice_obj):
arr = block[c.name]
arr[..., : ] = values
def smooth():
for block in dh.iterate(ghost_layers=True):
c_arr = block[c.name]
for i in range(n):
gaussian_filter(c_arr[..., i], sigma=2, output=c_arr[..., i])
def time_loop(steps):
dh.all_to_gpu()
for t in range(steps):
c_sync()
dh.run_kernel(μ_kernel)
μ_sync()
dh.run_kernel(c_kernel)
dh.swap(c.name, c_tmp.name)
#simplex_projection_2d(dh.cpu_arrays[c.name])
dh.all_to_cpu()
```
%% Cell type:code id: tags:
``` python
set_c(make_slice[:, :], [0, 0, 0, 0])
set_c(make_slice[:, 0.5:], [1, 0, 0, 0])
set_c(make_slice[:, :0.5], [0, 1, 0, 0])
set_c(make_slice[0.3:0.7, 0.3:0.7], [0, 0, 1, 0])
smooth()
```
%% Cell type:code id: tags:
``` python
#dh.load_all('n_phases_state_size200_stab10.npz')
```
%% Cell type:code id: tags:
``` python
plt.phase_plot(dh.gather_array(c.name))
```
%% Cell type:code id: tags:
``` python
neumann_angles_from_surface_tensions(lambda i, j: float(σ[i, j]))
```
%% Cell type:code id: tags:
``` python
import time
for i in range(10):
start = time.perf_counter()
time_loop(1_000)
end = time.perf_counter()
try:
print(i, end - start, liquid_lens_neumann_angles(dh.gather_array(c.name)))
except Exception:
print(i, end - start, "none found")
```
%% Cell type:code id: tags:
``` python
plt.subplot(1,3,1)
t = dh.gather_array(c.name, make_slice[25, :]).squeeze()
plt.plot(t);
plt.subplot(1,3,2)
plt.phase_plot(dh.gather_array(c.name), linewidth=1)
plt.subplot(1,3,3)
plt.scalar_field(dh.gather_array(μ.name)[:, :, 2])
plt.colorbar();
```
%% Cell type:code id: tags:
``` python
assert not np.isnan(dh.max(c.name))
```
%% Cell type:code id: tags:
``` python
t = dh.gather_array(c.name, make_slice[25, 55:90]).squeeze()
plt.hlines(0.5, 0, 30)
plt.plot(t);
```
......