olb::MixedScaleBoussinesqCouplingPostProcessor2D< T, DESCRIPTOR > Class Template Reference

#include <navierStokesAdvectionDiffusionCouplingPostProcessor2D.h>

Inheritance diagram for olb::MixedScaleBoussinesqCouplingPostProcessor2D< T, DESCRIPTOR >:

Collaboration diagram for olb::MixedScaleBoussinesqCouplingPostProcessor2D< T, DESCRIPTOR >:

Public Member Functions
	MixedScaleBoussinesqCouplingPostProcessor2D (int x0_, int x1_, int y0_, int y1_, T gravity_, T T0_, T deltaTemp_, std::vector< T > dir_, T PrTurb_, std::vector< BlockStructureD< 2 > * > partners_)
int	extent () const override
	Extent of application area (0 for purely local operations)
int	extent (int whichDirection) const override
	Extent of application area along a direction (0 or 1)
void	process (BlockLattice< T, DESCRIPTOR > &blockLattice) override
	Execute post-processing step.
void	processSubDomain (BlockLattice< T, DESCRIPTOR > &blockLattice, int x0_, int x1_, int y0_, int y1_) override
	Execute post-processing step on a sublattice.
Public Member Functions inherited from olb::PostProcessor2D< T, DESCRIPTOR >
	PostProcessor2D ()
virtual	~PostProcessor2D ()
std::string &	getName ()
	read and write access to name
std::string const &	getName () const
	read only access to name
int	getPriority () const
	read only access to priority

Additional Inherited Members
Protected Attributes inherited from olb::PostProcessor2D< T, DESCRIPTOR >
int	_priority

Detailed Description

template<typename T, typename DESCRIPTOR>
class olb::MixedScaleBoussinesqCouplingPostProcessor2D< T, DESCRIPTOR >

Definition at line 182 of file navierStokesAdvectionDiffusionCouplingPostProcessor2D.h.

Constructor & Destructor Documentation

MixedScaleBoussinesqCouplingPostProcessor2D()

template<typename T , typename DESCRIPTOR >

olb::MixedScaleBoussinesqCouplingPostProcessor2D< T, DESCRIPTOR >::MixedScaleBoussinesqCouplingPostProcessor2D	(	int	x0_,
		int	x1_,
		int	y0_,
		int	y1_,
		T	gravity_,
		T	T0_,
		T	deltaTemp_,
		std::vector< T >	dir_,
		T	PrTurb_,
		std::vector< BlockStructureD< 2 > * >	partners_ )

Definition at line 421 of file navierStokesAdvectionDiffusionCouplingPostProcessor2D.hh.

  :  x0(x0_), x1(x1_), y0(y0_), y1(y1_),
     gravity(gravity_), T0(T0_), deltaTemp(deltaTemp_),
     dir(dir_), PrTurb(PrTurb_), partners(partners_)
{
  // we normalize the direction of force vector
  T normDir = T();
  for (unsigned iD = 0; iD < dir.size(); ++iD) {
    normDir += dir[iD]*dir[iD];
  }
  normDir = util::sqrt(normDir);
  for (unsigned iD = 0; iD < dir.size(); ++iD) {
    dir[iD] /= normDir;
  }
 
  for (unsigned iD = 0; iD < dir.size(); ++iD) {
    forcePrefactor[iD] = gravity * dir[iD];
  }
 
  tauTurbADPrefactor = descriptors::invCs2<T,descriptors::D2Q5<descriptors::VELOCITY,descriptors::TAU_EFF,descriptors::CUTOFF_HEAT_FLUX>>() / descriptors::invCs2<T,DESCRIPTOR>() / PrTurb;
  tPartner = static_cast<BlockLattice<T,descriptors::D2Q5<descriptors::VELOCITY,descriptors::TAU_EFF,descriptors::CUTOFF_HEAT_FLUX>> *>(partners[0]);
}

References olb::util::sqrt().

Here is the call graph for this function:

Member Function Documentation

extent() [1/2]

template<typename T , typename DESCRIPTOR >

int olb::MixedScaleBoussinesqCouplingPostProcessor2D< T, DESCRIPTOR >::extent ( ) const

inlineoverridevirtual

Extent of application area (0 for purely local operations)

Implements olb::PostProcessor2D< T, DESCRIPTOR >.

Definition at line 187 of file navierStokesAdvectionDiffusionCouplingPostProcessor2D.h.

  {
    return 0;
  }

extent() [2/2]

template<typename T , typename DESCRIPTOR >

int olb::MixedScaleBoussinesqCouplingPostProcessor2D< T, DESCRIPTOR >::extent ( int direction ) const

inlineoverridevirtual

