Reply To: other non-Newtonian constitutive models

[Adrian]

OpenLB offers quite a lot of flexibility for extending and recombining already implemented models. From the software perspective it is straight forward to modify the local relaxation time. For illustration we can take a look at the declaration of the PowerLawBGKdynamics:

template <typename T, typename DESCRIPTOR, typename MOMENTA=momenta::BulkTuple>
using PowerLawBGKdynamics = dynamics::Tuple<
  T, DESCRIPTOR,
  MOMENTA,
  equilibria::SecondOrder,
  powerlaw::OmegaFromCell<collision::BGK>
>;

These dynamics consist of given momenta (bulk by default but changeable via template parameter), second order equilibrium and a BGK collision locally modified to compute the relaxation frequency omega from the cell via powerlaw::OmegaFromCell.

So in order to implement other non-Newtonian consitutive models you only need to implement the powerlaw::OmegaFromCell part and drop it into the declaration of your new dynamics. I would suggest to look at the code in src/dynamics/powerLawBGKdynamics.h for reference (also the various Smagorinsky-type LES dynamics where the same facilities are used).

For a very reduced overview:

template <typename COLLISION>
struct CustomNonNewtonianOmegaFromCell {
  // [...]
  template <typename DESCRIPTOR, typename MOMENTA, typename EQUILIBRIUM>
  struct type {
    using MomentaF = typename MOMENTA::template type<DESCRIPTOR>;
    using CollisionO = typename COLLISION::template type<DESCRIPTOR, MOMENTA, EQUILIBRIUM>;

    template <typename CELL, typename PARAMETERS, typename V=typename CELL::value_t>
    CellStatistic<V> apply(CELL& cell, PARAMETERS& parameters) any_platform
    {
      V rho, u[DESCRIPTOR::d], pi[util::TensorVal<DESCRIPTOR>::n] { };
      MomentaF().computeAllMomenta(cell, rho, u, pi); // momenta are computed here for usage in omega computation
      const V omega = // compute new omega here
      parameters.template set<descriptors::OMEGA>(omega); // default omega value is changed here
      return CollisionO().apply(cell, parameters); // arbitrary omega-accepting collision operator is called here
    }
  };
};

Once this works our automatic code generator may be used to transparently transform this into performance optimized code.

I hope that this helps, feel free to ask further questions.