Panzer_BlockedEpetraLinearObjFactory.hpp

Go to the documentation of this file.

// @HEADER
// ***********************************************************************
//
//           Panzer: A partial differential equation assembly
//       engine for strongly coupled complex multiphysics systems
//                 Copyright (2011) Sandia Corporation
//
// Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
// the U.S. Government retains certain rights in this software.
//
// Redistribution and use in source and binary forms, with or without
// modification, are permitted provided that the following conditions are
// met:
//
// 1. Redistributions of source code must retain the above copyright
// notice, this list of conditions and the following disclaimer.
//
// 2. Redistributions in binary form must reproduce the above copyright
// notice, this list of conditions and the following disclaimer in the
// documentation and/or other materials provided with the distribution.
//
// 3. Neither the name of the Corporation nor the names of the
// contributors may be used to endorse or promote products derived from
// this software without specific prior written permission.
//
// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
//
// Questions? Contact Roger P. Pawlowski (rppawlo@sandia.gov) and
// Eric C. Cyr (eccyr@sandia.gov)
// ***********************************************************************
// @HEADER

#ifndef   __Panzer_BlockedEpetraLinearObjFactory_hpp__
#define   __Panzer_BlockedEpetraLinearObjFactory_hpp__

#include <map>

// Epetra includes
#include "Epetra_Map.h"
#include "Epetra_CrsGraph.h"
#include "Epetra_Import.h"
#include "Epetra_Export.h"

#include "PanzerDiscFE_config.hpp"
#include "Panzer_GlobalIndexer.hpp"
#include "Panzer_LinearObjFactory.hpp"
#include "Panzer_EpetraLinearObjContainer.hpp"
#include "Panzer_BlockedEpetraLinearObjContainer.hpp"
#include "Panzer_BlockedDOFManager.hpp"
#include "Panzer_CloneableEvaluator.hpp"

#include "Panzer_GatherOrientation.hpp"
#include "Panzer_GatherSolution_BlockedEpetra.hpp"
#include "Panzer_GatherTangent_BlockedEpetra.hpp"
#include "Panzer_ScatterResidual_BlockedEpetra.hpp"
#include "Panzer_ScatterDirichletResidual_BlockedEpetra.hpp"
#include "Panzer_ScatterResidual_Epetra.hpp"
#include "Panzer_ScatterDirichletResidual_Epetra.hpp"
#include "Panzer_GatherSolution_Epetra.hpp"
#include "Panzer_GatherTangent_Epetra.hpp"
#include "Panzer_GatherOrientation.hpp"
#include "Panzer_ThyraObjFactory.hpp"
#include "Panzer_HashUtils.hpp"

#include "Thyra_BlockedLinearOpBase.hpp"
#include "Thyra_ProductVectorBase.hpp"

#include "Teuchos_RCP.hpp"
#include "Teuchos_DefaultMpiComm.hpp"
#include "Teuchos_OpaqueWrapper.hpp"

