-
Quentin Khan authoredQuentin Khan authored
FAdaptChebKernel.hpp 8.67 KiB
#ifndef FADAPTCHEBKERNEL_HPP
#define FADAPTCHEBKERNEL_HPP
#include "Kernels/Chebyshev/FChebSymKernel.hpp"
#include "FComputeClassDescriptor.hpp"
#ifdef _OPENMP
#include <omp.h>
#endif
template<
class FReal,
class CellClass,
class ContainerClass,
class MatrixKernelClass,
int ORDER,
int NVALS = 1
>
class FAdaptChebKernel : public FChebSymKernel<FReal, CellClass, ContainerClass, MatrixKernelClass, ORDER, NVALS>{
using FBase = FChebSymKernel<FReal, CellClass, ContainerClass, MatrixKernelClass, ORDER, NVALS>;
public:
using FBase::FBase;
using FBase::M2M;
using FBase::M2L;
using FBase::L2L;
using FBase::P2P;
void setup() {
#ifdef _OPENMP
if(omp_get_thread_num() == 0)
#endif
{
std::cout << "Symetric Chebyshev adaptive kernel" << '\n';
std::cout << "Compute class: "
<< typeid(typename ComputeClassDescriptor<FReal>::type).name()
<< '\n';
}
}
template<class SymbolicData>
void P2M(typename CellClass::multipole_t* const LeafCell,
const SymbolicData * const LeafSymbData,
const ContainerClass* const SourceParticles)
{
FReal leafBoxWidth = FBase::BoxWidth / FReal(1 << LeafSymbData->getLevel());
const FPoint<FReal> LeafCellCenter =
FBase::getCellCenter(LeafSymbData->getCoordinate(), static_cast<int>(LeafSymbData->getLevel()));
FBase::Interpolator->applyP2M(LeafCellCenter, leafBoxWidth,
LeafCell->get(0), SourceParticles);
}
template<class SymbolicData>
void L2P(const typename CellClass::local_expansion_t* const LeafCell,
const SymbolicData * const LeafSymbData,
ContainerClass* const TargetParticles)
{
FReal leafBoxWidth = FBase::BoxWidth / FReal(1 << LeafSymbData->getLevel());
const FPoint<FReal> LeafCellCenter =
FBase::getCellCenter(LeafSymbData->getCoordinate(), static_cast<int>(LeafSymbData->getLevel()));
// apply Sx and Px (grad Sx)
FBase::Interpolator->applyL2PTotal(LeafCellCenter, leafBoxWidth,
LeafCell->get(0), TargetParticles);
}
template<class SymbolicData>
void P2L(typename CellClass::local_expansion_t* const local,
const SymbolicData * const symb,
const ContainerClass* const particles,
const SymbolicData * const /*source_symb*/)
{
using ComputeClass = typename ComputeClassDescriptor<FReal>::type;
constexpr std::size_t FRealCount = ComputeClassDescriptor<FReal>::count;
// Target cell: local
const FReal localCellWidth(FBase::BoxWidth / FReal(1 << symb->getLevel()));
const FPoint<FReal> localCellCenter(
FBase::getCellCenter(symb->getCoordinate(),static_cast<int>(symb->getLevel())));
// interpolation points of target (X) cell
FPoint<FReal> X[FBase::nnodes];
FChebTensor<FReal,ORDER>::setRoots(localCellCenter, localCellWidth, X);
// Particles attributes
const ComputeClass * const posX = (const ComputeClass * const) particles->getPositions()[0];
const ComputeClass * const posY = (const ComputeClass * const) particles->getPositions()[1];
const ComputeClass * const posZ = (const ComputeClass * const) particles->getPositions()[2];
const ComputeClass * const physicalValues = (const ComputeClass * const) particles->getPhysicalValues();
const FReal* pX = particles->getPositions()[0];
const FReal* pY = particles->getPositions()[1];
const FReal* pZ = particles->getPositions()[2];
const FReal* pV = particles->getPhysicalValues();
// apply P2L
for(int idxRhs = 0 ; idxRhs < NVALS ; ++idxRhs){
for (unsigned int m = 0; m<FBase::nnodes; ++m) {
ComputeClass XX = FMath::ConvertTo<ComputeClass>(X[m].getX());
ComputeClass XY = FMath::ConvertTo<ComputeClass>(X[m].getY());
ComputeClass XZ = FMath::ConvertTo<ComputeClass>(X[m].getZ());
std::size_t idxPart = 0;
// Compute using vectorization for all but the last array elements
ComputeClass tmpLocalExp = FMath::Zero<ComputeClass>();
for (;
idxPart < ((particles->getNbParticles())
/ FRealCount);
++idxPart)
{
tmpLocalExp +=
FBase::MatrixKernel->evaluate(
XX, XY, XZ,
posX[idxPart], posY[idxPart], posZ[idxPart])
* physicalValues[idxPart];
}
local->get(idxRhs)[m] += FMath::ConvertTo<FReal>(tmpLocalExp);
// Compute the last array elements one by one if they exist
if(idxPart < ((particles->getNbParticles() + FRealCount - 1) / FRealCount)) {
auto Xx = X[m].getX();
auto Xy = X[m].getY();
auto Xz = X[m].getZ();
for(idxPart = FRealCount * (particles->getNbParticles() / FRealCount);
idxPart < static_cast<std::size_t>(particles->getNbParticles());
++idxPart)
{
local->get(idxRhs)[m] +=
FBase::MatrixKernel->evaluate(
Xx, Xy, Xz,
pX[idxPart], pY[idxPart], pZ[idxPart])
* pV[idxPart];
}
}
}
}// NVALS
}
template<class SymbolicData>
void M2P(const typename CellClass::multipole_t* const pole,
const SymbolicData * const symb,
ContainerClass* const particles,
const SymbolicData * const /*source_symb*/)
{
using ComputeClass = typename ComputeClassDescriptor<FReal>::type;
constexpr std::size_t FRealCount = ComputeClassDescriptor<FReal>::count;
// Source cell: pole
const FReal poleCellWidth(FBase::BoxWidth / FReal(1 << symb->getLevel()));
const FPoint<FReal> poleCellCenter(
FBase::getCellCenter(symb->getCoordinate(),static_cast<int>(symb->getLevel())));
// interpolation points of source (Y) cell
FPoint<FReal> Y[FBase::nnodes];
FChebTensor<FReal,ORDER>::setRoots(poleCellCenter, poleCellWidth, Y);
// read positions
const ComputeClass* const posX = (const ComputeClass* const)(particles->getPositions()[0]);
const ComputeClass* const posY = (const ComputeClass* const)(particles->getPositions()[1]);
const ComputeClass* const posZ = (const ComputeClass* const)(particles->getPositions()[2]);
// get potential
ComputeClass* const physVal = (ComputeClass* const)(particles->getPhysicalValues());
ComputeClass* const potentials = (ComputeClass* const)(particles->getPotentials());
ComputeClass* const fx = (ComputeClass* const)(particles->getForcesX());
ComputeClass* const fy = (ComputeClass* const)(particles->getForcesY());
ComputeClass* const fz = (ComputeClass* const)(particles->getForcesZ());
for(int idxRhs = 0 ; idxRhs < NVALS ; ++idxRhs){
// apply M2P
for (unsigned int n=0; n<FBase::nnodes; ++n){
ComputeClass MultipoleExpansion =
FMath::ConvertTo<ComputeClass, FReal>(pole->get(idxRhs)[n]);
ComputeClass YX = FMath::ConvertTo<ComputeClass, FReal>(Y[n].getX());
ComputeClass YY = FMath::ConvertTo<ComputeClass, FReal>(Y[n].getY());
ComputeClass YZ = FMath::ConvertTo<ComputeClass, FReal>(Y[n].getZ());
for(std::size_t idxPart = 0;
idxPart < ( (particles->getNbParticles() + FRealCount - 1)
/ FRealCount);
++idxPart)
{
ComputeClass Kxy[1];
ComputeClass dKxy[3];
FBase::MatrixKernel->evaluateBlockAndDerivative(
posX[idxPart], posY[idxPart], posZ[idxPart],
YX, YY, YZ,
Kxy,dKxy);
potentials[idxPart] += Kxy[0] * MultipoleExpansion;
fx[idxPart] += dKxy[0] * physVal[idxPart] * MultipoleExpansion;
fy[idxPart] += dKxy[1] * physVal[idxPart] * MultipoleExpansion;
fz[idxPart] += dKxy[2] * physVal[idxPart] * MultipoleExpansion;
}
}// Particles
}// NVALS
}
};
#endif /* FADAPTCHEBKERNEL_HPP */