entropic_eq_srt.py 3.7 KB
Newer Older
Martin Bauer's avatar
Martin Bauer committed
1
import sympy as sp
2
from pystencils import Assignment
Martin Bauer's avatar
Martin Bauer committed
3
from lbmpy.maxwellian_equilibrium import get_weights
Martin Bauer's avatar
Martin Bauer committed
4
5
6
7
8
from lbmpy.methods.abstractlbmethod import AbstractLbMethod, LbmCollisionRule
from lbmpy.methods.conservedquantitycomputation import DensityVelocityComputation


class EntropicEquilibriumSRT(AbstractLbMethod):
Martin Bauer's avatar
Martin Bauer committed
9
    def __init__(self, stencil, relaxation_rate, force_model, conserved_quantity_calculation):
Martin Bauer's avatar
Martin Bauer committed
10
11
        super(EntropicEquilibriumSRT, self).__init__(stencil)

Martin Bauer's avatar
Martin Bauer committed
12
13
14
15
16
        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
        self.shearRelaxationRate = relaxation_rate
Martin Bauer's avatar
Martin Bauer committed
17
18

    @property
Martin Bauer's avatar
Martin Bauer committed
19
    def conserved_quantity_computation(self):
Martin Bauer's avatar
Martin Bauer committed
20
21
22
23
24
25
26
        return self._cqc

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

    @property
Martin Bauer's avatar
Martin Bauer committed
27
28
    def zeroth_order_equilibrium_moment_symbol(self, ):
        return self._cqc.zeroth_order_moment_symbol
Martin Bauer's avatar
Martin Bauer committed
29
30

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

Martin Bauer's avatar
Martin Bauer committed
34
35
36
37
    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
38

Martin Bauer's avatar
Martin Bauer committed
39
40
41
42
43
    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
44

Martin Bauer's avatar
Martin Bauer committed
45
46
        if conserved_quantity_equations is None:
            conserved_quantity_equations = self._cqc.equilibrium_input_equations_from_pdfs(f)
47
        all_subexpressions = conserved_quantity_equations.all_assignments
Martin Bauer's avatar
Martin Bauer committed
48
49
50
51
52

        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):
53
54
                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
55
56
            eq.append(f_i)

Martin Bauer's avatar
Martin Bauer committed
57
58
        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
59

Martin Bauer's avatar
Martin Bauer committed
60
61
62
63
        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, )))
            force_subexpressions = [Assignment(sym, forceModelTerm)
64
                                    for sym, forceModelTerm in zip(force_term_symbols, force_model_terms)]
Martin Bauer's avatar
Martin Bauer committed
65
66
67
            all_subexpressions += force_subexpressions
            collision_eqs = [Assignment(eq.lhs, eq.rhs + forceTermSymbol)
                             for eq, forceTermSymbol in zip(collision_eqs, force_term_symbols)]
Martin Bauer's avatar
Martin Bauer committed
68

Martin Bauer's avatar
Martin Bauer committed
69
        return LbmCollisionRule(self, collision_eqs, all_subexpressions)
Martin Bauer's avatar
Martin Bauer committed
70

Martin Bauer's avatar
Martin Bauer committed
71
72
    def get_collision_rule(self):
        return self._get_collision_rule_with_relaxation_rate(self._relaxationRate)
Martin Bauer's avatar
Martin Bauer committed
73
74


Martin Bauer's avatar
Martin Bauer committed
75
def create_srt_entropic(stencil, relaxation_rate, force_model, compressible):
Martin Bauer's avatar
Martin Bauer committed
76
77
    if not compressible:
        raise NotImplementedError("entropic-srt only implemented for compressible models")
Martin Bauer's avatar
Martin Bauer committed
78
79
    density_velocity_computation = DensityVelocityComputation(stencil, compressible, force_model)
    return EntropicEquilibriumSRT(stencil, relaxation_rate, force_model, density_velocity_computation)