Nous avons procédé ce jeudi matin 08 avril 2021 à une MAJ de sécurité urgente. Nous sommes passé de la version 13.9.3 à la version 13.9.5 les releases notes correspondantes sont ici:
https://about.gitlab.com/releases/2021/03/17/security-release-gitlab-13-9-4-released/
https://about.gitlab.com/releases/2021/03/31/security-release-gitlab-13-10-1-released/

FUnifTensorialKernel.hpp 10.6 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15
// ===================================================================================
// Copyright ScalFmm 2011 INRIA, Olivier Coulaud, Bérenger Bramas, Matthias Messner
// olivier.coulaud@inria.fr, berenger.bramas@inria.fr
// This software is a computer program whose purpose is to compute the FMM.
//
// This software is governed by the CeCILL-C and LGPL licenses and
// abiding by the rules of distribution of free software.
//
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
// GNU General Public and CeCILL-C Licenses for more details.
// "http://www.cecill.info".
// "http://www.gnu.org/licenses".
// ===================================================================================
16 17 18
// Keep in private GIT
// @SCALFMM_PRIVATE

19 20 21 22
#ifndef FUNIFTENSORIALKERNEL_HPP
#define FUNIFTENSORIALKERNEL_HPP

#include "../../Utils/FGlobal.hpp"
23

24 25 26
#include "../../Utils/FSmartPointer.hpp"

#include "./FAbstractUnifKernel.hpp"
27 28
#include "./FUnifM2LHandler.hpp"
#include "./FUnifTensorialM2LHandler.hpp" //PB: temporary version
29 30 31 32 33 34 35 36 37 38 39 40 41 42 43

class FTreeCoordinate;

/**
 * @author Pierre Blanchard (pierre.blanchard@inria.fr)
 * @class FUnifTensorialKernel
 * @brief
 * Please read the license
 *
 * This kernels implement the Lagrange interpolation based FMM operators. It
 * implements all interfaces (P2P,P2M,M2M,M2L,L2L,L2P) which are required by
 * the FFmmAlgorithm and FFmmAlgorithmThread.
 *
 * PB: 3 IMPORTANT remarks !!!
 *
44 45
 * 1) Handling tensorial kernels (DIM,NRHS,NLHS) and having multiple rhs 
 * (NVALS) are considered 2 distinct features and are currently combined.
46
 *
47 48 49 50 51 52
 * 2) When it comes to applying M2L it is NOT much faster to loop over 
 * NRHSxNLHS inside applyM2L (at least for the Lagrange case).
 * 2-bis) During precomputation the tensorial matrix kernels are evaluated 
 * blockwise, but this is not always possible. 
 * In fact, in the ChebyshevSym variant the matrix kernel needs to be 
 * evaluated compo-by-compo since we currently use a scalar ACA.
53
 *
54 55 56
 * 3) We currently use multiple 1D FFT instead of multidim FFT since embedding
 * is circulant. Multidim FFT could be used if embedding were block circulant.
 * TODO investigate possibility of block circulant embedding
57 58 59 60 61 62 63 64 65 66
 *
 * @tparam CellClass Type of cell
 * @tparam ContainerClass Type of container to store particles
 * @tparam MatrixKernelClass Type of matrix kernel function
 * @tparam ORDER Lagrange interpolation order
 */
