CD_EBHelmholtzOp.H

Go to the documentation of this file.

/*
 * SPDX-FileCopyrightText: 2021-2026 SINTEF Energy Research
 *
 * SPDX-License-Identifier: GPL-3.0-or-later
 */
 
#ifndef CD_EBHELMHOLTZOP_H
#define CD_EBHELMHOLTZOP_H
 
// Std includes
#include <map>
 
// Chombo includes
#include <BaseEBBC.H>
#include <LevelTGA.H>
#include <BaseIVFAB.H>
#include <VCAggStencil.H>
 
// Our includes
#include <CD_Timer.H>
#include <CD_Location.H>
#include <CD_EBMultigridInterpolator.H>
#include <CD_EBCoarAve.H>
#include <CD_EBReflux.H>
#include <CD_EBMGRestrict.H>
#include <CD_EBMGProlong.H>
#include <CD_EBHelmholtzEBBC.H>
#include <CD_EBHelmholtzDomainBC.H>
#include <CD_NamespaceHeader.H>
 
class EBHelmholtzOp : public LevelTGAHelmOp<LevelData<EBCellFAB>, EBFluxFAB>
{
public:
  enum class Smoother
  {
    NoRelax,
    PointJacobi,
    GauSaiRedBlack,
    GauSaiMultiColor,
  };
 
  EBHelmholtzOp() = delete;
 
  EBHelmholtzOp(const EBHelmholtzOp& a_other) = delete;
 
  EBHelmholtzOp(const EBHelmholtzOp&& a_other) = delete;
 
  EBHelmholtzOp(const Location::Cell                             a_dataLocation,
                const EBLevelGrid&                               a_eblgFine,
                const EBLevelGrid&                               a_eblg,
                const EBLevelGrid&                               a_eblgCoFi,
                const EBLevelGrid&                               a_eblgCoar,
                const EBLevelGrid&                               a_eblgCoarMG,
                const RefCountedPtr<LevelData<BaseFab<bool>>>&   a_validCells,
                const RefCountedPtr<EBMultigridInterpolator>&    a_interpolator,
                const RefCountedPtr<EBReflux>&                   a_fluxReg,
                const RefCountedPtr<EBCoarAve>&                  a_coarAve,
                const RefCountedPtr<EBHelmholtzDomainBC>&        a_domainBC,
                const RefCountedPtr<EBHelmholtzEBBC>&            a_ebBC,
                const RealVect&                                  a_probLo,
                const Real&                                      a_dx,
                const int&                                       a_refToFine,
                const int&                                       a_refToCoar,
                const bool&                                      a_hasFine,
                const bool&                                      a_hasCoar,
                const bool&                                      a_hasMGObjects,
                const Real&                                      a_alpha,
                const Real&                                      a_beta,
                const RefCountedPtr<LevelData<EBCellFAB>>&       a_Acoef,
                const RefCountedPtr<LevelData<EBFluxFAB>>&       a_Bcoef,
                const RefCountedPtr<LevelData<BaseIVFAB<Real>>>& a_BcoefIrreg,
                const IntVect&                                   a_ghostPhi,
                const IntVect&                                   a_ghostRhs,
                const Smoother&                                  a_smoother,
                const Real&                                      a_relaxFactor);
 
  virtual ~EBHelmholtzOp();
 
  EBHelmholtzOp&
  operator=(const EBHelmholtzOp& a_oper) = delete;
 
  EBHelmholtzOp&
  operator=(const EBHelmholtzOp&& a_oper) = delete;
 
  void
  turnOffCFInterp();
 
  void
  turnOnCFInterp();
 
  void
  turnOffExchange();
 
  void
  turnOnExchange();
 
  void
  turnOffCoarsening();
 
  void
  turnOnCoarsening();
 
  void
  setAcoAndBco(const RefCountedPtr<LevelData<EBCellFAB>>&       a_Acoef,
               const RefCountedPtr<LevelData<EBFluxFAB>>&       a_Bcoef,
               const RefCountedPtr<LevelData<BaseIVFAB<Real>>>& a_BcoefIrreg);
 
  const RefCountedPtr<LevelData<EBCellFAB>>&
  getAcoef();
 
  const RefCountedPtr<LevelData<EBFluxFAB>>&
  getBcoef();
 
  const RefCountedPtr<LevelData<BaseIVFAB<Real>>>&
  getBcoefIrreg();
 
  void
  coarsenCell(LevelData<EBCellFAB>& a_phi, const LevelData<EBCellFAB>& a_phiFine);
 
