Interface Walberla LB

[RudolfWeeber]

Current status:

• All LB features except pressure tensor are in
• Pressure tensor can be merged as soon as corresponding PR in PbmPy and Walberla are merged
• Electrokinetics is WIP and tracked in a separate ticket
• Documentation is not yet adapted
• GPu support and support for switching accuracy to single precision on the CPu is not yet done.

The current plan is to merge this after the release of 4.2

[codecov]

Codecov Report

Merging #2701 into python will decrease coverage by 7%.
The diff coverage is 38%.

@@           Coverage Diff            @@
##           python   #2701     +/-   ##
========================================
- Coverage      89%     81%     -8%     
========================================
  Files         557     559      +2     
  Lines       24326   25963   +1637     
========================================
- Hits        21698   21100    -598     
- Misses       2628    4863   +2235     

Impacted Files	Coverage Δ
src/core/MpiCallbacks.hpp	97% <ø> (ø)
src/core/communication.hpp	100% <ø> (ø)
src/core/electrostatics_magnetostatics/coulomb.cpp	79% <ø> (-1%)	⬇️
...d_based_algorithms/FluctuatingMRT_LatticeModel.cpp	0% <0%> (ø)
...rid_based_algorithms/FluctuatingMRT_LatticeModel.h	0% <0%> (ø)
...based_algorithms/LbWalberlaD3Q19FluctuatingMRT.hpp	0% <0%> (ø)
...ore/grid_based_algorithms/lb_particle_coupling.hpp	100% <ø> (ø)
.../grid_based_algorithms/lbboundaries/LBBoundary.hpp	73% <ø> (-27%)	⬇️
src/core/grid_based_algorithms/philox_rand.h	0% <0%> (ø)
src/core/grid_based_algorithms/lb_interface.cpp	36% <45%> (-35%)	⬇️
... and 54 more

---

Continue to review full report at Codecov.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update 95a9464...379b1d3. Read the comment docs.

[RudolfWeeber]

Implementation status of Espresso's lb interface is here:
https://github.com/espressomd/espresso/wiki/Walberla_Integartion

[RudolfWeeber]

@fweik, this could be a starting point for integrating boundary support.
Particle coupling only works on a single core for now, if skin>0

The LbWalberla class has

get_node_velocity_at_boundary(const Utils::Vector3i &node) const;
bool set_node_velocity_at_boundary(const Utils::Vector3i node,
const Utils::Vector3d v);
bool remove_node_from_boundary(const Utils::Vector3i &node);
all of which are parallel calls.
Velocities are in LB units, indices are global.

[RudolfWeeber]

@mkuron, there still seems to be an issue with the ghost communication.
Could you please take a look at the constructor of LbWalberla?
The intention was to have a ghost layer of size
int(skin/agrid +1)
and to ghost-communicate the full set of populations. Then, velocity interpolation of particles outside the box domain and on ghosts should work.
Unfortunately, that is not what happens. The lb_momentum_conservation test cannot be run on more than one node.

[mkuron]

The communication looks fine. Are you sure it's not the domain decomposition where the error is coming from? If PBCs work, then the communication is also working.

[RudolfWeeber]

@mkuron, I investigated the Nan-proliferation from Walberla's UBB() further.
It only appears, when nodes are marked as boundary at one of the domain boundaires. As far as I can tell, this case is not covered by Walberla's testing of the UBB. There, alsways SimpleUBB is used for channels.

I extracted the channel test to
src/core/unit_tests/LbWalberla.cpp (boundary_test_shear) to get rid of potential Espresso interference.
Could you please take a look at the constructor of the
LbWalberla class in src/core/grid_based_algorithms/LbWalberla.cpp and at the
LbWalberla::LB_Boundary_handling class in src/core/grid_based_algorithms/LbWalberla.hpp
Maybe a boundary related communication or the like is still missing?
The UBB uses ghost-layer fields. Do we have to add that to the communication in the LbWalberla constructor?

[RudolfWeeber]

https://github.com/RudolfWeeber/espresso/tree/walberla/src/core/grid_based_algorithms

[mkuron]

when nodes are marked as boundary at one of the domain boundaires

You need to consistently mark cells on both sides of the PBC yourself. Communication cannot take care of that for you. So when you flag a boundary at x=0, you need to do the same at x=L.

[RudolfWeeber]

On Thu, Aug 15, 2019 at 04:21:07AM -0700, Michael Kuron wrote: > when nodes are marked as boundary at one of the domain boundaires You need to consistently mark cells on both sides of the PBC yourself. Communication cannot take care of that for you. So when you flag a boundary at x=0, you need to do the same at x=L.

