Skip to content
Snippets Groups Projects
Commit e7823f08 authored by Helen Schottenhamml's avatar Helen Schottenhamml
Browse files

Merge branch 'FixRR' into 'master'

Relaxation rates should be floats

See merge request !105
parents c49363ff ff68674d
1 merge request!105Relaxation rates should be floats
Pipeline #35385 passed with stages
in 53 minutes and 33 seconds
......@@ -2,6 +2,7 @@ import abc
from collections import namedtuple
import sympy as sp
from sympy.core.numbers import Zero
from pystencils import Assignment, AssignmentCollection
......@@ -52,6 +53,7 @@ class AbstractLbMethod(abc.ABC):
"""Returns a qxq diagonal matrix which contains the relaxation rate for each moment on the diagonal"""
d = sp.zeros(len(self.relaxation_rates))
for i in range(0, len(self.relaxation_rates)):
# note that 0.0 is converted to sp.Zero here. It is not possible to prevent this.
d[i, i] = self.relaxation_rates[i]
return d
......@@ -104,6 +106,9 @@ class AbstractLbMethod(abc.ABC):
for relaxation_rate in rr:
if relaxation_rate not in unique_relaxation_rates:
relaxation_rate = sp.sympify(relaxation_rate)
# special treatment for zero, sp.Zero would be an integer ..
if isinstance(relaxation_rate, Zero):
relaxation_rate = 0.0
if not isinstance(relaxation_rate, sp.Symbol):
rt_symbol = sp.Symbol(f"rr_{len(subexpressions)}")
subexpressions[relaxation_rate] = rt_symbol
......
......@@ -605,11 +605,11 @@ def _get_relaxation_info_dict(relaxation_rates, nested_moments, dim):
for group in nested_moments:
for moment in group:
if get_order(moment) <= 1:
result[moment] = 0
result[moment] = 0.0
elif is_shear_moment(moment, dim):
result[moment] = relaxation_rates[0]
else:
result[moment] = 1
result[moment] = 1.0
# if relaxation rate for each moment is specified they are all used
if len(relaxation_rates) == number_of_moments:
......@@ -634,7 +634,7 @@ def _get_relaxation_info_dict(relaxation_rates, nested_moments, dim):
next_rr = False
for moment in group:
if get_order(moment) <= 1:
result[moment] = 0
result[moment] = 0.0
elif is_shear_moment(moment, dim):
result[moment] = shear_rr
elif is_bulk_moment(moment, dim):
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment