Possible error in MovingWall example

Palabos Bug Reports
#1

Hi everyone,

I have a question about the movingWall show case, but first let me add a quick disclaimer: I just started using LBM and hence may be a bit off about the following.

In the show case [2] the InamuroIteration DataProcessor is called with the tau parameter equal to the relaxation time of the fluid. However following Inamuro [1] this is only true for a single iteration. In the paper the corrected Eulerian velocity u_l(x) at iteration l is calculated as u_0(x) + \Delta t g_l(x). This corresponds to setting the tau parameter of the InamuroIteration DataProcessor to one.

However to incorporate the body force into the LBM method the fluid velocity used to calculate the equilibrium should be given by u_0 + \tau g_l(x) before the collision. Below I added a DataProcessor which modifies the Eulerian velocity given by the InamuroIteration to include the additional relaxation time.

This is different from using the relaxation time in the InamuroIteration as for example the vertex forces will differ starting from the second iteration.

So to summarize: using the DataProcessor below I think the fluid-object coupling in the show case [2] should read

[code=“cpp”]
instantiateImmersedWallData(vertices, areas, container);
for (int i = 0; i < param.ibIter; i++) {
inamuroIteration(SurfaceVelocity(timeLB),
*rhoBar, *j, container, 1.0);
}

    inamuroFixRelaxationVelocity(*lattice, *j, 1.0 / params.omega);




[1]: Inamuro, "Lattice Boltzmann methods for moving boundary flows", Fluid Dyn. Res 44, 024001 (2012).
[2]: examples/showCases/movingWall/movingWall.cpp

[b]Edit:[/b] fix getTypeOfModification error

[code="cpp"]
// header
template<typename T, template<typename U> class Descriptor>
class InamuroFixRelaxationVelocity : public plb::BoxProcessingFunctional3D
{
public:
    InamuroFixRelaxationVelocity(T relaxationTime);

    virtual void processGenericBlocks(plb::Box3D domain, std::vector<plb::AtomicBlock3D*> fields);
    virtual InamuroFixRelaxationVelocity<T,Descriptor>* clone() const;
    virtual void getTypeOfModification(std::vector<plb::modif::ModifT>& modified) const;
    virtual plb::BlockDomain::DomainT appliesTo() const;

private:
    T relaxationTime;
};

template<typename T, template<typename U> class Descriptor>
void inamuroFixRelaxationVelocity(plb::MultiBlockLattice3D<T, Descriptor>& lattice,
		plb::MultiTensorField3D<T, 3>& current, T relaxationTime)
{
	std::vector<plb::MultiBlock3D*> args;
	args.push_back(&lattice);
	args.push_back(&current);

    plb::applyProcessingFunctional (new InamuroFixRelaxationVelocity<T, Descriptor>(relaxationTime),
    		lattice.getBoundingBox(), args);
}

[code=“cpp”]
// implementation
template<typename T, template class Descriptor>
InamuroFixRelaxationVelocity<T,Descriptor>::InamuroFixRelaxationVelocity(T relaxationTime)
: relaxationTime(relaxationTime)
{

}

template<typename T, template class Descriptor>
void InamuroFixRelaxationVelocity<T,Descriptor>::processGenericBlocks (
plb::Box3D domain, std::vectorplb::AtomicBlock3D* fields )
{
plb::BlockLattice3D<T,Descriptor>& lattice = dynamic_cast<plb::BlockLattice3D<T,Descriptor>>(fields[0]);
plb::TensorField3D<T,3>& jField = dynamic_cast<plb::TensorField3D<T,3>>(fields[1]);

plb::Dot3D off = plb::computeRelativeDisplacement(lattice, jField);

for (plb::plint iX=domain.x0; iX<=domain.x1; ++iX) {
   for (plb::plint iY=domain.y0; iY<=domain.y1; ++iY) {
        for (plb::plint iZ=domain.z0; iZ<=domain.z1; ++iZ) {
            plb::Cell<T,Descriptor> const& cell = lattice.get(iX,iY,iZ);

            plb::Array<T, 3> current;
            plb::momentTemplates<T,Descriptor>::get_j(cell, current);

            plb::Array<T, 3>& relaxationCurrent = jField.get(iX + off.x,iY + off.y, iZ + off.z);
            relaxationCurrent = (1 - relaxationTime) * current + relaxationTime * relaxationCurrent;
        }
    }
}

}

template<typename T, template class Descriptor>
InamuroFixRelaxationVelocity<T,Descriptor>* InamuroFixRelaxationVelocity<T,Descriptor>::clone() const
{
return new InamuroFixRelaxationVelocity<T,Descriptor>(*this);
}

template<typename T, template class Descriptor>
void InamuroFixRelaxationVelocity<T,Descriptor>::getTypeOfModification(std::vectorplb::modif::ModifT& modified) const {
modified[0] = plb::modif::nothing; // lattice
modified[1] = plb::modif::staticVariables; // j
}

template<typename T, template class Descriptor>
plb::BlockDomain::DomainT InamuroFixRelaxationVelocity<T,Descriptor>::appliesTo() const {
return plb::BlockDomain::bulk;
}