lbmpy issues https://i10git.cs.fau.de/pycodegen/lbmpy/-/issues 2021-04-30T13:22:18+02:00 https://i10git.cs.fau.de/pycodegen/lbmpy/-/issues/14 Wrong momentum density calculation for compressible methods with forces. 2021-04-30T13:22:18+02:00 Helen Schottenhamml Wrong momentum density calculation for compressible methods with forces. In _conservedquantitycomputation.py_ (line 210f) the calculation of the momentum density uses the _macroscopic velocity shift_ python mom_density_eq_coll = apply_force_model_shift(&#39;macroscopic_velocity_shift&#39;, dim, mom_density_eq_coll... In _conservedquantitycomputation.py_ (line 210f) the calculation of the momentum density uses the _macroscopic velocity shift_ python mom_density_eq_coll = apply_force_model_shift('macroscopic_velocity_shift', dim, mom_density_eq_coll, self._forceModel, self._compressible)  which defaults to <img src="https://latex.codecogs.com/svg.latex? &space;\frac{\mathbf{F}}{2\rho}" />. When using this momentum density to calculcate the velocity (as it is done in all generated lattice models in waLBerla), we divide by the density again, resulting in the overall velocity equation <img src="https://latex.codecogs.com/svg.latex? &space;\mathbf{u} = \frac{1}{\rho} \sum_i f_i \mathbf{c}_i + \frac{\mathbf{F}}{2\rho^2}" />, which is obviously wrong. Affected are all _lbmpy_ simulations that work directly on the momentum denisty, and all _waLBerla_ simulations that use generated lattice models and calculate the velocity at some point (in particular, the VTK output will be slightly off). Helen Schottenhamml Helen Schottenhamml https://i10git.cs.fau.de/pycodegen/lbmpy/-/issues/12 Apply force in moment space 2020-11-17T19:31:43+01:00 Michael Kuron mkuron@icp.uni-stuttgart.de Apply force in moment space As discussed with @bauer and @holzer, we should try applying the force in moment space. This should be equivalent to the Guo force model, but save significant FLOPs. To support this effort, the automatic Chapman-Enskog analysis should b... As discussed with @bauer and @holzer, we should try applying the force in moment space. This should be equivalent to the Guo force model, but save significant FLOPs. To support this effort, the automatic Chapman-Enskog analysis should be extended to support forces. https://i10git.cs.fau.de/pycodegen/lbmpy/-/issues/11 Guo and Buick force is wrong when applied to non-SRT LBs 2020-09-16T22:34:25+02:00 Michael Kuron mkuron@icp.uni-stuttgart.de Guo and Buick force is wrong when applied to non-SRT LBs As a simple test, we apply a constant force to a fluid and measure the resulting total momentum. It should be force * number of cells * number of time steps. We do this for different relaxation and force models. python from pystencil... As a simple test, we apply a constant force to a fluid and measure the resulting total momentum. It should be force * number of cells * number of time steps. We do this for different relaxation and force models. python from pystencils.session import * from lbmpy.session import * import lbmpy from lbmpy.macroscopic_value_kernels import macroscopic_values_setter L = (40, 40) stencil = get_stencil("D2Q9") eta = 0.15 omega = lbmpy.relaxationrates.relaxation_rate_from_lattice_viscosity(eta) F = [2e-4, -3e-4] dh = ps.create_data_handling(L, periodicity=True, default_target='cpu') src = dh.add_array('src', values_per_cell=len(stencil)) dst = dh.add_array_like('dst', 'src') ρ = dh.add_array('rho') u = dh.add_array('u', values_per_cell=dh.dim) collision = create_lb_update_rule(method="trt", stencil=stencil, relaxation_rate=omega, compressible=True, force_model='guo', force=F, kernel_type='collide_only', optimization={'symbolic_field': src}) stream = create_stream_pull_with_output_kernel(collision.method, src, dst, {'density': ρ, 'velocity': u}) opts = {'cpu_openmp': True, 'cpu_vectorize_info': None, 'target': dh.default_target} stream_kernel = ps.create_kernel(stream, **opts).compile() collision_kernel = ps.create_kernel(collision, **opts).compile() def init(): dh.fill(ρ.name, 1) dh.fill(u.name, 0) setter = macroscopic_values_setter(collision.method, velocity=(0,)*dh.dim, pdfs=src.center_vector, density=ρ.center) kernel = ps.create_kernel(setter, ghost_layers=0).compile() dh.run_kernel(kernel) sync_pdfs = dh.synchronization_function([src.name]) # needed before stream, but after collision def time_loop(steps): dh.all_to_gpu() for i in range(steps): dh.run_kernel(collision_kernel) sync_pdfs() dh.run_kernel(stream_kernel) dh.swap(src.name, dst.name) dh.all_to_cpu() t = 100 init() time_loop(t) total = np.sum(dh.gather_array(u.name), axis=(0,1)) print(total/np.prod(L)/F/t)  We see that for SRT or Luo, the script always prints 1.0, so all force applied has been converted into momentum. For TRT with Guo, it produces an omega-dependent number, which means that force has not been applied correctly. | Method | Momentum per Force | | ------ | ------------------ | | SRT Guo | 1.0 | | SRT Luo | 1.0 | | SRT Buick | 1.0 | | TRT Guo | * | | TRT Luo | 1.0 | | TRT Buick | * | | omega | * | | ----- | ----- | | 2 | 0 | | 1.8 | 0.22903226 | | 1.5 | 0.55769231 | | 1.2 | 0.87058824 | | 1 | 1.07142857 | | 0.8 | 1.26666667 | | 0.5 | 1.55 | | 0.25 | 1.77822581 | | 0.1 | 1.92 | | 0 | 2 | Michael Kuron mkuron@icp.uni-stuttgart.de Michael Kuron mkuron@icp.uni-stuttgart.de https://i10git.cs.fau.de/pycodegen/lbmpy/-/issues/10 KBC LBMs do not work with the current MRT description 2020-06-15T20:59:48+02:00 Markus Holzer KBC LBMs do not work with the current MRT description Conserved quantities are relaxed with zero. If a KBC is created it looks at the number of relaxation rates which has to be 2. Since the conserved quantities are already fixed with 0 it causes problems. Possible fix would be to state th... Conserved quantities are relaxed with zero. If a KBC is created it looks at the number of relaxation rates which has to be 2. Since the conserved quantities are already fixed with 0 it causes problems. Possible fix would be to state the relaxation time of the conserved quantities as a parameter for the create functions. Markus Holzer Markus Holzer https://i10git.cs.fau.de/pycodegen/lbmpy/-/issues/9 Boundaries don&#39;t work with OpenCL 2020-03-12T20:21:59+01:00 Michael Kuron mkuron@icp.uni-stuttgart.de Boundaries don't work with OpenCL Boundary conditions work fine on CPU and CUDA, but they don&#39;t have any effect on OpenCL. Here is a Poiseuille flow test case that shows that, just change the target in the third block. [gpu_boundary.ipynb](/uploads/88e7dc4e1c9201013e831c... Boundary conditions work fine on CPU and CUDA, but they don't have any effect on OpenCL. Here is a Poiseuille flow test case that shows that, just change the target in the third block. [gpu_boundary.ipynb](/uploads/88e7dc4e1c9201013e831c6f15de5832/gpu_boundary.ipynb) https://i10git.cs.fau.de/pycodegen/lbmpy/-/issues/8 test_phasefieldstep_direct fails after upgrade of scikit-image from 0.15.0 to... 2019-12-04T11:15:01+01:00 Michael Kuron mkuron@icp.uni-stuttgart.de test_phasefieldstep_direct fails after upgrade of scikit-image from 0.15.0 to 0.16.1 The *pycodegen-integration* job in pystencils now reports a failing test_phasefieldstep_direct (https://i10git.cs.fau.de/pycodegen/pystencils/-/jobs/327659). After a local test, I can confirm that this happens with scikit-image 0.16.1 an... The *pycodegen-integration* job in pystencils now reports a failing test_phasefieldstep_direct (https://i10git.cs.fau.de/pycodegen/pystencils/-/jobs/327659). After a local test, I can confirm that this happens with scikit-image 0.16.1 and did not happen with scikit-image 0.15.0. Martin Bauer Martin Bauer https://i10git.cs.fau.de/pycodegen/lbmpy/-/issues/7 Lees-Edwards LB 2021-01-22T16:47:22+01:00 Michael Kuron mkuron@icp.uni-stuttgart.de Lees-Edwards LB Will be implemented by @Bindgen according to https://doi.org/10.1023/A:1014595628808. It requires a modified equilibrium distribution and a communication scheme that interpolates populations across the periodic boundary. Will be implemented by @Bindgen according to https://doi.org/10.1023/A:1014595628808. It requires a modified equilibrium distribution and a communication scheme that interpolates populations across the periodic boundary. Sebastian Bindgen Sebastian Bindgen https://i10git.cs.fau.de/pycodegen/lbmpy/-/issues/6 The fluctuating feature should throw if the basis vectors of the underlying L... 2019-11-27T16:47:43+01:00 RudolfWeeber The fluctuating feature should throw if the basis vectors of the underlying LB aren't orthogonal Michael Kuron mkuron@icp.uni-stuttgart.de Michael Kuron mkuron@icp.uni-stuttgart.de https://i10git.cs.fau.de/pycodegen/lbmpy/-/issues/5 create_mrt_orthogonal&#39;s definition of orthogonality is inconsistent 2019-11-28T20:34:47+01:00 Michael Kuron mkuron@icp.uni-stuttgart.de create_mrt_orthogonal's definition of orthogonality is inconsistent Orthogonality can either be unweighted (as in a scalar product) or weighted (as in a weighted scalar product). lbmpy.methods.creationfunctions.create_mrt_orthogonal mixes the two. The D2Q9 is unweighted-orthogonal, the D3Q19 is weighte... Orthogonality can either be unweighted (as in a scalar product) or weighted (as in a weighted scalar product). lbmpy.methods.creationfunctions.create_mrt_orthogonal mixes the two. The D2Q9 is unweighted-orthogonal, the D3Q19 is weighted-orthogonal, and the D3Q27 is unweighted-orthogonal. ~~The D3Q15 is almost weighted-orthogonal, but contains a minor error somewhere that leads to two off-diagonal elements when checking for weighted-orthogonality~~(solved by !15). The Krüger book (and the original d'Humieres et al. paper) give the moments for an unweighted-orthogonal D3Q15 and D3Q19 and I don't know where the weighted D3Q15 version in the code came from. The weighted D3Q19 clearly comes from Schiller/Dünweg/Ladd. The D2Q9 can be made weighted-orthogonal by this simple change: diff --- a/lbmpy/methods/creationfunctions.py +++ b/lbmpy/methods/creationfunctions.py @@ -10,7 +10,7 @@ from lbmpy.maxwellian_equilibrium import ( compressible_to_incompressible_moment_value, get_cumulants_of_continuous_maxwellian_equilibrium, get_cumulants_of_discrete_maxwellian_equilibrium, get_moments_of_continuous_maxwellian_equilibrium, - get_moments_of_discrete_maxwellian_equilibrium) + get_moments_of_discrete_maxwellian_equilibrium, get_weights) from lbmpy.methods.abstractlbmethod import RelaxationInfo from lbmpy.methods.conservedquantitycomputation import DensityVelocityComputation from lbmpy.methods.cumulantbased import CumulantBasedLbMethod @@ -388,9 +388,10 @@ def create_mrt_orthogonal(stencil, relaxation_rate_getter=None, maxwellian_momen moment_to_relaxation_rate_dict = OrderedDict() if have_same_entries(stencil, get_stencil("D2Q9")): moments = get_default_moment_set_for_stencil(stencil) - orthogonal_moments = gram_schmidt(moments, stencil) + weights = get_weights(stencil, sp.Rational(1,3)) + orthogonal_moments = gram_schmidt(moments, stencil, weights) orthogonal_moments_scaled = [e * common_denominator(e) for e in orthogonal_moments] nested_moments = list(sort_moments_into_groups_of_same_order(orthogonal_moments_scaled).values()) elif have_same_entries(stencil, get_stencil("D3Q15")): sq = x ** 2 + y ** 2 + z ** 2 nested_moments = [  It can also be made to match the version from the Schiller/Dünweg/Ladd paper, which chooses slightly different moments, by this change: diff --- a/lbmpy/methods/creationfunctions.py +++ b/lbmpy/methods/creationfunctions.py @@ -10,7 +10,7 @@ from lbmpy.maxwellian_equilibrium import ( compressible_to_incompressible_moment_value, get_cumulants_of_continuous_maxwellian_equilibrium, get_cumulants_of_discrete_maxwellian_equilibrium, get_moments_of_continuous_maxwellian_equilibrium, - get_moments_of_discrete_maxwellian_equilibrium) + get_moments_of_discrete_maxwellian_equilibrium, get_weights) from lbmpy.methods.abstractlbmethod import RelaxationInfo from lbmpy.methods.conservedquantitycomputation import DensityVelocityComputation from lbmpy.methods.cumulantbased import CumulantBasedLbMethod @@ -388,9 +388,16 @@ def create_mrt_orthogonal(stencil, relaxation_rate_getter=None, maxwellian_momen moment_to_relaxation_rate_dict = OrderedDict() if have_same_entries(stencil, get_stencil("D2Q9")): moments = get_default_moment_set_for_stencil(stencil) - orthogonal_moments = gram_schmidt(moments, stencil) + weights = get_weights(stencil, sp.Rational(1,3)) + orthogonal_moments = gram_schmidt(moments, stencil, weights) orthogonal_moments_scaled = [e * common_denominator(e) for e in orthogonal_moments] nested_moments = list(sort_moments_into_groups_of_same_order(orthogonal_moments_scaled).values()) + sq = x ** 2 + y ** 2 + nested_moments = 3 * sq - 2 + nested_moments = 2 * x ** 2 - sq + nested_moments = (3 * sq - 4) * x + nested_moments = (3 * sq - 4) * y + nested_moments = 9 * sq ** 2 - 15 * sq + 2 elif have_same_entries(stencil, get_stencil("D3Q15")): sq = x ** 2 + y ** 2 + z ** 2 nested_moments = [  ---------------- How to check for unweighted orthogonality: python m = create_mrt_orthogonal(get_stencil("D3Q27"), maxwellian_moments=True) assert (m.moment_matrix * m.moment_matrix.T).is_diagonal()  How to check for weighted orthogonality: python m = create_mrt_orthogonal(get_stencil("D3Q19"), maxwellian_moments=True) w = get_weights(m.stencil, sp.Rational(1,3)) assert (sp.matrix_multiply_elementwise(m.moment_matrix, sp.Matrix([w]*19)) * m.moment_matrix.T).is_diagonal()  ----------------- The test case test_mrt_orthogonal in test_momentbased_methods_equilibrium.py currently only checks the two unweighted-orthogonal ones. ------------------ To resolve the inconsistency, all variants should be automatically constructed via Gram-Schmidt. An optional argument should be provided that allows the choice of weighted vs. unweighted. Also, some tricks may be necessary during orthogonalization to get those specific moments that certain people in literature use; these should also be selectable via an optional argument. Michael Kuron mkuron@icp.uni-stuttgart.de Michael Kuron mkuron@icp.uni-stuttgart.de https://i10git.cs.fau.de/pycodegen/lbmpy/-/issues/4 Tests for fluctuating MRT LB 2019-11-21T18:02:15+01:00 Michael Kuron mkuron@icp.uni-stuttgart.de Tests for fluctuating MRT LB RudolfWeeber RudolfWeeber https://i10git.cs.fau.de/pycodegen/lbmpy/-/issues/3 Implementation of FreeSlip Boundaries 2021-07-09T10:14:41+02:00 Markus Holzer Implementation of FreeSlip Boundaries FreeSlip Boundaries as in waLBerla should be added Link to the waLBerla implementation: [waLBerla FreeSlip](https://i10git.cs.fau.de/walberla/walberla/blob/master/src/lbm/boundary/FreeSlip.h) FreeSlip Boundaries as in waLBerla should be added Link to the waLBerla implementation: [waLBerla FreeSlip](https://i10git.cs.fau.de/walberla/walberla/blob/master/src/lbm/boundary/FreeSlip.h) https://i10git.cs.fau.de/pycodegen/lbmpy/-/issues/1 Relace update_with_default_parameters 2021-09-21T12:35:16+02:00 Stephan Seitz Relace update_with_default_parameters LBM methods and kernels can be constructed in a fabricator pattern with a lot of magic parameters bundled as dict. Users need to read the documentation or the source code to know about all the option. Some lazy users are to lazy to read.... LBM methods and kernels can be constructed in a fabricator pattern with a lot of magic parameters bundled as dict. Users need to read the documentation or the source code to know about all the option. Some lazy users are to lazy to read. Wouldn't it be better to provide a dict-like type for that with a __init__ that replaces update_with_default_parameters to provide the user with autocomplete and documentation (member-wise) about which fixed members this dict has? The type could still duck type or subclass a dict and would be explicitly convertible via dict(vars(type)) https://i10git.cs.fau.de/pycodegen/lbmpy/-/issues/2 Fluctuating (Randomized, Thermalized) LBM 2019-08-12T16:53:58+02:00 Martin Bauer Fluctuating (Randomized, Thermalized) LBM - [ ] get random number from external field, add kernel to test suite @winterhalter - [ ] physical test case to check random number correlation (auto correlation function?) @kuron - [x] final version should compute random numbers in th... - [ ] get random number from external field, add kernel to test suite @winterhalter - [ ] physical test case to check random number correlation (auto correlation function?) @kuron - [x] final version should compute random numbers in the kernel itself. 'mulhi' / 'mullo' functions are required - philox RNG - write a pystencils function similar to bitoperations for these operations - on x86 there are intrinsics for mulhi/mullo, alternatively one could also use AES instead of philox since there is also a intrinsic for that - on GPU there should also be an intrinsic - C fallback as one-liner - two versions required for drawing float & double random numbers - fast versions produce 2 or more RN per step, how to handle this? Fixes https://i10git.cs.fau.de/walberla/walberla/issues/80 Martin Bauer Martin Bauer