This class computes the inamuro BC with general dynamics. More...

#include <inamuroNewtonRaphsonDynamics.h>

Inheritance diagram for olb::InamuroNewtonRaphsonDynamics< T, DESCRIPTOR, Dynamics, MOMENTA, direction, orientation >:

Collaboration diagram for olb::InamuroNewtonRaphsonDynamics< T, DESCRIPTOR, Dynamics, MOMENTA, direction, orientation >:

Public Member Functions
	InamuroNewtonRaphsonDynamics (T omega)
	Constructor.

T	computeEquilibrium (int iPop, T rho, const T u[DESCRIPTOR::d], T uSqr) const override
	Compute equilibrium distribution function.

CellStatistic< T >	collide (Cell< T, DESCRIPTOR > &cell, LatticeStatistics< T > &statistics) override
	Collision step.

T	getOmega () const
	Get local relaxation parameter of the dynamics.

void	setOmega (T omega)
	Set local relaxation parameter of the dynamics.

Public Member Functions inherited from olb::legacy::BasicDynamics< T, DESCRIPTOR, MOMENTA >
T	computeEquilibrium (int iPop, T rho, const T u[DESCRIPTOR::d]) const override
	Return iPop equilibrium for given first and second momenta.

std::type_index	id () override
	Expose unique type-identifier for RTTI.

AbstractParameters< T, DESCRIPTOR > &	getParameters (BlockLattice< T, DESCRIPTOR > &block) override
	Parameters access for legacy post processors.

Public Member Functions inherited from olb::dynamics::CustomCollision< T, DESCRIPTOR, MOMENTA >
void	initialize (Cell< T, DESCRIPTOR > &cell) override
	Initialize dynamics-specific data for cell.

T	computeRho (ConstCell< T, DESCRIPTOR > &cell) const override
	Compute particle density.

void	computeU (ConstCell< T, DESCRIPTOR > &cell, T u[DESCRIPTOR::d]) const override
	Compute fluid velocity.

void	computeJ (ConstCell< T, DESCRIPTOR > &cell, T j[DESCRIPTOR::d]) const override
	Compute fluid momentum.

void	computeStress (ConstCell< T, DESCRIPTOR > &cell, T rho, const T u[DESCRIPTOR::d], T pi[util::TensorVal< DESCRIPTOR >::n]) const override
	Compute stress tensor.

void	computeRhoU (ConstCell< T, DESCRIPTOR > &cell, T &rho, T u[DESCRIPTOR::d]) const override
	Compute fluid velocity and particle density.

void	computeAllMomenta (ConstCell< T, DESCRIPTOR > &cell, T &rho, T u[DESCRIPTOR::d], T pi[util::TensorVal< DESCRIPTOR >::n]) const override
	Compute all momenta up to second order.

void	defineRho (Cell< T, DESCRIPTOR > &cell, T rho) override
	Set particle density.

void	defineU (Cell< T, DESCRIPTOR > &cell, const T u[DESCRIPTOR::d]) override
	Set fluid velocity.

void	defineRhoU (Cell< T, DESCRIPTOR > &cell, T rho, const T u[DESCRIPTOR::d]) override
	Define fluid velocity and particle density.

void	defineAllMomenta (Cell< T, DESCRIPTOR > &cell, T rho, const T u[DESCRIPTOR::d], const T pi[util::TensorVal< DESCRIPTOR >::n]) override
	Define all momenta up to second order.

void	inverseShiftRhoU (ConstCell< T, DESCRIPTOR > &cell, T &rho, T u[DESCRIPTOR::d]) const override
	Calculate population momenta s.t. the physical momenta are reproduced by the computeRhoU.

Public Member Functions inherited from olb::Dynamics< T, DESCRIPTOR >
virtual	~Dynamics () any_platform

virtual std::string	getName () const
	Return human-readable name.

virtual CellStatistic< T >	collide (Cell< T, DESCRIPTOR > &cell)
	Perform purely-local collision step on Cell interface (legacy, to be deprecated)