namespace panzer
{

template <typename Traits,typename LocalOrdinalT>
class BlockedEpetraLinearObjFactory : public LinearObjFactory<Traits>
                                    , public ThyraObjFactory<double> {
public:

   BlockedEpetraLinearObjFactory(const Teuchos::RCP<const Teuchos::MpiComm<int> > & comm,
                                 const Teuchos::RCP<const GlobalIndexer> & gidProvider,
                                 bool useDiscreteAdjoint=false);

   BlockedEpetraLinearObjFactory(const Teuchos::RCP<const Teuchos::MpiComm<int> > & comm,
                                 const Teuchos::RCP<const GlobalIndexer> & gidProvider,
                                 const Teuchos::RCP<const GlobalIndexer> & colGidProvider,
                                 bool useDiscreteAdjoint=false);

   virtual ~BlockedEpetraLinearObjFactory();

/*************** Linear object factory methods *******************/

   virtual void readVector(const std::string & identifier,LinearObjContainer & loc,int id) const;

   virtual void writeVector(const std::string & identifier,const LinearObjContainer & loc,int id) const;

   virtual Teuchos::RCP<LinearObjContainer> buildLinearObjContainer() const;

   virtual Teuchos::RCP<LinearObjContainer> buildPrimitiveLinearObjContainer() const 
   { return buildLinearObjContainer(); }

   virtual Teuchos::RCP<LinearObjContainer> buildGhostedLinearObjContainer() const;

   virtual Teuchos::RCP<LinearObjContainer> buildPrimitiveGhostedLinearObjContainer() const 
   { return buildGhostedLinearObjContainer(); }

   virtual void globalToGhostContainer(const LinearObjContainer & container,
                                       LinearObjContainer & ghostContainer,int) const;
   virtual void ghostToGlobalContainer(const LinearObjContainer & ghostContainer,
                                       LinearObjContainer & container,int) const;

   virtual void adjustForDirichletConditions(const LinearObjContainer & localBCRows,
                                             const LinearObjContainer & globalBCRows,
                                             LinearObjContainer & ghostedObjs,
                                             bool zeroVectorRows=false, bool adjustX=false) const;

   virtual void applyDirichletBCs(const LinearObjContainer & counter,
                                  LinearObjContainer & result) const;

   virtual Teuchos::RCP<ReadOnlyVector_GlobalEvaluationData> buildReadOnlyDomainContainer() const;

   virtual Teuchos::RCP<WriteVector_GlobalEvaluationData> buildWriteDomainContainer() const;

   virtual Teuchos::MpiComm<int> getComm() const;

   template <typename EvalT>
   Teuchos::RCP<panzer::CloneableEvaluator> buildScatter() const
   { 
     if(!colDOFManagerContainer_->containsBlockedDOFManager() &&
        !rowDOFManagerContainer_->containsBlockedDOFManager())
       return Teuchos::rcp(new ScatterResidual_Epetra<EvalT,Traits,LocalOrdinalT,int>(rowDOFManagerContainer_->getFieldDOFManagers()[0],
                                                                                      colDOFManagerContainer_->getFieldDOFManagers()[0],
                                                                                      useDiscreteAdjoint_));

     return Teuchos::rcp(new ScatterResidual_BlockedEpetra<EvalT,Traits,LocalOrdinalT,int>(rowDOFManagerContainer_->getFieldDOFManagers(),
                                                                                           colDOFManagerContainer_->getFieldDOFManagers(),
                                                                                           useDiscreteAdjoint_)); 
   }

   template <typename EvalT>
   Teuchos::RCP<panzer::CloneableEvaluator > buildGather() const
   { 
     if(!colDOFManagerContainer_->containsBlockedDOFManager() &&
        !rowDOFManagerContainer_->containsBlockedDOFManager())
       return Teuchos::rcp(new GatherSolution_Epetra<EvalT,Traits,LocalOrdinalT,int>(rowDOFManagerContainer_->getFieldDOFManagers()[0]));
     return Teuchos::rcp(new GatherSolution_BlockedEpetra<EvalT,Traits,LocalOrdinalT,int>(rowDOFManagerContainer_->getFieldDOFManagers())); 
   }

   template <typename EvalT>
   Teuchos::RCP<panzer::CloneableEvaluator > buildGatherTangent() const
   { 
     if(!colDOFManagerContainer_->containsBlockedDOFManager() &&
        !rowDOFManagerContainer_->containsBlockedDOFManager())
       return Teuchos::rcp(new GatherTangent_Epetra<EvalT,Traits,LocalOrdinalT,int>(rowDOFManagerContainer_->getFieldDOFManagers()[0]));
     return Teuchos::rcp(new GatherTangent_BlockedEpetra<EvalT,Traits,LocalOrdinalT,int>(rowDOFManagerContainer_->getFieldDOFManagers())); 
   }

   template <typename EvalT>
   Teuchos::RCP<panzer::CloneableEvaluator > buildGatherDomain() const
   { 
     if(!colDOFManagerContainer_->containsBlockedDOFManager())
       return Teuchos::rcp(new GatherSolution_Epetra<EvalT,Traits,LocalOrdinalT,int>(colDOFManagerContainer_->getFieldDOFManagers()[0]));
     return Teuchos::rcp(new GatherSolution_BlockedEpetra<EvalT,Traits,LocalOrdinalT,int>(colDOFManagerContainer_->getFieldDOFManagers())); 
   }

   template <typename EvalT>
   Teuchos::RCP<panzer::CloneableEvaluator > buildGatherOrientation() const
  { return Teuchos::rcp(new GatherOrientation<EvalT,Traits,LocalOrdinalT,panzer::GlobalOrdinal>(rowDOFManagerContainer_->getFieldDOFManagers())); }

   template <typename EvalT>
   Teuchos::RCP<panzer::CloneableEvaluator> buildScatterDirichlet() const
   { 
     if(!colDOFManagerContainer_->containsBlockedDOFManager() &&
        !rowDOFManagerContainer_->containsBlockedDOFManager())
       return Teuchos::rcp(new ScatterDirichletResidual_Epetra<EvalT,Traits,LocalOrdinalT,int>(rowDOFManagerContainer_->getFieldDOFManagers()[0],
                                                                                               colDOFManagerContainer_->getFieldDOFManagers()[0])); 
     return Teuchos::rcp(new ScatterDirichletResidual_BlockedEpetra<EvalT,Traits,LocalOrdinalT,int>(rowDOFManagerContainer_->getFieldDOFManagers(),
                                                                                                    colDOFManagerContainer_->getFieldDOFManagers())); 
   }

/*************** Generic helper functions for container setup *******************/
   
   void initializeContainer(int,LinearObjContainer & loc) const;

   void initializeGhostedContainer(int,LinearObjContainer & loc) const;

/*************** Thyra based methods *******************/

   Teuchos::RCP<const Thyra::VectorSpaceBase<double> > getThyraDomainSpace() const;

   Teuchos::RCP<const Thyra::VectorSpaceBase<double> > getThyraRangeSpace() const;

   Teuchos::RCP<Thyra::VectorBase<double> > getThyraDomainVector() const;

   Teuchos::RCP<Thyra::VectorBase<double> > getThyraRangeVector() const;

   Teuchos::RCP<Thyra::LinearOpBase<double> > getThyraMatrix() const;

   // and now the ghosted versions

   Teuchos::RCP<const Thyra::VectorSpaceBase<double> > getGhostedThyraDomainSpace() const;

   Teuchos::RCP<const Thyra::VectorSpaceBase<double>>
   getGhostedThyraDomainSpace2() const;

   Teuchos::RCP<const Thyra::VectorSpaceBase<double> > getGhostedThyraRangeSpace() const;

   Teuchos::RCP<Thyra::VectorBase<double> > getGhostedThyraDomainVector() const;

   Teuchos::RCP<Thyra::VectorBase<double> > getGhostedThyraRangeVector() const;

   Teuchos::RCP<Thyra::LinearOpBase<double> > getGhostedThyraMatrix() const;

/*************** Epetra based methods *******************/

   virtual const Teuchos::RCP<Epetra_Map> getMap(int i) const;

   virtual const Teuchos::RCP<Epetra_Map> getColMap(int i) const;

   virtual const Teuchos::RCP<Epetra_Map> getGhostedMap(int i) const;

   virtual const Teuchos::RCP<Epetra_Map>
   getGhostedMap2(
     int i) const;

   virtual const Teuchos::RCP<Epetra_Map> getGhostedColMap(int i) const;

   virtual const Teuchos::RCP<Epetra_Map>
   getGhostedColMap2(
     int i) const;

   virtual const Teuchos::RCP<Epetra_CrsGraph> getGraph(int i,int j) const;

   virtual const Teuchos::RCP<Epetra_CrsGraph> getGhostedGraph(int i,int j) const;

   virtual const Teuchos::RCP<Epetra_Import> getGhostedImport(int i) const;

   virtual const Teuchos::RCP<Epetra_Import>
   getGhostedImport2(
     int i) const;

   virtual const Teuchos::RCP<Epetra_Import> getGhostedColImport(int i) const;

   virtual const Teuchos::RCP<Epetra_Import>
   getGhostedColImport2(
     int i) const;

   virtual const Teuchos::RCP<Epetra_Export> getGhostedExport(int j) const;

   virtual const Teuchos::RCP<Epetra_Export>
   getGhostedExport2(
     int i) const;

   virtual const Teuchos::RCP<Epetra_Export> getGhostedColExport(int j) const;

   virtual const Teuchos::RCP<Epetra_Export>
   getGhostedColExport2(
     int i) const;

   virtual const Teuchos::RCP<const Epetra_Comm> getEpetraComm() const;

   Teuchos::RCP<Epetra_CrsMatrix> getEpetraMatrix(int i,int j) const;
   Teuchos::RCP<Epetra_CrsMatrix> getGhostedEpetraMatrix(int i,int j) const;

   int getBlockRowCount() const;

   int getBlockColCount() const;

   Teuchos::RCP<const panzer::BlockedDOFManager> getGlobalIndexer() const
   { return rowDOFManagerContainer_->getBlockedIndexer(); }

   Teuchos::RCP<const panzer::GlobalIndexer> getRangeGlobalIndexer() const
   { return rowDOFManagerContainer_->getGlobalIndexer(); }

   Teuchos::RCP<const panzer::GlobalIndexer> getDomainGlobalIndexer() const
   { return colDOFManagerContainer_->getGlobalIndexer(); }

   const std::vector<Teuchos::RCP<const GlobalIndexer> > & getRangeGlobalIndexers() const
   { return rowDOFManagerContainer_->getFieldDOFManagers(); }

   const std::vector<Teuchos::RCP<const GlobalIndexer> > & getDomainGlobalIndexers() const
   { return colDOFManagerContainer_->getFieldDOFManagers(); }

   void addExcludedPair(int rowBlock,int colBlock);

   void addExcludedPairs(const std::vector<std::pair<int,int> > & exPairs);

protected:

   void initializeContainer_internal(int mem,ThyraObjContainer<double> & loc) const;

   void initializeGhostedContainer_internal(int mem,ThyraObjContainer<double> & loc) const;

/*************** Utility class for handling blocked and nonblocked DOF managers *******************/

  class DOFManagerContainer {
  public:
    DOFManagerContainer() {} 
    DOFManagerContainer(const Teuchos::RCP<const GlobalIndexer> & ugi) 
    { setGlobalIndexer(ugi); } 

    void setGlobalIndexer(const Teuchos::RCP<const GlobalIndexer> & ugi)
    {
      using Teuchos::RCP;
      using Teuchos::rcp_dynamic_cast;

      auto blockedDOFManager = rcp_dynamic_cast<const BlockedDOFManager>(ugi);
      auto flatDOFManager    = rcp_dynamic_cast<const GlobalIndexer>(ugi);

      if(blockedDOFManager!=Teuchos::null) {
        // set BlockedDOFManager
        blockedDOFManager_ = blockedDOFManager;
 
        // get all GID providers
        auto dofManagers =  blockedDOFManager_->getFieldDOFManagers();
        for(auto itr=dofManagers.begin();itr!=dofManagers.end();++itr)
          gidProviders_.push_back(*itr);
      }
      else if(flatDOFManager!=Teuchos::null) {
        // for absolute clarity, nullify the blockedDOFManager_
        blockedDOFManager_ = Teuchos::null;

        // you have only a single GID provider
        gidProviders_.push_back(flatDOFManager);
      }
      else {
        TEUCHOS_ASSERT(false);
      }
    }

    int getFieldBlocks() const
    { return Teuchos::as<int>(gidProviders_.size()); }

    bool containsBlockedDOFManager() const
    { return blockedDOFManager_ !=Teuchos::null; }

    Teuchos::RCP<const BlockedDOFManager> getBlockedIndexer() const
    {
      TEUCHOS_ASSERT(containsBlockedDOFManager());
      return blockedDOFManager_;
    }
 
    Teuchos::RCP<const GlobalIndexer> getGlobalIndexer() const
    {
      if(blockedDOFManager_!=Teuchos::null) 
        return blockedDOFManager_;

      TEUCHOS_ASSERT(gidProviders_.size()==1);
      return gidProviders_[0];
    }

    const std::vector<Teuchos::RCP<const GlobalIndexer> > & getFieldDOFManagers() const
    { return gidProviders_; }

  private:
    Teuchos::RCP<const BlockedDOFManager> blockedDOFManager_;
    std::vector<Teuchos::RCP<const GlobalIndexer> > gidProviders_;
  };

/*************** Generic methods/members *******************/

   // Get the global indexer associated with a particular block
   Teuchos::RCP<const GlobalIndexer> getGlobalIndexer(int i) const;

   Teuchos::RCP<const GlobalIndexer> getColGlobalIndexer(int i) const;

   void makeRoomForBlocks(std::size_t blockCnt,std::size_t colBlockCnt=0);

   Teuchos::RCP<const DOFManagerContainer> rowDOFManagerContainer_;
   Teuchos::RCP<const DOFManagerContainer> colDOFManagerContainer_;

   bool useColGidProviders_;

   // which block entries are ignored
   std::unordered_set<std::pair<int,int>,panzer::pair_hash> excludedPairs_;
  
/*************** Thyra based methods/members *******************/

   void ghostToGlobalThyraVector(const Teuchos::RCP<const Thyra::VectorBase<double> > & in,
                                 const Teuchos::RCP<Thyra::VectorBase<double> > & out,bool col) const;
   void ghostToGlobalThyraMatrix(const Thyra::LinearOpBase<double> & in,Thyra::LinearOpBase<double> & out) const;
   void globalToGhostThyraVector(const Teuchos::RCP<const Thyra::VectorBase<double> > & in,
                                 const Teuchos::RCP<Thyra::VectorBase<double> > & out,bool col) const;

   mutable Teuchos::RCP<const Thyra::VectorSpaceBase<double> > rangeSpace_;
   mutable Teuchos::RCP<const Thyra::VectorSpaceBase<double> > domainSpace_;

   mutable Teuchos::RCP<const Thyra::VectorSpaceBase<double> > ghostedRangeSpace_;
   mutable Teuchos::RCP<const Thyra::VectorSpaceBase<double> > ghostedDomainSpace_;

/*************** Epetra based methods/members *******************/

   void adjustForDirichletConditions(const Epetra_Vector & local_bcs,
                                     const Epetra_Vector & global_bcs,
                                     const Teuchos::Ptr<Epetra_Vector> & f,
                                     const Teuchos::Ptr<Epetra_CrsMatrix> & A,
                                     bool zeroVectorRows) const;

   void ghostToGlobalEpetraVector(int i,const Epetra_Vector & in,Epetra_Vector & out,bool col) const;
   void globalToGhostEpetraVector(int i,const Epetra_Vector & in,Epetra_Vector & out,bool col) const;
   void ghostToGlobalEpetraMatrix(int blockRow,const Epetra_CrsMatrix & in,Epetra_CrsMatrix & out) const;

   // get the map from the matrix

   virtual const Teuchos::RCP<Epetra_Map>
   buildMap(
     int i) const;

   virtual const Teuchos::RCP<Epetra_Map>
   buildGhostedMap(
     int i) const;

   virtual const Teuchos::RCP<Epetra_Map>
   buildGhostedMap2(
     int i) const;

   // get the map from the matrix

   virtual const Teuchos::RCP<Epetra_Map>
   buildColMap(
     int i) const;

   virtual const Teuchos::RCP<Epetra_Map>
   buildColGhostedMap(
     int i) const;

   virtual const Teuchos::RCP<Epetra_Map>
   buildColGhostedMap2(
     int i) const;

   // get the graph of the crs matrix
   virtual const Teuchos::RCP<Epetra_CrsGraph> buildGraph(int i,int j) const;
   virtual const Teuchos::RCP<Epetra_CrsGraph> buildGhostedGraph(int i,int j,bool optimizeStorage) const;
   virtual const Teuchos::RCP<Epetra_CrsGraph> buildFilteredGhostedGraph(int i,int j) const;

   // storage for Epetra graphs and maps
   Teuchos::RCP<const Epetra_Comm> eComm_;
   Teuchos::RCP<const Teuchos::OpaqueWrapper<MPI_Comm> > rawMpiComm_;
   Teuchos::RCP<Teuchos::MpiComm<int> > tComm_;

   mutable std::vector<Teuchos::RCP<Epetra_Map>> maps_;

   mutable std::vector<Teuchos::RCP<Epetra_Map>> ghostedMaps_;

   mutable std::vector<Teuchos::RCP<Epetra_Map>> ghostedMaps2_;

   mutable std::vector<Teuchos::RCP<Epetra_Import>> importers_;

   mutable std::vector<Teuchos::RCP<Epetra_Import>> importers2_;

   mutable std::vector<Teuchos::RCP<Epetra_Export>> exporters_;

   mutable std::vector<Teuchos::RCP<Epetra_Map>> colMaps_;

   mutable std::vector<Teuchos::RCP<Epetra_Map>> colGhostedMaps_;

   mutable std::vector<Teuchos::RCP<Epetra_Map>> colGhostedMaps2_;

   mutable std::vector<Teuchos::RCP<Epetra_Import>> colImporters_;

   mutable std::vector<Teuchos::RCP<Epetra_Import>> colImporters2_;

   mutable std::vector<Teuchos::RCP<Epetra_Export>> colExporters_;

   mutable std::unordered_map<std::pair<int,int>,Teuchos::RCP<Epetra_CrsGraph>,panzer::pair_hash> graphs_ ;
   mutable std::unordered_map<std::pair<int,int>,Teuchos::RCP<Epetra_CrsGraph>,panzer::pair_hash> ghostedGraphs_;

   bool useDiscreteAdjoint_;
};

} // end of namespace panzer

#endif // __Panzer_BlockedEpetraLinearObjFactory_hpp__