Skip to content
Snippets Groups Projects
Commit 63d44cc9 authored by Christoph Rettinger's avatar Christoph Rettinger
Browse files

Merge branch 'fluctuating' into 'master'

Fluctuating MRT: introduce separate relaxation rates

See merge request !275
parents 31a01773 b686edda
No related merge requests found
...@@ -70,7 +70,7 @@ int main( int argc, char ** argv ) ...@@ -70,7 +70,7 @@ int main( int argc, char ** argv )
// create fields // create fields
BlockDataID forceFieldId = field::addToStorage<VectorField_T>( blocks, "Force", real_t( 0.0 )); BlockDataID forceFieldId = field::addToStorage<VectorField_T>( blocks, "Force", real_t( 0.0 ));
LatticeModel_T latticeModel = LatticeModel_T( forceFieldId, omega, seed, temperature, uint_t(0) ); LatticeModel_T latticeModel = LatticeModel_T( forceFieldId, omega, 0, omega, omega, seed, temperature, uint_t(0) );
BlockDataID pdfFieldId = lbm::addPdfFieldToStorage( blocks, "pdf field", latticeModel, initialVelocity, real_t(1) ); BlockDataID pdfFieldId = lbm::addPdfFieldToStorage( blocks, "pdf field", latticeModel, initialVelocity, real_t(1) );
BlockDataID flagFieldId = field::addFlagFieldToStorage< FlagField_T >( blocks, "flag field" ); BlockDataID flagFieldId = field::addFlagFieldToStorage< FlagField_T >( blocks, "flag field" );
......
import sympy as sp import sympy as sp
import pystencils as ps import pystencils as ps
from lbmpy.creationfunctions import create_lb_collision_rule from lbmpy.creationfunctions import create_lb_collision_rule, create_mrt_orthogonal, force_model_from_string
from lbmpy.moments import is_bulk_moment, is_shear_moment, get_order
from lbmpy.stencils import get_stencil
from pystencils_walberla import CodeGeneration from pystencils_walberla import CodeGeneration
from lbmpy_walberla import generate_lattice_model from lbmpy_walberla import generate_lattice_model
...@@ -9,13 +11,40 @@ with CodeGeneration() as ctx: ...@@ -9,13 +11,40 @@ with CodeGeneration() as ctx:
temperature = sp.symbols("temperature") temperature = sp.symbols("temperature")
force_field, vel_field = ps.fields("force(3), velocity(3): [3D]", layout='fzyx') force_field, vel_field = ps.fields("force(3), velocity(3): [3D]", layout='fzyx')
def rr_getter(moment_group):
is_shear = [is_shear_moment(m, 3) for m in moment_group]
is_bulk = [is_bulk_moment(m, 3) for m in moment_group]
order = [get_order(m) for m in moment_group]
assert min(order) == max(order)
order = order[0]
if order < 2:
return 0
elif any(is_bulk):
assert all(is_bulk)
return sp.Symbol("omega_bulk")
elif any(is_shear):
assert all(is_shear)
return sp.Symbol("omega_shear")
elif order % 2 == 0:
assert order > 2
return sp.Symbol("omega_even")
else:
return sp.Symbol("omega_odd")
method = create_mrt_orthogonal(
stencil=get_stencil('D3Q19'),
compressible=True,
weighted=True,
relaxation_rate_getter=rr_getter,
force_model=force_model_from_string('guo', force_field.center_vector)
)
collision_rule = create_lb_collision_rule( collision_rule = create_lb_collision_rule(
stencil='D3Q19', compressible=True, fluctuating={ method,
fluctuating={
'temperature' : temperature, 'temperature' : temperature,
'block_offsets' : 'walberla', 'block_offsets' : 'walberla',
}, },
method='mrt', relaxation_rates=[omega_shear]*19,
force_model='guo', force=force_field.center_vector,
optimization={'cse_global': True} optimization={'cse_global': True}
) )
......
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