void	iniEquilibrium (Cell< T, DESCRIPTOR > &cell, T rho, const T u[DESCRIPTOR::d])
	Initialize to equilibrium distribution.

void	iniRegularized (Cell< T, DESCRIPTOR > &cell, T rho, const T u[DESCRIPTOR::d], const T pi[util::TensorVal< DESCRIPTOR >::n])
	Initialize cell to equilibrium and non-equilibrum part.

Additional Inherited Members
Public Types inherited from olb::dynamics::CustomCollision< T, DESCRIPTOR, MOMENTA >
using	value_t = T

using	descriptor_t = DESCRIPTOR

using	MomentaF = typename MOMENTA::template type<DESCRIPTOR>

Public Types inherited from olb::Dynamics< T, DESCRIPTOR >
using	value_t = T

using	descriptor_t = DESCRIPTOR

Detailed Description

template<typename T, typename DESCRIPTOR, typename Dynamics, typename MOMENTA, int direction, int orientation>
class olb::InamuroNewtonRaphsonDynamics< T, DESCRIPTOR, Dynamics, MOMENTA, direction, orientation >

This class computes the inamuro BC with general dynamics.

It uses the formula from the paper by Inamuro et al. but since there is no explict solution for a lattice different from the D2Q9 and for a speed of sound c_s=q/util::sqrt(3), we have to use a Newton-Raphson algorithm to implement these boundary conditions.

Definition at line 40 of file inamuroNewtonRaphsonDynamics.h.

Constructor & Destructor Documentation

◆ InamuroNewtonRaphsonDynamics()

template<typename T , typename DESCRIPTOR , typename Dynamics , typename MOMENTA , int direction, int orientation>

olb::InamuroNewtonRaphsonDynamics< T, DESCRIPTOR, Dynamics, MOMENTA, direction, orientation >::InamuroNewtonRaphsonDynamics ( T omega )

Constructor.

Definition at line 39 of file inamuroNewtonRaphsonDynamics.hh.

  : legacy::BasicDynamics<T,DESCRIPTOR,MOMENTA>(),
    _boundaryDynamics(omega),
    clout(std::cout,"InamuroNewtonRaphsonDynamics")
{
  this->getName() = "InamuroNewtonRaphsonDynamics";
  _xi[0] = (T)1;
  for (int iDim = 1; iDim < DESCRIPTOR::d; ++iDim) {
    _xi[iDim] = T();
  }
}

References olb::Dynamics< T, DESCRIPTOR >::getName().

Here is the call graph for this function:

Member Function Documentation

◆ collide()

template<typename T , typename DESCRIPTOR , typename Dynamics , typename MOMENTA , int direction, int orientation>

CellStatistic< T > olb::InamuroNewtonRaphsonDynamics< T, DESCRIPTOR, Dynamics, MOMENTA, direction, orientation >::collide	(	Cell< T, DESCRIPTOR > &	cell,
		LatticeStatistics< T > &	statistics )

override

Collision step.

Definition at line 59 of file inamuroNewtonRaphsonDynamics.hh.

