entropic_eq_srt.py 4.24 KB
Newer Older
Martin Bauer's avatar
Martin Bauer committed
1
import sympy as sp
Martin Bauer's avatar
Martin Bauer committed
2

Martin Bauer's avatar
Martin Bauer committed
3
from lbmpy.maxwellian_equilibrium import get_weights
Martin Bauer's avatar
Martin Bauer committed
4
5
from lbmpy.methods.abstractlbmethod import AbstractLbMethod, LbmCollisionRule
from lbmpy.methods.conservedquantitycomputation import DensityVelocityComputation
Martin Bauer's avatar
Martin Bauer committed
6
from pystencils import Assignment
Martin Bauer's avatar
Martin Bauer committed
7
8
9


class EntropicEquilibriumSRT(AbstractLbMethod):
Martin Bauer's avatar
Martin Bauer committed
10
11
12
    """Equilibrium from 'Minimal entropic kinetic models for hydrodynamics'
    Ansumali, S. ; Karlin, I. V;  Öttinger, H. C, (2003)
    """
Martin Bauer's avatar
Martin Bauer committed
13
    def __init__(self, stencil, relaxation_rate, force_model, conserved_quantity_calculation):
Martin Bauer's avatar
Martin Bauer committed
14
15
        super(EntropicEquilibriumSRT, self).__init__(stencil)

Martin Bauer's avatar
Martin Bauer committed
16
17
18
19
        self._cqc = conserved_quantity_calculation
        self._weights = get_weights(stencil, c_s_sq=sp.Rational(1, 3))
        self._relaxationRate = relaxation_rate
        self._forceModel = force_model
Martin Bauer's avatar
Martin Bauer committed
20
        self.shear_relaxation_rate = relaxation_rate
Martin Bauer's avatar
Martin Bauer committed
21
22

    @property
Martin Bauer's avatar
Martin Bauer committed
23
    def conserved_quantity_computation(self):
Martin Bauer's avatar
Martin Bauer committed
24
25
26
27
28
29
30
        return self._cqc

    @property
    def weights(self):
        return self._weights

    @property
Martin Bauer's avatar
Martin Bauer committed
31
32
    def zeroth_order_equilibrium_moment_symbol(self, ):
        return self._cqc.zeroth_order_moment_symbol
Martin Bauer's avatar
Martin Bauer committed
33
34

    @property
Martin Bauer's avatar
Martin Bauer committed
35
36
    def first_order_equilibrium_moment_symbols(self, ):
        return self._cqc.first_order_moment_symbols
Martin Bauer's avatar
Martin Bauer committed
37

Martin Bauer's avatar
Martin Bauer committed
38
39
40
41
    def get_equilibrium(self, conserved_quantity_equations=None, include_force_terms=False):
        return self._get_collision_rule_with_relaxation_rate(1,
                                                             conserved_quantity_equations=conserved_quantity_equations,
                                                             include_force_terms=include_force_terms)
Martin Bauer's avatar
Martin Bauer committed
42

43
44
45
46
    def get_equilibrium_terms(self):
        equilibrium = self.get_equilibrium()
        return sp.Matrix([eq.rhs for eq in equilibrium.main_assignments])

Martin Bauer's avatar
Martin Bauer committed
47
48
49
50
51
    def _get_collision_rule_with_relaxation_rate(self, relaxation_rate, include_force_terms=True,
                                                 conserved_quantity_equations=None):
        f = sp.Matrix(self.pre_collision_pdf_symbols)
        rho = self._cqc.zeroth_order_moment_symbol
        u = self._cqc.first_order_moment_symbols
Martin Bauer's avatar
Martin Bauer committed
52

Martin Bauer's avatar
Martin Bauer committed
53
54
        if conserved_quantity_equations is None:
            conserved_quantity_equations = self._cqc.equilibrium_input_equations_from_pdfs(f)
55
        all_subexpressions = conserved_quantity_equations.all_assignments
Martin Bauer's avatar
Martin Bauer committed
56
57
58
59
60

        eq = []
        for w_i, direction in zip(self.weights, self.stencil):
            f_i = rho * w_i
            for u_a, e_ia in zip(u, direction):
61
62
                b = sp.sqrt(1 + 3 * u_a ** 2)
                f_i *= (2 - b) * ((2 * u_a + b) / (1 - u_a)) ** e_ia
Martin Bauer's avatar
Martin Bauer committed
63
64
            eq.append(f_i)

Martin Bauer's avatar
Martin Bauer committed
65
66
        collision_eqs = [Assignment(lhs, (1 - relaxation_rate) * f_i + relaxation_rate * eq_i)
                         for lhs, f_i, eq_i in zip(self.post_collision_pdf_symbols, self.pre_collision_pdf_symbols, eq)]
Martin Bauer's avatar
Martin Bauer committed
67

Martin Bauer's avatar
Martin Bauer committed
68
69
70
        if (self._forceModel is not None) and include_force_terms:
            force_model_terms = self._forceModel(self)
            force_term_symbols = sp.symbols("forceTerm_:%d" % (len(force_model_terms, )))
Martin Bauer's avatar
Martin Bauer committed
71
72
            force_subexpressions = [Assignment(sym, force_model_term)
                                    for sym, force_model_term in zip(force_term_symbols, force_model_terms)]
Martin Bauer's avatar
Martin Bauer committed
73
            all_subexpressions += force_subexpressions
Martin Bauer's avatar
Martin Bauer committed
74
75
            collision_eqs = [Assignment(eq.lhs, eq.rhs + force_term_symbol)
                             for eq, force_term_symbol in zip(collision_eqs, force_term_symbols)]
76
77
78
        cr = LbmCollisionRule(self, collision_eqs, all_subexpressions)
        cr.simplification_hints['relaxation_rates'] = []
        return cr
Martin Bauer's avatar
Martin Bauer committed
79

Martin Bauer's avatar
Martin Bauer committed
80
    def get_collision_rule(self, conserved_quantity_equations=None, keep_rrs_symbolic=None):
81
82
        return self._get_collision_rule_with_relaxation_rate(self._relaxationRate,
                                                             conserved_quantity_equations=conserved_quantity_equations)
Martin Bauer's avatar
Martin Bauer committed
83
84


Martin Bauer's avatar
Martin Bauer committed
85
def create_srt_entropic(stencil, relaxation_rate, force_model, compressible):
Martin Bauer's avatar
Martin Bauer committed
86
87
    if not compressible:
        raise NotImplementedError("entropic-srt only implemented for compressible models")
Martin Bauer's avatar
Martin Bauer committed
88
89
    density_velocity_computation = DensityVelocityComputation(stencil, compressible, force_model)
    return EntropicEquilibriumSRT(stencil, relaxation_rate, force_model, density_velocity_computation)