Merge pull request #3738 from crdohrm/bddc

Bddc

packages/shylu/shylu_dd/bddc/CMakeLists.txt

@@ -5,7 +5,6 @@ ADD_SUBDIRECTORY(src)
 
 TRIBITS_ADD_TEST_DIRECTORIES(test)
 
-
-#TRIBITS_ADD_EXAMPLE_DIRECTORIES(example)
+TRIBITS_ADD_EXAMPLE_DIRECTORIES(example)
 
 TRIBITS_SUBPACKAGE_POSTPROCESS()

packages/shylu/shylu_dd/bddc/cmake/Dependencies.cmake

@@ -1,9 +1,9 @@
 SET(LIB_REQUIRED_DEP_PACKAGES Tpetra Teuchos Zoltan)
-SET(LIB_OPTIONAL_DEP_PACKAGES ShyLU_NodeTacho)
+SET(LIB_OPTIONAL_DEP_PACKAGES Amesos2 MueLu ShyLU_NodeTacho)
 SET(TEST_REQUIRED_DEP_PACKAGES Gtest)
 SET(TEST_OPTIONAL_DEP_PACKAGES)
-SET(LIB_REQUIRED_DEP_TPLS MPI METIS)
-SET(LIB_OPTIONAL_DEP_TPLS MKL PARDISO_MKL SuperLU)
+SET(LIB_REQUIRED_DEP_TPLS MPI BLAS LAPACK)
+SET(LIB_OPTIONAL_DEP_TPLS MKL PARDISO_MKL SuperLU METIS ParMETIS)
 SET(TEST_REQUIRED_DEP_TPLS)
 SET(TEST_OPTIONAL_DEP_TPLS)
 

packages/shylu/shylu_dd/bddc/cmake/ShyLU_DDBDDC_config.h.in

@@ -4,10 +4,7 @@
  ShyLU BDDC configuration.
 */
 
-/* define if we want to use MPI */
-#cmakedefine HAVE_SHYLU_DDBDDC_MPI
-
-#cmakedefine HAVE_SHYLU_DDBDDC_TPETRA
+#cmakedefine HAVE_SHYLU_DDBDDC_PARMETIS
 
 #cmakedefine HAVE_SHYLU_DDBDDC_METIS
 
@@ -19,6 +16,10 @@
 
 #cmakedefine HAVE_SHYLU_DDBDDC_SUPERLU
 
+#cmakedefine HAVE_SHYLU_DDBDDC_MUELU
+
 #cmakedefine HAVE_SHYLU_DDBDDC_SHYLU_NODETACHO
 
+#cmakedefine HAVE_SHYLU_DDBDDC_AMESOS2
+
 #endif // SHYLU_DDBDDCCONFIG_H

packages/shylu/shylu_dd/bddc/example/CMakeLists.txt

@@ -0,0 +1,27 @@
+
+INCLUDE_DIRECTORIES(${CMAKE_CURRENT_SOURCE_DIR})
+INCLUDE_DIRECTORIES(../test)
+INCLUDE_DIRECTORIES(${CMAKE_CURRENT_BINARY_DIR})
+
+TRIBITS_ADD_EXECUTABLE_AND_TEST(
+ example1
+ SOURCES example1.cpp
+ DEPLIBS
+ ARGS ""
+ COMM mpi
+ NUM_MPI_PROCS 8
+)
+
+TRIBITS_ADD_EXECUTABLE_AND_TEST(
+ simpleExample1
+ SOURCES simpleExample1.cpp
+ DEPLIBS
+ ARGS ""
+ COMM mpi
+ NUM_MPI_PROCS 8
+)
+
+TRIBITS_COPY_FILES_TO_BINARY_DIR(ShyLUBDDCExampleFiles
+ SOURCE_FILES problemMaker.xml problemMakerSimple.xml bddc.xml
+)
+

packages/shylu/shylu_dd/bddc/example/bddc.xml

