PoissonSolver.h

Go to the documentation of this file.

// Copyright 2019-2020 CERN and copyright holders of ALICE O2.
// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders.
// All rights not expressly granted are reserved.
//
// This software is distributed under the terms of the GNU General Public
// License v3 (GPL Version 3), copied verbatim in the file "COPYING".
//
// In applying this license CERN does not waive the privileges and immunities
// granted to it by virtue of its status as an Intergovernmental Organization
// or submit itself to any jurisdiction.
 
 
#ifndef ALICEO2_TPC_POISSONSOLVER_H_
#define ALICEO2_TPC_POISSONSOLVER_H_
 
#include "TPCSpaceCharge/RegularGrid3D.h"
#include "CommonConstants/MathConstants.h"
#include "TPCSpaceCharge/SpaceChargeParameter.h"
#include <vector>
 
namespace o2
{
namespace tpc
{
 
template <typename DataT>
class Vector3D;
 
template <typename DataT>
class DataContainer3D;
 
 
template <typename DataT = double>
class PoissonSolver
{
 public:
  using RegularGrid = RegularGrid3D<DataT>;
  using DataContainer = DataContainer3D<DataT>;
  using Vector = Vector3D<DataT>;
 
  PoissonSolver(const RegularGrid& gridProperties) : mGrid3D{gridProperties} {};
 
  void poissonSolver3D(DataContainer& matricesV, const DataContainer& matricesCharge, const int symmetry);
 
  void poissonSolver2D(DataContainer& matricesV, const DataContainer& matricesCharge);
 
  DataT getSpacingZ() const { return mGrid3D.getSpacingZ(); }
  DataT getSpacingR() const { return mGrid3D.getSpacingR(); }
  DataT getSpacingPhi() const { return mGrid3D.getSpacingPhi(); }
 
  static void setConvergenceError(const DataT error) { sConvergenceError = error; }
 
  static DataT getConvergenceError() { return sConvergenceError; }
 
  static int getNThreads() { return sNThreads; }
 
  static void setNThreads(int nThreads) { sNThreads = nThreads; }
 
 private:
  const RegularGrid& mGrid3D{};                                      
  const ParamSpaceCharge mParamGrid{mGrid3D.getParamSC()};           
  inline static DataT sConvergenceError{1e-6};                       
  static constexpr DataT INVTWOPI = 1. / o2::constants::math::TwoPI; 
  inline static int sNThreads{4};                                    
 
  static DataT getGridSizePhiInv();
 
  DataT getConvergenceError(const Vector& matricesCurrentV, Vector& prevArrayV) const;
 
  //
  void poissonMultiGrid3D2D(DataContainer& matricesV, const DataContainer& matricesCharge, const int symmetry);
 
  void poissonMultiGrid3D(DataContainer& matricesV, const DataContainer& matricesCharge, const int symmetry);
 
  void poissonMultiGrid2D(DataContainer& matricesV, const DataContainer& matricesCharge, const int iPhi = 0);
 
  void restrict2D(Vector& matricesCurrentCharge, const Vector& residue, const int tnRRow, const int tnZColumn, const int iphi) const;
 
  void restrict3D(Vector& matricesCurrentCharge, const Vector& residue, const int tnRRow, const int tnZColumn, const int newPhiSlice, const int oldPhiSlice) const;
 
  void restrictBoundary3D(Vector& matricesCurrentCharge, const Vector& residue, const int tnRRow, const int tnZColumn, const int newPhiSlice, const int oldPhiSlice) const;
 
  void relax3D(Vector& matricesCurrentV, const Vector& matricesCurrentCharge, const int tnRRow, const int tnZColumn, const int iPhi, const int symmetry, const DataT h2, const DataT tempRatioZ,
               const std::vector<DataT>& coefficient1, const std::vector<DataT>& coefficient2, const std::vector<DataT>& coefficient3, const std::vector<DataT>& coefficient4) const;
 
  void relax2D(Vector& matricesCurrentV, const Vector& matricesCurrentCharge, const int tnRRow, const int tnZColumn, const DataT h2, const DataT tempFourth, const DataT tempRatio,
               std::vector<DataT>& coefficient1, std::vector<DataT>& coefficient2);
 