  void
  coarsenFlux(LevelData<EBFluxFAB>& a_flux, const LevelData<EBFluxFAB>& a_fineFlux);
 
  void
  pointJacobiKernel(EBCellFAB&             a_Lcorr,
                    EBCellFAB&             a_corr,
                    const EBCellFAB&       a_resid,
                    const EBCellFAB&       a_Acoef,
                    const EBFluxFAB&       a_Bcoef,
                    const BaseIVFAB<Real>& a_BcoefIrreg,
                    const Box&             a_cellBox,
                    const DataIndex&       a_dit) const noexcept;
 
  void
  gauSaiRedBlackKernel(EBCellFAB&             a_Lcorr,
                       EBCellFAB&             a_corr,
                       const EBCellFAB&       a_resid,
                       const EBCellFAB&       a_Acoef,
                       const EBFluxFAB&       a_Bcoef,
                       const BaseIVFAB<Real>& a_BcoefIrreg,
                       const Box&             a_cellBox,
                       const DataIndex&       a_dit,
                       const int&             a_redBlack) const noexcept;
 
  void
  gauSaiMultiColorKernel(EBCellFAB&             a_Lcorr,
                         EBCellFAB&             a_corr,
                         const EBCellFAB&       a_resid,
                         const EBCellFAB&       a_Acoef,
                         const EBFluxFAB&       a_Bcoef,
                         const BaseIVFAB<Real>& a_BcoefIrreg,
                         const Box&             a_cellBox,
                         const DataIndex&       a_dit,
                         const IntVect&         a_color) const noexcept;
 
  void
  residual(LevelData<EBCellFAB>&       a_residual,
           const LevelData<EBCellFAB>& a_phi,
           const LevelData<EBCellFAB>& a_rhs,
           const bool                  a_homogeneousPhysBc) override;
 
  void
  preCond(LevelData<EBCellFAB>& a_corr, const LevelData<EBCellFAB>& a_residual) override final;
 
  void
  interpolateCF(LevelData<EBCellFAB>& a_phiFine, const LevelData<EBCellFAB>* a_phiCoar, const bool a_homogeneousCFBC);
 
  void
  homogeneousCFInterp(LevelData<EBCellFAB>& a_phi);
 
  void
  inhomogeneousCFInterp(LevelData<EBCellFAB>& a_phi, const LevelData<EBCellFAB>& a_phiCoar);
 
  void
  gauSaiMultiColor(LevelData<EBCellFAB>&       a_phi,
                   const LevelData<EBCellFAB>& a_Lphi,
                   const LevelData<EBCellFAB>& a_rhs,
                   const IntVect&              a_color) const;
 
  void
  applyOp(LevelData<EBCellFAB>&             a_Lphi,
          const LevelData<EBCellFAB>&       a_phi,
          const LevelData<EBCellFAB>* const a_phiCoar,
          const bool                        a_homogeneousPhysBC,
          const bool                        a_homogeneousCFBC);
 
  void
  applyOp(LevelData<EBCellFAB>& a_Lphi, const LevelData<EBCellFAB>& a_phi, bool a_homogeneousPhysBc) override final;
 
  void
  applyOp(EBCellFAB&             a_Lphi,
          EBCellFAB&             a_phi,
          const EBCellFAB&       a_Acoef,
          const EBFluxFAB&       a_Bcoef,
          const BaseIVFAB<Real>& a_BcoefIrreg,
          const Box&             a_cellBox,
          const DataIndex&       a_dit,
          const bool             a_homogeneousPhysBC) const noexcept;
 
  void
  applyOpRegular(EBCellFAB&             a_Lphi,
                 EBCellFAB&             a_phi,
                 const EBCellFAB&       a_Acoef,
                 const EBFluxFAB&       a_Bcoef,
                 const BaseIVFAB<Real>& a_BcoefIrreg,
                 const Box&             a_cellBox,
                 const DataIndex&       a_dit,
                 const bool             a_homogeneousPhysBC) const noexcept;
 
  void
  applyDomainFlux(EBCellFAB&       a_phi,
                  const EBFluxFAB& a_Bcoef,
                  const Box&       a_cellBox,
                  const DataIndex& a_dit,
                  const bool       a_homogeneousPhysBc) const noexcept;
 
  void
  fillDomainFlux(EBFluxFAB& a_flux, const EBCellFAB& a_phi, const Box& a_cellBox, const DataIndex& a_dit);
 
