Working version

CMakeLists.txt

@@ -114,10 +114,6 @@ CMAKE_DEPENDENT_OPTION(NEKTAR_USE_DIRECT_BLAS_CALLS
     "NEKTAR_USE_BLAS_LAPACK" OFF)
 MARK_AS_ADVANCED(NEKTAR_USE_DIRECT_BLAS_CALLS)
 
-# Support for PETSc solvers
-OPTION(NEKTAR_USE_PETSC
-    "Enable PETSc parallel matrix solver support." OFF)
-
 OPTION(NEKTAR_USE_MEMORY_POOLS
     "Use memory pools to accelerate memory allocation." ON)
 MARK_AS_ADVANCED(NEKTAR_USE_MEMORY_POOLS)
@@ -172,6 +168,7 @@ INCLUDE (ThirdPartyZlib)
 INCLUDE (ThirdPartyFFTW)
 INCLUDE (ThirdPartyArpack)
 INCLUDE (ThirdPartyMPI)
+INCLUDE (ThirdPartyPETSc)
 INCLUDE (ThirdPartyVTK)
 INCLUDE (ThirdPartyQT4)
 INCLUDE (ThirdPartySMV)
@@ -227,10 +224,6 @@ IF( NEKTAR_USE_BLAS_LAPACK )
     ADD_DEFINITIONS( -DNEKTAR_USING_LAPACK -DNEKTAR_USING_BLAS )
 ENDIF( NEKTAR_USE_BLAS_LAPACK )
 
-IF( NEKTAR_USE_PETSC )
-    INCLUDE (FindPETSc)
-ENDIF( NEKTAR_USE_PETSC )
-
 IF( NEKTAR_USE_TINYXML_STL )
     ADD_DEFINITIONS( -DTIXML_USE_STL)
 ENDIF( NEKTAR_USE_TINYXML_STL )

cmake/NektarCommon.cmake

@@ -163,6 +163,9 @@ MACRO(ADD_NEKTAR_EXECUTABLE name component sources)
                 PROPERTIES COMPILE_FLAGS "${THE_COMPILE_FLAGS} -DMKL_ILP64")
     ENDIF( NEKTAR_USE_MKL AND MKL_FOUND )
         
+    IF( NEKTAR_USE_PETSC )
+        TARGET_LINK_LIBRARIES(${name} ${PETSC_LIBRARIES} )
+    ENDIF( NEKTAR_USE_PETSC )
 
     TARGET_LINK_LIBRARIES(${name}
         optimized LibUtilities debug LibUtilities-g

library/LibUtilities/BasicUtils/SessionReader.cpp

@@ -1654,7 +1654,12 @@ namespace Nektar
                         "IterativeMultiLevelStaticCond"                    ||
                     m_solverInfo["GLOBALSYSSOLN"] == "XxtFull"             ||
                     m_solverInfo["GLOBALSYSSOLN"] == "XxtStaticCond"       ||
-                    m_solverInfo["GLOBALSYSSOLN"] == "XxtMultiLevelStaticCond",
+                    m_solverInfo["GLOBALSYSSOLN"] ==
+                        "XxtMultiLevelStaticCond"                          ||
+                    m_solverInfo["GLOBALSYSSOLN"] == "PETScFull"           ||
+                    m_solverInfo["GLOBALSYSSOLN"] == "PETScStaticCond"     ||
+                    m_solverInfo["GLOBALSYSSOLN"] ==
+                        "PETScMultiLevelStaticCond",
                     "A parallel solver must be used when run in parallel.");
             }
             

library/MultiRegions/AssemblyMap/AssemblyMapCG2D.cpp

@@ -1286,6 +1286,7 @@ namespace Nektar
                 case eIterativeStaticCond:
                 case eXxtFullMatrix:
                 case eXxtStaticCond:
+                case ePETScStaticCond:
                     {
                         NoReordering(boostGraphObj,perm,iperm);
                     }
@@ -1298,6 +1299,7 @@ namespace Nektar
                 case eDirectMultiLevelStaticCond:
                 case eIterativeMultiLevelStaticCond:
                 case eXxtMultiLevelStaticCond:
+                case ePETScMultiLevelStaticCond:
                     {
                         MultiLevelBisectionReordering(boostGraphObj,perm,iperm,bottomUpGraph,partVerts,mdswitch);
                     }

library/MultiRegions/CMakeLists.txt

@@ -116,6 +116,18 @@ IF(NEKTAR_USE_MPI)
     )
 ENDIF(NEKTAR_USE_MPI)
 
