Reply To: Manipulating populations with postprocessors on GPU

[jmctighe]

Hi Adrian,

You are definitely correct, the interpolation point is causing the problem. Although, instead of it being a point outside of the domain, the value of x_m is being calculated as nan. It seems that all of the cell or parameter values used by the post processor are zeros. It might be important to note that the error is occurring during the lattice preparation/boundary setting procedure, before setting the initial condition. Interested to know your thoughts on this.

I’ve included the prepareLattice and setShearPeriodic functions below. I updated the post processor to use DOMAIN_EXTENT and LATTICE_SHEAR as parameters, but this had no effect on the error.

prepareLattice:

// Set up the geometry of the simulation
void prepareLattice(SuperLattice<T,DESCRIPTOR>& sLattice,
                    const UnitConverter<T,DESCRIPTOR>& converter,
                    SuperGeometry<T,3>& superGeometry,
                    ShearPeriodicParameters<T>& shearPeriodicParameters)
{
  OstreamManager clout(std::cout,"prepareLattice");
  clout << "Prepare Lattice ..." << std::endl;

  /// Material=1 -->bulk dynamics
  sLattice.defineDynamics<BulkDynamics>(superGeometry, 1);
  sLattice.addPostProcessor<stage::PostStream>(superGeometry.getMaterialIndicator({1,2}),meta::id<TrackCovarVelocity>{});
  /// Material = 2 --> boundary node
  sLattice.defineDynamics<BulkDynamics>(superGeometry, 2);
  /// === Set Initial Conditions == ///
  sLattice.setParameter<LATTICE_SHEAR>(shearPeriodicParameters.latticeShear);
  sLattice.setParameter<DOMAIN_EXTENT>(shearPeriodicParameters.extents);

  setshearPeriodic(sLattice, superGeometry, 2, shearPeriodicParameters);
  setInitialConditions(sLattice, converter, superGeometry);

  sLattice.setParameter<descriptors::OMEGA>( converter.getLatticeRelaxationFrequency() );
  //sLattice.setParameter<collision::LES::SMAGORINSKY>(T(0.12));

  // Make the lattice ready for simulation
  sLattice.initialize();

  clout << "Prepare Lattice ... OK" << std::endl;
}

setShearPeriodic.h:

#ifndef SET_SHEAR_PERIODIC_H
#define SET_SHEAR_PERIODIC_H

#include "shearPeriodicPostProcessor.h"