@@ -0,0 +1,48 @@
+<ParameterList name="BDDC">
+ <Parameter name="Use Corners" type="bool" value="false"/>
+ <Parameter name="Use Edges" type="bool" value="true"/>
+ <Parameter name="Use Faces" type="bool" value="false"/>
+ <Parameter name="Corner Weight Type" type="string" value="stiffness"/>
+ <Parameter name="Edge Weight Type" type="string" value="stiffness"/>
+ <Parameter name="Face Weight Type" type="string" value="stiffness"/>
+ <Parameter name="Use Vertex Coarse Space" type="bool" value="false"/>
+ <Parameter name="Economic Vertex Coarse Space" type="bool" value="false"/>
+ <Parameter name="Interface Preconditioner" type="bool" value="false"/>
+ <Parameter name="Use PComm" type="bool" value="true"/>
+ <!-- Krylov Solver options: -->
+ <!-- PCG: preconditioned conjugate gradient -->
+ <!-- GCR: generalized conjugate residual -->
+ <Parameter name="Krylov Solver" type="string" value="PCG"/>
+ <Parameter name="Maximum Iterations" type="int" value="20"/>
+ <Parameter name="Convergence Tolerance" type="double" value="1e-6"/>
+ <Parameter name="Convergence Check Interval" type="int" value="100"/>
+ <Parameter name="Maximum Stored Directions" type="int" value="200"/>
+ <!-- Solver options: -->
+ <!-- LAPACK: uses LAPACK banded or dense solver (not best option) -->
+ <!-- SuperLU: requires enabling SuperLU TPL -->
+ <!-- Pardiso: requires enabling Intel/MKL Pardiso TPL -->
+ <!-- KLU2: requires enabling Amesos2 package -->
+ <!-- Tacho: requires enabling Tacho in Shylu_node subpackage -->
+ <!-- MueLu: requires enabling MueLu package -->
+ <Parameter name="Dirichlet Solver" type="string" value="LAPACK"/>
+ <Parameter name="Neumann Solver" type="string" value="LAPACK"/>
+ <Parameter name="Coarse Solver" type="string" value="LAPACK"/>
+ <Parameter name="Estimate Condition Number" type="bool" value="true"/>
+ <Parameter name="Print Summary" type="int" value="3"/>
+ <Parameter name="numSubdomainsPerCoarseSubdomain" type="int" value="1000"/>
+ <Parameter name="numCoarseSubdomainsPerMpiRank" type="int" value="1"/>
+ <!-- coarsening options: -->
+ <!-- Graph: graph-based -->
+ <!-- RCB: recursive coordinate bisection -->
+ <!-- RIB: recursive inetial bisection -->
+ <Parameter name="Coarsening Option" type="string" value="RCB"/>
+ <!-- graph package options = PHG or ParMETIS -->
+ <Parameter name="Graph Package" type="string" value="PHG"/>
+ <Parameter name="Output Timing" type="bool" value="false"/>
+ <Parameter name="Operator Name" type="string" value="none"/> 
+ <Parameter name="Print Coarsest Matrix" type="bool" value="false"/>
+ <Parameter name="Print Vertex Matrix" type="bool" value="false"/>
+ <Parameter name="Print Psi Matrix" type="bool" value="false"/>
+ <Parameter name="Print Operator Matrix" type="bool" value="false"/>
+ <Parameter name="Print Preconditioner Matrix" type="bool" value="false"/>
+</ParameterList>

packages/shylu/shylu_dd/bddc/example/example1.cpp

