Skip to content
Snippets Groups Projects
Select Git revision
  • f0e5d9571000f745ead2c476eeed75a91f957c78
  • master default protected
  • rgpu
  • oras_vs_osm
  • refactor_coupled
  • lumi-stable
  • fix-compile-without-mpi
  • clean_multirhs
  • oras_comp
  • hpddm_integration
  • blockProduct
  • multiSrcs
  • splitPrePro
  • reuseGCR
  • helmholtz_2d_ddm
  • fix-template-instanciantion-clang-macos
  • customSchwarz
  • hp-convergence-test
  • fix_krylov
  • solverCorrection
  • boris-martin-master-patch-52103
  • gmshddm_1_0_0
22 results

Formulation.h

Blame
  • Formulation.h 7.16 KiB
    // GmshDDM - Copyright (C) 2019-2022, A. Royer, C. Geuzaine, Université de Liège
    //
    // See the LICENSE.txt file for license information. Please report all
    // issues on https://gitlab.onelab.info/gmsh/ddm/issues
    
    #ifndef H_GMSHDDM_FORMULATION
    #define H_GMSHDDM_FORMULATION
    
    #include "InterfaceCompoundField.h"
    #include "InterfaceField.h"
    #include "SubdomainCompoundField.h"
    #include "SubdomainField.h"
    
    #include <gmshfem/Formulation.h>
    #include <unordered_map>
    #include <vector>
    
    typedef struct _p_Mat *Mat;
    typedef struct _p_Vec *Vec;
    typedef struct _p_KSP *KSP;
    
    namespace gmshddm
    {
    
    
      namespace problem
      {
    
        template< class T_Scalar >
        int MatVectProductA(Mat A, Vec X, Vec Y);
        template< class T_Scalar >
        int MatVectProductI_A(Mat A, Vec X, Vec Y);
        template< class T_Scalar >
        int MatVectProductImpl(Mat A, Vec X, Vec Y, bool IA);
        template< class T_Scalar >
        int MatMatProductImpl(Mat A, Mat X, Mat Y, bool IA);
        template< class T_Scalar >
        int MatMatProductA(Mat A, Mat X, Mat Y);
        template< class T_Scalar >
        int MatMatProductI_A(Mat A, Mat X, Mat Y);
    
        template< class T_Scalar >
        class AbstractIterativeSolver;
    
        template< class T_Scalar >
        class Formulation
        {
         private:
          const std::string _name;
          std::vector< gmshfem::problem::Formulation< T_Scalar > * > _volume;
          std::vector< std::unordered_map< unsigned int, gmshfem::problem::Formulation< T_Scalar > * > > _surface;
          std::unordered_map< field::InterfaceFieldInterface< T_Scalar > *, field::InterfaceFieldInterface< T_Scalar > * > _interfaceFields;
          std::unordered_map< unsigned int, std::vector< unsigned long long > > _interfaceMapping;
          gmshfem::algebra::Vector< T_Scalar > _rhs;
          std::vector<gmshfem::algebra::Vector< T_Scalar >> _multi_rhs; // For many shots
    
          unsigned int _numberOfIterations;
          bool _physicalCommutator;
          bool _artificialCommutator;
          bool _IA;
          std::set< unsigned int > _physicalSourceTerms;
          std::set< unsigned int > _artificialSourceTerms;
          std::vector< std::set< unsigned int > > _numberedPhysicalSourceTerms;
          std::vector< bool > _activeShots;
    
    
          enum class sourceStatus {
            ACTIVE,
            INACTIVE,
            NOT_IN_NUMBERED_TERMS
          };
          sourceStatus indexStatusInShots(unsigned id) const;
    
    
    
         protected:
          void _runIterationOperations(const unsigned int iteration, const gmshfem::scalar::Precision< T_Scalar > rnorm);
          gmshfem::algebra::Vector< T_Scalar > _sortDofValues(const gmshfem::algebra::Vector< T_Scalar > &vector, const unsigned int fieldTag);
          void _fillG(const T_Scalar *array);
          void _fromPetscToInterfaceFields(Vec g);
    
          void _preProVolume();
          void _infoNumberOfDofs() const;
          void _assembleAllVolume();
          void _solveAllVolume();
          void _assembleAndSolveSurface();
          void _extractRHS();
    
          void _activateBilinear (const unsigned long long idom);
          void _deactivateBilinear (const unsigned long long idom);
    
          void _buildPetscRHS(Vec* RHS);
          void _buildPetscAllRHS(Mat* RHS);
    
         public:
          Formulation(const std::string &name, const std::vector< std::vector< unsigned int > > &topology);
          Formulation(const std::string &name, const std::vector< gmshfem::problem::Formulation< T_Scalar > * > &volume, const std::vector< std::unordered_map< unsigned int, gmshfem::problem::Formulation< T_Scalar > * > > &surface);
          ~Formulation();
    
          unsigned int size() const;
          unsigned int size(unsigned int i) const;
    
          const gmshfem::function::ScalarFunction< T_Scalar > physicalSource(const gmshfem::function::ScalarFunction< T_Scalar > &f);
          const gmshfem::function::VectorFunction< T_Scalar > physicalSource(const gmshfem::function::VectorFunction< T_Scalar > &f);
          void physicalSourceTerm(const unsigned int termTag);
          const gmshfem::function::ScalarFunction< T_Scalar > artificialSource(const gmshfem::function::ScalarFunction< T_Scalar > &f);
          const gmshfem::function::VectorFunction< T_Scalar > artificialSource(const gmshfem::function::VectorFunction< T_Scalar > &f);
          void artificialSourceTerm(const unsigned int termTag);
          void togglePhysicalAndArtificialSourceTerms();
    
          // Multi source options
          void setShotNumber(unsigned N);
          void enableShot(unsigned i);
          void disableShot(unsigned i);
          void toggleShot(unsigned i);
          bool isShotEnabled(unsigned i) const;
          void disableAllShots();
          void numberedPhysicalSourceTerm(const unsigned int id, const unsigned int termTag);
    
          // TMP
          void compareRHSexport();
    
    
          void addInterfaceField(field::InterfaceFieldInterface< T_Scalar > &interfaceField);
          void addInterfaceField(field::InterfaceFieldInterface< T_Scalar > &interfaceField, field::InterfaceFieldInterface< T_Scalar > &interfaceFieldMapping);
          field::InterfaceFieldInterface< T_Scalar > *getInterfaceField(const std::string &name) const;
    
          gmshfem::problem::Formulation< T_Scalar > &operator()(unsigned int i);
          gmshfem::problem::Formulation< T_Scalar > &operator()(unsigned int i, unsigned int j);
    
          gmshfem::common::Timer pre(bool mustAssemble = true);
          gmshfem::common::Timer assemble(int storeID = -1);
    
          gmshfem::common::Timer solve(const std::string &solver, const double tolerance = 1e-6, const int iterMax = 1000, const bool sameMatrixWithArtificialAndPhysicalSources = false, const bool skipFinalSolutionComputation = false);
          gmshfem::common::Timer solve(AbstractIterativeSolver<T_Scalar> &solver, const double tolerance = 1e-6, const int iterMax = 1000, const bool sameMatrixWithArtificialAndPhysicalSources = false, const bool skipFinalSolutionComputation = false);
          gmshfem::algebra::MatrixCCS< T_Scalar > computeMatrix();
          gmshfem::algebra::MatrixCCS< T_Scalar > computeMatrixBlock();
    
          void setSolutionIntoInterfaceFields(const std::vector< T_Scalar > &solution);
          const gmshfem::algebra::Vector< T_Scalar > &getRHS() const;
          gmshfem::algebra::Vector< T_Scalar > getGlobalRHS(); //Built on the fly, so not a reference. Modifies internally the Petsc representation of the local RHS.
    
    
          unsigned int numberOfIterations() const;
    
          friend int MatVectProductA< T_Scalar >(Mat A, Vec X, Vec Y);
          friend int MatVectProductI_A< T_Scalar >(Mat A, Vec X, Vec Y);
          friend int MatMatProductI_A< T_Scalar >(Mat A, Mat X, Mat Y);
          friend int MatMatProductA< T_Scalar >(Mat A, Mat X, Mat Y);
          friend int MatVectProductImpl< T_Scalar >(Mat A, Vec X, Vec Y, bool IA);
          friend int MatMatProductImpl< T_Scalar >(Mat A, Mat X, Mat Y, bool IA);
    
    
          // T_Scalar, PetscScalar, T_Interger
          template< class T_Scalar1, class T_Scalar2, class T_Integer >
          struct PetscInterface {
    
            static const T_Scalar1 *arrayInterface(const T_Scalar2 *array, const T_Integer size)
            {
              T_Scalar1 *ret = new T_Scalar1[size];
    #pragma omp parallel for
              for(T_Integer i = 0; i < size; ++i) {
                ret[i] = T_Scalar1(array[i]);
              }
              return ret;
            }
    
            static void freeArray(const T_Scalar1 *array)
            {
              delete[] array;
            }
          };
    
        };
    
    
      } // namespace problem
    
    
    } // namespace gmshddm
    
    #endif // H_GMSHDDM_FORMULATION