namespace olb {

//======================================================================
// ======== Shear Periodic parameters ======//
//======================================================================

template <typename T>
struct ShearPeriodicParameters{
  ///  distance from cell to velocity sampling point in lattice units
  Vector<T,3> extents;
  T latticeWallDistance;
  T latticeShear;
};

template<typename T, typename DESCRIPTOR>
void setshearPeriodic(SuperLattice<T, DESCRIPTOR>& sLattice,
                           FunctorPtr<SuperIndicatorF<T,DESCRIPTOR::d>>&& boundaryIndicator,
                           FunctorPtr<SuperIndicatorF<T,DESCRIPTOR::d>>&& bulkIndicator,
                           ShearPeriodicParameters<T>& shearPeriodicParameters,
                           IndicatorF<T,DESCRIPTOR::d>* indicatorAnalyticalBoundary = nullptr)
{
  int _overlap = 3;
  OstreamManager clout(std::cout, "shearPeriodicSetter");

  AnalyticalConst<DESCRIPTOR::d, T, T> zeroVector(T(0.), T(0.));

  sLattice.template defineField<descriptors::Y1>(std::move(boundaryIndicator), zeroVector);

  auto& load = sLattice.getLoadBalancer();
  for (int iC=0; iC < load.size(); ++iC) {
    setshearPeriodic<T,DESCRIPTOR>(sLattice.getBlock(iC),
                          (bulkIndicator->getBlockIndicatorF(iC)).getBlockGeometry(),
                          boundaryIndicator->getBlockIndicatorF(iC),
                          bulkIndicator->getBlockIndicatorF(iC),
                          shearPeriodicParameters,
                          indicatorAnalyticalBoundary);
  }
  addPoints2CommBC<T,DESCRIPTOR>(sLattice, std::forward<decltype(boundaryIndicator)>(boundaryIndicator), _overlap);
  clout.setMultiOutput(false);
  sLattice.template setParameter<descriptors::SAMPLING_DISTANCE>( shearPeriodicParameters.extents[1] );
}

template<typename T, typename DESCRIPTOR>
void setshearPeriodic(SuperLattice<T, DESCRIPTOR>& sLattice,
                           SuperGeometry<T,DESCRIPTOR::d>& superGeometry,
                           int materialOfSolidObstacle,
                           ShearPeriodicParameters<T>& shearPeriodicParameters,
                           IndicatorF<T,DESCRIPTOR::d>* indicatorAnalyticalBoundary = nullptr,
                           std::vector<int> bulkMaterials = std::vector<int>(1,1))
{
  // Getting the indicators by material numbers and calling the superLattice method via the indicators:
  setshearPeriodic<T,DESCRIPTOR>(sLattice,
                      FunctorPtr<SuperIndicatorF<T,DESCRIPTOR::d>>(superGeometry.getMaterialIndicator(materialOfSolidObstacle)),
                      FunctorPtr<SuperIndicatorF<T,DESCRIPTOR::d>>(superGeometry.getMaterialIndicator(std::move(bulkMaterials))),
                      shearPeriodicParameters,
                      indicatorAnalyticalBoundary);
}

template<typename T, typename DESCRIPTOR>
void setshearPeriodic(BlockLattice<T,DESCRIPTOR>& block,
                           BlockGeometry<T,DESCRIPTOR::d>& blockGeometry,
                           BlockIndicatorF<T,DESCRIPTOR::d>& boundaryIndicator,
                           BlockIndicatorF<T,DESCRIPTOR::d>& bulkIndicator,
                           ShearPeriodicParameters<T>& shearPeriodicParameters,
                           IndicatorF<T,DESCRIPTOR::d>* indicatorAnalyticalBoundary = nullptr)
{
  OstreamManager clout(std::cout, "shearPeriodicSetter");
  clout.setMultiOutput(true);

  // Defining boundary distance y1 and distance to exchange location y2 in wall normal direction for every boundary cell
  const T deltaR = blockGeometry.getDeltaR();
  // for each solid cell:
  int kMax = 2;
  for(int k=1; k < kMax; k++){
    block.forSpatialLocations([&](LatticeR<DESCRIPTOR::d> solidLatticeR) {
      // Check if cell is solid cell
      if (boundaryIndicator(solidLatticeR)) {
        for (int iPop=1; iPop < DESCRIPTOR::q; ++iPop) {
          Vector<T,DESCRIPTOR::d> boundaryLatticeR(solidLatticeR + k*descriptors::c<DESCRIPTOR>(iPop));
          if (blockGeometry.getNeighborhoodRadius(boundaryLatticeR) >= 1) {
            if (blockGeometry.isInside(boundaryLatticeR)) {
              Vector<T,DESCRIPTOR::d> boundaryPhysR{ };
              blockGeometry.getPhysR(boundaryPhysR,boundaryLatticeR);
              // check if neighbor is fluid cell
              if (bulkIndicator(boundaryLatticeR)) {
                if (indicatorAnalyticalBoundary) {
                  // Calculate surface normal
                  Vector<T,DESCRIPTOR::d> normal = indicatorAnalyticalBoundary->surfaceNormal(boundaryPhysR, deltaR);
                  T y1 = 0.;
                  // Calculate boundary distance y1 in surface normal direction
                  indicatorAnalyticalBoundary->distance(y1, boundaryPhysR, normal, blockGeometry.getIcGlob());
                  if(y1 == T(0)) {
                    for ( int iD = 0; iD < DESCRIPTOR::d; iD++ ) {
                      normal[iD] *= T(-1);
                    }
                    indicatorAnalyticalBoundary->distance(y1, boundaryPhysR, normal, blockGeometry.getIcGlob());
                  }
                  y1 /= deltaR; // y1 in lattice units
                  if ( util::abs(y1) > T(k) ) {
                    y1 = util::sign(y1) * T(k);
                  }
                  block.get(boundaryLatticeR).template setField<descriptors::Y1>(-y1*normal);
                } else {
                  Vector<T,DESCRIPTOR::d> normal;
                  for(int jPop = 0; jPop < DESCRIPTOR::q; jPop++){
                    Vector<T,DESCRIPTOR::d> boundaryLatticeR2(boundaryLatticeR + k*descriptors::c<DESCRIPTOR>(jPop));
                    if(boundaryIndicator(boundaryLatticeR2)){
                      for ( int iD = 0; iD < DESCRIPTOR::d; iD++ ) {
                        normal[iD] += k*descriptors::c<DESCRIPTOR>(jPop,iD);
                      }
                    }
                  }
                  T normalNorm = util::norm<DESCRIPTOR::d>(normal);
                  if(normalNorm != T(0)) {
                    for ( int iD = 0; iD < DESCRIPTOR::d; iD++ ) {
                      normal[iD] /= normalNorm;
                    }
                  }
                  for ( int iD = 0; iD < DESCRIPTOR::d; iD++ ) {
                    if(util::abs(normal[iD]) < T(0.4)) {
                      normal[iD] = T(0);
                    }
                  }
                  normalNorm = util::norm<DESCRIPTOR::d>(normal);
                  for ( int iD = 0; iD < DESCRIPTOR::d; iD++ ) {
                    normal[iD] /= normalNorm;
                  }
                  if(normal[0] != T(0) && normal[0] > T(0)) {
                    normal[0] /= normal[0];
                  }
                  else if(normal[0] != T(0) && normal[0] < T(0)) {
                    normal[0] /= normal[0];
                    normal[0] *= T(-1.);
                  }
                  if(normal[1] != T(0) && normal[1] > T(0)) {
                    normal[1] /= normal[1];
                  }
                  else if(normal[1] != T(0) && normal[1] < T(0)) {
                    normal[1] /= normal[1];
                    normal[1] *= T(-1.);
                  }
                  if ( DESCRIPTOR::d == 3) {
                    if(normal[2] != T(0) && normal[2] > T(0)) {
                      normal[2] /= normal[2];
                    }
                    else if(normal[2] != T(0) && normal[2] < T(0)) {
                      normal[2] /= normal[2];
                      normal[2] *= T(-1.);
                    }
                  }
                  if(normalNorm != T(0)) {
                    auto field = block.get(boundaryLatticeR).template getField<descriptors::Y1>();
                    if( util::norm<DESCRIPTOR::d>(field) == T(0)) {
                      T y1 = (k-1) + shearPeriodicParameters.latticeWallDistance;
                      block.get(boundaryLatticeR).template setField<descriptors::Y1>(-y1*normal);
                    }
                  }
                }
                //T pi[util::TensorVal<DESCRIPTOR>::n] {T(0)};
                //block.get(boundaryLatticeR).template setField<descriptors::TENSOR>(pi);
                // Setting shear periodic post processor
                block.get(boundaryLatticeR).template setField<descriptors::LOCATION>(boundaryLatticeR);
                block.addPostProcessor(typeid(stage::PreCollide),
                        boundaryLatticeR,
                        meta::id<shearPeriodicPostProcessor>{});
              }
            }
          }
        }
      }
    });
  }
  
}

}

#endif