Extent of application area along a direction (0 or 1)

Implements olb::PostProcessor2D< T, DESCRIPTOR >.

Definition at line 191 of file navierStokesAdvectionDiffusionCouplingPostProcessor2D.h.

  {
    return 0;
  }

process()

template<typename T , typename DESCRIPTOR >

void olb::MixedScaleBoussinesqCouplingPostProcessor2D< T, DESCRIPTOR >::process ( BlockLattice< T, DESCRIPTOR > & blockLattice )

overridevirtual

Execute post-processing step.

Implements olb::PostProcessor2D< T, DESCRIPTOR >.

Definition at line 618 of file navierStokesAdvectionDiffusionCouplingPostProcessor2D.hh.

{
  processSubDomain(blockLattice, x0, x1, y0, y1);
}

processSubDomain()

template<typename T , typename DESCRIPTOR >

void olb::MixedScaleBoussinesqCouplingPostProcessor2D< T, DESCRIPTOR >::processSubDomain ( BlockLattice< T, DESCRIPTOR > & blockLattice, int x0_, int x1_, int y0_, int y1_ )

overridevirtual

Execute post-processing step on a sublattice.

Molecular realaxation time

Turbulent realaxation time

Effective realaxation time

Implements olb::PostProcessor2D< T, DESCRIPTOR >.

Definition at line 448 of file navierStokesAdvectionDiffusionCouplingPostProcessor2D.hh.

