Reply To: Overlap for the FD BC and the simultaneous use of Bouzidi and a Thermal lattice

[vaillant]

For the second bug, let’s take one boundary cell. For initializing Bouzidi, setBouzidiVelocityBoundary is called, which add a VelocityBouzidiLinearPostProcessorGenerator3D to the neighbour fluid cell. This Post Processor is doing the Bouzidi algorithm directly on the fluid node so I understand NoDynamics is correct for the boundary node since it just store the population.

This works well with one lattice. But if you add another lattice to resolve the Advection-diffusion Equation, there will be a problem.

If we still take the same boundary cell, we initialize it with setAdvectionDiffusionTemperatureBoundary, and if we assume for example that the boundary cell is an edge, so discreteNormal[0] == 3, the dynamic is then AdvectionDiffusionEdgesDynamics and the momenta is EquilibriumBM.
When arrive the moment to collide this boundary cell, collide of AdvectionDiffusionEdgesDynamics need the velocity of the fluid to compute the equilibrium, so it calls :
auto u = cell.template getField<descriptors::VELOCITY>();

We go look at what have been stored in this field. It’s done in NavierStokesAdvectionDiffusionCouplingPostProcessor3D in processSubDomain with :

auto u = tPartner->get(iX,iY,iZ).template getField<descriptors::VELOCITY>();
blockLattice.get(iX,iY,iZ).computeU(u.data());
tPartner->get(iX,iY,iZ).template setField<descriptors::VELOCITY>(u);

The blockLattice here is one of the first lattice. So on the boundary cell, the dynamic is NoDynamics and so computeU(u.data()) will do :

  for (int iD=0; iD<DESCRIPTOR::d; ++iD) {
    u[iD] = T();
  }

.
So even if there is a velocity applied on the boundary cell, the value store in the field of the second lattice will be 0. So in AdvectionDiffusionEdgesDynamics, when you compute the equilibrium with lbH::equilibriumFirstOrder(unknownIndexes[iPop], temperature, u.data()), u comes from the field and so it is 0, and you will have a value of the equilibrium accordingly. In the meantime, the neighbour fluid cell retrieves the correct velocity from his field since computeU(u.data()) of for example BGKdynamics would have give it the right value.
So to resume, for the thermal lattice, the boundary cell ‘see’ no velocity while its fluid neighbour ‘see’ a velocity. This is why, I think, it crash when I try it. And with correcting the value of the velocity in the field, it’s working fine.