  void
  applyOpIrregular(EBCellFAB&             a_Lphi,
                   const EBCellFAB&       a_phi,
                   const EBCellFAB&       a_Acoef,
                   const EBFluxFAB&       a_Bcoef,
                   const BaseIVFAB<Real>& a_BcoefIrreg,
                   const BaseIVFAB<Real>& a_alphaDiagWeight,
                   const Box&             a_cellBox,
                   const DataIndex&       a_dit,
                   const bool             a_homogeneousPhysBC) const noexcept;
 
  void
  create(LevelData<EBCellFAB>& a_lhs, const LevelData<EBCellFAB>& a_rhs) override final;
 
  void
  assign(LevelData<EBCellFAB>& a_lhs, const LevelData<EBCellFAB>& a_rhs) override final;
 
  void
  assignCopier(LevelData<EBCellFAB>& a_lhs, const LevelData<EBCellFAB>& a_rhs, const Copier& a_copier) override final;
 
  void
  assignLocal(LevelData<EBCellFAB>& a_lhs, const LevelData<EBCellFAB>& a_rhs) override final;
 
  void
  buildCopier(Copier& a_copier, const LevelData<EBCellFAB>& a_lhs, const LevelData<EBCellFAB>& a_rhs) override;
 
  Real
  dotProduct(const LevelData<EBCellFAB>& a_lhs, const LevelData<EBCellFAB>& a_rhs) override final;
 
  void
  incr(LevelData<EBCellFAB>& a_lhs, const LevelData<EBCellFAB>& a_rhs, const Real a_scale) override final;
 
  void
  axby(LevelData<EBCellFAB>&       a_lhs,
       const LevelData<EBCellFAB>& a_x,
       const LevelData<EBCellFAB>& a_y,
       const Real                  a_a,
       const Real                  a_b) override final;
 
  void
  scale(LevelData<EBCellFAB>& a_lhs, const Real& a_scale) override final;
 
  static void
  scaleLocal(LevelData<EBCellFAB>& a_lhs, const LevelData<EBCellFAB>& a_rhs) noexcept;
 
  Real
  norm(const LevelData<EBCellFAB>& a_rhs, const int a_order) override final;
 
  void
  setToZero(LevelData<EBCellFAB>& a_lhs) override final;
 
  void
  createCoarser(LevelData<EBCellFAB>& a_coarse, const LevelData<EBCellFAB>& a_fine, bool a_ghosted) override final;
 
  void
  relax(LevelData<EBCellFAB>& a_correction, const LevelData<EBCellFAB>& a_residual, int a_iterations) override final;
 
  void
  restrictResidual(LevelData<EBCellFAB>&       a_resCoar,
                   LevelData<EBCellFAB>&       a_phi,
                   const LevelData<EBCellFAB>& a_rhs) override final;
 
  void
  prolongIncrement(LevelData<EBCellFAB>& a_phi, const LevelData<EBCellFAB>& a_correctCoarse) override final;
 
  int
  refToCoarser() override final;
 
  void
  AMROperator(LevelData<EBCellFAB>&             a_Lphi,
              const LevelData<EBCellFAB>&       a_phiFine,
              const LevelData<EBCellFAB>&       a_phi,
              const LevelData<EBCellFAB>&       a_phiCoar,
              const bool                        a_homogeneousPhysBC,
              AMRLevelOp<LevelData<EBCellFAB>>* a_finerOp) override final;
 
  void
  refluxFreeAMROperator(LevelData<EBCellFAB>&             a_Lphi,
                        const LevelData<EBCellFAB>&       a_phiFine,
                        const LevelData<EBCellFAB>&       a_phi,
                        const LevelData<EBCellFAB>&       a_phiCoar,
                        const bool                        a_homogeneousPhysBC,
                        AMRLevelOp<LevelData<EBCellFAB>>* a_finerOp);
 
  void
  AMRResidual(LevelData<EBCellFAB>&             a_residual,
              const LevelData<EBCellFAB>&       a_phiFine,
              const LevelData<EBCellFAB>&       a_phi,
              const LevelData<EBCellFAB>&       a_phiCoar,
              const LevelData<EBCellFAB>&       a_rhs,
              bool                              a_homogeneousPhysBC,
              AMRLevelOp<LevelData<EBCellFAB>>* a_finerOp) override final;
 
  void
  AMRResidualNF(LevelData<EBCellFAB>&       a_residual,
                const LevelData<EBCellFAB>& a_phi,
                const LevelData<EBCellFAB>& a_phiCoar,
                const LevelData<EBCellFAB>& a_rhs,
                bool                        a_homogeneousPhysBC) override final;
 
