#include <shanChenForcedPostProcessor.h>

Collaboration diagram for olb::ShanChenForcedPostProcessor< POTENTIAL >:

Classes
struct	G

struct	OMEGA_A

struct	OMEGA_B

struct	RHO0

Public Types
using	parameters

Public Member Functions
template<typename CELLS , typename PARAMETERS >
void	apply (CELLS &cells, PARAMETERS &parameters) any_platform

Static Public Attributes
static constexpr OperatorScope	scope = OperatorScope::PerCellWithParameters

Detailed Description

template<typename POTENTIAL>
struct olb::ShanChenForcedPostProcessor< POTENTIAL >

Definition at line 61 of file shanChenForcedPostProcessor.h.

Member Typedef Documentation

◆ parameters

template<typename POTENTIAL >

using olb::ShanChenForcedPostProcessor< POTENTIAL >::parameters

Initial value:

 typename POTENTIAL::parameters::template decompose_into<
    meta::list<RHO0,G,OMEGA_A,OMEGA_B>::template include
  >

Definition at line 69 of file shanChenForcedPostProcessor.h.

Member Function Documentation

◆ apply()

template<typename POTENTIAL >

template<typename CELLS , typename PARAMETERS >

void olb::ShanChenForcedPostProcessor< POTENTIAL >::apply	(	CELLS &	cells,
		PARAMETERS &	parameters )

inline

Definition at line 74 of file shanChenForcedPostProcessor.h.

  {
    using V = typename CELLS::template value_t<names::A>::value_t;
    using DESCRIPTOR = typename CELLS::template value_t<names::A>::descriptor_t;
 
    auto& cellA = cells.template get<names::A>();
    auto& cellB = cells.template get<names::B>();
 
    auto rho0 = parameters.template get<RHO0>();
 
    Vector<V,2> rhoField{
      cellA.template getFieldComponent<descriptors::STATISTIC>(0) * rho0[0],
      cellB.template getFieldComponent<descriptors::STATISTIC>(0) * rho0[1]
    };
    Vector<V,DESCRIPTOR::d> uA{};
    lbm<DESCRIPTOR>::computeJ(cellA, uA);
    Vector<V,DESCRIPTOR::d> uB{};
    lbm<DESCRIPTOR>::computeJ(cellB, uB);
 
    V omegaA = parameters.template get<OMEGA_A>();
    V omegaB = parameters.template get<OMEGA_B>();
    // Computation of the common velocity, shared among the two populations
    V rhoTot = rhoField[0]*omegaA
             + rhoField[1]*omegaB;
 
    Vector<V,DESCRIPTOR::d> uTot = (uA*rho0[0]*omegaA + uB*rho0[1]*omegaB) / rhoTot;
 
    // Computation of the interaction potential
    Vector<V,DESCRIPTOR::d> rhoBlockContribution{};
    Vector<V,DESCRIPTOR::d> rhoPartnerContribution{};
    V psi1 = POTENTIAL().compute(rhoField[0], parameters);
    V psi2 = POTENTIAL().compute(rhoField[1], parameters);
 
    for (int iPop = 0; iPop < DESCRIPTOR::q; ++iPop) {
      auto nextCellA = cellA.neighbor(descriptors::c<DESCRIPTOR>(iPop));
      auto nextCellB = cellB.neighbor(descriptors::c<DESCRIPTOR>(iPop));
      V rhoA = POTENTIAL().compute(nextCellA.template getFieldComponent<descriptors::STATISTIC>(0)*rho0[0], parameters);
      V rhoB = POTENTIAL().compute(nextCellB.template getFieldComponent<descriptors::STATISTIC>(0)*rho0[1], parameters);
      rhoBlockContribution   += psi2 * rhoA * descriptors::c<DESCRIPTOR>(iPop) * descriptors::t<V,DESCRIPTOR>(iPop);
      rhoPartnerContribution += psi1 * rhoB * descriptors::c<DESCRIPTOR>(iPop) * descriptors::t<V,DESCRIPTOR>(iPop);
    }
 
    // Computation and storage of the final velocity, consisting
    //   of u and the momentum difference due to interaction
    //   potential plus external force
    auto externalBlockForce   = cellA.template getField<descriptors::EXTERNAL_FORCE>();
    auto externalPartnerForce = cellB.template getField<descriptors::EXTERNAL_FORCE>();
 
    auto g = parameters.template get<G>();
 
    cellA.template setField<descriptors::VELOCITY>(uTot);
    cellB.template setField<descriptors::VELOCITY>(uTot);
    cellA.template setField<descriptors::FORCE>(externalBlockForce
        - g*rhoPartnerContribution/rhoField[0]);
    cellB.template setField<descriptors::FORCE>(externalPartnerForce
        - g*rhoBlockContribution/rhoField[1]);
  }