Skip to content

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

OpenLB – Open Source Lattice Boltzmann Code Forums on OpenLB General Topics Method to generate effective relaxation time as output in VTK Reply To: Method to generate effective relaxation time as output in VTK

#10060
mwkim
Participant

SORRY FOR LATE UPDATE.
I resolved the issue.

Basically, Yuji’s suggestion is on the right way.
However, We should also define variables in CellStatistic.

I attached example codes.


template <typename COLLISION, typename DESCRIPTOR, typename MOMENTA, typename EQUILIBRIUM>
struct SmagorinskyEffectiveOmegaTauEffOut {
  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>
  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;

    V tauEff_vtk = cell.template getField<descriptors::TAU_EFF>();
    tauEff_vtk = tauMol + tauTurb;
    cell.template setField<descriptors::TAU_EFF>(tauEff_vtk);

    return  V{1} / tauEff;
  }

  template <typename CELL, typename PARAMETERS, typename V=typename CELL::value_t>
  CellStatistic<V> apply(CELL& cell, PARAMETERS& parameters) any_platform {
    parameters.template set<descriptors::OMEGA>(
      computeEffectiveOmega(cell, parameters));
    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 tauEff_vtk = cell.template getField<descriptors::TAU_EFF>();
    V piNeqNorm = util::sqrt(piNeqNormSqr);
    V preFactor = smagorinsky*smagorinsky
                * descriptors::invCs2<V,DESCRIPTOR>()*descriptors::invCs2<V,DESCRIPTOR>()
                * 2 * util::sqrt(2);
    V tauMol = V{1} / omega;
    V tauTurb = V{0.5} * (util::sqrt(tauMol*tauMol + preFactor / rho * piNeqNorm) - tauMol);
    tauEff_vtk = tauMol + tauTurb;
    cell.template setField<descriptors::TAU_EFF>(tauEff_vtk);
    return CollisionO().apply(cell, parameters);
  }
};

I hope this would be helpful for some others.

Best,
Minwoo Kim