  void
  AMRResidualNC(LevelData<EBCellFAB>&             a_residual,
                const LevelData<EBCellFAB>&       a_phiFine,
                const LevelData<EBCellFAB>&       a_phi,
                const LevelData<EBCellFAB>&       a_rhs,
                bool                              a_homogeneousPhysBC,
                AMRLevelOp<LevelData<EBCellFAB>>* a_finerOp) override final;
 
  void
  AMROperatorNF(LevelData<EBCellFAB>&       a_Lphi,
                const LevelData<EBCellFAB>& a_phi,
                const LevelData<EBCellFAB>& a_phiCoar,
                bool                        a_homogeneousPhysBC) override final;
 
  void
  AMROperatorNC(LevelData<EBCellFAB>&             a_Lphi,
                const LevelData<EBCellFAB>&       a_phiFine,
                const LevelData<EBCellFAB>&       a_phi,
                bool                              a_homogeneousPhysBC,
                AMRLevelOp<LevelData<EBCellFAB>>* a_finerOp) override final;
 
  void
  AMRRestrict(LevelData<EBCellFAB>&       a_residualCoarse,
              const LevelData<EBCellFAB>& a_residual,
              const LevelData<EBCellFAB>& a_correction,
              const LevelData<EBCellFAB>& a_coarseCorrection,
              bool                        a_skip_res) override final;
 
  void
  AMRProlong(LevelData<EBCellFAB>& a_correction, const LevelData<EBCellFAB>& a_coarseCorrection) override final;
 
  void
  AMRUpdateResidual(LevelData<EBCellFAB>&       a_residual,
                    const LevelData<EBCellFAB>& a_correction,
                    const LevelData<EBCellFAB>& a_coarseCorrection) override final;
 
  void
  setAlphaAndBeta(const Real& a_alpha, const Real& a_beta) override final;
 
  void
  createCoarsened(LevelData<EBCellFAB>& a_lhs, const LevelData<EBCellFAB>& a_rhs, const int& a_refRat) override final;
 
  void
  diagonalScale(LevelData<EBCellFAB>& a_rhs, bool a_kappaWeighted) override final;
 
  void
  divideByIdentityCoef(LevelData<EBCellFAB>& a_rhs) override final;
 
  void
  applyOpNoBoundary(LevelData<EBCellFAB>& a_Lphi, const LevelData<EBCellFAB>& a_phi) override final;
 
  void
  fillGrad(const LevelData<EBCellFAB>& a_phi) override final;
 
  void
  getFlux(EBFluxFAB&                  a_flux,
          const LevelData<EBCellFAB>& a_data,
          const Box&                  a_grid,
          const DataIndex&            a_dit,
          Real                        a_scale) override final;
 
  void
  allocateFlux() const noexcept;
 
  void
  deallocateFlux() const noexcept;
 
  LevelData<EBFluxFAB>&
  getFlux() const;
 
  const LevelData<EBCellFAB>&
  getRelaxationCoeff() const noexcept;
 
  LayoutData<VoFIterator>&
  getVofIterIrreg() const noexcept;
 
protected:
  static constexpr int m_nComp = 1;
 
  static constexpr int m_comp = 0;
 
  Timer m_timer;
 
  Interval m_interval;
 
  Smoother m_smoother;
 
  Location::Cell m_dataLocation;
 
  bool m_refluxFree;
 
  bool m_profile;
 
  bool m_hasMGObjects;
 
  bool m_hasFine;
 
  bool m_hasCoar;
 
  int m_refToCoar;
 
  int m_refToFine;
 
  int m_numSmoothPreCond;
 
  bool m_doInterpCF;
 
  bool m_doExchange;
 
  bool m_doCoarsen;
 
  IntVect m_ghostPhi;
 
  IntVect m_ghostRhs;
 
  Real m_alpha;
 
  Real m_beta;
 
  Real m_dx;
 
  Real m_relaxFactor;
 
  RealVect m_vecDx;
 
  RealVect m_probLo;
 
  Copier m_exchangeCopier;
 
  Copier m_exchangeCopierFine;
 
  EBLevelGrid m_eblgFine;
 
  EBLevelGrid m_eblg;
 
  EBLevelGrid m_eblgCoFi;
 
  EBLevelGrid m_eblgCoar;
 
  EBLevelGrid m_eblgCoarMG;
 
  EBMGRestrict m_restrictOp;
 
  EBMGRestrict m_restrictOpMG;
 
  EBMGProlong m_prolongOp;
 
  EBMGProlong m_prolongOpMG;
 
  RefCountedPtr<LevelData<BaseFab<bool>>> m_validCells;
 