@@ -0,0 +1,205 @@
+
+//@HEADER
+// ************************************************************************
+// 
+// ShyLU: Hybrid preconditioner package
+// Copyright 2012 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 Michael A. Heroux ([email protected]) 
+// 
+// ************************************************************************
+//@HEADER
+
+#include <mpi.h>
+#include "ShyLU_DDBDDC_config.h"
+#if defined(HAVE_SHYLU_DDBDDC_SHYLU_NODETACHO)
+#include "Kokkos_Core.hpp"
+#endif
+#include "shylu_PreconditionerBDDC.hpp"
+#include "shylu_KrylovSolver.hpp"
+#include "setupTest.hpp"
+#include "Teuchos_XMLParameterListHelpers.hpp"
+#include "Tpetra_Core.hpp"
+
+using Teuchos::RCP;
+
+typedef int LO; // Local Ordinal
+typedef Tpetra::Map<>::global_ordinal_type GO; // Global Ordinal
+typedef double SX; // floating point data type
+typedef double SM; // real (magnitude) for SX
+
+void readInputFiles(MPI_Comm Comm,
+ const std::string fileNamePM, 
+ const std::string fileNameBDDC, 
+ const std::string fileNameMueLu, 
+ const std::string fileNameNodalAMG, 
+ RCP<Teuchos::ParameterList> & parametersPM,
+ RCP<Teuchos::ParameterList> & parametersBDDC,
+ RCP<Teuchos::ParameterList> & parametersMueLu,
+ RCP<Teuchos::ParameterList> & parametersNodalAMG)
+{
+ RCP<const Teuchos::Comm<int> > TComm = 
+ Teuchos::rcp( new Teuchos::MpiComm<int>(Comm) );
+ Teuchos::ParameterList paramsPM, paramsBDDC, paramsMueLu, paramsNodalAMG;
+
+ Teuchos::updateParametersFromXmlFileAndBroadcast
+ (fileNamePM, Teuchos::Ptr<Teuchos::ParameterList>(&paramsPM), *TComm);
+ parametersPM = Teuchos::rcp( new Teuchos::ParameterList(paramsPM) );
+
+ Teuchos::updateParametersFromXmlFileAndBroadcast
+ (fileNameBDDC, Teuchos::Ptr<Teuchos::ParameterList>(&paramsBDDC), *TComm);
+ parametersBDDC = Teuchos::rcp( new Teuchos::ParameterList(paramsBDDC) );
+
+ if (fileNameMueLu != "") {
+ Teuchos::updateParametersFromXmlFileAndBroadcast
+ (fileNameMueLu, Teuchos::Ptr<Teuchos::ParameterList>(&paramsMueLu), 
+ *TComm);
+ parametersMueLu = Teuchos::rcp( new Teuchos::ParameterList(paramsMueLu) );
+ }
+
+ if (fileNameNodalAMG != "") {
+ Teuchos::updateParametersFromXmlFileAndBroadcast
+ (fileNameNodalAMG, Teuchos::Ptr<Teuchos::ParameterList>(&paramsNodalAMG),
+ *TComm);
+ parametersNodalAMG = 
+ Teuchos::rcp( new Teuchos::ParameterList(paramsNodalAMG) );
+ }
+}
+
+void runTest(RCP<Teuchos::ParameterList> & parametersPM,
+ RCP<Teuchos::ParameterList> & parametersBDDC,
+ RCP<Teuchos::ParameterList> & parametersMueLu,
+ RCP<Teuchos::ParameterList> & parametersNodalAMG,
+ const std::string meshDataFile,
+ const int myPID,
+ MPI_Comm Comm,
+ const bool resetFile,
+ int & numIterations)
+{
+ // setup test data
+ bddc::setupTest<LO,GO,SX,SM> 
+ test(parametersPM, parametersBDDC, parametersMueLu,
+ parametersNodalAMG, meshDataFile, Comm);
+ LO numNode, *nodeBegin(nullptr), *localDofs(nullptr);
+ const GO *nodeGlobalIDs(nullptr);
+ const SM *xCoord(nullptr), *yCoord(nullptr), *zCoord(nullptr);
+ std::vector<LO*> subRowBeginPtr, subColumnsPtr;
+ std::vector<SX*> subValuesPtr;
+ bddc::OperatorBase<SX>* Operator(nullptr);
+ // extract test data
+ test.getProblemData(numNode, nodeBegin, localDofs, nodeGlobalIDs,
+ xCoord, yCoord, zCoord, subRowBeginPtr,
+ subColumnsPtr, subValuesPtr, Operator);
+ const std::vector< std::vector<LO> > & subNodes = test.getSubNodes();
+ std::vector<int> nodeSend;
+ bddc::getNodeSend(numNode, nodeGlobalIDs, Comm, nodeSend);
+ // initialize preconditioner
+ int level(0);
+ MPI_Barrier(Comm);
+ double startTimeBDDCPre = test.getTime();
+ RCP< bddc::PreconditionerBDDC<SX,SM,LO,GO> > Preconditioner =
+ rcp( new bddc::PreconditionerBDDC<SX,SM,LO,GO>
+ (numNode, nodeBegin, localDofs, nodeGlobalIDs, xCoord, yCoord, zCoord, 
+ subNodes, subRowBeginPtr.data(), subColumnsPtr.data(), 
+ subValuesPtr.data(), parametersBDDC, Comm, level, Operator, 
+ &nodeSend) );
+ MPI_Barrier(Comm);
+ double stopTimeBDDCPre = test.getTime();
+ const LO numMyRows = Preconditioner->NumMyRows();
+ std::vector<SX> rhs;
+ test.getRhs(numMyRows, rhs);
+ double startTimeKrylovInit(0), startTimeKrylovSolve(0);
+ double stopTimeKrylovInit(0), stopTimeKrylovSolve(0);
+ test.openFiles(Preconditioner, resetFile);
+ MPI_Barrier(Comm);
+ // initialize Krylov solver
+ startTimeKrylovInit = test.getTime();
+ RCP< bddc::KrylovSolver<SX,SM,LO,GO> > Solver =
+ rcp ( new bddc::KrylovSolver<SX,SM,LO,GO>(Preconditioner, parametersBDDC) );
+ std::vector<SX> sol(numMyRows), Ax(numMyRows);
+ MPI_Barrier(Comm);
+ // solve equations
+ stopTimeKrylovInit = test.getTime();
+ startTimeKrylovSolve = test.getTime();
+ LO numberOfSolves = parametersPM->get("Number of Solves", 1);
+ for (int i=0; i<numberOfSolves; i++) {
+ if (i > 0) test.updateRhs(rhs.data(), numMyRows);
+ Solver->Solve(rhs.data(), sol.data());
+ }
+ MPI_Barrier(Comm);
+ stopTimeKrylovSolve = test.getTime();
+ Preconditioner->ApplyFullOperator(sol.data(), Ax.data());
+ for (LO i=0; i<numMyRows; i++) {
+ Ax[i] = rhs[i] - Ax[i];
+ }
+ // check solution
+ SM normError = Preconditioner->Norm2(Ax.data(), numMyRows);
+ SM normRhs = Preconditioner->Norm2(rhs.data(), numMyRows);
+ if (myPID == 0) {
+ std::cout << "BDDC Preconditioner initialization time = "
+ << stopTimeBDDCPre - startTimeBDDCPre << std::endl;
+ std::cout << "KrylovSolver initialization time = "
+ << stopTimeKrylovInit - startTimeKrylovInit << std::endl;
+ std::cout << "KrylovSolver solve solve time = "
+ << stopTimeKrylovSolve - startTimeKrylovSolve << std::endl;
+ std::cout << "relative residual = " << normError/normRhs << std::endl;
+ }
+ const bool outputTiming = parametersBDDC->get("Output Timing", false);
+ if (outputTiming) {
+ int krylovMethod = parametersBDDC->get("Krylov Method", 0);
+ test.printTimings(Preconditioner, Solver, krylovMethod);
+ Preconditioner->printTimings("BDDC_timers.dat");
+ }
+}
+
+int main (int argc, char *argv[]) {
+ // manage init and finalize of MPI and Kokkos
+ Tpetra::ScopeGuard tpetraScope(&argc, &argv);
+ const std::string fileNamePM = "problemMaker.xml";
+ const std::string fileNameBDDC = "bddc.xml";
+ const std::string fileNameMueLu = "";
+ const std::string fileNameNodalAMG = "";
+ const std::string meshDataFile = "";
+ RCP<Teuchos::ParameterList> parametersPM, parametersBDDC, parametersMueLu,
+ parametersNodalAMG;
+ MPI_Comm Comm = MPI_COMM_WORLD;
+ int myPID;
+ MPI_Comm_rank(Comm, &myPID); 
+ readInputFiles(Comm, fileNamePM, fileNameBDDC, fileNameMueLu, 
+ fileNameNodalAMG, parametersPM, parametersBDDC,
+ parametersMueLu, parametersNodalAMG);
+ int numIterations;
+ bool resetFile = true;
+ runTest(parametersPM, parametersBDDC, parametersMueLu, parametersNodalAMG, 
+ meshDataFile, myPID, Comm, resetFile, numIterations);
+}