{
 
  const T C_nu = 0.04;
  const T C_alpha = 0.5;
  const T deltaT = 1.0;
 
  const T invCs2_g = descriptors::invCs2<T,descriptors::D2Q5<descriptors::VELOCITY,descriptors::TAU_EFF,descriptors::CUTOFF_HEAT_FLUX>>();
 
  int newX0, newX1, newY0, newY1;
  if ( util::intersect (
         x0, x1, y0, y1,
         x0_, x1_, y0_, y1_,
         newX0, newX1, newY0, newY1 ) ) {
 
    auto& heatFluxCache = blockLattice.template getField<HEAT_FLUX_CACHE>()[0];
    for (int iX=newX0-1; iX<=newX1+1; ++iX) {
      for (int iY=newY0-1; iY<=newY1+1; ++iY) {
        const T temperature = tPartner->get(iX,iY).computeRho();
        heatFluxCache[blockLattice.getCellId(iX, iY)] = temperature;
 
        // computation of the bousinessq force
        T temperatureDifference = temperature - T0;
        blockLattice.get(iX,iY).template setField<descriptors::FORCE>(temperatureDifference*forcePrefactor);
 
        FieldD<T,DESCRIPTOR,descriptors::VELOCITY> u;
        blockLattice.get(iX,iY).computeU(u.data());
        tPartner->get(iX,iY).template setField<descriptors::VELOCITY>(u);
      }
    }
 
    for (int iX=newX0; iX<=newX1; ++iX) {
      for (int iY=newY0; iY<=newY1; ++iY) {
 
        auto u_pp = tPartner->get(iX+1, iY+1).template getField<descriptors::VELOCITY>();
        auto u_0p = tPartner->get(iX, iY+1).template getField<descriptors::VELOCITY>();
        auto u_np = tPartner->get(iX-1, iY+1).template getField<descriptors::VELOCITY>();
        auto u_p0 = tPartner->get(iX+1, iY  ).template getField<descriptors::VELOCITY>();
        auto u_pn = tPartner->get(iX+1, iY-1).template getField<descriptors::VELOCITY>();
        auto u_00 = tPartner->get(iX, iY  ).template getField<descriptors::VELOCITY>();
        auto u_0n = tPartner->get(iX, iY-1).template getField<descriptors::VELOCITY>();
        auto u_n0 = tPartner->get(iX-1, iY  ).template getField<descriptors::VELOCITY>();
        auto u_nn = tPartner->get(iX-1, iY-1).template getField<descriptors::VELOCITY>();
 
        const T *h_pp = & heatFluxCache[blockLattice.getCellId(iX+1, iY+1)];
        const T *h_0p = & heatFluxCache[blockLattice.getCellId(iX, iY+1)];
        const T *h_np = & heatFluxCache[blockLattice.getCellId(iX-1, iY+1)];
        const T *h_p0 = & heatFluxCache[blockLattice.getCellId(iX+1, iY  )];
        const T *h_pn = & heatFluxCache[blockLattice.getCellId(iX+1, iY-1)];
        const T *h_00 = & heatFluxCache[blockLattice.getCellId(iX, iY  )];
        const T *h_0n = & heatFluxCache[blockLattice.getCellId(iX, iY-1)];
        const T *h_n0 = & heatFluxCache[blockLattice.getCellId(iX-1, iY  )];
        const T *h_nn = & heatFluxCache[blockLattice.getCellId(iX-1, iY-1)];
 
        Vector<T, 2> filtered_u;
        T filtered_h;
        filtered_h =((h_pp[0] + 2.*h_0p[0] + h_np[0])
                     + 2.*(h_p0[0] + 2.*h_00[0] + h_n0[0])
                     + (h_pn[0] + 2.*h_0n[0] + h_nn[0]))*0.25*0.25;
 
        filtered_u =((u_pp + 2.*u_0p + u_np)
                     + 2.*(u_p0 + 2.*u_00 + u_n0)
                     + (u_pn + 2.*u_0n + u_nn))*0.25*0.25;
 
        Vector<T,2> filtered_u_reduced = u_00 - filtered_u;
        Vector<T,2> filtered_heatFlux_reduced = h_00[0]*u_00 - filtered_h*filtered_u;
 
        T cutoffKinEnergy = filtered_u_reduced[0]*filtered_u_reduced[0]
                            + filtered_u_reduced[1]*filtered_u_reduced[1];
        T cutoffHeatFlux = filtered_heatFlux_reduced[0]*filtered_heatFlux_reduced[0]
                           + filtered_heatFlux_reduced[1]*filtered_heatFlux_reduced[1];
 
        blockLattice.get(iX,iY).template setField<descriptors::CUTOFF_KIN_ENERGY>(util::pow(0.5*cutoffKinEnergy, 0.25));
        tPartner->get(iX,iY).template setField<descriptors::CUTOFF_HEAT_FLUX>(util::pow(0.5*cutoffHeatFlux, 0.25));
        // cout << cutoffKinEnergy << " " << u_00[0] << " " << u_00[1] << " " << cutoffKinEnergy_14[0] << std::endl;
 
        // tau coupling
        auto tauNS = blockLattice.get(iX,iY).template getField<descriptors::TAU_EFF>();
        auto tauAD = tPartner->get(iX,iY).template getField<descriptors::TAU_EFF>();
 
        T tau_mol_NS = 1. / blockLattice.get(iX,iY).getDynamics()->getParameters(blockLattice).template getOrFallback<descriptors::OMEGA>(0);
        T tau_mol_AD = 1. / tPartner->get(iX,iY).getDynamics()->getParameters(*tPartner).template getOrFallback<descriptors::OMEGA>(0);
 
        const T temperature = tPartner->get(iX,iY).computeRho();
 
        // computation of the bousinessq force
        T temperatureDifference = temperature - T0;
        blockLattice.get(iX,iY).template
        setField<descriptors::FORCE>(temperatureDifference*forcePrefactor);
 
        auto u = tPartner->get(iX,iY).template getField<descriptors::VELOCITY>();
        T rho, pi[util::TensorVal<DESCRIPTOR>::n], j[DESCRIPTOR::d];
        // blockLattice.get(iX,iY).computeAllMomenta(rho, u, pi);
        rho = blockLattice.get(iX,iY).computeRho();
        blockLattice.get(iX,iY).computeStress(pi);
 
        auto force = blockLattice.get(iX,iY).template getField<descriptors::FORCE>();
 
        int iPi = 0;
        for (int Alpha=0; Alpha<DESCRIPTOR::d; ++Alpha) {
          for (int Beta=Alpha; Beta<DESCRIPTOR::d; ++Beta) {
            pi[iPi] += rho/2.*(force[Alpha]*u[Beta] + u[Alpha]*force[Beta]);
            ++iPi;
          }
        }
        const Vector<T,3> piSqr = {pi[0]*pi[0], pi[1]*pi[1], pi[2]*pi[2]};
        const T PiNeqNormSqr = piSqr[0] + 2.0*piSqr[1] + piSqr[2];
        const T PiNeqNorm    = util::sqrt(PiNeqNormSqr);
 
        tPartner->get(iX,iY).computeJ(j);
        const T tmp_preFactor = invCs2_g / rho / tauAD;
        const Vector<T,2> jNeq = {(j[0] - temperature * u[0]), (j[1] - temperature * u[1])};
        const Vector<T,2> jNeqSqr = {jNeq[0]*jNeq[0], jNeq[1]*jNeq[1]};
        const T jNeqSqr_prefacor = 2. * 0.25 * (jNeq[0] + jNeq[1]) * (jNeq[0] + jNeq[1]);
 
        const T TnormSqr = jNeqSqr_prefacor*PiNeqNormSqr;
        const T Tnorm    = util::sqrt(TnormSqr);
 
        // T tau_turb_NS = 0.5*(util::sqrt(tau_mol_NS*tau_mol_NS + dynamic_cast<SmagorinskyDynamics<T,DESCRIPTOR>*>(blockLattice.get(iX,iY).getDynamics())->getPreFactor()/rho*PiNeqNorm) - tau_mol_NS);
 
        // const T tmp_A = C_nu * util::sqrt(util::sqrt(2.)/2.) * descriptors::invCs2<T,DESCRIPTOR>() * descriptors::invCs2<T,DESCRIPTOR>() * util::sqrt(PiNeqNorm / rho) * cutoffKinEnergy_14[0];
        // const T tmp_A_2 = tmp_A * tmp_A;
        // const T tmp_A_4 = tmp_A_2 * tmp_A_2;
 
        // const T tau_mol_NS_2 = tau_mol_NS * tau_mol_NS;
        // const T tau_mol_NS_3 = tau_mol_NS_2 * tau_mol_NS;
 
        // const T tmp_1_3 = 1./3.;
        // const T tmp_2_13 = util::pow(2., tmp_1_3);
        // const T tmp_3_3_12 = 3. * util::sqrt(3.);
 
        // const T tmp_sqrtA = util::sqrt(27.*tmp_A_4-4.*tmp_A_2*tau_mol_NS_3);
 
        //   // T tau_turb_NS = 1/3 ((27 A^2 + 3 util::sqrt(3) util::sqrt(27 A^4 - 4 A^2 b^3) - 2 b^3)^(1/3)/2^(1/3) + (2^(1/3) b^2)/(27 A^2 + 3 util::sqrt(3) util::sqrt(27 A^4 - 4 A^2 b^3) - 2 b^3)^(1/3) - b)
        // T tau_turb_NS = ( util::pow(27.*tmp_A_2 + tmp_3_3_12*util::sqrt(27.*tmp_A_4-4.*tmp_A_2*tau_mol_NS_3)-2.*tau_mol_NS_3, tmp_1_3) / tmp_2_13
        //                    + (tmp_2_13*tau_mol_NS_2) / util::pow(27.*tmp_A_2+tmp_3_3_12*util::sqrt(27.*tmp_A_4-4.*tmp_A_2*tau_mol_NS_3) - 2.*tau_mol_NS_3, tmp_1_3)
        //                    - tau_mol_NS
        //                    ) * tmp_1_3;
 
        // if ( tau_turb_NS != tau_turb_NS )
        //   tau_turb_NS = 0.;
 
        //cout << tau_turb_NS << " " << 27. * tmp_A_2 << " " << 4. * tau_mol_NS_3 << " " << PiNeqNorm << " " << " " << rho << std::endl;
 
        auto cutoffKinEnergy_14 = blockLattice.get(iX,iY).template getField<descriptors::CUTOFF_KIN_ENERGY>();
        auto cutoffHeatFlux_14 = tPartner->get(iX,iY).template getField<descriptors::CUTOFF_HEAT_FLUX>();
 
        const T tmp_A = C_nu * util::sqrt(util::sqrt(2.)/2.) * descriptors::invCs2<T,DESCRIPTOR>() * descriptors::invCs2<T,DESCRIPTOR>() * util::sqrt(PiNeqNorm / rho / tauNS) * cutoffKinEnergy_14;
        const T tau_turb_NS = tmp_A;
 
        // T tau_turb_AD = tau_turb_NS * tauTurbADPrefactor;
        const T tmp_B = C_alpha * descriptors::invCs2<T,DESCRIPTOR>() / rho * util::sqrt(2.0 * Tnorm * invCs2_g / tauNS / tauAD) * cutoffHeatFlux_14;
        const T tau_turb_AD = tmp_B;
        // cout << jNeq[0] << " " << jNeq[1] << " " << util::sqrt(Tnorm * invCs2_g / tauNS / tauAD) << " " << TnormSqr << std::endl;
 
        blockLattice.get(iX,iY).template setField<descriptors::TAU_EFF>(tau_mol_NS+tau_turb_NS);
        tPartner->get(iX,iY).template setField<descriptors::TAU_EFF>(tau_mol_AD+tau_turb_AD);
 
      }
    }
  }
 
}

References olb::Vector< T, D >::data(), olb::BlockLattice< T, DESCRIPTOR >::get(), olb::BlockStructureD< D >::getCellId(), olb::util::intersect(), olb::util::pow(), and olb::util::sqrt().

Here is the call graph for this function:

---

The documentation for this class was generated from the following files:

• src/dynamics/navierStokesAdvectionDiffusionCouplingPostProcessor2D.h
• src/dynamics/navierStokesAdvectionDiffusionCouplingPostProcessor2D.hh