  RefCountedPtr<EBHelmholtzDomainBC> m_domainBc;
 
  RefCountedPtr<EBHelmholtzEBBC> m_ebBc;
 
  RefCountedPtr<EBMultigridInterpolator> m_interpolator;
 
  RefCountedPtr<EBReflux> m_fluxReg;
 
  RefCountedPtr<EBCoarAve> m_coarAve;
 
  RefCountedPtr<LevelData<EBCellFAB>> m_Acoef;
 
  RefCountedPtr<LevelData<EBFluxFAB>> m_Bcoef;
 
  RefCountedPtr<LevelData<BaseIVFAB<Real>>> m_BcoefIrreg;
 
  LevelData<EBCellFAB> m_relCoef;
 
  mutable LevelData<EBFluxFAB>* m_flux;
 
  LayoutData<BaseIFFAB<Real>> m_interpolant[SpaceDim];
 
  LayoutData<BaseIFFAB<FaceStencil>> m_interpStencil[SpaceDim];
 
  LayoutData<BaseIFFAB<VoFStencil>> m_centroidFluxStencil[SpaceDim];
 
  LayoutData<BaseIVFAB<VoFStencil>> m_relaxStencils;
 
  LayoutData<RefCountedPtr<VCAggStencil>> m_aggRelaxStencil;
 
  LayoutData<BaseIVFAB<Real>> m_alphaDiagWeight;
 
  LayoutData<BaseIVFAB<Real>> m_betaDiagWeight;
 
  mutable LayoutData<VoFIterator> m_vofIterIrreg;
 
  mutable LayoutData<VoFIterator> m_vofIterMulti;
 
  mutable LayoutData<VoFIterator> m_vofIterStenc;
 
  mutable LayoutData<VoFIterator> m_vofIterDomLo[SpaceDim];
 
  mutable LayoutData<VoFIterator> m_vofIterDomHi[SpaceDim];
 
  std::map<std::pair<int, Side::LoHiSide>, Box> m_sideBox;
 
  struct DomainBndryBoxes
  {
    bool touchesDomain = false;
 
    int hasLo[SpaceDim] = {D_DECL(0, 0, 0)};
 
    int hasHi[SpaceDim] = {D_DECL(0, 0, 0)};
 
    Box loBox[SpaceDim];
 
    Box hiBox[SpaceDim];
  };
 
  LayoutData<DomainBndryBoxes> m_domainBndryBoxes;
 
  Vector<IntVect> m_colors;
 
  void
  relaxPointJacobi(LevelData<EBCellFAB>& a_correction, const LevelData<EBCellFAB>& a_residual, const int a_iterations);
 
  void
  relaxGSRedBlack(LevelData<EBCellFAB>& a_correction, const LevelData<EBCellFAB>& a_residual, const int a_iterations);
 
  void
  relaxGSMultiColor(LevelData<EBCellFAB>& a_correction, const LevelData<EBCellFAB>& a_residual, const int a_iterations);
 
  void
  computeDiagWeight();
 
  void
  computeRelaxationCoefficient();
 
  void
  makeAggStencil();
 
  void
  defineStencils();
 
  VoFStencil
  getFaceCenterFluxStencil(const FaceIndex& a_face, const DataIndex& a_dit) const;
 
  VoFStencil
  getFaceCentroidFluxStencil(const FaceIndex& a_face, const DataIndex& a_dit) const;
 
  void
  computeFlux(EBFaceFAB&       a_fluxCentroid,
              const EBCellFAB& a_phi,
              const Box&       a_cellBox,
              const DataIndex& a_dit,
              const int        a_dir);
 
  void
  computeFaceCenteredFlux(EBFaceFAB&       a_fluxCenter,
                          const EBCellFAB& a_phi,
                          const Box&       a_cellBox,
                          const DataIndex& a_dit,
                          const int        a_dir);
 
  void
  computeFaceCentroidFlux(EBFaceFAB&       a_flux,
                          const EBCellFAB& a_phi,
                          const Box&       a_cellBox,
                          const DataIndex& a_dit,
                          const int        a_dir);
 
  void
  computeFlux(const LevelData<EBCellFAB>& a_phi);
 
  void
  reflux(LevelData<EBCellFAB>&             a_Lphi,
         const LevelData<EBCellFAB>&       a_phiFine,
         const LevelData<EBCellFAB>&       a_phi,
         AMRLevelOp<LevelData<EBCellFAB>>& a_finerOp);
};
 
#include <CD_NamespaceFooter.H>
 
#endif