Where L is the grid dimension and the index runs from 0 to L-1?

[mkuron]

Exactly. So your ghost cells would be at x=-1 and X=L.

[RudolfWeeber]

OK, so it's not pretty, but works.
Since the opposing ghost cells can be on a different MPI rank, all ranks now go through the full LB grid and mark boundary cells shifted by a full lattice size in all combinations of coordiantes as well now.

In the steady state case, the shear profiles now match (testsuite/python/lb_shear.py with 500 integration steps).
The time evolution is not yet correct.
Maybe, the viscosity used in LbWalberla::LbWalberla has to be converted to lattice units.

[RudolfWeeber]

With the viscosity converted to lattice units, the time-dependent shear profile is now reproduced

[RudolfWeeber]

So, the boundary force from Walberla for the couette flow along the shear direction matches the expected result. (lb_shear.py)

@mkuron, the UBB also reports a force perpendicular to the wall. The Walberla test for BoundaryForcecouette does not check this.
What is the expected value? It is of similar magnitude to the equilibrium pressure times surface area, but off by a factor of not quite 2-3 depending on parameters.

Equilibrium pressure is:
p_eq = DENS * AGRID2 / TIME_STEP2 / 3

agrid, density, and viscosity are all !=1, but setting them to 1 doesn't help, so not unit conversions, I think.

[KaiSzuttor]

@RudolfWeeber you seem to have forgotten to commit grid_based_algorithms/lbboundaries/LBBoundary.cpp

[RudolfWeeber]

Added it.

[RudolfWeeber]

Current state (to my understanding)

Working:

• Basic LB setup and most property getters/setters
• parallelization up to 2 cores (Walberla enumerates mpi nodes differently)
• Velocity boundary conditions validated via onset of Couette flow
• partial support for forces on boundaries (tangential force on boundaires for Couette flow)
• standard particle coupling (there are slight differences between what Walberla is doing and what Espresso is doing with regards to forces applied to the LB)
• inertialess tracer coupling
• lb electrodydrodynamics
  Not working:
• ENGINE (probably affected by the same bug as ENGINE for Espresso's cpu lb)
• setting bulk viscosity independently
• non-equal time steps for MD and Walberla LB
• Thermalization

[KaiSzuttor]

why is this the expected momentum?

[RudolfWeeber]

This is the current state as per the test suite.
Thermalization is waiting for the type erasure work. That will come tomorrow, if all goes well.

Walberla is currently limited to 2 mpi ranks.

Working:

• lb_buoyancy_force.py Passes on 2 cores
• lb_poiseuille.py passes on 2 cores
• lb_thermo_virtual.py Working on 2 cores
• lb_walberla.py passes on 2 cores, but probably obsolete by now
• linear_momentum_lb.py Passing on 2 cores
• lb_electrohydrodynamics.py

Woprking on 1 core:

• lb_boundary.py On 2 nodes, the head node doesn't see correct boundary flags on boundaries stored on the 2nd node. Probably Mpi ranks should only answer, if the site is in their local domain
• lb_boundary_volume_force.py passes on 1 core. Wrong values on 2 core. Issue with boundary force reduction?

Partial:

• lb.py Passes on 2 cores with thermalization tests disabled
• lb_momentum_conservation.py Worse results on 2 nodes than on 1. Liekly bug. Generally unclear, what tolerance is acceptable. Unit-test for viscous coupling needed
• lb_shear.py the velocity profile is correct on 2 cores, but the stress tensor is wrong on 1 and 2 cores (Walberla: Figure out unit issues with stress tensor #3464)

Broken:

• engine_lb.py broken. Slight mismatch. This should be looked at after the unit-testing for particle coupling is done. Note that the total momentum used in the test now contains dt * F/2 for all point forces applied to lb in the previous time step. This is a behaviour change.
• lb_boundary_velocity.py Broken. Not entirely clear, why it eer worked, though
• lb_interpolation.py Broken (not examined)
• lb_poiseuille_cylinder.py Broken, but values appear to be proportional
• lb_stokes_sphere.py Results off by some 10 percent. Either issue with different boundary handling or with different LB model
• virtual_sites_tracers_walberla.py Broken. Liekly, the case of intertialess tracers with no lb active is not handled correctly

Waiting for thermalziation:

• lb_thermostat.py Waiting for thermalization
• lb_density.py Needs thermalization. Waiting for type reasure work
• lb.py thermalization part

Other:

• lb_streaming.py Not applicable until the same LB model as in Espresso is used.