Commit 57e8751a authored by berenger-bramas's avatar berenger-bramas

start to develop a clean kernel and to put data in separate files between

real computation and useful calculus like spherical harmonic.
This is just a draft.

git-svn-id: svn+ssh://scm.gforge.inria.fr/svn/scalfmm/scalfmm/trunk@292 2616d619-271b-44dc-8df4-d4a8f33a7222
parent fb8920ef
#ifndef FELECBASICKERNELS_HPP
#define FELECBASICKERNELS_HPP
// [--License--]
#include "../Components/FAbstractKernels.hpp"
#include "../Utils/FGlobal.hpp"
#include "../Utils/FTrace.hpp"
#include "FHarmonic.hpp"
/**
* @author Berenger Bramas (berenger.bramas@inria.fr)
* @brief
* Please read the license
*
*/
template< class ParticleClass, class CellClass, class ContainerClass>
class FElecForcesKernels : public FAbstractKernels<ParticleClass,CellClass,ContainerClass> {
int devP;
FHarmonic harmonic;
public:
FElecForcesKernels(const int inDevP)
: devP(inDevP), harmonic(inDevP) {
}
/** Default destructor */
~FElecForcesKernels(){
}
/** P2M with a cell and all its particles */
void P2M(CellClass* const inPole, const ContainerClass* const ) {
FComplexe* const cellMultiPole = inPole->getMultipole();
// Copying the position is faster than using cell position
const F3DPosition polePosition = inPole->getPosition();
// For all particles in the leaf box
typename ContainerClass::ConstBasicIterator iterParticle(*inParticles);
while( iterParticle.hasNotFinished()){
// P2M
particleToMultiPole(cellMultiPole, polePosition, iterParticle.data());
iterParticle.gotoNext();
}
}
/** Do nothing */
void M2M(CellClass* const FRestrict , const CellClass*const FRestrict *const FRestrict , const int ) {
}
/** Do nothing */
virtual void M2L(CellClass* const FRestrict , const CellClass* [], const int , const int ) {
}
/** Do nothing */
void L2L(const CellClass* const FRestrict , CellClass* FRestrict *const FRestrict , const int ) {
}
/** L2P with a cell and all its particles */
void L2P(const CellClass* const local, ContainerClass* const particles){
const FComplexe* const cellLocal = local->getLocal();
// Copying the position is faster than using cell position
const F3DPosition localPosition = local->getPosition();
// For all particles in the leaf box
typename ContainerClass::BasicIterator iterTarget(*particles);
while( iterTarget.hasNotFinished() ){
// L2P
localToParticle(&iterTarget.data(), localPosition, cellLocal);
iterTarget.gotoNext();
}
}
/** This P2P has to be used when target != sources
* It will proceed an direct interation no mutual
*
* It takes all the target particles from the current leaf,
* then it computes the sources/targets interaction in this leaf,
* then it computes the sources/targets inteactions between this leaf and the
* neighbors.
*/
void P2P(ContainerClass* const FRestrict targets, const ContainerClass* const FRestrict sources,
const ContainerClass* const directNeighbors[26], const int size) {
{ // Compute interaction in this leaf
typename ContainerClass::BasicIterator iterTarget(*targets);
while( iterTarget.hasNotFinished() ){
// We copy the target particle to work with a particle in the heap
ParticleClass target( iterTarget.data() );
// For all the source particles in the same leaf
typename ContainerClass::ConstBasicIterator iterSameBox(*sources);
while( iterSameBox.hasNotFinished() ){
//(&iterSameBox.data() != &iterTarget.data())
directInteraction(target, iterSameBox.data());
iterSameBox.gotoNext();
}
// Set data and progress
iterTarget.setData(target);
iterTarget.gotoNext();
}
}
{ // Compute interactions with other leaves
// For all the neigbors leaves
for(int idxDirectNeighbors = 0 ; idxDirectNeighbors < size ; ++idxDirectNeighbors){
// For all particles in current leaf
typename ContainerClass::BasicIterator iterTarget(*targets);
while( iterTarget.hasNotFinished() ){
// For all the particles in the other leaf
typename ContainerClass::ConstBasicIterator iterSource(*directNeighbors[idxDirectNeighbors]);
while( iterSource.hasNotFinished() ){
directInteraction(target, iterSource.data());
iterSource.gotoNext();
}
// Set data and progress
iterTarget.setData(target);
iterTarget.gotoNext();
}
}
}
}
/** This P2P has to be used when target == sources
* It will proceed a direct interation >> mutual
*
* It takes all the particles from the current leaf,
* then it computes the interactions in this leaf,
* then it computes the inteactions between this leaf and the
* neighbors.
*/
void P2P(const MortonIndex inCurrentIndex,
ContainerClass* const FRestrict targets, const ContainerClass* const FRestrict /*sources*/,
ContainerClass* const directNeighbors[26], const MortonIndex inNeighborsIndex[26], const int size){
{ // Compute interaction in this leaf
typename ContainerClass::BasicIterator iterTarget(*targets);
while( iterTarget.hasNotFinished() ){
// We copy the target particle to work with a particle in the heap
ParticleClass target( iterTarget.data() );
// For all particles after the current one
typename ContainerClass::ConstBasicIterator iterSameBox = iterTarget;
iterSameBox.gotoNext();
while( iterSameBox.hasNotFinished() ){
directInteractionMutual(&target, &iterSameBox.data());
iterSameBox.gotoNext();
}
// Set data and progress
iterTarget.setData(target);
iterTarget.gotoNext();
}
}
{ // Compute interactions with other leaves
// For all the neigbors leaves
for(int idxDirectNeighbors = 0 ; idxDirectNeighbors < size ; ++idxDirectNeighbors){
if(inCurrentIndex < inNeighborsIndex[idxDirectNeighbors] ){
// For all particles in current leaf
typename ContainerClass::BasicIterator iterTarget(*targets);
while( iterTarget.hasNotFinished() ){
ParticleClass target( iterTarget.data() );
// For all the particles in the other leaf
typename ContainerClass::ConstBasicIterator iterSource(*directNeighbors[idxDirectNeighbors]);
while( iterSource.hasNotFinished() ){
directInteractionMutual(&target, &iterSource.data());
iterSource.gotoNext();
}
// Set data and progress
iterTarget.setData(target);
iterTarget.gotoNext();
}
}
}
}
}
///////////////////////////////////////////////////////////////////////////////
// Periodic
///////////////////////////////////////////////////////////////////////////////
/** Before Downward */
void M2L(CellClass* const FRestrict , const CellClass* [189], FTreeCoordinate [189], const int , const int ) {
}
/** After Downward */
void P2P(const MortonIndex ,
ContainerClass* const FRestrict , const ContainerClass* const FRestrict ,
ContainerClass* const [26], const FTreeCoordinate [26], const int ) {
}
///////////////////////////////////////////////////////////////////////////////
// Computation
///////////////////////////////////////////////////////////////////////////////
/** P2M computation
*/
void particleToMultiPole(FComplexe*const cellMultiPole, const F3DPosition& inPolePosition ,
const CellClass& particle){
harmonic.computeInner( FSpherical(particle.getPosition() - inPolePosition) );
FReal minus_one_pow_j = 1.0;//(-1)^j
const FReal valueParticle = particle.getPhysicalValue();
int p_exp_term = 0; // p_Y_term
for(int jIdx = 0 ; jIdx <= devP ; ++jIdx){
for(int kIdx = 0 ; kIdx <= jIdx ; ++kIdx, ++p_Y_term, ++p_exp_term){
harmonic.result(p_exp_term).mulRealAndImag( valueParticle * minus_one_pow_j );
cellMultiPole[p_exp_term] += harmonic.result(p_exp_term);
}
minus_one_pow_j = -minus_one_pow_j;
}
}
void localToParticle(ParticleClass*const particle, const F3DPosition& local_position,
const FComplexe*const local_exp){
FReal force_vector_in_local_base_x = 0;
FReal force_vector_in_local_base_y = 0;
FReal force_vector_in_local_base_z = 0;
const FSpherical spherical(particle->getPosition() - local_position);
harmonic.computeInnerTheta( spherical );
// The maximum degree used here will be P.
const FComplexe* p_Y_term = current_thread_Y+1;
const FComplexe* p_Y_theta_derivated_term = current_thread_Y_theta_derivated+1;
const FComplexe* p_local_exp_term = local_exp+1;
for (int j = 1 ; j <= devP ; ++j ){
FReal exp_term_aux_real = 0.0;
FReal exp_term_aux_imag = 0.0;
// k=0:
// F_r:
exp_term_aux_real = ( (p_Y_term->getReal() * p_local_exp_term->getReal()) - (p_Y_term->getImag() * p_local_exp_term->getImag()) );
exp_term_aux_imag = ( (p_Y_term->getReal() * p_local_exp_term->getImag()) + (p_Y_term->getImag() * p_local_exp_term->getReal()) );
force_vector_in_local_base_x = ( force_vector_in_local_base_x + FReal(j) * exp_term_aux_real );
// F_phi: k=0 => nothing to do for F_phi
// F_theta:
exp_term_aux_real = ( (p_Y_theta_derivated_term->getReal() * p_local_exp_term->getReal()) - (p_Y_theta_derivated_term->getImag() * p_local_exp_term->getImag()) );
exp_term_aux_imag = ( (p_Y_theta_derivated_term->getReal() * p_local_exp_term->getImag()) + (p_Y_theta_derivated_term->getImag() * p_local_exp_term->getReal()) );
force_vector_in_local_base_y = ( force_vector_in_local_base_y + exp_term_aux_real );
++p_local_exp_term;
++p_Y_term;
++p_Y_theta_derivated_term;
// k>0:
for (int k=1; k<=j ;++k, ++p_local_exp_term, ++p_Y_term, ++p_Y_theta_derivated_term){
// F_r:
exp_term_aux_real = ( (p_Y_term->getReal() * p_local_exp_term->getReal()) - (p_Y_term->getImag() * p_local_exp_term->getImag()) );
exp_term_aux_imag = ( (p_Y_term->getReal() * p_local_exp_term->getImag()) + (p_Y_term->getImag() * p_local_exp_term->getReal()) );
force_vector_in_local_base_x = (force_vector_in_local_base_x + FReal(2 * j) * exp_term_aux_real );
// F_phi:
force_vector_in_local_base_z = ( force_vector_in_local_base_z - FReal(2 * k) * exp_term_aux_imag);
// F_theta:
exp_term_aux_real = ( (p_Y_theta_derivated_term->getReal() * p_local_exp_term->getReal()) - (p_Y_theta_derivated_term->getImag() * p_local_exp_term->getImag()) );
exp_term_aux_imag = ( (p_Y_theta_derivated_term->getReal() * p_local_exp_term->getImag()) + (p_Y_theta_derivated_term->getImag() * p_local_exp_term->getReal()) );
force_vector_in_local_base_y = (force_vector_in_local_base_y + FReal(2.0) * exp_term_aux_real );
}
}
// We want: - gradient(POTENTIAL_SIGN potential).
// The -(- 1.0) computing is not the most efficient programming ...
force_vector_in_local_base_x = ( force_vector_in_local_base_x * FReal(-1.0) / spherical.getR());
force_vector_in_local_base_y = ( force_vector_in_local_base_y * FReal(-1.0) / spherical.getR());
force_vector_in_local_base_z = ( force_vector_in_local_base_z * FReal(-1.0) / (spherical.getR() * spherical.getSinTheta()));
/////////////////////////////////////////////////////////////////////
//spherical_position_Set_ph
//FMB_INLINE COORDINATES_T angle_Convert_in_MinusPi_Pi(COORDINATES_T a){
FReal ph = FMath::Fmod(spherical.getPhi(), FReal(2)*FMath::FPi);
if (ph > M_PI) ph -= FReal(2) * FMath::FPi;
if (ph < -M_PI + FMath::Epsilon) ph += FReal(2) * FMath::Epsilon;
//spherical_position_Set_th
FReal th = FMath::Fmod(FMath::ACos(spherical.getCosTheta()), FReal(2) * FMath::FPi);
if (th < 0.0) th += 2*FMath::FPi;
if (th > FMath::FPi){
th = 2*FMath::FPi - th;
//spherical_position_Set_ph(p, spherical_position_Get_ph(p) + M_PI);
ph = FMath::Fmod(ph + FMath::FPi, 2*FMath::FPi);
if (ph > M_PI) ph -= 2*FMath::FPi;
if (ph < -M_PI + FMath::Epsilon) ph += 2 * FMath::Epsilon;
th = FMath::Fmod(th, 2*FMath::FPi);
if (th > M_PI) th -= 2*FMath::FPi;
if (th < -M_PI + FMath::Epsilon) th += 2 * FMath::Epsilon;
}
//spherical_position_Set_r
//FReal rh = spherical.r;
if (spherical.r < 0){
//rh = -spherical.r;
//spherical_position_Set_ph(p, M_PI - spherical_position_Get_th(p));
ph = FMath::Fmod(FMath::FPi - th, 2*FMath::FPi);
if (ph > M_PI) ph -= 2*FMath::FPi;
if (ph < -M_PI + FMath::Epsilon) ph += 2 * FMath::Epsilon;
//spherical_position_Set_th(p, spherical_position_Get_th(p) + M_PI);
th = FMath::Fmod(th + FMath::FPi, 2*FMath::FPi);
if (th < 0.0) th += 2*FMath::FPi;
if (th > FMath::FPi){
th = 2*FMath::FPi - th;
//spherical_position_Set_ph(p, spherical_position_Get_ph(p) + M_PI);
ph = FMath::Fmod(ph + FMath::FPi, 2*FMath::FPi);
if (ph > M_PI) ph -= 2*FMath::FPi;
if (ph < -M_PI + FMath::Epsilon) ph += 2 * FMath::Epsilon;
th = FMath::Fmod(th, 2*FMath::FPi);
if (th > M_PI) th -= 2*FMath::FPi;
if (th < -M_PI + FMath::Epsilon) th += 2 * FMath::Epsilon;
}
}
const FReal cos_theta = FMath::Cos(th);
const FReal cos_phi = FMath::Cos(ph);
const FReal sin_theta = FMath::Sin(th);
const FReal sin_phi = FMath::Sin(ph);
FReal force_vector_tmp_x = (
cos_phi * sin_theta * force_vector_in_local_base_x +
cos_phi * cos_theta * force_vector_in_local_base_y +
(-sin_phi) * force_vector_in_local_base_z);
FReal force_vector_tmp_y = (
sin_phi * sin_theta * force_vector_in_local_base_x +
sin_phi * cos_theta * force_vector_in_local_base_y +
cos_phi * force_vector_in_local_base_z);
FReal force_vector_tmp_z = (
cos_theta * force_vector_in_local_base_x +
(-sin_theta) * force_vector_in_local_base_y);
const FReal physicalValue = particle->getPhysicalValue();
force_vector_tmp_x *= physicalValue;
force_vector_tmp_y *= physicalValue;
force_vector_tmp_z *= physicalValue;
particle->incForces( force_vector_tmp_x, force_vector_tmp_y, force_vector_tmp_z );
particle->incPotential(expansion_Evaluate_local_with_Y_already_computed(local_exp));
}
/** P2P mutual interaction
* F = q * q' / r²
*/
void directInteractionMutual(ParticleClass*const FRestrict target, ParticleClass*const FRestrict source){
FReal dx = -(target->getPosition().getX() - source->getPosition().getX());
FReal dy = -(target->getPosition().getY() - source->getPosition().getY());
FReal dz = -(target->getPosition().getZ() - source->getPosition().getZ());
FReal inv_square_distance = FReal(1.0) / (dx*dx + dy*dy + dz*dz);
FReal inv_distance = FMath::Sqrt(inv_square_distance);
inv_square_distance *= inv_distance;
inv_square_distance *= target->getPhysicalValue() * source->getPhysicalValue();
dx *= inv_square_distance;
dy *= inv_square_distance;
dz *= inv_square_distance;
target->incForces( dx, dy, dz);
target->incPotential( inv_distance * source->getPhysicalValue() );
source.incForces( (-dx), (-dy), (-dz));
source.incPotential( inv_distance * target->getPhysicalValue() );
}
/** P2P NO mutual interaction
* F = q * q' / r²
*/
void directInteraction(ParticleClass*const FRestrict target, const ParticleClass& source){
FReal dx = -(target->getPosition().getX() - source.getPosition().getX());
FReal dy = -(target->getPosition().getY() - source.getPosition().getY());
FReal dz = -(target->getPosition().getZ() - source.getPosition().getZ());
FReal inv_square_distance = FReal(1.0) / (dx*dx + dy*dy + dz*dz);
FReal inv_distance = FMath::Sqrt(inv_square_distance);
inv_square_distance *= inv_distance;
inv_square_distance *= target->getPhysicalValue() * source.getPhysicalValue();
dx *= inv_square_distance;
dy *= inv_square_distance;
dz *= inv_square_distance;
target->incForces( dx, dy, dz);
target->incPotential( inv_distance * source.getPhysicalValue() );
}
};
#endif //FELECBASICKERNELS_HPP
// [--END--]
#ifndef FHARMONIC_HPP
#define FHARMONIC_HPP
#include "../Utils/FComplexe.hpp"
/** Todo move this class outside when kernels will be developed */
class FSpherical {
FReal r, cosTheta, sinTheta, phi;
public:
FSpherical(): r(0), cosTheta(0), sinTheta(0), phi(0){
}
explicit FSpherical(const F3DPosition& inVector){
const FReal x2y2 = (inVector.getX() * inVector.getX()) + (inVector.getY() * inVector.getY());
this->r = FMath::Sqrt( x2y2 + (inVector.getZ() * inVector.getZ()));
this->phi = FMath::Atan2(inVector.getY(),inVector.getX());
this->cosTheta = inVector.getZ() / outSphere->r; // cos_th = z/r
this->sinTheta = FMath::Sqrt(x2y2) / outSphere->r; // sin_th = sqrt(x^2 + y^2)/r
}
FReal getR() const{
return r;
}
FReal getCosTheta() const{
return cosTheta;
}
FReal getSinTheta() const{
return sinTheta;
}
FReal getPhi() const{
return phi;
}
};
class FHarmonic {
const int devP;
const int expSize;
FComplexe* harmonic;
FComplexe* cosSin;
FReal* legendre;
FComplexe* thetaDerivatedResult;
FReal* sphereHarmoInnerCoef;
FReal* sphereHarmoOuterCoef;
void sphericalHarmonicInitialize(){
FReal factOuter = 1.0;
for(int idxP = 0 ; idxP <= devP; factOuter *= FReal(++idxP) ){
sphereHarmoOuterCoef[idxP] = factOuter;
}
FReal* currentInner = sphereHarmoInnerCoef;
FReal factInner = 1.0;
FReal powN1idxP = 1.0;
for(int idxP = 0 ; idxP <= this->devP ; factInner *= FReal(++idxP), powN1idxP = -powN1idxP){
for(int idxMP = 0, fact_l_m = int(factInner); idxMP <= idxP ; fact_l_m *= idxP+(++idxMP), ++currentInner){
*currentInner = powN1idxP / FReal(fact_l_m);
}
}
}
void computeCosSin(const FReal inPhi){
for(int idxl = 0, idxlMod4 = 0; idxl <= devP ; ++idxl, ++idxlMod4){
if(idxlMod4 == 4) idxlMod4 = 0;
const FReal angle = FReal(idxl) * inPhi + this->PiArrayOuter[idxlMod4];
cosSin[idxl].setReal( FMath::Sin(angle + FMath::FPiDiv2) );
cosSin[idxl].setImag( FMath::Sin(angle) );
}
}
void computeLegendre(const FReal inCosTheta, const FReal inSinTheta){
int idxCurrent = 0;
legendre[idxCurrent++] = 1.0; // P_0^0(cosTheta) = 1
legendre[idxCurrent++] = inCosTheta; // P_1^{0} using (3)
// Compute P_1^1 using (2 bis) and store it into results_array
const FReal invSinTheta = -inSinTheta; // somx2 => -sinTheta
legendre[idxCurrent++] = invSinTheta;
// l>1:
int idxCurrent1m = 1; //pointer on P_{l-1}^m P_1^0
int idxCurrent2m = 0; //pointer on P_{l-2}^m P_0^0
FReal fact = 3.0;
// Remark: p_results_array_l_minus_1_m and p_results_array_l_minus_2_m
// just need to be incremented at each iteration.
for(int idxl = 2; idxl <= devP ; ++idxl ){
for( int idxm = 0; idxm <= idxl - 2 ; ++idxm , ++idxCurrent , ++idxCurrent1m , ++idxCurrent2m ){
// Compute P_l^m, l >= m+2, using (1) and store it into results_array:
legendre[idxCurrent] = (inCosTheta * FReal( 2 * idxl - 1 ) * legendre[idxCurrent1m] - FReal( idxl + idxm - 1 )
* legendre[idxCurrent2m] ) / FReal( idxl - idxm );
}
// p_results_array_l_minus_1_m now points on P_{l-1}^{l-1}
// Compute P_l^{l-1} using (3) and store it into ptrResults:
legendre[idxCurrent++] = inCosTheta * FReal( 2 * idxl - 1 ) * legendre[idxCurrent1m];
// Compute P_l^l using (2 bis) and store it into results_array:
legendre[idxCurrent++] = fact * invSinTheta * legendre[idxCurrent1m];
fact += FReal(2.0);
++idxCurrent1m;
}
}
public:
explicit FHarmonic(const int inDevP)
: devP(inDevP),expSize(int(((inDevP)+1) * ((inDevP)+2) * 0.5)),
harmonic(0), cosSin(0), legendre(0), thetaDerivatedResult(0),
sphereHarmoInnerCoef(0), sphereHarmoOuterCoef(0) {
harmonic = new FComplexe[expSize];
cosSin = new FComplexe[2 * devP + 1];
legendre = new FReal[expSize];
thetaDerivatedResult = new FComplexe[expSize];
sphereHarmoInnerCoef = new FReal[expSize];
sphereHarmoOuterCoef = new FReal[devP + 1];
sphericalHarmonicInitialize();
}
~FHarmonic(){
delete[] harmonic;
delete[] cosSin;
delete[] legendre;
delete[] thetaDerivatedResult;
delete[] sphereHarmoInnerCoef;
delete[] sphereHarmoOuterCoef;
}
FComplexe* data(){
return harmonic;
}
const FComplexe* data() const {
return harmonic;
}
FComplexe& result(const int index){
return harmonic[index];
}
const FComplexe& result(const int index) const{
return harmonic[index];
}
void computeInner(const FSpherical& inSphere){
computeCosSin(inSphere.getPhi());
// p_associated_Legendre_function_Array
computeLegendre(inSphere.getCosTheta(), inSphere.getSinTheta());
FComplexe* currentResult = harmonic;
int idxLegendre = 0;//ptr_associated_Legendre_function_Array
int idxSphereHarmoCoef = 0;
FReal idxRl = 1.0 ;
for(int idxl = 0; idxl <= devP ; ++idxl, idxRl *= inSphere.getR()){
for(int idxm = 0 ; idxm <= idxl ; ++idxm, ++currentResult, ++idxSphereHarmoCoef, ++idxLegendre){
const FReal magnitude = sphereHarmoInnerCoef[idxSphereHarmoCoef] * idxRl * legendre[idxLegendre];
currentResult->setReal( magnitude * cosSin[idxm].getReal() );
currentResult->setImag( magnitude * cosSin[idxm].getImag() );
}
}
}
void computeOuter(const FSpherical& inSphere){
computeCosSin(inSphere.getPhi());
// p_associated_Legendre_function_Array
computeLegendre(inSphere.getCosTheta(), inSphere.getSinTheta());
int idxLegendre = 0;
FComplexe* currentResult = harmonic;
const FReal invR = 1/inSphere.getR();
FReal idxRl1 = invR;
for(int idxl = 0 ; idxl <= devP ; ++idxl, idxRl1 *= invR){
for(int idxm = 0 ; idxm <= idxl ; ++idxm, ++currentResult, ++idxLegendre){
const FReal magnitude = this->sphereHarmoOuterCoef[idxl-idxm] * idxRl1 * legendre[idxLegendre];
currentResult->setReal( magnitude * cosSin[idxm].getReal() );
currentResult->setImag( magnitude * cosSin[idxm].getImag() );
}
}
}