  void interp2D(Vector& matricesCurrentV, const Vector& matricesCurrentVC, const int tnRRow, const int tnZColumn, const int iphi) const;
 
  void interp3D(Vector& matricesCurrentV, const Vector& matricesCurrentVC, const int tnRRow, const int tnZColumn, const int newPhiSlice, const int oldPhiSlice) const;
 
  void addInterp3D(Vector& matricesCurrentV, const Vector& matricesCurrentVC, const int tnRRow, const int tnZColumn, const int newPhiSlice, const int oldPhiSlice) const;
 
  void addInterp2D(Vector& matricesCurrentV, const Vector& matricesCurrentVC, const int tnRRow, const int tnZColumn, const int tnPhi) const;
 
  void vCycle3D2D(const int symmetry, const int gridFrom, const int gridTo, const int nPre, const int nPost, const DataT ratioZ, const DataT ratioPhi, std::vector<Vector>& tvArrayV,
                  std::vector<Vector>& tvCharge, std::vector<Vector>& tvResidue, std::vector<DataT>& coefficient1, std::vector<DataT>& coefficient2, std::vector<DataT>& coefficient3,
                  std::vector<DataT>& coefficient4, std::vector<DataT>& inverseCoefficient4) const;
 
  void vCycle3D(const int symmetry, const int gridFrom, const int gridTo, const int nPre, const int nPost, const DataT ratioZ, std::vector<Vector>& tvArrayV, std::vector<Vector>& tvCharge,
                std::vector<Vector>& tvResidue, std::vector<DataT>& coefficient1, std::vector<DataT>& coefficient2, std::vector<DataT>& coefficient3,
                std::vector<DataT>& coefficient4, std::vector<DataT>& inverseCoefficient4) const;
 
  void vCycle2D(const int gridFrom, const int gridTo, const int nPre, const int nPost, const DataT gridSizeR, const DataT ratio, std::vector<Vector>& tvArrayV,
                std::vector<Vector>& tvCharge, std::vector<Vector>& tvResidue);
 
  void wCycle2D(const int gridFrom, const int gridTo, const int gamma, const int nPre, const int nPost, const DataT gridSizeR, const DataT ratio,
                std::vector<Vector>& tvArrayV, std::vector<Vector>& tvCharge, std::vector<Vector>& tvResidue);
 
  void residue3D(Vector& residue, const Vector& matricesCurrentV, const Vector& matricesCurrentCharge, const int tnRRow, const int tnZColumn, const int tnPhi, const int symmetry, const DataT ih2, const DataT tempRatioZ,
                 const std::vector<DataT>& coefficient1, const std::vector<DataT>& coefficient2, const std::vector<DataT>& coefficient3, const std::vector<DataT>& inverseCoefficient4) const;
 
  void residue2D(Vector& residue, const Vector& matricesCurrentV, const Vector& matricesCurrentCharge, const int tnRRow, const int tnZColumn, const DataT ih2, const DataT inverseTempFourth,
                 const DataT tempRatio, std::vector<DataT>& coefficient1, std::vector<DataT>& coefficient2);
 
  void restrictBoundary2D(Vector& matricesCurrentCharge, const Vector& residue, const int tnRRow, const int tnZColumn, const int tnPhi) const;
 
  // calculate coefficients
  void calcCoefficients(unsigned int from, unsigned int to, const DataT h, const DataT tempRatioZ, const DataT tempRatioPhi, std::vector<DataT>& coefficient1,
                        std::vector<DataT>& coefficient2, std::vector<DataT>& coefficient3, std::vector<DataT>& coefficient4) const;
 
  // calculate coefficients for 2D poisson solver
  void calcCoefficients2D(unsigned int from, unsigned int to, const DataT h, std::vector<DataT>& coefficient1, std::vector<DataT>& coefficient2) const;
 
  bool isPowerOfTwo(const int i) const
  {
    return ((i > 0) && !(i & (i - 1)));
  };
};
 
} // namespace tpc
} // namespace o2
 
#endif