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
  • accessor_choice
  • csebug
  • fluct_zero_centered
  • fluctuating
  • fluctuating_lb_test
  • fluctuation_test
  • improved_comm
  • master
  • poiseuille
  • test_martin
  • release/0.2.1
  • release/0.2.2
  • release/0.2.3
  • release/0.2.4
  • release/0.2.5
  • release/0.2.6
16 results
Show changes
Showing
with 3781 additions and 70 deletions
Source diff could not be displayed: it is too large. Options to address this: view the blob.
......@@ -7,34 +7,31 @@ This setup is used to illustrate the difference of D3Q19 and D3Q27 for high Reyn
The deficiencies of the D3Q19 model can be resolved by choosing a better equilibrium (D3Q19_new)
"""
import os
import sys
import numpy as np
import sympy as sp
from lbmpy.boundaries import FixedDensity, NoSlip
from lbmpy.creationfunctions import LBMConfig
from lbmpy.enums import Method, Stencil
from lbmpy.geometry import add_pipe_inflow_boundary, add_pipe_walls
from lbmpy.lbstep import LatticeBoltzmannStep
from lbmpy.relaxationrates import relaxation_rate_from_lattice_viscosity
from lbmpy.stencils import LBStencil
from pystencils import create_data_handling
from pystencils.slicing import make_slice, slice_from_direction
def rod_simulation(stencil_name, diameter, parallel=False, entropic=False, re=150,
time_to_simulate=17.0, eval_interval=0.5, write_vtk=True, write_numpy=True,
optimization_params=None):
if optimization_params is None:
optimization_params = {}
time_to_simulate=17.0, eval_interval=0.5, write_vtk=True, write_numpy=True):
if stencil_name == 'D3Q19_new':
stencil = 'D3Q19'
maxwellian_moments = True
stencil = Stencil.D3Q19
continuous_equilibrium = True
elif stencil_name == 'D3Q19_old':
stencil = 'D3Q19'
maxwellian_moments = False
stencil = Stencil.D3Q19
continuous_equilibrium = False
elif stencil_name == 'D3Q27':
stencil = 'D3Q27'
maxwellian_moments = False
stencil = Stencil.D3Q27
continuous_equilibrium = False
else:
raise ValueError("Unknown stencil name " + stencil_name)
d = diameter
......@@ -53,23 +50,15 @@ def rod_simulation(stencil_name, diameter, parallel=False, entropic=False, re=15
lattice_viscosity = u_mean_at_throat * constriction_diameter / re
omega = relaxation_rate_from_lattice_viscosity(lattice_viscosity)
method_parameters = {'stencil': stencil,
'compressible': True,
'method': 'srt',
'relaxation_rates': [omega],
'entropic': entropic,
'maxwellian_moments': maxwellian_moments, }
lbm_config = LBMConfig(stencil=LBStencil(stencil), compressible=True, method=Method.SRT, relaxation_rate=omega,
entropic=entropic, continuous_equilibrium=continuous_equilibrium)
kernel_params = {}
if entropic:
method_parameters['method'] = 'mrt3'
method_parameters['relaxation_rates'] = sp.symbols("omega_fix omega_fix omega_free")
kernel_params['omega_fix'] = omega
print("ω=", omega)
dh = create_data_handling(domain_size=(length, d, d), parallel=parallel)
sc = LatticeBoltzmannStep(data_handling=dh, kernel_params=kernel_params, optimization=optimization_params,
name=stencil_name, **method_parameters)
sc = LatticeBoltzmannStep(data_handling=dh, kernel_params=kernel_params,
name=stencil_name, lbm_config=lbm_config)
# ----------------- Boundary Setup --------------------------------
......@@ -100,8 +89,6 @@ def rod_simulation(stencil_name, diameter, parallel=False, entropic=False, re=15
# ----------------- Run --------------------------------------------
scenario_name = stencil_name
if entropic:
scenario_name += "_entropic"
if not os.path.exists(scenario_name):
os.mkdir(scenario_name)
......@@ -124,15 +111,15 @@ def rod_simulation(stencil_name, diameter, parallel=False, entropic=False, re=15
if np.isnan(max_vel):
raise ValueError("Simulation went unstable")
if write_numpy:
np.save(scenario_name + '/velocity_%06d' % (sc.time_steps_run,), sc.velocity_slice().filled(0.0))
np.save(scenario_name + f'/velocity_{sc.time_steps_run}', sc.velocity_slice().filled(0.0))
def test_rod_scenario_simple():
rod_simulation("D3Q19_new", re=100, diameter=20, parallel=False, entropic=False,
time_to_simulate=0.2, eval_interval=0.1, write_vtk=False, write_numpy=False)
if __name__ == '__main__':
# High Re (Entropic method)
rod_simulation(stencil_name=sys.argv[1], re=500, diameter=80, entropic=True, time_to_simulate=17,
parallel=False, optimization_params={'target': 'gpu'})
# depricated usage
# if __name__ == '__main__':
# # High Re (Entropic method)
# rod_simulation(stencil_name=sys.argv[1], re=500, diameter=80, entropic=True, time_to_simulate=17,
# parallel=False, optimization_params={'target': Target.GPU})
......@@ -8,7 +8,9 @@ a cylinder. In Flow simulation with high-performance computers II (pp. 547-566).
import warnings
import numpy as np
import pytest
from pystencils.backends.simd_instruction_sets import get_supported_instruction_sets
from lbmpy.boundaries.boundaryconditions import NoSlip
from lbmpy.geometry import get_pipe_velocity_field
from lbmpy.relaxationrates import relaxation_rate_from_lattice_viscosity
......@@ -122,7 +124,7 @@ def schaefer_turek_2d(cells_per_diameter, u_max=0.3, max_lattice_velocity=0.05,
(u_max, domain_size, dx, dt, omega, re_lattice))
initial_velocity = get_pipe_velocity_field(domain_size, max_lattice_velocity)
scenario = create_channel(domain_size=domain_size, u_max=max_lattice_velocity, kernel_params={'omega_0': omega},
scenario = create_channel(domain_size=domain_size, u_max=max_lattice_velocity, relaxation_rate=omega,
initial_velocity=initial_velocity, **kwargs)
scenario.boundary_handling.set_boundary(NoSlip('obstacle'), mask_callback=geometry_callback)
parameter_info['u_bar_l'] = u_bar_l
......@@ -146,8 +148,9 @@ def long_run(steady=True, **kwargs):
plt.show()
@pytest.mark.skipif(not get_supported_instruction_sets(), reason='cannot detect CPU instruction set')
def test_schaefer_turek():
opt = {'vectorization': {'instruction_set': 'avx', 'assume_aligned': True}, 'openmp': 2}
opt = {'vectorization': {'instruction_set': get_supported_instruction_sets()[-1], 'assume_aligned': True}, 'openmp': 2}
sc_2d_1 = schaefer_turek_2d(30, max_lattice_velocity=0.08, optimization=opt)
sc_2d_1.run(30000)
result = evaluate_static_quantities(sc_2d_1)
......@@ -156,5 +159,5 @@ def test_schaefer_turek():
if __name__ == '__main__':
long_run(entropic=True, method='trt-kbc-n1', compressible=True,
optimization={'target': 'gpu', 'gpuIndexingParams': {'blockSize': (16, 8, 2)}})
long_run(entropic=True, method='trt_kbc_n1', compressible=True,
optimization={'target': Target.GPU, 'gpuIndexingParams': {'blockSize': (16, 8, 2)}})
......@@ -2,15 +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 pytest
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):
......@@ -50,27 +58,26 @@ 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)
def fit_and_get_slope(x_values, y_values):
matrix = np.vstack([x_values, np.ones(len(x_values))]).T
m, _ = np.linalg.lstsq(matrix, y_values, rcond=None)[0]
m, _ = np.linalg.lstsq(matrix, y_values, rcond=1e-14)[0]
return m
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):
......@@ -89,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 = []
......@@ -146,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,
'maxwellian_moments': mbEq,
'optimization': {'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
'maxwellian_moments': True, 'equilibrium_order': 3,
'optimization': {'cse_global': True}},
{'method': 'trt-kbc-n4', 'relaxation_rates': ["omega", "omega_f"], 'entropic': True, # entropic order 2
'maxwellian_moments': True, 'equilibrium_order': 2},
{'method': 'trt-kbc-n4', 'relaxation_rates': ["omega", "omega_f"], 'entropic': True, # entropic order 3
'maxwellian_moments': True, 'equilibrium_order': 3},
# entropic cumulant:
{'method': 'trt-kbc-n4', 'relaxation_rates': ["omega", "omega_f"], 'entropic': True,
'cumulant': True, 'maxwellian_moments': True, 'equilibrium_order': 3},
{'method': 'trt-kbc-n2', 'relaxation_rates': ["omega", "omega_f"], 'entropic': True,
'cumulant': True, 'maxwellian_moments': True, 'equilibrium_order': 3},
{'method': 'mrt3', 'relaxation_rates': ["omega", "omega_f", "omega_f"], 'entropic': True,
'cumulant': True, 'maxwellian_moments': 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)
......@@ -212,14 +252,25 @@ def create_full_parameter_study():
def test_shear_wave():
pytest.importorskip('cupy')
params = {
'y_size': 1,
'l_0': 32,
'u_0': 0.096,
'v_0': 0.1,
'stencil': 'D3Q19',
'compressible': True,
"optimization": {"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], maxwellian_moments=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()
......@@ -6,29 +6,29 @@ Silva, Semiao
python3 scenario_square_channel.py server
python3 scenario_square_channel.py client --host i10staff41 -P '{ "optimization" : { "target" : "gpu"} }'
python3 scenario_square_channel.py client --host i10staff41 -P '{ "optimization" : { "target" : Target.GPU} }'
"""
import numpy as np
import sympy as sp
from lbmpy import Stencil, Method, ForceModel, LBStencil
from lbmpy.methods.creationfunctions import relaxation_rate_from_magic_number
from lbmpy.scenarios import create_channel
from pystencils import make_slice
defaultParameters = {
'stencil': 'D3Q19',
'method': 'srt',
'stencil': LBStencil(Stencil.D3Q19),
'method': Method.SRT,
'lambda_plus_sq': 4 / 25,
'square_size': 15,
'quadratic': True,
're': 10,
'compressible': True,
'cumulant': False,
'maxwellian_moments': False,
'continuous_equilibrium': False,
'equilibrium_order': 2,
'force_model': 'guo',
'force_model': ForceModel.GUO,
'c_s_sq': 1 / 3,
'analytic_initial_velocity': False,
......@@ -119,7 +119,7 @@ def run(convergence_check_interval=1000, convergence_threshold=1e-12, plot_resul
params['relaxation_rates'] = [omega, relaxation_rate_from_magic_number(omega, 3 / 16)]
stencil = params['stencil']
viscosity_factor = 1 / 2 if stencil == 'D3Q15' and params['maxwellian_moments'] else 1 / 3
viscosity_factor = 1 / 2 if stencil == LBStencil(Stencil.D3Q15) and params['continuous_equilibrium'] else 1 / 3
print("Running size %d quadratic %d analyticInit %d " %
(square_size, quadratic, analytic_initial_velocity) + str(params))
......@@ -209,9 +209,9 @@ def run(convergence_check_interval=1000, convergence_threshold=1e-12, plot_resul
def parameter_filter(p):
if p.cumulant and p.compressible:
return None
if p.cumulant and not p.maxwellian_moments:
if p.cumulant and not p.continuous_equilibrium:
return None
if p.cumulant and p.stencil == 'D3Q15':
if p.cumulant and p.stencil == LBStencil(Stencil.D3Q15):
return None
if not p.quadratic and not p.reynolds_nr_accuracy:
# analytical formula not valid for rectangular channel
......@@ -232,14 +232,13 @@ def small_study():
common_degrees_of_freedom = [
('reynolds_nr_accuracy', [1e-8, None]),
('analytic_initial_velocity', [True]),
('force_model', ['luo']),
('cumulant', [False]),
('method', ['srt']),
('force_model', [ForceModel.LUO]),
('method', [Method.SRT]),
('equilibrium_order', [2]),
('stencil', ['D3Q19']),
('stencil', [LBStencil(Stencil.D3Q19)]),
('compressible', [True]),
('quadratic', [True, False]),
('maxwellian_moments', [False, True]),
('continuous_equilibrium', [False, True]),
]
convergence_study = common_degrees_of_freedom.copy()
convergence_study += [('square_size', [15, 25, 35, 45, 53])]
......@@ -258,13 +257,13 @@ def create_full_parameter_study():
('cumulant', [False]),
('compressible', [False, True]),
('reynolds_nr_accuracy', [None, 1e-8]),
('stencil', ['D3Q19', 'D3Q15']),
('stencil', [LBStencil(Stencil.D3Q19), LBStencil(Stencil.D3Q15)]),
('analytic_initial_velocity', [False]),
('force_model', ['guo', 'simple', 'silva', 'luo']),
('method', ['srt', 'trt']),
('force_model', [ForceModel.GUO, ForceModel.SIMPLE, ForceModel.SILVA, ForceModel.LUO]),
('method', [Method.SRT, Method.TRT]),
('equilibrium_order', [2, 3]),
('quadratic', [True, False]),
('maxwellian_moments', [False, True]),
('continuous_equilibrium', [False, True]),
('use_mean_for_reynolds', [True, False]),
]
......@@ -292,12 +291,12 @@ def d3q15_cs_sq_half_study():
('compressible', [False, True]),
('reynolds_nr_accuracy', [None, ]),
('analytic_initial_velocity', [False]),
('force_model', ['guo', 'silva']),
('method', ['srt', 'trt']),
('force_model', [ForceModel.GUO, ForceModel.SILVA]),
('method', [Method.SRT, Method.TRT]),
('equilibrium_order', [2, 3]),
('stencil', ['D3Q15']),
('stencil', [LBStencil(Stencil.D3Q15)]),
('quadratic', [True, ]),
('maxwellian_moments', [True, ]),
('continuous_equilibrium', [True, ]),
('c_s_sq', [1 / 3]),
('square_size', [45, 85]),
]
......@@ -313,11 +312,11 @@ def d3q27_study():
('compressible', [False]),
('reynolds_nr_accuracy', [None, ]),
('analytic_initial_velocity', [False]),
('force_model', ['guo', 'silva']),
('method', ['srt']),
('force_model', [ForceModel.GUO, ForceModel.SILVA]),
('method', [Method.SRT]),
('equilibrium_order', [2]),
('stencil', ['D3Q27']),
('maxwellian_moments', [True, ]),
('stencil', [LBStencil(Stencil.D3Q27)]),
('continuous_equilibrium', [True, ]),
('c_s_sq', [1 / 3]),
('square_size', [15, 25, 35, 45, 53, 85, 135]),
('use_mean_for_reynolds', [False]),
......@@ -329,13 +328,17 @@ def d3q27_study():
def test_square_channel():
res = run(convergence_check_interval=1000, convergence_threshold=1e-5, plot_result=False, lambda_plus_sq=4 / 25,
re=10, square_size=53, quadratic=True, analytic_initial_velocity=False, reynolds_nr_accuracy=None,
force_model='buick', stencil='D3Q19', maxwellian_moments=False, equilibrium_order=2, compressible=True)
force_model=ForceModel.BUICK, stencil=LBStencil(Stencil.D3Q19),
continuous_equilibrium=False, equilibrium_order=2, compressible=True)
assert 1e-5 < res['normalized_spurious_vel_max'] < 1.2e-5
# TODO test again if compressible works when !113 is merged
res = run(convergence_check_interval=1000, convergence_threshold=1e-5, plot_result=False, lambda_plus_sq=4 / 25,
re=10, square_size=53, quadratic=True, analytic_initial_velocity=False, reynolds_nr_accuracy=None,
force_model='buick', stencil='D3Q19', maxwellian_moments=True, equilibrium_order=2, compressible=True)
force_model=ForceModel.BUICK, stencil=LBStencil(Stencil.D3Q19),
continuous_equilibrium=True, equilibrium_order=2, compressible=False)
assert res['normalized_spurious_vel_max'] < 1e-14
......
%% Cell type:code id: tags:
``` python
from lbmpy.session import *
from lbmpy.phasefield.analytical import *
from pystencils.fd import evaluate_diffs
```
%% Cell type:markdown id: tags:
# Analytical checks for 3-Phase model
Here you can inspect the components of the free energy. View only bulk or interface contributions, before and after transformation from $U \rightarrow (\rho, \phi,\psi)$:
%% Cell type:code id: tags:
``` python
order_parameters = sp.symbols("rho phi psi")
rho, phi, psi = order_parameters
F, _ = free_energy_functional_3_phases(include_bulk=True,
include_interface=True,
transformed=True,
expand_derivatives=False)
F
```
%% Output
$\displaystyle \frac{\alpha^{2} \kappa_{0} {\partial (\frac{\phi}{2} - \frac{\psi}{2} + \frac{\rho}{2}) }^{2}}{2} + \frac{\alpha^{2} \kappa_{1} {\partial (- \frac{\phi}{2} - \frac{\psi}{2} + \frac{\rho}{2}) }^{2}}{2} + \frac{\alpha^{2} \kappa_{2} {\partial \psi}^{2}}{2} + \frac{\kappa_{0} \left(\frac{\phi}{2} - \frac{\psi}{2} + \frac{\rho}{2}\right)^{2} \left(- \frac{\phi}{2} + \frac{\psi}{2} - \frac{\rho}{2} + 1\right)^{2}}{2} + \frac{\kappa_{1} \left(- \frac{\phi}{2} - \frac{\psi}{2} + \frac{\rho}{2}\right)^{2} \left(\frac{\phi}{2} + \frac{\psi}{2} - \frac{\rho}{2} + 1\right)^{2}}{2} + \frac{\kappa_{2} \psi^{2} \left(1 - \psi\right)^{2}}{2}$
2 2 2 2 2
α ⋅κ₀⋅D(phi/2 - psi/2 + rho/2) α ⋅κ₁⋅D(-phi/2 - psi/2 + rho/2) α ⋅κ₂⋅D(p
─────────────────────────────── + ──────────────────────────────── + ─────────
2 2 2
2 2 2 2
⎛φ ψ ρ⎞ ⎛ φ ψ ρ ⎞ ⎛ φ ψ ρ⎞ ⎛φ ψ ρ ⎞
2 κ₀⋅⎜─ - ─ + ─⎟ ⋅⎜- ─ + ─ - ─ + 1⎟ κ₁⋅⎜- ─ - ─ + ─⎟ ⋅⎜─ + ─ - ─ + 1⎟
si) ⎝2 2 2⎠ ⎝ 2 2 2 ⎠ ⎝ 2 2 2⎠ ⎝2 2 2 ⎠
──── + ────────────────────────────────── + ──────────────────────────────────
2 2
2 2
κ₂⋅ψ ⋅(1 - ψ)
+ ──────────────
2
%% Cell type:markdown id: tags:
### Analytically checking the phase transition profile
%% Cell type:markdown id: tags:
Automatically deriving chemical potential as functional derivative of free energy
%% Cell type:code id: tags:
``` python
F, _ = free_energy_functional_3_phases(order_parameters)
mu_diff_eq = chemical_potentials_from_free_energy(F, order_parameters)
mu_diff_eq[0]
```
%% Output
$\displaystyle - \frac{\alpha^{2} \kappa_{0} {\partial {\partial \phi}}}{4} + \frac{\alpha^{2} \kappa_{0} {\partial {\partial \psi}}}{4} - \frac{\alpha^{2} \kappa_{0} {\partial {\partial \rho}}}{4} + \frac{\alpha^{2} \kappa_{1} {\partial {\partial \phi}}}{4} + \frac{\alpha^{2} \kappa_{1} {\partial {\partial \psi}}}{4} - \frac{\alpha^{2} \kappa_{1} {\partial {\partial \rho}}}{4} + \frac{\kappa_{0} \phi^{3}}{8} - \frac{3 \kappa_{0} \phi^{2} \psi}{8} + \frac{3 \kappa_{0} \phi^{2} \rho}{8} - \frac{3 \kappa_{0} \phi^{2}}{8} + \frac{3 \kappa_{0} \phi \psi^{2}}{8} - \frac{3 \kappa_{0} \phi \psi \rho}{4} + \frac{3 \kappa_{0} \phi \psi}{4} + \frac{3 \kappa_{0} \phi \rho^{2}}{8} - \frac{3 \kappa_{0} \phi \rho}{4} + \frac{\kappa_{0} \phi}{4} - \frac{\kappa_{0} \psi^{3}}{8} + \frac{3 \kappa_{0} \psi^{2} \rho}{8} - \frac{3 \kappa_{0} \psi^{2}}{8} - \frac{3 \kappa_{0} \psi \rho^{2}}{8} + \frac{3 \kappa_{0} \psi \rho}{4} - \frac{\kappa_{0} \psi}{4} + \frac{\kappa_{0} \rho^{3}}{8} - \frac{3 \kappa_{0} \rho^{2}}{8} + \frac{\kappa_{0} \rho}{4} - \frac{\kappa_{1} \phi^{3}}{8} - \frac{3 \kappa_{1} \phi^{2} \psi}{8} + \frac{3 \kappa_{1} \phi^{2} \rho}{8} - \frac{3 \kappa_{1} \phi^{2}}{8} - \frac{3 \kappa_{1} \phi \psi^{2}}{8} + \frac{3 \kappa_{1} \phi \psi \rho}{4} - \frac{3 \kappa_{1} \phi \psi}{4} - \frac{3 \kappa_{1} \phi \rho^{2}}{8} + \frac{3 \kappa_{1} \phi \rho}{4} - \frac{\kappa_{1} \phi}{4} - \frac{\kappa_{1} \psi^{3}}{8} + \frac{3 \kappa_{1} \psi^{2} \rho}{8} - \frac{3 \kappa_{1} \psi^{2}}{8} - \frac{3 \kappa_{1} \psi \rho^{2}}{8} + \frac{3 \kappa_{1} \psi \rho}{4} - \frac{\kappa_{1} \psi}{4} + \frac{\kappa_{1} \rho^{3}}{8} - \frac{3 \kappa_{1} \rho^{2}}{8} + \frac{\kappa_{1} \rho}{4}$
2 2 2 2 2
α ⋅κ₀⋅D(D(phi)) α ⋅κ₀⋅D(D(psi)) α ⋅κ₀⋅D(D(rho)) α ⋅κ₁⋅D(D(phi)) α ⋅κ
- ─────────────── + ─────────────── - ─────────────── + ─────────────── + ────
4 4 4 4
2 3 2 2 2
₁⋅D(D(psi)) α ⋅κ₁⋅D(D(rho)) κ₀⋅φ 3⋅κ₀⋅φ ⋅ψ 3⋅κ₀⋅φ ⋅ρ 3⋅κ₀⋅φ 3⋅κ₀
─────────── - ─────────────── + ───── - ───────── + ───────── - ─────── + ────
4 4 8 8 8 8
2 2 3 2
⋅φ⋅ψ 3⋅κ₀⋅φ⋅ψ⋅ρ 3⋅κ₀⋅φ⋅ψ 3⋅κ₀⋅φ⋅ρ 3⋅κ₀⋅φ⋅ρ κ₀⋅φ κ₀⋅ψ 3⋅κ₀⋅ψ ⋅
───── - ────────── + ──────── + ───────── - ──────── + ──── - ───── + ────────
8 4 4 8 4 4 8 8
2 2 3 2 3
ρ 3⋅κ₀⋅ψ 3⋅κ₀⋅ψ⋅ρ 3⋅κ₀⋅ψ⋅ρ κ₀⋅ψ κ₀⋅ρ 3⋅κ₀⋅ρ κ₀⋅ρ κ₁⋅φ 3
─ - ─────── - ───────── + ──────── - ──── + ───── - ─────── + ──── - ───── - ─
8 8 4 4 8 8 4 8
2 2 2 2 2
⋅κ₁⋅φ ⋅ψ 3⋅κ₁⋅φ ⋅ρ 3⋅κ₁⋅φ 3⋅κ₁⋅φ⋅ψ 3⋅κ₁⋅φ⋅ψ⋅ρ 3⋅κ₁⋅φ⋅ψ 3⋅κ₁⋅φ⋅ρ
──────── + ───────── - ─────── - ───────── + ────────── - ──────── - ─────────
8 8 8 8 4 4 8
3 2 2 2
3⋅κ₁⋅φ⋅ρ κ₁⋅φ κ₁⋅ψ 3⋅κ₁⋅ψ ⋅ρ 3⋅κ₁⋅ψ 3⋅κ₁⋅ψ⋅ρ 3⋅κ₁⋅ψ⋅ρ κ₁⋅ψ
+ ──────── - ──── - ───── + ───────── - ─────── - ───────── + ──────── - ────
4 4 8 8 8 8 4 4
3 2
κ₁⋅ρ 3⋅κ₁⋅ρ κ₁⋅ρ
+ ───── - ─────── + ────
8 8 4
%% Cell type:markdown id: tags:
Checking if expected profile is a solution of the differential equation:
%% Cell type:code id: tags:
``` python
x = sp.symbols("x")
expectedProfile = analytic_interface_profile(x)
expectedProfile
```
%% Output
$\displaystyle \frac{\tanh{\left(\frac{x}{2 \alpha} \right)}}{2} + \frac{1}{2}$
⎛ x ⎞
tanh⎜───⎟
⎝2⋅α⎠ 1
───────── + ─
2 2
%% Cell type:markdown id: tags:
Checking a phase transition from $C_0$ to $C_2$. This means that $\rho=1$ while $phi$ and $psi$ are the analytical profile or 1-analytical profile
%% Cell type:code id: tags:
``` python
for eq in mu_diff_eq:
eq = eq.subs({rho: 1,
phi: 1 - expectedProfile,
psi: expectedProfile})
eq = evaluate_diffs(eq, x).expand()
assert eq == 0
```
%% Cell type:markdown id: tags:
### Checking the surface tensions parameters
%% Cell type:code id: tags:
``` python
F, _ = free_energy_functional_3_phases(order_parameters)
F = expand_diff_linear(F, functions=order_parameters) # expand derivatives using product rule
two_phase_free_energy = F.subs({rho: 1,
phi: 1 - expectedProfile,
psi: expectedProfile})
two_phase_free_energy = sp.simplify(evaluate_diffs(two_phase_free_energy, x))
two_phase_free_energy
```
%% Output
$\displaystyle \frac{\kappa_{0} + \kappa_{2}}{16 \cosh^{4}{\left(\frac{x}{2 \alpha} \right)}}$
κ₀ + κ₂
─────────────
4⎛ x ⎞
16⋅cosh ⎜───⎟
⎝2⋅α⎠
%% Cell type:code id: tags:
``` python
gamma = cosh_integral(two_phase_free_energy, x)
gamma
```
%% Output
$\displaystyle \frac{\alpha \left(\kappa_{0} + \kappa_{2}\right)}{6}$
α⋅(κ₀ + κ₂)
───────────
6
%% Cell type:code id: tags:
``` python
alpha, k0, k2 = sp.symbols("alpha, kappa_0, kappa_2")
assert gamma == alpha/6 * (k0 + k2)
```
%% Cell type:code id: tags:
``` python
from lbmpy.session import *
from lbmpy.phasefield.analytical import *
from functools import partial
```
%% Cell type:markdown id: tags:
# Analytical checks for N-Phase model
%% Cell type:markdown id: tags:
### Formulation of free energy
%% Cell type:markdown id: tags:
In the next cell you can inspect the formulation of the free energy functional.
Bulk and interface terms can be viewed separately.
%% Cell type:code id: tags:
``` python
num_phases = 3
order_params = symbolic_order_parameters(num_symbols=num_phases-1)
f2 = partial(n_phases_correction_function, beta=1)
F = free_energy_functional_n_phases(order_parameters=order_params,
include_interface=False, f2=f2,
include_bulk=True )
F
```
%% Output
$\displaystyle \frac{3 \left(\tau_{0 1} \left(\phi_{0}^{2} \left(1 - \phi_{0}\right)^{2} + \phi_{1}^{2} \left(1 - \phi_{1}\right)^{2} - \begin{cases} - \left(\phi_{0} + \phi_{1}\right)^{2} & \text{for}\: \phi_{0} + \phi_{1} < 0 \\- \left(- \phi_{0} - \phi_{1} + 1\right)^{2} & \text{for}\: \phi_{0} + \phi_{1} > 1 \\\left(\phi_{0} + \phi_{1}\right)^{2} \left(- \phi_{0} - \phi_{1} + 1\right)^{2} & \text{otherwise} \end{cases}\right) + \tau_{0 2} \left(\phi_{0}^{2} \left(1 - \phi_{0}\right)^{2} + \left(\phi_{0} + \phi_{1}\right)^{2} \left(- \phi_{0} - \phi_{1} + 1\right)^{2} - \begin{cases} - \left(1 - \phi_{1}\right)^{2} & \text{for}\: \phi_{1} - 1 > 0 \\- \phi_{1}^{2} & \text{for}\: \phi_{1} - 1 < -1 \\\phi_{1}^{2} \left(1 - \phi_{1}\right)^{2} & \text{otherwise} \end{cases}\right) + \tau_{1 2} \left(\phi_{1}^{2} \left(1 - \phi_{1}\right)^{2} + \left(\phi_{0} + \phi_{1}\right)^{2} \left(- \phi_{0} - \phi_{1} + 1\right)^{2} - \begin{cases} - \left(1 - \phi_{0}\right)^{2} & \text{for}\: \phi_{0} - 1 > 0 \\- \phi_{0}^{2} & \text{for}\: \phi_{0} - 1 < -1 \\\phi_{0}^{2} \left(1 - \phi_{0}\right)^{2} & \text{otherwise} \end{cases}\right)\right)}{2 \alpha}$
⎛ ⎛ ⎛⎧ 2
⎜ ⎜ ⎜⎪ -(φ₀ + φ₁) for φ₀
⎜ ⎜ ⎜⎪
⎜ ⎜ 2 2 2 2 ⎜⎪ 2
3⋅⎜τ₀ ₁⋅⎜φ₀ ⋅(1 - φ₀) + φ₁ ⋅(1 - φ₁) - ⎜⎨ -(-φ₀ - φ₁ + 1) for φ₀
⎜ ⎜ ⎜⎪
⎜ ⎜ ⎜⎪ 2 2
⎜ ⎜ ⎜⎪(φ₀ + φ₁) ⋅(-φ₀ - φ₁ + 1) othe
⎝ ⎝ ⎝⎩
──────────────────────────────────────────────────────────────────────────────
⎞⎞ ⎛ ⎛⎧ 2
+ φ₁ < 0⎟⎟ ⎜ ⎜⎪ -(1 - φ₁)
⎟⎟ ⎜ ⎜⎪
⎟⎟ ⎜ 2 2 2 2 ⎜⎪ 2
+ φ₁ > 1⎟⎟ + τ₀ ₂⋅⎜φ₀ ⋅(1 - φ₀) + (φ₀ + φ₁) ⋅(-φ₀ - φ₁ + 1) - ⎜⎨ -φ₁
⎟⎟ ⎜ ⎜⎪
⎟⎟ ⎜ ⎜⎪ 2
rwise ⎟⎟ ⎜ ⎜⎪φ₁ ⋅(1 - φ₁)
⎠⎠ ⎝ ⎝⎩
──────────────────────────────────────────────────────────────────────────────
2⋅α
⎞⎞ ⎛ ⎛⎧
for φ₁ - 1 > 0 ⎟⎟ ⎜ ⎜⎪ -
⎟⎟ ⎜ ⎜⎪
⎟⎟ ⎜ 2 2 2 2 ⎜⎪
for φ₁ - 1 < -1⎟⎟ + τ₁ ₂⋅⎜φ₁ ⋅(1 - φ₁) + (φ₀ + φ₁) ⋅(-φ₀ - φ₁ + 1) - ⎜⎨
⎟⎟ ⎜ ⎜⎪
2 ⎟⎟ ⎜ ⎜⎪
otherwise ⎟⎟ ⎜ ⎜⎪φ₀
⎠⎠ ⎝ ⎝⎩
──────────────────────────────────────────────────────────────────────────────
2 ⎞⎞⎞
(1 - φ₀) for φ₀ - 1 > 0 ⎟⎟⎟
⎟⎟⎟
2 ⎟⎟⎟
-φ₀ for φ₀ - 1 < -1⎟⎟⎟
⎟⎟⎟
2 2 ⎟⎟⎟
⋅(1 - φ₀) otherwise ⎟⎟⎟
⎠⎠⎠
───────────────────────────────
%% Cell type:markdown id: tags:
### Analytically checking the phase transition profile
First we define the order parameters and free energy, including bulk and interface terms:
%% Cell type:code id: tags:
``` python
F = free_energy_functional_n_phases(order_parameters=symbolic_order_parameters(num_symbols=num_phases-1))
```
%% Cell type:markdown id: tags:
Then we automatically derive the differential equations for the chemial potential $\mu$
%% Cell type:code id: tags:
``` python
mu_diff_eq = chemical_potentials_from_free_energy(F, order_params)
# there is one equation less than phases
assert len(mu_diff_eq) == num_phases - 1
# show the first one
mu_diff_eq[0]
```
%% Output
$\displaystyle 3 \alpha \tau_{0 1} {\partial {\partial \phi_{1}}} - 6 \alpha \tau_{0 2} {\partial {\partial \phi_{0}}} - 3 \alpha \tau_{0 2} {\partial {\partial \phi_{1}}} - 3 \alpha \tau_{1 2} {\partial {\partial \phi_{1}}} + \frac{12 \phi_{0}^{3} \tau_{0 2}}{\alpha} - \frac{18 \phi_{0}^{2} \phi_{1} \tau_{0 1}}{\alpha} + \frac{18 \phi_{0}^{2} \phi_{1} \tau_{0 2}}{\alpha} + \frac{18 \phi_{0}^{2} \phi_{1} \tau_{1 2}}{\alpha} - \frac{18 \phi_{0}^{2} \tau_{0 2}}{\alpha} - \frac{18 \phi_{0} \phi_{1}^{2} \tau_{0 1}}{\alpha} + \frac{18 \phi_{0} \phi_{1}^{2} \tau_{0 2}}{\alpha} + \frac{18 \phi_{0} \phi_{1}^{2} \tau_{1 2}}{\alpha} + \frac{18 \phi_{0} \phi_{1} \tau_{0 1}}{\alpha} - \frac{18 \phi_{0} \phi_{1} \tau_{0 2}}{\alpha} - \frac{18 \phi_{0} \phi_{1} \tau_{1 2}}{\alpha} + \frac{6 \phi_{0} \tau_{0 2}}{\alpha} - \frac{6 \phi_{1}^{3} \tau_{0 1}}{\alpha} + \frac{6 \phi_{1}^{3} \tau_{0 2}}{\alpha} + \frac{6 \phi_{1}^{3} \tau_{1 2}}{\alpha} + \frac{9 \phi_{1}^{2} \tau_{0 1}}{\alpha} - \frac{9 \phi_{1}^{2} \tau_{0 2}}{\alpha} - \frac{9 \phi_{1}^{2} \tau_{1 2}}{\alpha} - \frac{3 \phi_{1} \tau_{0 1}}{\alpha} + \frac{3 \phi_{1} \tau_{0 2}}{\alpha} + \frac{3 \phi_{1} \tau_{1 2}}{\alpha}$
3⋅α⋅τ₀ ₁⋅D(D(phi_1)) - 6⋅α⋅τ₀ ₂⋅D(D(phi_0)) - 3⋅α⋅τ₀ ₂⋅D(D(phi_1)) - 3⋅α⋅τ₁ ₂⋅
3 2 2 2
12⋅φ₀ ⋅τ₀ ₂ 18⋅φ₀ ⋅φ₁⋅τ₀ ₁ 18⋅φ₀ ⋅φ₁⋅τ₀ ₂ 18⋅φ₀ ⋅φ₁⋅τ₁ ₂
D(D(phi_1)) + ─────────── - ────────────── + ────────────── + ────────────── -
α α α α
2 2 2 2
18⋅φ₀ ⋅τ₀ ₂ 18⋅φ₀⋅φ₁ ⋅τ₀ ₁ 18⋅φ₀⋅φ₁ ⋅τ₀ ₂ 18⋅φ₀⋅φ₁ ⋅τ₁ ₂ 18⋅φ₀⋅φ₁⋅τ₀
─────────── - ────────────── + ────────────── + ────────────── + ────────────
α α α α α
3 3
₁ 18⋅φ₀⋅φ₁⋅τ₀ ₂ 18⋅φ₀⋅φ₁⋅τ₁ ₂ 6⋅φ₀⋅τ₀ ₂ 6⋅φ₁ ⋅τ₀ ₁ 6⋅φ₁ ⋅τ₀ ₂ 6⋅φ₁
─ - ───────────── - ───────────── + ───────── - ────────── + ────────── + ────
α α α α α
3 2 2 2
⋅τ₁ ₂ 9⋅φ₁ ⋅τ₀ ₁ 9⋅φ₁ ⋅τ₀ ₂ 9⋅φ₁ ⋅τ₁ ₂ 3⋅φ₁⋅τ₀ ₁ 3⋅φ₁⋅τ₀ ₂ 3⋅φ₁⋅τ
────── + ────────── - ────────── - ────────── - ───────── + ───────── + ──────
α α α α α α α
₁ ₂
───
%% Cell type:markdown id: tags:
Next we check, that the $tanh$ is indeed a solution of this equation...
%% Cell type:code id: tags:
``` python
x = sp.symbols("x")
expected_profile = analytic_interface_profile(x)
expected_profile
```
%% Output
$\displaystyle \frac{\tanh{\left(\frac{x}{2 \alpha} \right)}}{2} + \frac{1}{2}$
⎛ x ⎞
tanh⎜───⎟
⎝2⋅α⎠ 1
───────── + ─
2 2
%% Cell type:markdown id: tags:
... by inserting it for $\phi_0$, and setting all other order parameters to zero
%% Cell type:code id: tags:
``` python
# zero other parameters
diff_eq = mu_diff_eq[0].subs({p: 0 for p in order_params[1:]})
# insert analytical solution
diff_eq = diff_eq.subs(order_params[0], expected_profile)
diff_eq.factor()
```
%% Output
$\displaystyle - \frac{3 \tau_{0 2} \left(4 \alpha^{2} {\partial {\partial (\frac{\tanh{\left(\frac{x}{2 \alpha} \right)}}{2} + \frac{1}{2}) }} - \tanh^{3}{\left(\frac{x}{2 \alpha} \right)} + \tanh{\left(\frac{x}{2 \alpha} \right)}\right)}{2 \alpha}$
⎛ 2 3⎛ x ⎞ ⎛ x ⎞⎞
-3⋅τ₀ ₂⋅⎜4⋅α ⋅D(D(tanh(x/(2*alpha))/2 + 1/2)) - tanh ⎜───⎟ + tanh⎜───⎟⎟
⎝ ⎝2⋅α⎠ ⎝2⋅α⎠⎠
────────────────────────────────────────────────────────────────────────
2⋅α
%% Cell type:markdown id: tags:
finally the differentials have to be evaluated...
%% Cell type:code id: tags:
``` python
from pystencils.fd import evaluate_diffs
diff_eq = evaluate_diffs(diff_eq, x)
diff_eq
```
%% Output
$\displaystyle \frac{3 \tau_{0 2} \left(1 - \tanh^{2}{\left(\frac{x}{2 \alpha} \right)}\right) \tanh{\left(\frac{x}{2 \alpha} \right)}}{2 \alpha} + \frac{12 \tau_{0 2} \left(\frac{\tanh{\left(\frac{x}{2 \alpha} \right)}}{2} + \frac{1}{2}\right)^{3}}{\alpha} - \frac{18 \tau_{0 2} \left(\frac{\tanh{\left(\frac{x}{2 \alpha} \right)}}{2} + \frac{1}{2}\right)^{2}}{\alpha} + \frac{6 \tau_{0 2} \left(\frac{\tanh{\left(\frac{x}{2 \alpha} \right)}}{2} + \frac{1}{2}\right)}{\alpha}$
3
⎛ ⎛ x ⎞ ⎞ ⎛ ⎛
⎜tanh⎜───⎟ ⎟ ⎜tanh⎜─
⎛ 2⎛ x ⎞⎞ ⎛ x ⎞ ⎜ ⎝2⋅α⎠ 1⎟ ⎜ ⎝2
3⋅τ₀ ₂⋅⎜1 - tanh ⎜───⎟⎟⋅tanh⎜───⎟ 12⋅τ₀ ₂⋅⎜───────── + ─⎟ 18⋅τ₀ ₂⋅⎜──────
⎝ ⎝2⋅α⎠⎠ ⎝2⋅α⎠ ⎝ 2 2⎠ ⎝ 2
───────────────────────────────── + ──────────────────────── - ───────────────
2⋅α α α
2
x ⎞ ⎞ ⎛ ⎛ x ⎞ ⎞
──⎟ ⎟ ⎜tanh⎜───⎟ ⎟
⋅α⎠ 1⎟ ⎜ ⎝2⋅α⎠ 1⎟
─── + ─⎟ 6⋅τ₀ ₂⋅⎜───────── + ─⎟
2⎠ ⎝ 2 2⎠
───────── + ──────────────────────
α
%% Cell type:markdown id: tags:
.. and the result simplified...
%% Cell type:code id: tags:
``` python
assert diff_eq.expand() == 0
diff_eq.expand()
```
%% Output
$\displaystyle 0$
0
%% Cell type:markdown id: tags:
...and indeed the expected tanh profile satisfies this differential equation.
Next lets check for the interface between phase 0 and phase 1:
%% Cell type:code id: tags:
``` python
for diff_eq in mu_diff_eq:
eq = diff_eq.subs({order_params[0]: expected_profile,
order_params[1]: 1 - expected_profile})
assert evaluate_diffs(eq, x).expand() == 0
```
%% Cell type:markdown id: tags:
### Checking the surface tensions parameters
Computing the excess free energy per unit area of an interface between two phases.
This should be exactly the surface tension parameter.
%% Cell type:code id: tags:
``` python
order_params = symbolic_order_parameters(num_symbols=num_phases-1)
F = free_energy_functional_n_phases(order_parameters=order_params)
```
%% Cell type:code id: tags:
``` python
two_phase_free_energy = F.subs({order_params[0]: expected_profile,
order_params[1]: 1 - expected_profile})
# Evaluate differentials and simplification
two_phase_free_energy = sp.simplify(evaluate_diffs(two_phase_free_energy, x))
```
%% Cell type:code id: tags:
``` python
excess_free_energy = sp.Integral(two_phase_free_energy, x)
excess_free_energy
```
%% Output
$\displaystyle \int \frac{3 \tau_{0 1}}{8 \alpha \cosh^{4}{\left(\frac{x}{2 \alpha} \right)}}\, dx$
⎮ 3⋅τ₀ ₁
⎮ ────────────── dx
⎮ 4⎛ x ⎞
⎮ 8⋅α⋅cosh ⎜───⎟
⎮ ⎝2⋅α⎠
%% Cell type:markdown id: tags:
Sympy cannot integrate this automatically - help with a manual substitution $\frac{x}{2\alpha} \rightarrow u$
%% Cell type:code id: tags:
``` python
coshTerm = list(excess_free_energy.atoms(sp.cosh))[0]
transformed_int = excess_free_energy.transform(coshTerm.args[0], sp.Symbol("u", real=True))
transformed_int
```
%% Output
$\displaystyle \int \frac{3 \tau_{0 1}}{4 \cosh^{4}{\left(u \right)}}\, du$
⎮ 3⋅τ₀ ₁
⎮ ────────── du
⎮ 4
⎮ 4⋅cosh (u)
%% Cell type:markdown id: tags:
Now the integral can be done:
%% Cell type:code id: tags:
``` python
result = sp.integrate(transformed_int.args[0], (transformed_int.args[1][0], -sp.oo, sp.oo))
assert result == symmetric_symbolic_surface_tension(0,1)
result
```
%% Output
$\displaystyle \tau_{0 1}$
τ₀ ₁
......@@ -46,7 +46,7 @@
},
{
"cell_type": "code",
"execution_count": 5,
"execution_count": 3,
"metadata": {},
"outputs": [
{
......@@ -113,7 +113,7 @@
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.6.6"
"version": "3.8.2"
}
},
"nbformat": 4,
......
......@@ -144,7 +144,7 @@
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.6.4"
"version": "3.8.2"
}
},
"nbformat": 4,
......
Source diff could not be displayed: it is too large. Options to address this: view the blob.
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.6.7"
"version": "3.11.0rc1"
}
},
"nbformat": 4,
......
%% Cell type:code id: tags:
``` python
import pytest
pytest.importorskip('cupy')
```
%% Output
<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 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);
```
%% Cell type:code id: tags:
``` python
from lbmpy.session import *
from lbmpy.phasefield.analytical import *
from lbmpy.phasefield.n_phase_boyer import *
from lbmpy.phasefield.kerneleqs import *
from pystencils.fd.spatial import discretize_spatial
from pystencils.simp import sympy_cse_on_assignment_list
from lbmpy.phasefield.cahn_hilliard_lbm import *
from pystencils.fd.derivation import *
one = sp.sympify(1)
```
%% Cell type:markdown id: tags:
# Chemical potential
Current state:
- not stable (yet)
- without LB coupling the model is stable
%% Cell type:code id: tags:
``` python
domain_size = (100, 100)
n = 4
c = sp.symbols(f"c_:{n}")
simple_potential = False
omega_h = 1.3
if simple_potential:
f = free_energy_functional_n_phases_penalty_term(c, interface_width=1, kappa=(0.05, 0.05/2, 0.05/4))
else:
ε = one * 4
mobility = one * 2 / 1000
κ = (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 )
α, _ = diffusion_coefficients(σ)
f_b, f_if, mu_b, mu_if = chemical_potential_n_phase_boyer(c, ε, σ, one,
assume_nonnegative=True, zero_threshold=1e-10)
```
%% Cell type:markdown id: tags:
# Data Setup
%% Cell type:code id: tags:
``` python
dh = create_data_handling(domain_size, periodicity=True, default_ghost_layers=2)
p_linearization = symmetric_tensor_linearization(dh.dim)
c_field = dh.add_array("c", values_per_cell=n)
rho_field = dh.add_array("rho")
mu_field = dh.add_array("mu", values_per_cell=n, latex_name="\\mu")
pressure_field = dh.add_array("P", values_per_cell=len(p_linearization))
force_field = dh.add_array("F", values_per_cell=dh.dim)
u_field = dh.add_array("u", values_per_cell=dh.dim)
# Distribution functions for each order parameter
pdf_field = []
pdf_dst_field = []
for i in range(n):
pdf_field_local = dh.add_array(f"f{i}", values_per_cell=9) # 9 for D2Q9
pdf_dst_field_local = dh.add_array(f"f{i}_dst", values_per_cell=9, latex_name="f%d_{dst}" % i)
pdf_field.append(pdf_field_local)
pdf_dst_field.append(pdf_dst_field_local)
# Distribution functions for the hydrodynamics
pdf_hydro_field = dh.add_array("fh", values_per_cell=9)
pdf_hydro_dst_field = dh.add_array("fh_dst", values_per_cell=9, latex_name="fh_{dst}")
```
%% Cell type:markdown id: tags:
### μ-Kernel
%% Cell type:code id: tags:
``` python
if simple_potential:
mu_assignments = mu_kernel(f, c, c_field, mu_field, discretization='isotropic')
else:
mu_subs = {a: b for a, b in zip(c, c_field.center_vector)}
mu_if_discrete = [discretize_spatial(e.subs(mu_subs), dx=1, stencil='isotropic') for e in mu_if]
mu_b_discrete = [e.subs(mu_subs) for e in mu_b]
mu_assignments = [Assignment(mu_field(i),
mu_if_discrete[i] + mu_b_discrete[i]) for i in range(n)]
mu_assignments = sympy_cse_on_assignment_list(mu_assignments)
μ_kernel = create_kernel(mu_assignments).compile()
```
%% Cell type:code id: tags:
``` python
second_neighbor_stencil = [(i, j)
for i in (-2, -1, 0, 1, 2)
for j in (-2, -1, 0, 1, 2)
]
x_diff = FiniteDifferenceStencilDerivation((0,), second_neighbor_stencil)
x_diff.set_weight((2, 0), sp.Rational(1, 10))
x_diff.assume_symmetric(0, anti_symmetric=True)
x_diff.assume_symmetric(1)
x_diff_stencil = x_diff.get_stencil(isotropic=True)
y_diff = FiniteDifferenceStencilDerivation((1,), second_neighbor_stencil)
y_diff.set_weight((0, 2), sp.Rational(1, 10))
y_diff.assume_symmetric(1, anti_symmetric=True)
y_diff.assume_symmetric(0)
y_diff_stencil = y_diff.get_stencil(isotropic=True)
```
%% Cell type:code id: tags:
``` python
μ = mu_field
μ_discretization = {}
for i in range(n):
μ_discretization.update({Diff(μ(i), 0): x_diff_stencil.apply(μ(i)),
Diff(μ(i), 1): y_diff_stencil.apply(μ(i))})
force_rhs = force_from_phi_and_mu(order_parameters=c_field.center_vector, dim=dh.dim, mu=mu_field.center_vector)
force_rhs = force_rhs.subs(μ_discretization)
force_assignments = [Assignment(force_field(i), force_rhs[i]) for i in range(dh.dim)]
force_kernel = create_kernel(force_assignments).compile()
```
%% Cell type:markdown id: tags:
## Lattice Boltzmann kernels
For each order parameter a Cahn-Hilliard LB method computes the time evolution.
%% Cell type:code id: tags:
``` python
if simple_potential:
mu_alpha = mu_field.center_vector
else:
mu_alpha = mobility * α * mu_field.center_vector
```
%% Cell type:code id: tags:
``` python
# Defining the Cahn-Hilliard Collision assignments
ch_kernels = []
for i in range(n):
stencil = LBStencil(Stencil.D2Q9)
ch_method = cahn_hilliard_lb_method(stencil, mu_alpha[i],
relaxation_rate=1.0, gamma=1.0)
lbm_config = LBMConfig(lb_method=ch_method, velocity_input=u_field.center_vector,
compressible=True, output={'density': c_field(i)})
lbm_opt = LBMOptimisation(symbolic_field=pdf_field[i],
symbolic_temporary_field=pdf_dst_field[i])
kernel = create_lb_function(lbm_config=lbm_config, lbm_optimisation=lbm_opt)
ch_kernels.append(kernel)
```
%% Cell type:code id: tags:
``` python
lbm_config = LBMConfig(relaxation_rate=omega_h, force=force_field, compressible=True,
output={'velocity': u_field})
lbm_opt = LBMOptimisation(symbolic_field=pdf_hydro_field,
symbolic_temporary_field=pdf_hydro_dst_field)
hydro_lbm = create_lb_function(lbm_config=lbm_config, lbm_optimisation=lbm_opt)
```
%% Cell type:markdown id: tags:
# Initialization
%% Cell type:code id: tags:
``` python
def set_c(slice_obj, values):
for block in dh.iterate(slice_obj):
arr = block[c_field.name]
arr[..., : ] = values
def init():
dh.fill(u_field.name, 0)
dh.fill(mu_field.name, 0)
dh.fill(force_field.name, 0)
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])
pdf_sync_fns = dh.synchronization_function([f.name for f in pdf_field])
hydro_sync_fn=dh.synchronization_function([pdf_hydro_field.name])
c_sync_fn = dh.synchronization_function([c_field.name])
mu_sync_fn = dh.synchronization_function([mu_field.name])
def time_loop(steps):
for i in range(steps):
c_sync_fn()
dh.run_kernel(μ_kernel)
mu_sync_fn()
dh.run_kernel(force_kernel)
hydro_sync_fn()
dh.run_kernel(hydro_lbm)
dh.swap(pdf_hydro_field.name, pdf_hydro_dst_field.name)
pdf_sync_fns()
for i in range(n):
dh.run_kernel(ch_kernels[i])
dh.swap(pdf_field[i].name,pdf_dst_field[i].name)
```
%% Cell type:code id: tags:
``` python
init()
plt.phase_plot(dh.gather_array(c_field.name))
```
%% Output
%% Cell type:code id: tags:
``` python
time_loop(10)
```
%% Cell type:code id: tags:
``` python
#plt.scalar_field(dh.cpu_arrays[mu_field.name][..., 0])
plt.scalar_field(dh.cpu_arrays[u_field.name][..., 0])
plt.colorbar();
```
%% Output
Source diff could not be displayed: it is too large. Options to address this: view the blob.