From 36029a9ac9531fb1e3a5155ff9ede760f9a24284 Mon Sep 17 00:00:00 2001
From: Felix Winterhalter <felix.winterhalter@fau.de>
Date: Fri, 2 Aug 2019 16:04:36 +0200
Subject: [PATCH] Variance equations from temperature
---
lbmpy/creationfunctions.py | 21 ++++++++++++++++-----
lbmpy/fluctuatinglb.py | 14 ++++++++++++++
lbmpy_tests/test_fluctuation.py | 3 ++-
3 files changed, 32 insertions(+), 6 deletions(-)
diff --git a/lbmpy/creationfunctions.py b/lbmpy/creationfunctions.py
index 6c18b64..936ec0d 100644
--- a/lbmpy/creationfunctions.py
+++ b/lbmpy/creationfunctions.py
@@ -178,7 +178,7 @@ from lbmpy.fieldaccess import (
AAEvenTimeStepAccessor, AAOddTimeStepAccessor, CollideOnlyInplaceAccessor, EsoTwistEvenTimeStepAccessor,
EsoTwistOddTimeStepAccessor, PdfFieldAccessor, PeriodicTwoFieldsAccessor, StreamPullTwoFieldsAccessor,
StreamPushTwoFieldsAccessor)
-from lbmpy.fluctuatinglb import fluctuation_correction, method_with_rescaled_equilibrium_values
+from lbmpy.fluctuatinglb import fluctuation_correction, fluctuating_variance_equations
from lbmpy.methods import create_mrt3, create_mrt_orthogonal, create_mrt_raw, create_srt, create_trt, create_trt_kbc
from lbmpy.methods.creationfunctions import create_generic_mrt
from lbmpy.methods.cumulantbased import CumulantBasedLbMethod
@@ -309,9 +309,6 @@ def create_lb_collision_rule(lb_method=None, optimization={}, **kwargs):
if rho_in is not None and isinstance(rho_in, Field):
rho_in = rho_in.center
- if params['fluctuating']:
- lb_method = method_with_rescaled_equilibrium_values(lb_method)
-
if u_in is not None:
density_rhs = sum(lb_method.pre_collision_pdf_symbols) if rho_in is None else rho_in
eqs = [Assignment(cqc.zeroth_order_moment_symbol, density_rhs)]
@@ -326,13 +323,26 @@ def create_lb_collision_rule(lb_method=None, optimization={}, **kwargs):
collision_rule = simplification(collision_rule)
if params['fluctuating']:
- variances = SymbolGen("variance") if params['fluctuating'] is True else params['fluctuating']
+ variances = params["fluctuating"]
+ if params["fluctuating"] is True:
+ variances = [v for v, _ in zip(iter(SymbolGen("variance")), lb_method.moments)]
+
correction = fluctuation_correction(lb_method, random_symbol(collision_rule.subexpressions, dim=lb_method.dim),
variances)
+
for i, corr in enumerate(correction):
collision_rule.main_assignments[i] = Assignment(collision_rule.main_assignments[i].lhs,
collision_rule.main_assignments[i].rhs + corr)
+ if params["temperature"] is not None:
+ temperature = sp.Symbol("temperature") if params["temperature"] is True else params["temperature"]
+ variance_equations = fluctuating_variance_equations(method=lb_method,
+ temperature=temperature,
+ c_s_sq=params["c_s_sq"],
+ variances=variances)
+ collision_rule.subexpressions += variance_equations
+ collision_rule.topological_sort(sort_subexpressions=True, sort_main_assignments=False)
+
if params['entropic']:
if params['smagorinsky']:
raise ValueError("Choose either entropic or smagorinsky")
@@ -528,6 +538,7 @@ def update_with_default_parameters(params, opt_params=None, fail_on_unknown_para
'omega_output_field': None,
'smagorinsky': False,
'fluctuating': False,
+ 'temperature': None,
'output': {},
'velocity_input': None,
diff --git a/lbmpy/fluctuatinglb.py b/lbmpy/fluctuatinglb.py
index c0d4ab3..076e063 100644
--- a/lbmpy/fluctuatinglb.py
+++ b/lbmpy/fluctuatinglb.py
@@ -23,6 +23,20 @@ import sympy as sp
from lbmpy.moments import MOMENT_SYMBOLS
from pystencils.simp.assignment_collection import SymbolGen
+from pystencils import Assignment
+
+
+def fluctuating_variance_equations(method, temperature=sp.Symbol("fluctuating_temperature"),
+ c_s_sq=sp.Symbol("c_s") ** 2,
+ variances=SymbolGen("variance")):
+ """Produces variance equations according to (3.54) in Schiller08"""
+ normalization_factors = abs(method.moment_matrix) * sp.Matrix(method.weights)
+ relaxation_rates_sqr = [r ** 2 for r in method.relaxation_rates]
+ density = method.zeroth_order_equilibrium_moment_symbol
+ mu = temperature * density / c_s_sq
+
+ return [Assignment(v, sp.sqrt(mu * norm * (1 - rr)))
+ for v, norm, rr in zip(iter(variances), normalization_factors, relaxation_rates_sqr)]
def method_with_rescaled_equilibrium_values(base_method):
diff --git a/lbmpy_tests/test_fluctuation.py b/lbmpy_tests/test_fluctuation.py
index f3e45f5..88b01ca 100644
--- a/lbmpy_tests/test_fluctuation.py
+++ b/lbmpy_tests/test_fluctuation.py
@@ -2,7 +2,8 @@ from lbmpy.scenarios import create_channel
def test_fluctuating_generation_pipeline():
ch = create_channel((40, 10), method='mrt3', relaxation_rates=[1.5, 1, 1], force=1e-5,
- fluctuating=[1e-3] * 9,
+ fluctuating=True,
+ temperature=1e-9,
kernel_params={'time_step': 1})
ch.run(10)
--
GitLab