+IF(NEKTAR_USE_PETSC)
+    SET(MULTI_REGIONS_HEADERS ${MULTI_REGIONS_HEADERS}
+        GlobalLinSysPETSc.h
+        GlobalLinSysPETScStaticCond.h
+    )
+    SET(MULTI_REGIONS_SOURCES ${MULTI_REGIONS_SOURCES}
+        GlobalLinSysPETSc.cpp
+        GlobalLinSysPETScStaticCond.cpp
+    )
+ENDIF(NEKTAR_USE_PETSC)
+
+
 ADD_DEFINITIONS(-DMULTI_REGIONS_EXPORTS)
 
 ADD_NEKTAR_LIBRARY(MultiRegions lib ${NEKTAR_LIBRARY_TYPE} ${MULTI_REGIONS_SOURCES} ${MULTI_REGIONS_HEADERS} ${ASSEMBLY_MAP_HEADERS})

library/MultiRegions/GlobalLinSys.cpp

@@ -45,7 +45,7 @@ namespace Nektar
 {
     namespace MultiRegions
     {
-        std::string GlobalLinSys::lookupIds[8] = {
+        std::string GlobalLinSys::lookupIds[10] = {
             LibUtilities::SessionReader::RegisterEnumValue(
                 "GlobalSysSoln", "DirectFull",
                 MultiRegions::eDirectFullMatrix),
@@ -69,7 +69,13 @@ namespace Nektar
                 MultiRegions::eXxtFullMatrix),
             LibUtilities::SessionReader::RegisterEnumValue(
                 "GlobalSysSoln", "XxtStaticCond",
-                MultiRegions::eXxtStaticCond)
+                MultiRegions::eXxtStaticCond),
+            LibUtilities::SessionReader::RegisterEnumValue(
+                "GlobalSysSoln", "PETScStaticCond",
+                MultiRegions::ePETScStaticCond),
+            LibUtilities::SessionReader::RegisterEnumValue(
+                "GlobalSysSoln", "PETScMultiLevelStaticCond",
+                MultiRegions::ePETScMultiLevelStaticCond)
         };
 
         std::string GlobalLinSys::def = LibUtilities::SessionReader::

library/MultiRegions/GlobalLinSysPETSc.cpp

+///////////////////////////////////////////////////////////////////////////////
+//
+// File GlobalLinSys.cpp
+//
+// For more information, please see: http://www.nektar.info
+//
+// The MIT License
+//
+// Copyright (c) 2006 Division of Applied Mathematics, Brown University (USA),
+// Department of Aeronautics, Imperial College London (UK), and Scientific
+// Computing and Imaging Institute, University of Utah (USA).
+//
+// License for the specific language governing rights and limitations under
+// Permission is hereby granted, free of charge, to any person obtaining a
+// copy of this software and associated documentation files (the "Software"),
+// to deal in the Software without restriction, including without limitation
+// the rights to use, copy, modify, merge, publish, distribute, sublicense,
+// and/or sell copies of the Software, and to permit persons to whom the
+// Software is furnished to do so, subject to the following conditions:
+//
+// The above copyright notice and this permission notice shall be included
+// in all copies or substantial portions of the Software.
+//
+// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
+// OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
+// THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
+// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
+// DEALINGS IN THE SOFTWARE.
+//
+// Description: GlobalLinSys definition
+//
+///////////////////////////////////////////////////////////////////////////////
+
+#include <MultiRegions/GlobalLinSysPETSc.h>
+
+namespace Nektar
+{
+    namespace MultiRegions
+    {
+        /**
+         * @class GlobalLinSysPETSc
+         *
+         * Solves a linear system using direct methods.
+         */
+
+
+        /// Constructor for full direct matrix solve.
+        GlobalLinSysPETSc::GlobalLinSysPETSc(
+            const GlobalLinSysKey                &pKey,
+            const boost::weak_ptr<ExpList>       &pExp,
+            const boost::shared_ptr<AssemblyMap> &pLocToGloMap)
+            : GlobalLinSys(pKey, pExp, pLocToGloMap)
+        {
+            // Initialise PETSc
+            PetscInitialize(0, NULL, NULL, NULL);
+
+            // Create matrix
+            MatCreate(PETSC_COMM_WORLD, &m_matrix);
+        }
+
+        GlobalLinSysPETSc::~GlobalLinSysPETSc()
+        {
+        }
+
+        void GlobalLinSysPETSc::v_SolveLinearSystem(
+            const int                          pNumRows,
+            const Array<OneD,const NekDouble> &pInput,
+                  Array<OneD,      NekDouble> &pOutput,
+            const AssemblyMapSharedPtr        &locToGloMap,
+            const int                          pNumDir)
+        {
+            int i;
+            for (i = pNumDir; i < pNumRows; ++i)
+            {
+                VecSetValue(m_b, i-pNumDir, pInput[i], INSERT_VALUES);
+            }
+
+            PetscErrorCode ierr = KSPSolve(m_ksp, m_b, m_x);
+            PetscScalar   *avec;
+            VecGetArray(m_x, &avec);
+            KSPView(m_ksp, PETSC_VIEWER_STDOUT_WORLD);
+
+            PetscInt its;
+            KSPGetIterationNumber(m_ksp,&its);
+
+            cout << "iteration = " << its << endl;
+
+            for (i = 0; i < pNumRows; ++i)
+            {
+                pOutput[i] = avec[i];
+            }
+        }
+    }
+}

library/MultiRegions/GlobalLinSysPETSc.h

+///////////////////////////////////////////////////////////////////////////////
+//
+// File GlobalLinSys.h
+//
+// For more information, please see: http://www.nektar.info
+//
+// The MIT License
+//
+// Copyright (c) 2006 Division of Applied Mathematics, Brown University (USA),
+// Department of Aeronautics, Imperial College London (UK), and Scientific
+// Computing and Imaging Institute, University of Utah (USA).
+//
+// License for the specific language governing rights and limitations under
+// Permission is hereby granted, free of charge, to any person obtaining a
+// copy of this software and associated documentation files (the "Software"),
+// to deal in the Software without restriction, including without limitation
+// the rights to use, copy, modify, merge, publish, distribute, sublicense,
+// and/or sell copies of the Software, and to permit persons to whom the
+// Software is furnished to do so, subject to the following conditions:
+//
+// The above copyright notice and this permission notice shall be included
+// in all copies or substantial portions of the Software.
+//
+// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
+// OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
+// THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
+// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
+// DEALINGS IN THE SOFTWARE.
+//
+// Description: GlobalLinSysPETSc header
+//
+///////////////////////////////////////////////////////////////////////////////
+#ifndef NEKTAR_LIB_MULTIREGIONS_GLOBALLINSYSPETSC_H
+#define NEKTAR_LIB_MULTIREGIONS_GLOBALLINSYSPETSC_H
+
+#include <MultiRegions/MultiRegionsDeclspec.h>
+#include <MultiRegions/GlobalLinSys.h>
+
+#include <petscmat.h>
+#include <petscksp.h>
+
+namespace Nektar
+{
+    namespace MultiRegions
+    {
+        // Forward declarations
+        class ExpList;
+
+        /// A global linear system.
+        class GlobalLinSysPETSc : public GlobalLinSys
+        {
+        public:
+            /// Constructor for full direct matrix solve.
+            MULTI_REGIONS_EXPORT GlobalLinSysPETSc(
+                const GlobalLinSysKey                &pKey,
+                const boost::weak_ptr<ExpList>       &pExp,
+                const boost::shared_ptr<AssemblyMap> &pLocToGloMap);
+            
+            MULTI_REGIONS_EXPORT virtual ~GlobalLinSysPETSc();
+
+            virtual void v_SolveLinearSystem(
+                const int                          pNumRows,
+                const Array<OneD,const NekDouble> &pInput,
+                      Array<OneD,      NekDouble> &pOutput,
+                const AssemblyMapSharedPtr        &locToGloMap,
+                const int                          pNumDir);
+
+        protected:
+            Mat m_matrix;
+            Vec m_x, m_b;
+            KSP m_ksp;
+        };
+    }
+}
+
+#endif

library/MultiRegions/GlobalLinSysPETScStaticCond.h

+///////////////////////////////////////////////////////////////////////////////
+//
+// File GlobalLinSys.h
+//
+// For more information, please see: http://www.nektar.info
+//
+// The MIT License
+//
+// Copyright (c) 2006 Division of Applied Mathematics, Brown University (USA),
+// Department of Aeronautics, Imperial College London (UK), and Scientific
+// Computing and Imaging Institute, University of Utah (USA).
+//
+// License for the specific language governing rights and limitations under
+// Permission is hereby granted, free of charge, to any person obtaining a
+// copy of this software and associated documentation files (the "Software"),
+// to deal in the Software without restriction, including without limitation
+// the rights to use, copy, modify, merge, publish, distribute, sublicense,
+// and/or sell copies of the Software, and to permit persons to whom the
+// Software is furnished to do so, subject to the following conditions:
+//
+// The above copyright notice and this permission notice shall be included
+// in all copies or substantial portions of the Software.
+//
+// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
+// OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
+// THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
+// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
+// DEALINGS IN THE SOFTWARE.
+//
+// Description: GlobalLinSysStaticCond header
+//
+///////////////////////////////////////////////////////////////////////////////
+#ifndef NEKTAR_LIB_MULTIREGIONS_GLOBALLINSYSDIRECTSTATICCOND_H
+#define NEKTAR_LIB_MULTIREGIONS_GLOBALLINSYSDIRECTSTATICCOND_H
+
+#include <MultiRegions/GlobalLinSysPETSc.h>
+#include <MultiRegions/MultiRegionsDeclspec.h>
+
+namespace Nektar
+{
+    namespace MultiRegions
+    {
+        // Forward declarations
+        class ExpList;
+        class GlobalLinSysPETScStaticCond;
+
+        typedef boost::shared_ptr<GlobalLinSysPETScStaticCond>
+                                        GlobalLinSysPETScStaticCondSharedPtr;
+
+        /// A global linear system.
+        class GlobalLinSysPETScStaticCond : public GlobalLinSysPETSc
+        {
+        public:
+            /// Creates an instance of this class
+            static GlobalLinSysSharedPtr create(
+                        const GlobalLinSysKey &pLinSysKey,
+                        const boost::weak_ptr<ExpList> &pExpList,
+                        const boost::shared_ptr<AssemblyMap>
+                                                               &pLocToGloMap)
+            {
+                return MemoryManager<GlobalLinSysPETScStaticCond>
+                    ::AllocateSharedPtr(pLinSysKey, pExpList, pLocToGloMap);
+            }
+
+            /// Name of class
+            MULTI_REGIONS_EXPORT static std::string className;
+            static std::string className2;
+
+            /// Constructor for full direct matrix solve.
+            MULTI_REGIONS_EXPORT GlobalLinSysPETScStaticCond(
+                        const GlobalLinSysKey &mkey,
+                        const boost::weak_ptr<ExpList> &pExpList,
+                        const boost::shared_ptr<AssemblyMap>
+                                                                &locToGloMap);
+
+            /// Constructor for full direct matrix solve.
+            MULTI_REGIONS_EXPORT GlobalLinSysPETScStaticCond(
+                        const GlobalLinSysKey &mkey,
+                        const boost::weak_ptr<ExpList> &pExpList,
+                        const DNekScalBlkMatSharedPtr pSchurCompl,
+                        const DNekScalBlkMatSharedPtr pBinvD,
+                        const DNekScalBlkMatSharedPtr pC,
+                        const DNekScalBlkMatSharedPtr pInvD,
+                        const boost::shared_ptr<AssemblyMap>
+                                                                &locToGloMap);
+
+            MULTI_REGIONS_EXPORT virtual ~GlobalLinSysPETScStaticCond();
+
+        private:
+            /// Schur complement for Direct Static Condensation.
+            GlobalLinSysPETScStaticCondSharedPtr m_recursiveSchurCompl;
+
+            /// Block matrices at this level
+            DNekScalBlkMatSharedPtr m_schurCompl;
+            DNekScalBlkMatSharedPtr m_BinvD;
+            DNekScalBlkMatSharedPtr m_C;
+            DNekScalBlkMatSharedPtr m_invD;
+
+            /// Solve the linear system for given input and output vectors
+            /// using a specified local to global map.
+            virtual void v_Solve(
+                        const Array<OneD, const NekDouble>  &in,
+                              Array<OneD,       NekDouble>  &out,
+                        const AssemblyMapSharedPtr &locToGloMap,
+                        const Array<OneD, const NekDouble>  &dirForcing
+                                                        = NullNekDouble1DArray);
+
+            /// Initialise this object
+            void Initialise(
+                const boost::shared_ptr<AssemblyMap>& locToGloMap);
+
+            /// Set up the storage for the Schur complement or the top level
+            /// of the multi-level Schur complement.
+            void SetupTopLevel(
+                    const boost::shared_ptr<AssemblyMap>& locToGloMap);
+
+            /// Assemble the Schur complement matrix.
+            void AssembleSchurComplement(
+                const boost::shared_ptr<AssemblyMap>& locToGloMap);
+
+            ///
+            void ConstructNextLevelCondensedSystem(
+                const boost::shared_ptr<AssemblyMap>& locToGloMap);
+         };
+    }
+}
+
+#endif

library/MultiRegions/MultiRegions.hpp

@@ -86,6 +86,9 @@ namespace Nektar
             eXxtFullMatrix,
             eXxtStaticCond,
             eXxtMultiLevelStaticCond,
+            ePETScFullMatrix,
+            ePETScStaticCond,
+            ePETScMultiLevelStaticCond,
             eSIZE_GlobalSysSolnType
         };
 
@@ -101,7 +104,10 @@ namespace Nektar
             "IterativeMultiLevelStaticCond",
             "XxtFull",
             "XxtStaticCond",
-            "XxtMultiLevelStaticCond"
+            "XxtMultiLevelStaticCond",
+            "PETScFull",
+            "PETScStaticCond",
+            "PETScMultiLevelStaticCond"
         };
 
         /// Type of Galerkin projection.