template < class CellClass,	class ContainerClass,	class MatrixKernelClass, int ORDER, int NVALS = 1>
class FUnifTensorialKernel
  : public FAbstractUnifKernel< CellClass, ContainerClass, MatrixKernelClass, ORDER, NVALS>
{
67
  enum {nRhs = MatrixKernelClass::NRHS,
68 69 70
        nLhs = MatrixKernelClass::NLHS,
        nPot = MatrixKernelClass::NPOT,
        nPV = MatrixKernelClass::NPV};
71 72 73

protected://PB: for OptiDis

74
  // private types
75
  typedef FUnifTensorialM2LHandler<ORDER,MatrixKernelClass,MatrixKernelClass::Type> M2LHandlerClass;
76 77 78 79 80

  // using from
  typedef FAbstractUnifKernel< CellClass, ContainerClass, MatrixKernelClass, ORDER, NVALS>
  AbstractBaseClass;

81 82 83
  /// Needed for P2P and M2L operators
  const MatrixKernelClass *const MatrixKernel;

84
  /// Needed for M2L operator
85
  const M2LHandlerClass M2LHandler;
86 87 88 89 90 91 92 93

public:
  /**
   * The constructor initializes all constant attributes and it reads the
   * precomputed and compressed M2L operators from a binary file (an
   * runtime_error is thrown if the required file is not valid).
   */
  FUnifTensorialKernel(const int inTreeHeight,
94
                       const FReal inBoxWidth,
95
                       const FPoint& inBoxCenter,
96 97 98 99 100
                       const MatrixKernelClass *const inMatrixKernel,
                       const FReal inBoxWidthExtension)
    : FAbstractUnifKernel< CellClass, ContainerClass, MatrixKernelClass, ORDER, NVALS>(inTreeHeight,inBoxWidth,inBoxCenter,inBoxWidthExtension),
      MatrixKernel(inMatrixKernel),
      M2LHandler(MatrixKernel,
101
                 inTreeHeight,
102 103
                 inBoxWidth,
                 inBoxWidthExtension) 
104
  { }
105 106 107 108 109 110


  void P2M(CellClass* const LeafCell,
           const ContainerClass* const SourceParticles)
  {
    const FPoint LeafCellCenter(AbstractBaseClass::getLeafCellCenter(LeafCell->getCoordinate())); 
111 112
    const FReal ExtendedLeafCellWidth(AbstractBaseClass::BoxWidthLeaf 
                                      + AbstractBaseClass::BoxWidthExtension);
113

114
    for(int idxV = 0 ; idxV < NVALS ; ++idxV){
115 116

      // 1) apply Sy
117
      AbstractBaseClass::Interpolator->applyP2M(LeafCellCenter, ExtendedLeafCellWidth,
118 119
                                                LeafCell->getMultipole(idxV*nRhs), SourceParticles);

120 121 122 123 124
      for(int idxRhs = 0 ; idxRhs < nRhs ; ++idxRhs){
        // update multipole index
        int idxMul = idxV*nRhs + idxRhs;

        // 2) apply Discrete Fourier Transform
125 126
        M2LHandler.applyZeroPaddingAndDFT(LeafCell->getMultipole(idxMul), 
                                          LeafCell->getTransformedMultipole(idxMul));
127 128 129

      }
    }// NVALS
130 131 132 133 134
  }


  void M2M(CellClass* const FRestrict ParentCell,
           const CellClass*const FRestrict *const FRestrict ChildCells,
135
           const int TreeLevel)
136
  {
137 138 139 140 141 142 143 144 145
    for(int idxV = 0 ; idxV < NVALS ; ++idxV){
      for(int idxRhs = 0 ; idxRhs < nRhs ; ++idxRhs){
        // update multipole index
        int idxMul = idxV*nRhs + idxRhs;

        // 1) apply Sy
        FBlas::scal(AbstractBaseClass::nnodes, FReal(0.), ParentCell->getMultipole(idxMul));
        for (unsigned int ChildIndex=0; ChildIndex < 8; ++ChildIndex){
          if (ChildCells[ChildIndex]){
146 147 148 149
            AbstractBaseClass::Interpolator->applyM2M(ChildIndex, 
                                                      ChildCells[ChildIndex]->getMultipole(idxMul),
                                                      ParentCell->getMultipole(idxMul), 
                                                      TreeLevel/*Cell width extension specific*/);
150
          }
151
        }
152
        // 2) Apply Discete Fourier Transform
153 154
        M2LHandler.applyZeroPaddingAndDFT(ParentCell->getMultipole(idxMul), 
                                          ParentCell->getTransformedMultipole(idxMul));
155
      }
156
    }// NVALS
157 158 159 160 161 162 163 164
  }


  void M2L(CellClass* const FRestrict TargetCell,
           const CellClass* SourceCells[343],
           const int /*NumSourceCells*/,
           const int TreeLevel)
  {
165
    const FReal CellWidth(AbstractBaseClass::BoxWidth / FReal(FMath::pow(2, TreeLevel)));
166
    const FReal ExtendedCellWidth(CellWidth + AbstractBaseClass::BoxWidthExtension);
167
    const FReal scale(MatrixKernel->getScaleFactor(ExtendedCellWidth));
168 169 170

    for(int idxV = 0 ; idxV < NVALS ; ++idxV){
      for (int idxLhs=0; idxLhs < nLhs; ++idxLhs){
171

172 173 174 175
          // update local index
          const int idxLoc = idxV*nLhs + idxLhs;

          // load transformed local expansion
176
          FComplex *const TransformedLocalExpansion = TargetCell->getTransformedLocal(idxLoc);
177

178 179
          // update idxRhs
          const int idxRhs = idxLhs % nPV; 
180

181
          // update multipole index
182 183
          const int idxMul = idxV*nRhs + idxRhs;

184
          // get index in matrix kernel
185
          const unsigned int d = MatrixKernel->getPosition(idxLhs);
186

187 188
          for (int idx=0; idx<343; ++idx){
            if (SourceCells[idx]){
189

190
              M2LHandler.applyFC(idx, TreeLevel, scale, d,
191
                                  SourceCells[idx]->getTransformedMultipole(idxMul),
192
                                  TransformedLocalExpansion);
193

194
            }
195
          }
196
      }// NLHS=NPOT*NPV
197
    }// NVALS
198 199 200 201 202
  }


  void L2L(const CellClass* const FRestrict ParentCell,
           CellClass* FRestrict *const FRestrict ChildCells,
203
           const int TreeLevel)
204
  {
205 206 207 208
    for(int idxV = 0 ; idxV < NVALS ; ++idxV){
      for(int idxLhs = 0 ; idxLhs < nLhs ; ++idxLhs){
        int idxLoc = idxV*nLhs + idxLhs;
        // 1) Apply Inverse Discete Fourier Transform
209 210
        M2LHandler.unapplyZeroPaddingAndDFT(ParentCell->getTransformedLocal(idxLoc),
                                            const_cast<CellClass*>(ParentCell)->getLocal(idxLoc));
211 212 213
        // 2) apply Sx
        for (unsigned int ChildIndex=0; ChildIndex < 8; ++ChildIndex){
          if (ChildCells[ChildIndex]){
214 215 216 217
            AbstractBaseClass::Interpolator->applyL2L(ChildIndex, 
                                                      ParentCell->getLocal(idxLoc), 
                                                      ChildCells[ChildIndex]->getLocal(idxLoc),
                                                      TreeLevel/*Cell width extension specific*/);
218
          }
219 220
        }
      }
221
    }// NVALS
222 223 224 225 226 227
  }

  void L2P(const CellClass* const LeafCell,
           ContainerClass* const TargetParticles)
  {
    const FPoint LeafCellCenter(AbstractBaseClass::getLeafCellCenter(LeafCell->getCoordinate()));
228 229
    const FReal ExtendedLeafCellWidth(AbstractBaseClass::BoxWidthLeaf 
                                      + AbstractBaseClass::BoxWidthExtension);
230

231 232 233 234
    for(int idxV = 0 ; idxV < NVALS ; ++idxV){
      for(int idxLhs = 0 ; idxLhs < nLhs ; ++idxLhs){
        int idxLoc = idxV*nLhs + idxLhs;
        // 1)  Apply Inverse Discete Fourier Transform
235 236
        M2LHandler.unapplyZeroPaddingAndDFT(LeafCell->getTransformedLocal(idxLoc), 
                                            const_cast<CellClass*>(LeafCell)->getLocal(idxLoc));
237

238 239 240
      }

      // 2.a) apply Sx
241
      AbstractBaseClass::Interpolator->applyL2P(LeafCellCenter, ExtendedLeafCellWidth,
242
                                                LeafCell->getLocal(idxV*nLhs), TargetParticles);
243

244
      // 2.b) apply Px (grad Sx)
245
      AbstractBaseClass::Interpolator->applyL2PGradient(LeafCellCenter, ExtendedLeafCellWidth,
246
                                                        LeafCell->getLocal(idxV*nLhs), TargetParticles);
247

248
    }// NVALS
249 250 251 252 253 254 255 256
  }

  void P2P(const FTreeCoordinate& /* LeafCellCoordinate */, // needed for periodic boundary conditions
           ContainerClass* const FRestrict TargetParticles,
           const ContainerClass* const FRestrict /*SourceParticles*/,
           ContainerClass* const NeighborSourceParticles[27],
           const int /* size */)
  {
257
    DirectInteractionComputer<MatrixKernelClass::NCMP, NVALS>::P2P(TargetParticles,NeighborSourceParticles,MatrixKernel);
258 259 260 261 262 263
  }


  void P2PRemote(const FTreeCoordinate& /*inPosition*/,
                 ContainerClass* const FRestrict inTargets, const ContainerClass* const FRestrict /*inSources*/,
                 ContainerClass* const inNeighbors[27], const int /*inSize*/){
264
    DirectInteractionComputer<MatrixKernelClass::NCMP, NVALS>::P2PRemote(inTargets,inNeighbors,27,MatrixKernel);
265 266 267 268 269
  }

};


270
#endif //FUNIFTENSORIALKERNEL_HPP
271 272

// [--END--]