{
  typedef DESCRIPTOR L;
 
  T rho, u[L::d];
  MOMENTA().computeRhoU(cell,rho,u);
 
  constexpr auto missingIndexes = util::subIndexOutgoing<L,direction,orientation>();
  std::vector<int> knownIndexes;
  bool test[L::q];
  for (int iPop = 0; iPop < L::q; ++iPop) {
    test[iPop] = true;
  }
 
  for (unsigned iPop = 0; iPop < missingIndexes.size(); ++iPop) {
    test[missingIndexes[iPop]] = false;
  }
  for (int iPop = 0; iPop < L::q; ++iPop) {
    if (test[iPop]) {
      knownIndexes.push_back(iPop);
    }
  }
 
  T approxMomentum[L::d];
 
  computeApproxMomentum(approxMomentum,cell,rho,u,_xi,knownIndexes,missingIndexes);
 
  T error = computeError(rho, u,approxMomentum);
  int counter = 0;
 
  while (error > 1.0e-15) {
    ++counter;
 
    T gradError[L::d], gradGradError[L::d][L::d];
    computeGradGradError(gradGradError,gradError,rho,u,_xi,approxMomentum,missingIndexes);
 
    bool everythingWentFine = newtonRaphson(_xi,gradError,gradGradError);
    if ((counter > 100000) || everythingWentFine == false) {
      // if we need more that 100000 iterations or
      // if we have a problem with the inversion of the
      // jacobian matrix, we stop the program and
      // print this error message on the screen.
      clout << "Failed to converge...." << std::endl;
      clout << "Error = " << error << std::endl;
      clout << "u = (" << rho*u[0];
      for (int iD=1; iD<DESCRIPTOR::d; ++iD) {
        clout << ", " << rho*u[iD];
      }
      clout << ")" << std::endl;
      clout << "uApprox = (" << approxMomentum[0];
      for (int iD=1; iD<DESCRIPTOR::d; ++iD) {
        clout << ", " << approxMomentum[iD];
      }
      clout << ")" << std::endl;
      clout << "xi = (" << _xi[0];
      for (int iD=1; iD<DESCRIPTOR::d; ++iD) {
        clout << ", " << _xi[iD];
      }
      clout << ")" << std::endl;
 
      exit(1);
    }
 
    computeApproxMomentum(approxMomentum,cell,rho,u,_xi,knownIndexes,missingIndexes);
    error = computeError(rho, u,approxMomentum);
 
  }
 
  T uCs[L::d];
  int counterDim = 0;
  for (int iDim = 0; iDim < L::d; ++iDim) {
    if (direction == iDim) {
      ++counterDim;
      uCs[iDim] = u[iDim];
    }
    else {
      uCs[iDim] = u[iDim] + _xi[iDim+1-counterDim];
    }
  }
 
  T uCsSqr = util::normSqr<T,L::d>(uCs);
 
  for (unsigned iPop = 0; iPop < missingIndexes.size(); ++iPop) {
    cell[missingIndexes[iPop]] = computeEquilibrium(missingIndexes[iPop],_xi[0],uCs,uCsSqr);
  }
 
  _boundaryDynamics.collide(cell, statistics);
}

◆ computeEquilibrium()

template<typename T , typename DESCRIPTOR , typename Dynamics , typename MOMENTA , int direction, int orientation>

T olb::InamuroNewtonRaphsonDynamics< T, DESCRIPTOR, Dynamics, MOMENTA, direction, orientation >::computeEquilibrium	(	int	iPop,
		T	rho,
		const T	u[DESCRIPTOR::d],
		T	uSqr ) const

override

Compute equilibrium distribution function.

Definition at line 52 of file inamuroNewtonRaphsonDynamics.hh.

{
  return _boundaryDynamics.computeEquilibrium(iPop, rho, u, uSqr);
}

◆ getOmega()

template<typename T , typename DESCRIPTOR , typename Dynamics , typename MOMENTA , int direction, int orientation>

T olb::InamuroNewtonRaphsonDynamics< T, DESCRIPTOR, Dynamics, MOMENTA, direction, orientation >::getOmega ( ) const

Get local relaxation parameter of the dynamics.

Definition at line 151 of file inamuroNewtonRaphsonDynamics.hh.

{
  return _boundaryDynamics.getOmega();
}

◆ setOmega()

template<typename T , typename DESCRIPTOR , typename Dynamics , typename MOMENTA , int direction, int orientation>

void olb::InamuroNewtonRaphsonDynamics< T, DESCRIPTOR, Dynamics, MOMENTA, direction, orientation >::setOmega ( T omega )

Set local relaxation parameter of the dynamics.

Definition at line 157 of file inamuroNewtonRaphsonDynamics.hh.

{
  _boundaryDynamics.setOmega(omega);
}

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

src/boundary/inamuroNewtonRaphsonDynamics.h
src/boundary/inamuroNewtonRaphsonDynamics.hh

Public Member Functions

Additional Inherited Members

Detailed Description

Constructor & Destructor Documentation

◆ InamuroNewtonRaphsonDynamics()

Member Function Documentation

◆ collide()

◆ computeEquilibrium()

◆ getOmega()

◆ setOmega()