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
  • improved_comm
  • master
  • 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
42 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 4648 additions and 70 deletions
from lbmpy.creationfunctions import LBMConfig
from lbmpy.enums import Method, Stencil, CollisionSpace
from lbmpy.maxwellian_equilibrium import generate_equilibrium_by_matching_moments
from lbmpy.moments import extract_monomials
from lbmpy.stencils import LBStencil
import pytest
import sympy as sp
from pystencils.simp import AssignmentCollection
from pystencils import Assignment
from lbmpy.creationfunctions import create_lb_method
from lbmpy.methods import CollisionSpaceInfo
from lbmpy.moment_transforms import (
FastCentralMomentTransform,
PdfsToCentralMomentsByMatrix,
PdfsToCentralMomentsByShiftMatrix,
BinomialChimeraTransform)
sympy_numeric_version = [int(x, 10) for x in sp.__version__.split('.')]
if len(sympy_numeric_version) < 3:
sympy_numeric_version.append(0)
sympy_numeric_version.reverse()
sympy_version = sum(x * (100 ** i) for i, x in enumerate(sympy_numeric_version))
reference_equilibria = dict()
@pytest.mark.skipif(sympy_version < 10200,
reason="Old Sympy Versions take too long to form the inverse moment matrix")
@pytest.mark.parametrize('stencil_name', [Stencil.D2Q9, Stencil.D3Q19, Stencil.D3Q27])
@pytest.mark.parametrize('cm_transform', [PdfsToCentralMomentsByMatrix,
FastCentralMomentTransform,
PdfsToCentralMomentsByShiftMatrix,
BinomialChimeraTransform])
def test_equilibrium_pdfs(stencil_name, cm_transform):
stencil = LBStencil(stencil_name)
cspace = CollisionSpaceInfo(CollisionSpace.CUMULANTS, central_moment_transform_class=cm_transform)
lbm_config = LBMConfig(stencil=stencil, method=Method.CUMULANT, compressible=True, zero_centered=False,
collision_space_info=cspace)
c_lb_method = create_lb_method(lbm_config=lbm_config)
rho = c_lb_method.zeroth_order_equilibrium_moment_symbol
u = c_lb_method.first_order_equilibrium_moment_symbols
pdfs = c_lb_method.post_collision_pdf_symbols
cqe_assignments = [Assignment(sym, sym) for sym in u] + [Assignment(rho, rho)]
cqe = AssignmentCollection(main_assignments=cqe_assignments)
method_equilibrium_eqs = c_lb_method.get_equilibrium(cqe, False, False).main_assignments_dict
# Reference Equations
ref_equilibrium = reference_equilibria.get(stencil_name, None)
if ref_equilibrium is None:
raw_moments = list(extract_monomials(c_lb_method.cumulants, dim=stencil.D))
ref_equilibrium = generate_equilibrium_by_matching_moments(
stencil, tuple(raw_moments), rho=rho, u=u, c_s_sq=sp.Rational(1, 3), order=2*stencil.D)
reference_equilibria[stencil_name] = ref_equilibrium
for i in range(stencil.Q):
method_eq = method_equilibrium_eqs[pdfs[i]]
ref_eq = ref_equilibrium[i]
assert method_eq.expand() - ref_eq.expand() == sp.sympify(0), \
f"Equilibrium equation for pdf index {i} did not match."
import pytest
import numpy as np
from dataclasses import replace
from lbmpy.advanced_streaming import BetweenTimestepsIndexing, Timestep, get_timesteps
from lbmpy.creationfunctions import create_lb_function, create_lb_collision_rule,\
create_lb_method, LBMConfig, LBMOptimisation
from lbmpy.boundaries import LatticeBoltzmannBoundaryHandling, NoSlip, UBB
from lbmpy.enums import Method, Stencil
from lbmpy.relaxationrates import relaxation_rate_from_lattice_viscosity
from lbmpy.macroscopic_value_kernels import pdf_initialization_assignments
from lbmpy.stencils import LBStencil
from pystencils import create_kernel, create_data_handling, Assignment, Target, CreateKernelConfig
from pystencils.slicing import slice_from_direction, get_slice_before_ghost_layer
def flow_around_sphere(stencil, galilean_correction, fourth_order_correction, L_LU, total_steps):
if galilean_correction and stencil.Q != 27:
pytest.skip()
target = Target.GPU
streaming_pattern = 'aa'
timesteps = get_timesteps(streaming_pattern)
u_max = 0.05
Re = 500000
kinematic_viscosity = (L_LU * u_max) / Re
initial_velocity = (u_max,) + (0,) * (stencil.D - 1)
omega_v = relaxation_rate_from_lattice_viscosity(kinematic_viscosity)
channel_size = (10 * L_LU,) + (5 * L_LU,) * (stencil.D - 1)
sphere_position = (channel_size[0] // 3,) + (channel_size[1] // 2,) * (stencil.D - 1)
sphere_radius = L_LU // 2
lbm_config = LBMConfig(stencil=stencil, method=Method.CUMULANT, relaxation_rate=omega_v,
compressible=True,
galilean_correction=galilean_correction, fourth_order_correction=fourth_order_correction)
lbm_opt = LBMOptimisation(pre_simplification=True)
config = CreateKernelConfig(target=target)
lb_method = create_lb_method(lbm_config=lbm_config)
def get_extrapolation_kernel(timestep):
boundary_assignments = []
indexing = BetweenTimestepsIndexing(
pdf_field, stencil, streaming_pattern=streaming_pattern, prev_timestep=timestep)
f_out, _ = indexing.proxy_fields
for i, d in enumerate(stencil):
if d[0] == -1:
asm = Assignment(f_out.neighbor(0, 1)(i), f_out.center(i))
boundary_assignments.append(asm)
boundary_assignments = indexing.substitute_proxies(
boundary_assignments)
iter_slice = get_slice_before_ghost_layer((1,) + (0,) * (stencil.D - 1))
extrapolation_ast = create_kernel(
boundary_assignments,
config=CreateKernelConfig(
iteration_slice=iter_slice,
target=target,
skip_independence_check=True)
)
return extrapolation_ast.compile()
dh = create_data_handling(channel_size, periodicity=False, default_layout='fzyx', default_target=target)
u_field = dh.add_array('u', stencil.D)
rho_field = dh.add_array('rho', 1)
pdf_field = dh.add_array('pdfs', stencil.Q)
dh.fill(u_field.name, 0.0, ghost_layers=True)
dh.fill(rho_field.name, 0.0, ghost_layers=True)
dh.to_gpu(u_field.name)
dh.to_gpu(rho_field.name)
lbm_opt = replace(lbm_opt, symbolic_field=pdf_field)
bh = LatticeBoltzmannBoundaryHandling(lb_method, dh, pdf_field.name,
streaming_pattern=streaming_pattern, target=target)
wall = NoSlip()
inflow = UBB(initial_velocity)
bh.set_boundary(inflow, slice_from_direction('W', stencil.D))
directions = ('N', 'S', 'T', 'B') if stencil.D == 3 else ('N', 'S')
for direction in directions:
bh.set_boundary(wall, slice_from_direction(direction, stencil.D))
outflow_kernels = [get_extrapolation_kernel(Timestep.EVEN), get_extrapolation_kernel(Timestep.ODD)]
def sphere_boundary_callback(x, y, z=None):
x = x - sphere_position[0]
y = y - sphere_position[1]
z = z - sphere_position[2] if z is not None else 0
return np.sqrt(x ** 2 + y ** 2 + z ** 2) <= sphere_radius
bh.set_boundary(wall, mask_callback=sphere_boundary_callback)
init_eqs = pdf_initialization_assignments(lb_method, 1.0, initial_velocity, pdf_field,
streaming_pattern=streaming_pattern,
previous_timestep=timesteps[0])
init_kernel = create_kernel(init_eqs, config=config).compile()
output = {
'density': rho_field,
'velocity': u_field
}
lbm_config = replace(lbm_config, output=output)
lb_collision_rule = create_lb_collision_rule(lb_method=lb_method, lbm_config=lbm_config, lbm_optimisation=lbm_opt)
lb_kernels = []
for t in timesteps:
lbm_config = replace(lbm_config, timestep=t)
lbm_config = replace(lbm_config, streaming_pattern=streaming_pattern)
lb_kernels.append(create_lb_function(collision_rule=lb_collision_rule,
lbm_config=lbm_config, lbm_optimisation=lbm_opt, config=config))
timestep = timesteps[0]
dh.run_kernel(init_kernel)
stability_check_frequency = 1000
for i in range(total_steps):
bh(prev_timestep=timestep)
dh.run_kernel(outflow_kernels[timestep.idx])
timestep = timestep.next()
dh.run_kernel(lb_kernels[timestep.idx])
if i % stability_check_frequency == 0:
dh.to_cpu(u_field.name)
assert np.isfinite(dh.cpu_arrays[u_field.name]).all()
@pytest.mark.parametrize('stencil', [Stencil.D2Q9, Stencil.D3Q19, Stencil.D3Q27])
@pytest.mark.parametrize('galilean_correction', [False, True])
@pytest.mark.parametrize('fourth_order_correction', [False, True])
def test_flow_around_sphere_short(stencil, galilean_correction, fourth_order_correction):
pytest.importorskip('cupy')
stencil = LBStencil(stencil)
if fourth_order_correction and stencil.Q != 27:
pytest.skip("Fourth-order correction only defined for D3Q27 stencil.")
flow_around_sphere(stencil, galilean_correction, fourth_order_correction, 5, 200)
@pytest.mark.parametrize('stencil', [Stencil.D2Q9, Stencil.D3Q19, Stencil.D3Q27])
@pytest.mark.parametrize('galilean_correction', [False, True])
@pytest.mark.parametrize('fourth_order_correction', [False, 0.01])
@pytest.mark.longrun
def test_flow_around_sphere_long(stencil, galilean_correction, fourth_order_correction):
pytest.importorskip('cupy')
stencil = LBStencil(stencil)
if fourth_order_correction and stencil.Q != 27:
pytest.skip("Fourth-order correction only defined for D3Q27 stencil.")
flow_around_sphere(stencil, galilean_correction, fourth_order_correction, 20, 3000)
%% Cell type:code id: tags:
``` python
%load_ext autoreload
%autoreload 2
```
%% Cell type:code id: tags:
``` python
from dataclasses import replace
import sympy as sp
import numpy as np
from numpy.testing import assert_allclose
from pystencils.session import *
from pystencils import Target
from lbmpy.session import *
from lbmpy.macroscopic_value_kernels import macroscopic_values_getter, macroscopic_values_setter
```
%% Cell type:code id: tags:
``` python
try:
import cupy
except ImportError:
cupy = None
target = ps.Target.CPU
print('No cupy installed')
if cupy:
target = ps.Target.GPU
```
%% Cell type:markdown id: tags:
## Pipe Flow Scenario
%% Cell type:code id: tags:
``` python
class PeriodicPipeFlow:
def __init__(self, lbm_config, lbm_optimisation, config):
wall_boundary = NoSlip()
self.stencil = lbm_config.stencil
self.streaming_pattern = lbm_config.streaming_pattern
self.target = config.target
self.gpu = self.target == Target.GPU
# Stencil
self.q = stencil.Q
self.dim = stencil.D
# Streaming
self.inplace = is_inplace(self.streaming_pattern)
self.timesteps = get_timesteps(self.streaming_pattern)
self.zeroth_timestep = self.timesteps[0]
# Domain, Data Handling and PDF fields
self.pipe_length = 60
self.pipe_radius = 15
self.domain_size = (self.pipe_length, ) + (2 * self.pipe_radius,) * (self.dim - 1)
self.periodicity = (True, ) + (False, ) * (self.dim - 1)
force = (0.0001, ) + (0.0,) * (self.dim - 1)
self.force = lbm_config.force
self.dh = create_data_handling(domain_size=self.domain_size,
periodicity=self.periodicity, default_target=self.target)
self.pdfs = self.dh.add_array('pdfs', self.q)
if not self.inplace:
self.pdfs_tmp = self.dh.add_array_like('pdfs_tmp', self.pdfs.name)
lbm_optimisation = replace(lbm_optimisation, symbolic_field=self.pdfs)
if not self.inplace:
lbm_optimisation = replace(lbm_optimisation, symbolic_temporary_field=self.pdfs_tmp)
self.lb_collision = create_lb_collision_rule(lbm_config=lbm_config, lbm_optimisation=lbm_optimisation)
self.lb_method = self.lb_collision.method
self.lb_kernels = []
for t in self.timesteps:
lbm_config = replace(lbm_config, timestep=t, collision_rule=self.lb_collision)
self.lb_kernels.append(create_lb_function(lbm_config=lbm_config,
lbm_optimisation=lbm_optimisation,
config=config))
# Macroscopic Values
self.density = 1.0
self.density_field = self.dh.add_array('rho', 1)
u_x = 0.0
self.velocity = (u_x,) * self.dim
self.velocity_field = self.dh.add_array('u', self.dim)
single_gl_config = replace(config, ghost_layers=1)
setter = macroscopic_values_setter(
self.lb_method, self.density, self.velocity, self.pdfs,
streaming_pattern=self.streaming_pattern, previous_timestep=self.zeroth_timestep)
self.init_kernel = create_kernel(setter, config=single_gl_config).compile()
self.getter_kernels = []
for t in self.timesteps:
getter = macroscopic_values_getter(
self.lb_method, self.density_field, self.velocity_field, self.pdfs,
streaming_pattern=self.streaming_pattern, previous_timestep=t)
self.getter_kernels.append(create_kernel(getter, config=single_gl_config).compile())
# Periodicity
self.periodicity_handler = LBMPeriodicityHandling(
self.stencil, self.dh, self.pdfs.name, streaming_pattern=self.streaming_pattern)
# Boundary Handling
self.wall = wall_boundary
self.bh = LatticeBoltzmannBoundaryHandling(
self.lb_method, self.dh, self.pdfs.name,
streaming_pattern=self.streaming_pattern, target=self.target)
self.bh.set_boundary(boundary_obj=self.wall, mask_callback=self.mask_callback)
self.current_timestep = self.zeroth_timestep
def mask_callback(self, x, y, z=None):
y = y - self.pipe_radius
z = z - self.pipe_radius if z is not None else 0
return np.sqrt(y**2 + z**2) >= self.pipe_radius
def init(self):
self.current_timestep = self.zeroth_timestep
self.dh.run_kernel(self.init_kernel)
def step(self):
# Order matters! First communicate, then boundaries, otherwise
# periodicity handling overwrites reflected populations
# Periodicty
self.periodicity_handler(self.current_timestep)
# Boundaries
self.bh(prev_timestep=self.current_timestep)
# Here, the next time step begins
self.current_timestep = self.current_timestep.next()
# LBM Step
self.dh.run_kernel(self.lb_kernels[self.current_timestep.idx])
# Field Swaps
if not self.inplace:
self.dh.swap(self.pdfs.name, self.pdfs_tmp.name)
# Macroscopic Values
self.dh.run_kernel(self.getter_kernels[self.current_timestep.idx])
def run(self, iterations):
for _ in range(iterations):
self.step()
@property
def velocity_array(self):
if self.gpu:
self.dh.to_cpu(self.velocity_field.name)
return self.dh.gather_array(self.velocity_field.name)
def get_trimmed_velocity_array(self):
if self.gpu:
self.dh.to_cpu(self.velocity_field.name)
u = np.copy(self.dh.gather_array(self.velocity_field.name))
mask = self.bh.get_mask(None, self.wall)
for idx in np.ndindex(u.shape[:-1]):
if mask[idx] != 0:
u[idx] = np.full((self.dim, ), np.nan)
return u
```
%% Cell type:markdown id: tags:
## General Setup
%% Cell type:code id: tags:
``` python
stencil = LBStencil(Stencil.D3Q19)
dim = stencil.D
streaming_pattern = 'aa'
config = ps.CreateKernelConfig(target=target)
viscous_rr = 1.1
force = (0.0001, ) + (0.0,) * (dim - 1)
```
%% Cell type:markdown id: tags:
## 1. Reference: SRT Method
%% Cell type:code id: tags:
``` python
srt_config = LBMConfig(stencil=stencil, method=Method.SRT, relaxation_rate=viscous_rr,
force_model=ForceModel.SIMPLE, force=force, streaming_pattern=streaming_pattern)
srt_flow = PeriodicPipeFlow(srt_config, LBMOptimisation(), config)
srt_flow.init()
srt_flow.run(400)
```
%% Output
---------------------------------------------------------------------------
AttributeError Traceback (most recent call last)
Cell In[6], line 4
1 srt_config = LBMConfig(stencil=stencil, method=Method.SRT, relaxation_rate=viscous_rr,
2 force_model=ForceModel.SIMPLE, force=force, streaming_pattern=streaming_pattern)
----> 4 srt_flow = PeriodicPipeFlow(srt_config, LBMOptimisation(), config)
5 srt_flow.init()
6 srt_flow.run(400)
Cell In[4], line 47, in PeriodicPipeFlow.__init__(self, lbm_config, lbm_optimisation, config)
45 for t in self.timesteps:
46 lbm_config = replace(lbm_config, timestep=t, collision_rule=self.lb_collision)
---> 47 self.lb_kernels.append(create_lb_function(lbm_config=lbm_config,
48 lbm_optimisation=lbm_optimisation,
49 config=config))
51 # Macroscopic Values
52 self.density = 1.0
File ~/pystencils/lbmpy/lbmpy/creationfunctions.py:505, in create_lb_function(ast, lbm_config, lbm_optimisation, config, optimization, **kwargs)
502 ast = lbm_config.ast
504 if ast is None:
--> 505 ast = create_lb_ast(lbm_config.update_rule, lbm_config=lbm_config,
506 lbm_optimisation=lbm_optimisation, config=config)
508 res = ast.compile()
510 res.method = ast.method
File ~/pystencils/lbmpy/lbmpy/creationfunctions.py:530, in create_lb_ast(update_rule, lbm_config, lbm_optimisation, config, optimization, **kwargs)
525 update_rule = create_lb_update_rule(lbm_config.collision_rule, lbm_config=lbm_config,
526 lbm_optimisation=lbm_optimisation, config=config)
528 field_types = set(fa.field.dtype for fa in update_rule.defined_symbols if isinstance(fa, Field.Access))
--> 530 config = replace(config, data_type=collate_types(field_types), ghost_layers=1)
531 ast = create_kernel(update_rule, config=config)
533 ast.method = update_rule.method
File /usr/lib/python3.11/dataclasses.py:1492, in replace(obj, **changes)
1485 changes[f.name] = getattr(obj, f.name)
1487 # Create the new object, which calls __init__() and
1488 # __post_init__() (if defined), using all of the init fields we've
1489 # added and/or left in 'changes'. If there are values supplied in
1490 # changes that aren't fields, this will correctly raise a
1491 # TypeError.
-> 1492 return obj.__class__(**changes)
File <string>:24, in __init__(self, target, backend, function_name, data_type, default_number_float, default_number_int, iteration_slice, ghost_layers, cpu_openmp, cpu_vectorize_info, cpu_blocking, omp_single_loop, gpu_indexing, gpu_indexing_params, default_assignment_simplifications, cpu_prepend_optimizations, use_auto_for_assignments, index_fields, coordinate_names, allow_double_writes, skip_independence_check)
File ~/pystencils/pystencils/pystencils/config.py:177, in CreateKernelConfig.__post_init__(self)
174 self._check_type(dtype)
176 if not isinstance(self.data_type, dict):
--> 177 dt = copy(self.data_type) # The copy is necessary because BasicType has sympy shinanigans
178 self.data_type = defaultdict(self.DataTypeFactory(dt))
180 if isinstance(self.data_type, dict) and not isinstance(self.data_type, defaultdict):
File /usr/lib/python3.11/copy.py:102, in copy(x)
100 if isinstance(rv, str):
101 return x
--> 102 return _reconstruct(x, None, *rv)
File /usr/lib/python3.11/copy.py:273, in _reconstruct(x, memo, func, args, state, listiter, dictiter, deepcopy)
271 state = deepcopy(state, memo)
272 if hasattr(y, '__setstate__'):
--> 273 y.__setstate__(state)
274 else:
275 if isinstance(state, tuple) and len(state) == 2:
File ~/.local/lib/python3.11/site-packages/sympy/core/basic.py:144, in Basic.__setstate__(self, state)
143 def __setstate__(self, state):
--> 144 for name, value in state.items():
145 setattr(self, name, value)
AttributeError: 'tuple' object has no attribute 'items'
%% Cell type:code id: tags:
``` python
srt_u = srt_flow.get_trimmed_velocity_array()
ps.plot.vector_field_magnitude(srt_u[30,:,:,:])
ps.plot.colorbar()
```
%% Cell type:markdown id: tags:
## 2. Centered Cumulant Method with Central Moment-Space Forcing
%% Cell type:code id: tags:
``` python
cm_config = LBMConfig(method=Method.MONOMIAL_CUMULANT, stencil=stencil, compressible=True,
relaxation_rate=viscous_rr,
force=force, streaming_pattern=streaming_pattern)
lbm_opt = LBMOptimisation(pre_simplification=True)
cm_impl_f_flow = PeriodicPipeFlow(cm_config, lbm_opt, config)
```
%% Cell type:code id: tags:
``` python
lb_method = cm_impl_f_flow.lb_method
assert all(rr == 0 for rr in lb_method.relaxation_rates[1:4])
assert all(rr == viscous_rr for rr in lb_method.relaxation_rates[4:9])
```
%% Cell type:code id: tags:
``` python
cm_impl_f_flow.init()
cm_impl_f_flow.run(400)
```
%% Cell type:code id: tags:
``` python
cm_impl_f_u = cm_impl_f_flow.get_trimmed_velocity_array()
ps.plot.vector_field_magnitude(cm_impl_f_u[30,:,:,:])
ps.plot.colorbar()
```
%% Cell type:code id: tags:
``` python
assert_allclose(cm_impl_f_u, srt_u, rtol=1, atol=1e-4)
```
%% Cell type:markdown id: tags:
## 3. Centered Cumulant Method with Simple force model
%% Cell type:code id: tags:
``` python
from lbmpy.forcemodels import Simple
cm_expl_force_config = LBMConfig(method=Method.CUMULANT, stencil=stencil, compressible=True,
relaxation_rate=viscous_rr, force_model=Simple(force),
force=force, streaming_pattern=streaming_pattern)
lbm_opt = LBMOptimisation(pre_simplification=True)
cm_expl_f_flow = PeriodicPipeFlow(cm_expl_force_config, lbm_opt, config)
```
%% Cell type:code id: tags:
``` python
lb_method = cm_expl_f_flow.lb_method
assert all(rr == 0 for rr in lb_method.relaxation_rates[1:4])
assert all(rr == viscous_rr for rr in lb_method.relaxation_rates[4:9])
```
%% Cell type:code id: tags:
``` python
cm_expl_f_flow.init()
cm_expl_f_flow.run(400)
```
%% Cell type:code id: tags:
``` python
cm_expl_f_u = cm_expl_f_flow.get_trimmed_velocity_array()
ps.plot.vector_field_magnitude(cm_expl_f_u[30,:,:,:])
ps.plot.colorbar()
```
%% Cell type:code id: tags:
``` python
assert_allclose(cm_expl_f_u, srt_u, rtol=1, atol=1e-4)
assert_allclose(cm_expl_f_u, cm_impl_f_u, rtol=1, atol=1e-4)
```
......@@ -250,8 +250,8 @@ def run(re=6000, eval_interval=0.05, total_time=3.0, domain_size=100, u_0=0.05,
def create_full_parameter_study(gpu=False):
"""Creates a parameter study that can run the Kida vortex flow with entropic, KBC, Smagorinsky and MRT methods."""
opt_cpu = {'target': 'cpu', 'openmp': 4}
opt_gpu = {'target': 'gpu'}
opt_cpu = {'target': ps.Target.CPU, 'openmp': 4}
opt_gpu = {'target': ps.Target.GPU}
mrt_one = [{'method': 'mrt3', 'relaxation_rates': ['viscosity', 1, 1], 'stencil': stencil}
for stencil in ('D3Q19', 'D3Q27')]
......
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.
......@@ -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,11 @@ a cylinder. In Flow simulation with high-performance computers II (pp. 547-566).
import warnings
import numpy as np
import pytest
from pystencils import Target
from lbmpy._compat 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 +126,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 +150,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 +161,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.