packages/shylu/shylu_dd/bddc/example/problemMaker.xml

@@ -0,0 +1,39 @@
+<ParameterList name="BDDC">
+ <!-- Problem Type options: -->
+ <!-- Poisson-3D -->
+ <!-- Elasticity-3D -->
+ <!-- Poisson-2D -->
+ <!-- Elasticity-3D -->
+ <Parameter name="Problem Type String" type="string" value="Poisson-3D"/>
+ <Parameter name="Length Direction 1" type="double" value="1"/>
+ <Parameter name="Length Direction 2" type="double" value="1"/>
+ <Parameter name="Length Direction 3" type="double" value="1"/>
+ <Parameter name="Number Subregions Direction 1" type="int" value="2"/>
+ <Parameter name="Number Subregions Direction 2" type="int" value="2"/>
+ <Parameter name="Number Subregions Direction 3" type="int" value="2"/>
+ <Parameter name="Number Subdomains Per Subregion Direction 1" type="int" value="2"/>
+ <Parameter name="Number Subdomains Per Subregion Direction 2" type="int" value="2"/>
+ <Parameter name="Number Subdomains Per Subregion Direction 3" type="int" value="2"/>
+ <Parameter name="Number of Elements Per Subdomain Direction 1" type="int" value="4"/>
+ <Parameter name="Number of Elements Per Subdomain Direction 2" type="int" value="4"/>
+ <Parameter name="Number of Elements Per Subdomain Direction 3" type="int" value="4"/>
+ <Parameter name="Sort Subdomain Nodes" type="bool" value="true"/>
+ <Parameter name="Apply Left Side Essential BCs" type="bool" value="true"/>
+ <Parameter name="Apply Right Side Essential BCs" type="bool" value="false"/>
+ <!-- Diagonal Scale Factor > 1 removes singularities for floating domains -->
+ <Parameter name="Diagonal Scale Factor" type="double" value="1.0000"/>
+ <Parameter name="Output Model" type="bool" value="false"/>
+ <Parameter name="Load Direction" type="int" value="0"/>
+ <Parameter name="Material Property Option" type="int" value="0"/>
+ <Parameter name="E1" type="double" value="1.0"/>
+ <Parameter name="nu1" type="double" value="0.3"/>
+ <Parameter name="E2" type="double" value="1.0"/>
+ <Parameter name="nu2" type="double" value="0.3"/>
+ <Parameter name="Artificial Foundation Stiffness" type="double" value="0"/>
+ <Parameter name="omega" type="double" value="0"/>
+ <Parameter name="Adjust Nodal Coordinates" type="bool" value="true"/>
+ <Parameter name="Add Asymmetry" type="bool" value="false"/>
+ <Parameter name="Number of Solves" type="int" value="2"/>
+ <Parameter name="Quad Type 2 Point" type="int" value="1"/>
+ <Parameter name="Quad Type TP" type="int" value="1"/>
+</ParameterList>