Reply To: Method to generate effective relaxation time as output in VTK

[Yuji]

OK. I can suggest another step No.3.

Step 3-1. Implemetate new struct SmagorinskyEffectiveOmegaToWatchSGSv in src/dynamics/collisionLES.h

struct SmagorinskyEffectiveOmegaToWatchSGSv {
  using MomentaF = typename MOMENTA::template type<DESCRIPTOR>;
  using CollisionO = typename COLLISION::template type<DESCRIPTOR, MOMENTA, EQUILIBRIUM>;

  template <concepts::Cell CELL, concepts::Parameters PARAMETERS, typename V=typename CELL::value_t>
  V computeEffectiveOmega(CELL& cell, PARAMETERS& parameters) any_platform {
    V piNeqNormSqr { };
    MomentaF().computePiNeqNormSqr(cell, piNeqNormSqr);
    const V rho = MomentaF().computeRho(cell);
    const V omega = parameters.template get<descriptors::OMEGA>();
    const V smagorinsky = parameters.template get<collision::LES::SMAGORINSKY>();
    V piNeqNorm = util::sqrt(piNeqNormSqr);
    V preFactor = smagorinsky*smagorinsky
                * descriptors::invCs2<V,DESCRIPTOR>()*descriptors::invCs2<V,DESCRIPTOR>()
                * 2 * util::sqrt(2);
    /// Molecular realaxation time
    V tauMol = V{1} / omega;
    /// Turbulent realaxation time
    V tauTurb = V{0.5} * (util::sqrt(tauMol*tauMol + preFactor / rho * piNeqNorm) - tauMol);
    /// Effective realaxation time
    V tauEff = tauMol + tauTurb;
    cell.template setField<descriptors::SCALAR>(tauEff);
    return  V{1} / tauEff;
  }
  
  template <concepts::Cell CELL, concepts::Parameters PARAMETERS, typename V=typename CELL::value_t>
  CellStatistic<V> apply(CELL& cell, PARAMETERS& parameters) any_platform {
    parameters.template set<descriptors::OMEGA>(
      computeEffectiveOmega(cell, parameters));

    return CollisionO().apply(cell, parameters);
  }
};

template <typename COLLISION>
struct SmagorinskyEffectiveOmegaToWatchSGSv {
  using parameters = typename COLLISION::parameters::template include<
    descriptors::OMEGA, LES::SMAGORINSKY
  >;

  static_assert(COLLISION::parameters::template contains<descriptors::OMEGA>(),
                "COLLISION must be parametrized using relaxation frequency OMEGA");

  static std::string getName() {
    return "SmagorinskyEffectiveOmegaToWatchSGSv<" + COLLISION::getName() + ">";
  }

  template <typename DESCRIPTOR, typename MOMENTA, typename EQUILIBRIUM>
  using type = detail::SmagorinskyEffectiveOmegaToWatchSGSv<COLLISION,DESCRIPTOR,MOMENTA,EQUILIBRIUM>;
};

Step 3-2. Define BulkDynamics in your cpp, for example
using BulkDynamics = SmagorinskyEffectiveOmegaToWatchSGSv<T,DESCRIPTOR>;

If you get something, let me know.

Best regards,
Yuji