diff --git a/packages/muelu/test/maxwell/Belos.xml b/packages/muelu/test/maxwell/Belos.xml new file mode 100644 index 000000000000..54bfc3273e69 --- /dev/null +++ b/packages/muelu/test/maxwell/Belos.xml @@ -0,0 +1,18 @@ + + + + + + + + + + + + + + + + + + diff --git a/packages/muelu/test/maxwell/CMakeLists.txt b/packages/muelu/test/maxwell/CMakeLists.txt index 51f83e006ad2..9f10108ecd1f 100644 --- a/packages/muelu/test/maxwell/CMakeLists.txt +++ b/packages/muelu/test/maxwell/CMakeLists.txt @@ -2,6 +2,7 @@ ASSERT_DEFINED( ${PACKAGE_NAME}_ENABLE_Amesos2 ${PACKAGE_NAME}_ENABLE_Belos + ${PACKAGE_NAME}_ENABLE_ML ) INCLUDE_DIRECTORIES(${CMAKE_CURRENT_SOURCE_DIR}) @@ -29,7 +30,7 @@ IF (${PACKAGE_NAME}_ENABLE_Belos AND ${PACKAGE_NAME}_ENABLE_Amesos2) NUM_MPI_PROCS 4 ) - IF (${PACKAGE_NAME}_ENABLE_ML AND ${PACKAGE_NAME}_ENABLE_EPETRA) + IF (${PACKAGE_NAME}_ENABLE_ML AND ${PACKAGE_NAME}_ENABLE_Epetra) TRIBITS_ADD_TEST( Maxwell3D NAME "Maxwell3D-Tpetra-ML-list" @@ -59,7 +60,7 @@ IF (${PACKAGE_NAME}_ENABLE_Belos AND ${PACKAGE_NAME}_ENABLE_Amesos2) NUM_MPI_PROCS 4 ) - IF (${PACKAGE_NAME}_ENABLE_ML AND ${PACKAGE_NAME}_HAVE_EPETRA_SOLVER_STACK AND (NOT Xpetra_INT_LONG_LONG)) + IF (${PACKAGE_NAME}_ENABLE_ML) TRIBITS_ADD_TEST( Maxwell3D NAME "Maxwell3D-ML" @@ -72,7 +73,7 @@ IF (${PACKAGE_NAME}_ENABLE_Belos AND ${PACKAGE_NAME}_ENABLE_Amesos2) TRIBITS_COPY_FILES_TO_BINARY_DIR(Maxwell_cp - SOURCE_FILES M0.mat M1.mat S.mat D0.mat coords.mat Maxwell.xml Maxwell_ML.xml Maxwell_ML_MueLu.xml Maxwell_ML1.xml + SOURCE_FILES M0.mat M1.mat S.mat D0.mat coords.mat Belos.xml Maxwell.xml Maxwell_ML.xml Maxwell_ML_MueLu.xml Maxwell_ML1.xml ) IF(HAVE_MUELU_COMPLEX) diff --git a/packages/muelu/test/maxwell/Maxwell3D.cpp b/packages/muelu/test/maxwell/Maxwell3D.cpp index 788ac12f26d1..9ffbf14e991a 100644 --- a/packages/muelu/test/maxwell/Maxwell3D.cpp +++ b/packages/muelu/test/maxwell/Maxwell3D.cpp @@ -44,6 +44,7 @@ // // @HEADER #include +#include // Teuchos #include @@ -67,6 +68,10 @@ #include #include +using Teuchos::RCP; +using Teuchos::rcp; +using Teuchos::TimeMonitor; + #ifdef HAVE_MUELU_TPETRA #include #endif @@ -78,7 +83,6 @@ #include #include #include // => This header defines Belos::XpetraOp -//#include // => This header defines Belos::MueLuOp #endif // Stratimikos @@ -101,7 +105,7 @@ // Helper functions for compilation purposes template -struct ML_Wrapper{ +struct EpetraSolvers_Wrapper{ static void Generate_ML_MaxwellPreconditioner(Teuchos::RCP >& SM, Teuchos::RCP >& D0, Teuchos::RCP >& Kn, @@ -128,7 +132,7 @@ struct ML_Wrapper{ template -struct ML_Wrapper { +struct EpetraSolvers_Wrapper { static void Generate_ML_MaxwellPreconditioner(Teuchos::RCP >& SM, Teuchos::RCP >& D0, Teuchos::RCP >& Kn, @@ -136,13 +140,12 @@ struct ML_Wrapper::coordinateType,int,GlobalOrdinal,Kokkos::Compat::KokkosSerialWrapperNode> >& coords, Teuchos::ParameterList & mueluList, Teuchos::RCP >& mlopX) { +#if defined(HAVE_MUELU_ML) and defined(HAVE_MUELU_EPETRA) typedef double SC; typedef int LO; typedef GlobalOrdinal GO; typedef Kokkos::Compat::KokkosSerialWrapperNode NO; typedef typename Teuchos::ScalarTraits::coordinateType coordinate_type; - using Teuchos::RCP; - using Teuchos::rcp; typedef typename Xpetra::Matrix Matrix; RCP epetraSM = Xpetra::Helpers::Op2EpetraCrs(SM); @@ -176,6 +179,10 @@ struct ML_Wrapper(mlop)); +#else + TEUCHOS_TEST_FOR_EXCEPTION(true, std::runtime_error, + "Need ML & Epetra support"); +#endif } static void Generate_ML_RefMaxwellPreconditioner(Teuchos::RCP >& SM, @@ -188,13 +195,12 @@ struct ML_Wrapper::coordinateType,int,GlobalOrdinal,Kokkos::Compat::KokkosSerialWrapperNode> >& coords, Teuchos::ParameterList & mueluList, Teuchos::RCP >& mlopX) { +#if defined(HAVE_MUELU_ML) and defined(HAVE_MUELU_EPETRA) typedef double SC; typedef int LO; typedef GlobalOrdinal GO; typedef Kokkos::Compat::KokkosSerialWrapperNode NO; typedef typename Teuchos::ScalarTraits::coordinateType coordinate_type; - using Teuchos::RCP; - using Teuchos::rcp; RCP epetraSM = Xpetra::Helpers::Op2EpetraCrs(SM); RCP epetraD0 = Xpetra::Helpers::Op2EpetraCrs(D0); @@ -233,258 +239,131 @@ struct ML_Wrapper(mlop)); +#else + TEUCHOS_TEST_FOR_EXCEPTION(true, std::runtime_error, + "Need ML & Epetra support"); +#endif } }; #endif -// Main wrappers struct +// Setup & solve wrappers struct // Because C++ doesn't support partial template specialization of functions. // By default, do not try to run Stratimikos, since that only works for Scalar=double. template -struct MainWrappers { - static int main_(Teuchos::CommandLineProcessor &clp, Xpetra::UnderlyingLib lib, int argc, char *argv[]); +struct SetupSolveWrappers { + static bool SetupSolve(std::map inputs); }; // Partial template specialization on SC=double // This code branch gives the option to run with Stratimikos. template -struct MainWrappers { - static int main_(Teuchos::CommandLineProcessor &clp, Xpetra::UnderlyingLib lib, int argc, char *argv[]); +struct SetupSolveWrappers { + static bool SetupSolve(std::map inputs); }; // By default, do not try to run Stratimikos, since that only works for Scalar=double. template -int MainWrappers::main_(Teuchos::CommandLineProcessor &clp, Xpetra::UnderlyingLib lib, int argc, char *argv[]) { +bool SetupSolveWrappers::SetupSolve(std::map inputs) { #include -#if defined(HAVE_MUELU_TPETRA) && defined(HAVE_MUELU_IFPACK2) - -#include - - using Teuchos::RCP; using Teuchos::rcp; - using Teuchos::TimeMonitor; - - bool success = false; - bool verbose = true; - try { - RCP< const Teuchos::Comm > comm = Teuchos::DefaultComm::getComm(); - - RCP out = Teuchos::fancyOStream(Teuchos::rcpFromRef(std::cout)); - out->setOutputToRootOnly(0); - - bool printTimings = true; clp.setOption("timings", "notimings", &printTimings, "print timings to screen"); - std::string timingsFormat = "table-fixed"; clp.setOption("time-format", &timingsFormat, "timings format (table-fixed | table-scientific | yaml)"); - double scaling = 1.0; clp.setOption("scaling", &scaling, "scale mass term"); - std::string solverName = "Belos"; clp.setOption("solverName", &solverName, "Name of iterative linear solver " - "to use for solving the linear system. " - "(\"Belos\")"); - std::string belosSolverType = "Block CG"; clp.setOption("belosSolverType", &belosSolverType, "Name of the Belos linear solver"); - std::string precType = "MueLu-RefMaxwell"; clp.setOption("precType", &precType, "preconditioner to use (MueLu-RefMaxwell|ML-RefMaxwell|none)"); - std::string xml; - if (!TYPE_EQUAL(SC, std::complex) && !TYPE_EQUAL(SC, std::complex)) - xml = "Maxwell.xml"; - else - xml = "Maxwell_complex.xml"; - clp.setOption("xml", &xml, "xml file with solver parameters"); - double tol = 1e-10; clp.setOption("tol", &tol, "solver convergence tolerance"); - int maxIts = 200; clp.setOption("its", &maxIts, "maximum number of solver iterations"); - bool use_stacked_timer = false; clp.setOption("stacked-timer", "no-stacked-timer", &use_stacked_timer, "use stacked timer"); - - std::string S_file, SM_file, M1_file, M0_file, M0inv_file, D0_file, coords_file, rhs_file="", nullspace_file="", material_file = "", Ms_file=""; - - if (!TYPE_EQUAL(SC, std::complex)) { - S_file = "S.mat"; - SM_file = ""; - M1_file = "M1.mat"; - M0_file = "M0.mat"; - M0inv_file = ""; - D0_file = "D0.mat"; - } else { - S_file = "S_complex.mat"; - SM_file = ""; - M1_file = "M1_complex.mat"; - M0_file = "M0_complex.mat"; - M0inv_file = ""; - D0_file = "D0_complex.mat"; - } - coords_file = "coords.mat"; - - clp.setOption("S", &S_file); - clp.setOption("SM", &SM_file); - clp.setOption("Ms", &Ms_file); - clp.setOption("M1", &M1_file); - clp.setOption("M0", &M0_file); - clp.setOption("M0inv", &M0inv_file); - clp.setOption("D0", &D0_file); - clp.setOption("coords", &coords_file); - clp.setOption("material", &material_file); - clp.setOption("nullspace", &nullspace_file); - clp.setOption("rhs", &rhs_file); - - clp.recogniseAllOptions(true); - switch (clp.parse(argc, argv)) { - case Teuchos::CommandLineProcessor::PARSE_HELP_PRINTED: return EXIT_SUCCESS; - case Teuchos::CommandLineProcessor::PARSE_ERROR: - case Teuchos::CommandLineProcessor::PARSE_UNRECOGNIZED_OPTION: return EXIT_FAILURE; - case Teuchos::CommandLineProcessor::PARSE_SUCCESSFUL: break; - } - - comm->barrier(); + typedef Xpetra::MultiVector::magnitudeType, LO, GO, NO> coordMV; - Teuchos::RCP stacked_timer; - if (use_stacked_timer) - stacked_timer = rcp(new Teuchos::StackedTimer("Maxwell Driver")); - Teuchos::TimeMonitor::setStackedTimer(stacked_timer); - - auto globalTimeMonitor = TimeMonitor::getNewTimer("Maxwell: S - Global Time"); - auto tm = TimeMonitor::getNewTimer("Maxwell: 1 - Read and Build Matrices"); - - // Read matrices in from files - // gradient matrix - RCP D0_Matrix = Xpetra::IO::Read(D0_file, lib, comm); - // maps for nodal and edge matrices - RCP node_map = D0_Matrix->getDomainMap(); - RCP edge_map = D0_Matrix->getRangeMap(); - // edge mass matrix - RCP M1_Matrix = Xpetra::IO::Read(M1_file, edge_map); - RCP Ms_Matrix; - if (Ms_file != "") - Ms_Matrix = Xpetra::IO::Read(Ms_file, edge_map); - else - Ms_Matrix = M1_Matrix; - // build stiffness plus mass matrix (SM_Matrix) - RCP SM_Matrix; - if (SM_file == "") { - // edge stiffness matrix - RCP S_Matrix = Xpetra::IO::Read(S_file, edge_map); - Xpetra::MatrixMatrix::TwoMatrixAdd(*S_Matrix,false,(SC)1.0,*M1_Matrix,false,scaling,SM_Matrix,*out); - SM_Matrix->fillComplete(); - } else { - SM_Matrix = Xpetra::IO::Read(SM_file, edge_map); - } - RCP M0inv_Matrix; - if (M0inv_file == "") { - // nodal mass matrix - RCP M0_Matrix = Xpetra::IO::Read(M0_file, node_map); - // build lumped mass matrix inverse (M0inv_Matrix) - RCP diag = Utilities::GetLumpedMatrixDiagonal(M0_Matrix); - RCP M0inv_MatrixWrap = Teuchos::rcp(new CrsMatrixWrap(node_map, node_map, 0)); - RCP M0inv_CrsMatrix = M0inv_MatrixWrap->getCrsMatrix(); - Teuchos::ArrayRCP rowPtr; - Teuchos::ArrayRCP colInd; - Teuchos::ArrayRCP values; - Teuchos::ArrayRCP diags = diag->getData(0); - size_t nodeNumElements = node_map->getNodeNumElements(); - M0inv_CrsMatrix->allocateAllValues(nodeNumElements, rowPtr, colInd, values); - SC ONE = (SC)1.0; - for (size_t i = 0; i < nodeNumElements; i++) { - rowPtr[i] = i; colInd[i] = i; values[i] = ONE / diags[i]; - } - rowPtr[nodeNumElements] = nodeNumElements; - M0inv_CrsMatrix->setAllValues(rowPtr, colInd, values); - M0inv_CrsMatrix->expertStaticFillComplete(node_map, node_map); - M0inv_Matrix = Teuchos::rcp_dynamic_cast(M0inv_MatrixWrap); - } else if (M0inv_file == "none") { - // pass - } else { - M0inv_Matrix = Xpetra::IO::Read(M0inv_file, node_map); - } - // coordinates - RCP::coordinateType, LO, GO, NO> > coords = Xpetra::IO::coordinateType, LO, GO, NO>::ReadMultiVector(coords_file, node_map); + RCP SM_Matrix = *static_cast*>(inputs["SM"]); + RCP D0_Matrix = *static_cast*>(inputs["D0"]); + RCP M1_Matrix = *static_cast*>(inputs["M1"]); + RCP Ms_Matrix = *static_cast*>(inputs["Ms"]); + RCP M0inv_Matrix = *static_cast*>(inputs["M0inv"]); + RCP Kn_Matrix = *static_cast*>(inputs["Kn"]); - RCP nullspace = Teuchos::null; - if (nullspace_file != "") - nullspace = Xpetra::IO::ReadMultiVector(nullspace_file, edge_map); + RCP coords = *static_cast*>(inputs["coordinates"]); + RCP nullspace = *static_cast*>(inputs["nullspace"]); + RCP material = *static_cast*>(inputs["material"]); - RCP material = Teuchos::null; - if (material_file != "") { - material = Xpetra::IO::ReadMultiVector(material_file, node_map); - } + RCP B = *static_cast*>(inputs["B"]); + RCP X = *static_cast*>(inputs["X"]); + RCP X0 = *static_cast*>(inputs["X0"]); - // set parameters - Teuchos::ParameterList params; - Teuchos::updateParametersFromXmlFileAndBroadcast(xml,Teuchos::Ptr(¶ms),*comm); - - // setup LHS, RHS - RCP B; - if (rhs_file == "") { - B = MultiVectorFactory::Build(edge_map,1); - RCP vec = MultiVectorFactory::Build(edge_map,1); - vec -> putScalar((SC)1.0); - SM_Matrix->apply(*vec,*B); - } else - B = Xpetra::IO::ReadMultiVector(rhs_file, edge_map); + Teuchos::ParameterList params = *static_cast(inputs["params"]); + Teuchos::ParameterList belosParams = *static_cast(inputs["belosParams"]); + std::string solverName = *static_cast(inputs["solverName"]); + std::string belosSolverType = *static_cast(inputs["belosSolverType"]); + std::string precType = *static_cast(inputs["precType"]); + int numResolves = *static_cast(inputs["numResolves"]); - RCP X = MultiVectorFactory::Build(edge_map,1); - X -> putScalar((SC)0.0); - - comm->barrier(); - tm = Teuchos::null; - - if (solverName == "Belos") { - auto tm2 = TimeMonitor::getNewTimer("Maxwell: 2 - Build Belos solver etc"); - - // construct preconditioner - RCP > preconditioner; - if (precType=="MueLu-RefMaxwell") - preconditioner - = rcp( new MueLu::RefMaxwell(SM_Matrix,D0_Matrix,Ms_Matrix,M0inv_Matrix, - M1_Matrix,nullspace,coords,params) ); + Xpetra::UnderlyingLib lib = *static_cast(inputs["lib"]); + RCP > comm = *static_cast >*>(inputs["comm"]); + RCP out = *static_cast*>(inputs["out"]); + bool success = false; + auto tm2 = TimeMonitor::getNewTimer("Maxwell: 2 - Build solver and preconditioner"); +#ifdef HAVE_MUELU_BELOS + if (solverName == "Belos") { + // construct preconditioner + RCP preconditioner; + if (precType=="MueLu-RefMaxwell") { + preconditioner = rcp( new MueLu::RefMaxwell(SM_Matrix,D0_Matrix,Ms_Matrix,M0inv_Matrix, + M1_Matrix,nullspace,coords,params) ); #ifdef HAVE_MUELU_TPETRA { // A test to make sure we can wrap this guy as a MueLu::TpetraOperator RCP precOp = Teuchos::rcp_dynamic_cast(preconditioner); MueLu::TpetraOperator OpT(precOp); } +#endif // HAVE_MUELU_TPETRA + } + else + *out << "Preconditioner not supported\n"; + + +#ifdef HAVE_MUELU_TPETRA + { + // A test to make sure we can wrap this guy as a MueLu::TpetraOperator + RCP precOp = Teuchos::rcp_dynamic_cast(preconditioner); + MueLu::TpetraOperator OpT(precOp); + } #endif - // Belos linear problem -#ifdef HAVE_MUELU_BELOS - typedef MultiVector MV; - typedef Belos::OperatorT OP; - Teuchos::RCP belosOp = Teuchos::rcp(new Belos::XpetraOp(SM_Matrix)); // Turns a Xpetra::Matrix object into a Belos operator - - RCP > problem = rcp( new Belos::LinearProblem() ); - problem -> setOperator( belosOp ); - Teuchos::RCP belosPrecOp; - if (precType!="none") { - belosPrecOp = Teuchos::rcp(new Belos::XpetraOp(preconditioner)); // Turns a Xpetra::Matrix object into a Belos operator - problem -> setRightPrec( belosPrecOp ); - } + // Belos linear problem + typedef MultiVector MV; + typedef Belos::OperatorT OP; + Teuchos::RCP belosOp = Teuchos::rcp(new Belos::XpetraOp(SM_Matrix)); // Turns a Xpetra::Matrix object into a Belos operator + + RCP > problem = rcp( new Belos::LinearProblem() ); + problem -> setOperator( belosOp ); + Teuchos::RCP belosPrecOp; + if (precType != "none") { + belosPrecOp = Teuchos::rcp(new Belos::XpetraOp(preconditioner)); // Turns a Xpetra::Matrix object into a Belos operator + problem -> setRightPrec( belosPrecOp ); + } + problem -> setProblem( X, B ); - problem -> setProblem( X, B ); + bool set = problem->setProblem(); + if (set == false) { + *out << "\nERROR: Belos::LinearProblem failed to set up correctly!" << std::endl; + return false; + } - bool set = problem->setProblem(); - if (set == false) { - *out << "\nERROR: Belos::LinearProblem failed to set up correctly!" << std::endl; - return EXIT_FAILURE; - } + // Belos solver + RCP< Belos::SolverManager > solver; + RCP< Belos::SolverFactory > factory = rcp( new Belos::SolverFactory() ); + solver = factory->create(belosSolverType,Teuchos::rcpFromRef(belosParams.sublist(belosSolverType))); - // Belos solver - RCP< Belos::SolverManager > solver; - RCP< Belos::SolverFactory > factory = rcp( new Belos::SolverFactory() ); - RCP belosParams - = rcp( new Teuchos::ParameterList() ); - belosParams->set("Maximum Iterations", maxIts); - belosParams->set("Convergence Tolerance",tol); - belosParams->set("Verbosity", Belos::Errors + Belos::Warnings + Belos::StatusTestDetails); - belosParams->set("Output Frequency",1); - belosParams->set("Output Style",Belos::Brief); - belosParams->set("Implicit Residual Scaling","None"); - belosParams->set("Explicit Residual Test",true); - belosParams->set("Explicit Residual Scaling","None"); - solver = factory->create(belosSolverType,belosParams); - - comm->barrier(); - tm2=Teuchos::null; - auto tm3 = TimeMonitor::getNewTimer("Maxwell: 3 - Solve"); - - // set problem and solve - solver -> setProblem( problem ); + comm->barrier(); + tm2 = Teuchos::null; + + auto tm3 = TimeMonitor::getNewTimer("Maxwell: 3 - Solve"); + + // set problem and solve + solver -> setProblem( problem ); + for(int solveno = 0; solveno<=numResolves; solveno++) { + if (X0.is_null()) + X->putScalar(Teuchos::ScalarTraits::zero()); + else + X = X0; Belos::ReturnType status = solver -> solve(); int iters = solver -> getNumIters(); success = (iters<50 && status == Belos::Converged); @@ -492,298 +371,115 @@ int MainWrappers::main_(Teuchos::Command *out << "SUCCESS! Belos converged in " << iters << " iterations." << std::endl; else *out << "FAILURE! Belos did not converge fast enough." << std::endl; - tm3 = Teuchos::null; - } - comm->barrier(); - globalTimeMonitor = Teuchos::null; - - if (printTimings) { - RCP reportParams = rcp(new Teuchos::ParameterList); - if (timingsFormat == "yaml") { - reportParams->set("Report format", "YAML"); // "Table" or "YAML" - reportParams->set("YAML style", "compact"); // "spacious" or "compact" - } - reportParams->set("How to merge timer sets", "Union"); - reportParams->set("alwaysWriteLocal", false); - reportParams->set("writeGlobalStats", true); - reportParams->set("writeZeroTimers", false); - // FIXME: no "ignoreZeroTimers" - - const std::string filter = ""; - - std::ios_base::fmtflags ff(out->flags()); - if (timingsFormat == "table-fixed") *out << std::fixed; - else * out << std::scientific; - if (use_stacked_timer) { - stacked_timer->stop("Maxwell Driver"); - Teuchos::StackedTimer::OutputOptions options; - options.output_fraction = options.output_histogram = options.output_minmax = true; - stacked_timer->report(*out, comm, options); - } else - TimeMonitor::report(comm.ptr(), *out, filter, reportParams); - *out << std::setiosflags(ff); } - - TimeMonitor::clearCounters(); -#endif // #ifdef HAVE_MUELU_BELOS } - TEUCHOS_STANDARD_CATCH_STATEMENTS(verbose, std::cerr, success); +#endif // HAVE_MUELU_BELOS + comm->barrier(); - return ( success ? EXIT_SUCCESS : EXIT_FAILURE ); -#else - return EXIT_SUCCESS; -#endif -} // main + return success; +} // SetupSolve // This code branch gives the option to run with Stratimikos. template -int MainWrappers::main_(Teuchos::CommandLineProcessor &clp, Xpetra::UnderlyingLib lib, int argc, char *argv[]) { +bool SetupSolveWrappers::SetupSolve(std::map inputs) { typedef double Scalar; #include -#if defined(HAVE_MUELU_TPETRA) && defined(HAVE_MUELU_IFPACK2) + typedef Xpetra::MultiVector::magnitudeType, LO, GO, NO> coordMV; -#include + RCP SM_Matrix = *static_cast*>(inputs["SM"]); + RCP D0_Matrix = *static_cast*>(inputs["D0"]); + RCP M1_Matrix = *static_cast*>(inputs["M1"]); + RCP Ms_Matrix = *static_cast*>(inputs["Ms"]); + RCP M0inv_Matrix = *static_cast*>(inputs["M0inv"]); + RCP Kn_Matrix = *static_cast*>(inputs["Kn"]); - using Teuchos::RCP; using Teuchos::rcp; - using Teuchos::TimeMonitor; + RCP coords = *static_cast*>(inputs["coordinates"]); + RCP nullspace = *static_cast*>(inputs["nullspace"]); + RCP material = *static_cast*>(inputs["material"]); - bool success = false; - bool verbose = true; - try { - RCP< const Teuchos::Comm > comm = Teuchos::DefaultComm::getComm(); - - RCP out = Teuchos::fancyOStream(Teuchos::rcpFromRef(std::cout)); - out->setOutputToRootOnly(0); - - bool printTimings = true; clp.setOption("timings", "notimings", &printTimings, "print timings to screen"); - std::string timingsFormat = "table-fixed"; clp.setOption("time-format", &timingsFormat, "timings format (table-fixed | table-scientific | yaml)"); - double scaling = 1.0; clp.setOption("scaling", &scaling, "scale mass term"); - std::string solverName = "Belos"; clp.setOption("solverName", &solverName, "Name of iterative linear solver " - "to use for solving the linear system. " - "(\"Belos\" or \"Stratimikos\")"); - std::string belosSolverType = "Block CG"; clp.setOption("belosSolverType", &belosSolverType, "Name of the Belos linear solver"); - std::string precType = "MueLu-RefMaxwell"; clp.setOption("precType", &precType, "preconditioner to use (MueLu-RefMaxwell|ML-RefMaxwell|none)"); - std::string xml = ""; clp.setOption("xml", &xml, "xml file with solver parameters (default: \"Maxwell.xml\")"); - double tol = 1e-10; clp.setOption("tol", &tol, "solver convergence tolerance"); - int maxIts = 200; clp.setOption("its", &maxIts, "maximum number of solver iterations"); - bool use_stacked_timer = false; clp.setOption("stacked-timer", "no-stacked-timer", &use_stacked_timer, "use stacked timer"); - - - std::string S_file, SM_file, M1_file, M0_file, M0inv_file, D0_file, coords_file, rhs_file="", nullspace_file="", material_file = "", Ms_file="", sol_file="", Kn_file=""; - - S_file = "S.mat"; - SM_file = ""; - M1_file = "M1.mat"; - M0_file = "M0.mat"; - M0inv_file = ""; - D0_file = "D0.mat"; - coords_file = "coords.mat"; - - clp.setOption("S", &S_file); - clp.setOption("SM", &SM_file); - clp.setOption("Ms", &Ms_file); - clp.setOption("M1", &M1_file); - clp.setOption("M0", &M0_file); - clp.setOption("Kn", &Kn_file); - clp.setOption("M0inv", &M0inv_file); - clp.setOption("D0", &D0_file); - clp.setOption("coords", &coords_file); - clp.setOption("rhs", &rhs_file); - clp.setOption("nullspace", &nullspace_file); - clp.setOption("material", &material_file); - - clp.recogniseAllOptions(true); - switch (clp.parse(argc, argv)) { - case Teuchos::CommandLineProcessor::PARSE_HELP_PRINTED: return EXIT_SUCCESS; - case Teuchos::CommandLineProcessor::PARSE_ERROR: - case Teuchos::CommandLineProcessor::PARSE_UNRECOGNIZED_OPTION: return EXIT_FAILURE; - case Teuchos::CommandLineProcessor::PARSE_SUCCESSFUL: break; - } + RCP B = *static_cast*>(inputs["B"]); + RCP X = *static_cast*>(inputs["X"]); + RCP X0 = *static_cast*>(inputs["X0"]); - if (xml == ""){ - if (precType == "MueLu-RefMaxwell") - xml = "Maxwell.xml"; - else if (precType == "ML-RefMaxwell") - xml = "Maxwell_ML.xml"; - else if (precType == "ML-Maxwell") - xml = "Maxwell_ML1.xml"; - } + Teuchos::ParameterList params = *static_cast(inputs["params"]); + Teuchos::ParameterList belosParams = *static_cast(inputs["belosParams"]); + std::string solverName = *static_cast(inputs["solverName"]); + std::string belosSolverType = *static_cast(inputs["belosSolverType"]); + std::string precType = *static_cast(inputs["precType"]); + int numResolves = *static_cast(inputs["numResolves"]); - comm->barrier(); + Xpetra::UnderlyingLib lib = *static_cast(inputs["lib"]); + RCP > comm = *static_cast >*>(inputs["comm"]); + RCP out = *static_cast*>(inputs["out"]); - Teuchos::RCP stacked_timer; - if (use_stacked_timer) - stacked_timer = rcp(new Teuchos::StackedTimer("Maxwell Driver")); - Teuchos::TimeMonitor::setStackedTimer(stacked_timer); - - auto globalTimeMonitor = TimeMonitor::getNewTimer("Maxwell: S - Global Time"); - auto tm = TimeMonitor::getNewTimer("Maxwell: 1 - Read and Build Matrices"); - - // Read matrices in from files - // gradient matrix - RCP D0_Matrix = Xpetra::IO::Read(D0_file, lib, comm); - // maps for nodal and edge matrices - RCP node_map = D0_Matrix->getDomainMap(); - RCP edge_map = D0_Matrix->getRangeMap(); - // edge mass matrix - RCP M1_Matrix = Xpetra::IO::Read(M1_file, edge_map); - RCP Ms_Matrix; - if (Ms_file != "") - Ms_Matrix = Xpetra::IO::Read(Ms_file, edge_map); - else - Ms_Matrix = M1_Matrix; - // build stiffness plus mass matrix (SM_Matrix) - RCP SM_Matrix, S_Matrix; - if (SM_file == "") { - // edge stiffness matrix - S_Matrix = Xpetra::IO::Read(S_file, edge_map); - Xpetra::MatrixMatrix::TwoMatrixAdd(*S_Matrix,false,(SC)1.0,*M1_Matrix,false,scaling,SM_Matrix,*out); - SM_Matrix->fillComplete(); - } else { - SM_Matrix = Xpetra::IO::Read(SM_file, edge_map); - } - RCP M0inv_Matrix, M0_Matrix; - if (M0inv_file == "") { - // nodal mass matrix - M0_Matrix = Xpetra::IO::Read(M0_file, node_map); - // build lumped mass matrix inverse (M0inv_Matrix) - RCP diag = Utilities::GetLumpedMatrixDiagonal(M0_Matrix); - RCP M0inv_MatrixWrap = Teuchos::rcp(new CrsMatrixWrap(node_map, node_map, 0)); - RCP M0inv_CrsMatrix = M0inv_MatrixWrap->getCrsMatrix(); - Teuchos::ArrayRCP rowPtr; - Teuchos::ArrayRCP colInd; - Teuchos::ArrayRCP values; - Teuchos::ArrayRCP diags = diag->getData(0); - size_t nodeNumElements = node_map->getNodeNumElements(); - M0inv_CrsMatrix->allocateAllValues(nodeNumElements, rowPtr, colInd, values); - SC ONE = (SC)1.0; - for (size_t i = 0; i < nodeNumElements; i++) { - rowPtr[i] = i; colInd[i] = i; values[i] = ONE / diags[i]; - } - rowPtr[nodeNumElements] = nodeNumElements; - M0inv_CrsMatrix->setAllValues(rowPtr, colInd, values); - M0inv_CrsMatrix->expertStaticFillComplete(node_map, node_map); - M0inv_Matrix = Teuchos::rcp_dynamic_cast(M0inv_MatrixWrap); - } else if (M0inv_file == "none") { - // pass - } else { - M0inv_Matrix = Xpetra::IO::Read(M0inv_file, node_map); + bool success = false; + + auto tm2 = TimeMonitor::getNewTimer("Maxwell: 2 - Build solver and preconditioner"); +#ifdef HAVE_MUELU_BELOS + if (solverName == "Belos") { + // construct preconditioner + RCP preconditioner; + if (precType=="MueLu-RefMaxwell") { + preconditioner = rcp( new MueLu::RefMaxwell(SM_Matrix,D0_Matrix,Ms_Matrix,M0inv_Matrix, + M1_Matrix,nullspace,coords,params) ); } - // coordinates - RCP::magnitudeType, LO, GO, NO> > coords = Xpetra::IO::magnitudeType, LO, GO, NO>::ReadMultiVector(coords_file, node_map); - - RCP nullspace = Teuchos::null; - if (nullspace_file != "") - nullspace = Xpetra::IO::ReadMultiVector(nullspace_file, edge_map); - - RCP material = Teuchos::null; - if (material_file != "") { - material = Xpetra::IO::ReadMultiVector(material_file, node_map); + else if (precType=="ML-RefMaxwell") { + TEUCHOS_ASSERT(lib==Xpetra::UseEpetra); + EpetraSolvers_Wrapper::Generate_ML_RefMaxwellPreconditioner(SM_Matrix,D0_Matrix,Ms_Matrix,M0inv_Matrix, + M1_Matrix,nullspace,material,coords,params,preconditioner); + } else if (precType=="ML-Maxwell") { + TEUCHOS_ASSERT(lib==Xpetra::UseEpetra); + EpetraSolvers_Wrapper::Generate_ML_MaxwellPreconditioner(SM_Matrix,D0_Matrix,Kn_Matrix, + nullspace,coords,params,preconditioner); } - // set parameters - Teuchos::ParameterList params; - Teuchos::updateParametersFromXmlFileAndBroadcast(xml,Teuchos::Ptr(¶ms),*comm); - - // setup LHS, RHS - RCP B; - if (rhs_file == "") { - B = MultiVectorFactory::Build(edge_map,1); - RCP vec = MultiVectorFactory::Build(edge_map,1); - vec -> putScalar((SC)1.0); - SM_Matrix->apply(*vec,*B); - } else - B = Xpetra::IO::ReadMultiVector(rhs_file, edge_map); - RCP X = MultiVectorFactory::Build(edge_map,1); - X -> putScalar((SC)0.0); - - comm->barrier(); - tm = Teuchos::null; - - if (solverName == "Belos") { - auto tm2 = TimeMonitor::getNewTimer("Maxwell: 2 - Build Belos solver etc"); - - // construct preconditioner - RCP preconditioner; - if (precType=="MueLu-RefMaxwell") - preconditioner = rcp( new MueLu::RefMaxwell(SM_Matrix,D0_Matrix,Ms_Matrix,M0inv_Matrix, - M1_Matrix,nullspace,coords,params) ); - else if (precType=="ML-RefMaxwell") { -#if defined(HAVE_MUELU_ML) and defined(HAVE_MUELU_EPETRA) - TEUCHOS_ASSERT(lib==Xpetra::UseEpetra); - ML_Wrapper::Generate_ML_RefMaxwellPreconditioner(SM_Matrix,D0_Matrix,Ms_Matrix,M0inv_Matrix, - M1_Matrix,nullspace,material,coords,params,preconditioner); -#else - throw std::runtime_error("ML not available."); -#endif - } else if (precType=="ML-Maxwell") { -#if defined(HAVE_MUELU_ML) and defined(HAVE_MUELU_EPETRA) - TEUCHOS_ASSERT(lib==Xpetra::UseEpetra); - RCP Kn_Matrix; - if (Kn_file!="") - Kn_Matrix = Xpetra::IO::Read(Kn_file, node_map); - ML_Wrapper::Generate_ML_MaxwellPreconditioner(SM_Matrix,D0_Matrix,Kn_Matrix, - nullspace,coords,params,preconditioner); -#else - throw std::runtime_error("ML not available."); -#endif - } - #ifdef HAVE_MUELU_TPETRA - { - // A test to make sure we can wrap this guy as a MueLu::TpetraOperator - RCP precOp = Teuchos::rcp_dynamic_cast(preconditioner); - MueLu::TpetraOperator OpT(precOp); - } + { + // A test to make sure we can wrap this guy as a MueLu::TpetraOperator + RCP precOp = Teuchos::rcp_dynamic_cast(preconditioner); + MueLu::TpetraOperator OpT(precOp); + } #endif + // Belos linear problem + typedef MultiVector MV; + typedef Belos::OperatorT OP; + RCP belosOp = Teuchos::rcp(new Belos::XpetraOp(SM_Matrix)); // Turns a Xpetra::Matrix object into a Belos operator + + RCP > problem = rcp( new Belos::LinearProblem() ); + problem -> setOperator( belosOp ); + RCP belosPrecOp; + if (precType != "none") { + belosPrecOp = Teuchos::rcp(new Belos::XpetraOp(preconditioner)); // Turns a Xpetra::Matrix object into a Belos operator + problem -> setRightPrec( belosPrecOp ); + } + problem -> setProblem( X, B ); - // Belos linear problem -#ifdef HAVE_MUELU_BELOS - typedef MultiVector MV; - typedef Belos::OperatorT OP; - Teuchos::RCP belosOp = Teuchos::rcp(new Belos::XpetraOp(SM_Matrix)); // Turns a Xpetra::Matrix object into a Belos operator - - RCP > problem = rcp( new Belos::LinearProblem() ); - problem -> setOperator( belosOp ); - Teuchos::RCP belosPrecOp; - if (precType != "none") { - belosPrecOp = Teuchos::rcp(new Belos::XpetraOp(preconditioner)); // Turns a Xpetra::Matrix object into a Belos operator - problem -> setRightPrec( belosPrecOp ); - } - - problem -> setProblem( X, B ); - - bool set = problem->setProblem(); - if (set == false) { - *out << "\nERROR: Belos::LinearProblem failed to set up correctly!" << std::endl; - return EXIT_FAILURE; - } - - // Belos solver - RCP< Belos::SolverManager > solver; - RCP< Belos::SolverFactory > factory = rcp( new Belos::SolverFactory() ); - RCP belosParams - = rcp( new Teuchos::ParameterList() ); - belosParams->set("Maximum Iterations", maxIts); - belosParams->set("Convergence Tolerance",tol); - belosParams->set("Verbosity", Belos::Errors + Belos::Warnings + Belos::StatusTestDetails); - belosParams->set("Output Frequency",1); - belosParams->set("Output Style",Belos::Brief); - belosParams->set("Implicit Residual Scaling","None"); - belosParams->set("Explicit Residual Test",true); - belosParams->set("Explicit Residual Scaling","None"); + bool set = problem->setProblem(); + if (set == false) { + *out << "\nERROR: Belos::LinearProblem failed to set up correctly!" << std::endl; + return false; + } - solver = factory->create(belosSolverType,belosParams); + // Belos solver + RCP< Belos::SolverManager > solver; + RCP< Belos::SolverFactory > factory = rcp( new Belos::SolverFactory() ); + solver = factory->create(belosSolverType,Teuchos::rcpFromRef(belosParams.sublist(belosSolverType))); - comm->barrier(); + comm->barrier(); + tm2 = Teuchos::null; - auto tm3 = TimeMonitor::getNewTimer("Maxwell: 3 - Solve"); + auto tm3 = TimeMonitor::getNewTimer("Maxwell: 3 - Solve"); - // set problem and solve - solver -> setProblem( problem ); + // set problem and solve + solver -> setProblem( problem ); + for(int solveno = 0; solveno<=numResolves; solveno++) { + if (X0.is_null()) + X->putScalar(Teuchos::ScalarTraits::zero()); + else + X = X0; Belos::ReturnType status = solver -> solve(); int iters = solver -> getNumIters(); success = (iters<50 && status == Belos::Converged); @@ -791,69 +487,271 @@ int MainWrappers::main_(Teuchos::Command *out << "SUCCESS! Belos converged in " << iters << " iterations." << std::endl; else *out << "FAILURE! Belos did not converge fast enough." << std::endl; - tm3 = Teuchos::null; } + } +#endif // HAVE_MUELU_BELOS #ifdef HAVE_MUELU_STRATIMIKOS - if (solverName == "Stratimikos") { - auto tm4 = TimeMonitor::getNewTimer("Maxwell: 2 - Build Stratimikos solver"); - - // Build the rest of the Stratimikos list - Teuchos::ParameterList SList; - SList.set("Linear Solver Type","Belos"); - SList.sublist("Linear Solver Types").sublist("Belos").set("Solver Type", belosSolverType); - SList.sublist("Linear Solver Types").sublist("Belos").sublist("Solver Types").sublist(belosSolverType).set("Output Frequency",1); - SList.sublist("Linear Solver Types").sublist("Belos").sublist("Solver Types").sublist(belosSolverType).set("Maximum Iterations",maxIts); - SList.sublist("Linear Solver Types").sublist("Belos").sublist("Solver Types").sublist(belosSolverType).set("Convergence Tolerance",tol); - SList.sublist("Linear Solver Types").sublist("Belos").sublist("Solver Types").sublist(belosSolverType).set("Output Style",1); - SList.sublist("Linear Solver Types").sublist("Belos").sublist("Solver Types").sublist(belosSolverType).set("Verbosity",33); - SList.sublist("Linear Solver Types").sublist("Belos").sublist("Solver Types").sublist(belosSolverType).set("Implicit Residual Scaling","None"); - SList.sublist("Linear Solver Types").sublist("Belos").sublist("Solver Types").sublist(belosSolverType).set("Explicit Residual Test",true); - SList.sublist("Linear Solver Types").sublist("Belos").sublist("Solver Types").sublist(belosSolverType).set("Explicit Residual Scaling","None"); - SList.sublist("Linear Solver Types").sublist("Belos").sublist("VerboseObject").set("Verbosity Level", "medium"); - SList.set("Preconditioner Type","MueLuRefMaxwell"); - params.set("parameterlist: syntax","muelu"); - SList.sublist("Preconditioner Types").set("MueLuRefMaxwell",params); - // Add matrices to parameterlist - SList.sublist("Preconditioner Types").sublist("MueLuRefMaxwell").set("D0",D0_Matrix); - SList.sublist("Preconditioner Types").sublist("MueLuRefMaxwell").set("M0inv",M0inv_Matrix); - SList.sublist("Preconditioner Types").sublist("MueLuRefMaxwell").set("M1",M1_Matrix); - SList.sublist("Preconditioner Types").sublist("MueLuRefMaxwell").set("Ms",Ms_Matrix); - SList.sublist("Preconditioner Types").sublist("MueLuRefMaxwell").set("Coordinates",coords); - - // Build Thyra linear algebra objects - RCP > thyraA = Xpetra::ThyraUtils::toThyra(Teuchos::rcp_dynamic_cast(SM_Matrix)->getCrsMatrix()); - RCP< Thyra::VectorBase >thyraX = Teuchos::rcp_const_cast >(Xpetra::ThyraUtils::toThyraVector(X->getVectorNonConst(0))); - // TODO: Why do we loose a reference when running this with Epetra? - RCP >thyraB = Xpetra::ThyraUtils::toThyraVector(B->getVector(0)); - - // Build Stratimikos solver - Stratimikos::DefaultLinearSolverBuilder linearSolverBuilder; // This is the Stratimikos main class (= factory of solver factory). - Stratimikos::enableMueLuRefMaxwell(linearSolverBuilder); // Register MueLu as a Stratimikos preconditioner strategy. - linearSolverBuilder.setParameterList(rcp(&SList,false)); // Setup solver parameters using a Stratimikos parameter list. - - // Build a new "solver factory" according to the previously specified parameter list. - RCP > solverFactory = Thyra::createLinearSolveStrategy(linearSolverBuilder); - - // Build a Thyra operator corresponding to A^{-1} computed using the Stratimikos solver. - Teuchos::RCP > thyraInverseA = Thyra::linearOpWithSolve(*solverFactory, thyraA); - - comm->barrier(); - - tm4 = Teuchos::null; - auto tm5 = TimeMonitor::getNewTimer("Maxwell: 3 - Solve"); - - // Solve Ax = b. - Thyra::SolveStatus status = Thyra::solve(*thyraInverseA, Thyra::NOTRANS, *thyraB, thyraX.ptr()); - std::cout << status << std::endl; - - success = (status.solveStatus == Thyra::SOLVE_STATUS_CONVERGED); - tm5 = Teuchos::null; - } -#endif + if (solverName == "Stratimikos") { + // Build the rest of the Stratimikos list + Teuchos::ParameterList stratimikosParams; + stratimikosParams.set("Linear Solver Type","Belos"); + stratimikosParams.sublist("Linear Solver Types").sublist("Belos").set("Solver Type", belosSolverType); + stratimikosParams.sublist("Linear Solver Types").sublist("Belos").set("Solver Types", belosParams); + stratimikosParams.sublist("Linear Solver Types").sublist("Belos").sublist("VerboseObject").set("Verbosity Level", "medium"); + stratimikosParams.set("Preconditioner Type","MueLuRefMaxwell"); + params.set("parameterlist: syntax","muelu"); + stratimikosParams.sublist("Preconditioner Types").set("MueLuRefMaxwell",params); + // Add matrices to parameterlist + stratimikosParams.sublist("Preconditioner Types").sublist("MueLuRefMaxwell").set("D0", D0_Matrix); + stratimikosParams.sublist("Preconditioner Types").sublist("MueLuRefMaxwell").set("M0inv", M0inv_Matrix); + stratimikosParams.sublist("Preconditioner Types").sublist("MueLuRefMaxwell").set("M1", M1_Matrix); + stratimikosParams.sublist("Preconditioner Types").sublist("MueLuRefMaxwell").set("Ms", Ms_Matrix); + stratimikosParams.sublist("Preconditioner Types").sublist("MueLuRefMaxwell").set("Coordinates",coords); + + // Build Thyra linear algebra objects + RCP > thyraA = Xpetra::ThyraUtils::toThyra(Teuchos::rcp_dynamic_cast(SM_Matrix)->getCrsMatrix()); + RCP< Thyra::VectorBase >thyraX = Teuchos::rcp_const_cast >(Xpetra::ThyraUtils::toThyraVector(X->getVectorNonConst(0))); + // TODO: Why do we loose a reference when running this with Epetra? + RCP >thyraB = Xpetra::ThyraUtils::toThyraVector(B->getVector(0)); + + // Build Stratimikos solver + Stratimikos::DefaultLinearSolverBuilder linearSolverBuilder; // This is the Stratimikos main class (= factory of solver factory). + Stratimikos::enableMueLuRefMaxwell(linearSolverBuilder); // Register MueLu as a Stratimikos preconditioner strategy. + linearSolverBuilder.setParameterList(rcp(&stratimikosParams,false)); // Setup solver parameters using a Stratimikos parameter list. + + // Build a new "solver factory" according to the previously specified parameter list. + RCP > solverFactory = Thyra::createLinearSolveStrategy(linearSolverBuilder); + + // Build a Thyra operator corresponding to A^{-1} computed using the Stratimikos solver. + Teuchos::RCP > thyraInverseA = Thyra::linearOpWithSolve(*solverFactory, thyraA); + comm->barrier(); - globalTimeMonitor = Teuchos::null; + tm2 = Teuchos::null; + + auto tm5 = TimeMonitor::getNewTimer("Maxwell: 3 - Solve"); + + // Solve Ax = b. + Thyra::SolveStatus status = Thyra::solve(*thyraInverseA, Thyra::NOTRANS, *thyraB, thyraX.ptr()); + std::cout << status << std::endl; + + success = (status.solveStatus == Thyra::SOLVE_STATUS_CONVERGED); + } // HAVE_MUELU_STRATIMIKOS +#endif + comm->barrier(); + + return success; +} // SetupSolve + + +template +int main_(Teuchos::CommandLineProcessor &clp, Xpetra::UnderlyingLib lib, int argc, char *argv[]) { +#include + + RCP< const Teuchos::Comm > comm = Teuchos::DefaultComm::getComm(); + RCP out = Teuchos::fancyOStream(Teuchos::rcpFromRef(std::cout)); + out->setOutputToRootOnly(0); + + bool printTimings = true; clp.setOption("timings", "notimings", &printTimings, "print timings to screen"); + std::string timingsFormat = "table-fixed"; clp.setOption("time-format", &timingsFormat, "timings format (table-fixed | table-scientific | yaml)"); + double scaling = 1.0; clp.setOption("scaling", &scaling, "scale mass term"); + std::string solverName = "Belos"; clp.setOption("solverName", &solverName, "Name of iterative linear solver " + "to use for solving the linear system. " + "(\"Belos\" or \"Stratimikos\")"); + std::string belosSolverType = "Block CG"; clp.setOption("belosSolverType", &belosSolverType, "Name of the Belos linear solver"); + std::string precType = "MueLu-RefMaxwell"; clp.setOption("precType", &precType, "preconditioner to use (MueLu-RefMaxwell|ML-RefMaxwell|none)"); + std::string xml = ""; clp.setOption("xml", &xml, "xml file with solver parameters (default: \"Maxwell.xml\")"); + std::string belos_xml = "Belos.xml"; clp.setOption("belos_xml", &belos_xml, "xml file with Belos solver parameters (default: \"Belos.xml\")"); + std::string stratimikos_xml = "Stratimikos.xml"; clp.setOption("stratimikos_xml", &stratimikos_xml, "xml file with Stratimikos solver parameters (default: \"Stratimikos.xml\")"); + int numResolves = 0; clp.setOption("resolve", &numResolves, "#times to redo solve"); + double tol = 1e-10; clp.setOption("tol", &tol, "solver convergence tolerance"); + int maxIts = 200; clp.setOption("its", &maxIts, "maximum number of solver iterations"); + bool use_stacked_timer = false; clp.setOption("stacked-timer", + "no-stacked-timer", &use_stacked_timer, "use stacked timer"); + + std::string S_file, D0_file, M1_file, M0_file; + if (!Teuchos::ScalarTraits::isComplex) { + S_file = "S.mat"; + D0_file = "D0.mat"; + M1_file = "M1.mat"; + M0_file = "M0.mat"; + } else { + S_file = "S_complex.mat"; + D0_file = "D0_complex.mat"; + M1_file = "M1_complex.mat"; + M0_file = "M0_complex.mat"; + } + clp.setOption("S", &S_file); + std::string SM_file = ""; clp.setOption("SM", &SM_file); + std::string Kn_file = ""; clp.setOption("Kn", &Kn_file); + clp.setOption("D0", &D0_file); + clp.setOption("M1", &M1_file); + std::string Ms_file = ""; clp.setOption("Ms", &Ms_file); + clp.setOption("M0", &M0_file); + std::string M0inv_file = ""; clp.setOption("M0inv", &M0inv_file); + + std::string coords_file = "coords.mat"; clp.setOption("coords", &coords_file); + std::string nullspace_file = ""; clp.setOption("nullspace", &nullspace_file); + std::string material_file = ""; clp.setOption("material", &material_file); + + std::string rhs_file = ""; clp.setOption("rhs", &rhs_file); + std::string x0_file = ""; clp.setOption("x0", &x0_file); + + clp.recogniseAllOptions(true); + switch (clp.parse(argc, argv)) { + case Teuchos::CommandLineProcessor::PARSE_HELP_PRINTED: return EXIT_SUCCESS; + case Teuchos::CommandLineProcessor::PARSE_ERROR: + case Teuchos::CommandLineProcessor::PARSE_UNRECOGNIZED_OPTION: return EXIT_FAILURE; + case Teuchos::CommandLineProcessor::PARSE_SUCCESSFUL: break; + } + + if (xml == ""){ + if (precType == "MueLu-RefMaxwell") + xml = "Maxwell.xml"; + else if (precType == "ML-RefMaxwell") + xml = "Maxwell_ML.xml"; + else if (precType == "ML-Maxwell") + xml = "Maxwell_ML1.xml"; + else if (precType == "Hypre") + xml = "Hypre.xml"; + } - if (printTimings) { + RCP stacked_timer; + if (use_stacked_timer) + stacked_timer = rcp(new Teuchos::StackedTimer("Maxwell Driver")); + TimeMonitor::setStackedTimer(stacked_timer); + + auto globalTimeMonitor = TimeMonitor::getNewTimer("Maxwell: S - Global Time"); + auto tm = TimeMonitor::getNewTimer("Maxwell: 1 - Read and Build Matrices"); + + // Read matrices in from files + RCP D0_Matrix, SM_Matrix, M1_Matrix, Ms_Matrix, M0inv_Matrix, Kn_Matrix; + + // gradient matrix + D0_Matrix = Xpetra::IO::Read(D0_file, lib, comm); + + // maps for nodal and edge matrices + RCP node_map = D0_Matrix->getDomainMap(); + RCP edge_map = D0_Matrix->getRangeMap(); + + // build stiffness plus mass matrix (SM_Matrix) + if (SM_file == "") { + // edge stiffness matrix + RCP S_Matrix = Xpetra::IO::Read(S_file, edge_map); + M1_Matrix = Xpetra::IO::Read(M1_file, edge_map); + Xpetra::MatrixMatrix::TwoMatrixAdd(*S_Matrix,false,(SC)1.0,*M1_Matrix,false,scaling,SM_Matrix,*out); + SM_Matrix->fillComplete(); + } else { + SM_Matrix = Xpetra::IO::Read(SM_file, edge_map); + if (M1_file != "") + M1_Matrix = Xpetra::IO::Read(M1_file, edge_map); + } + if (Ms_file != "") + Ms_Matrix = Xpetra::IO::Read(Ms_file, edge_map); + else + Ms_Matrix = M1_Matrix; + if (Kn_file!="") + Kn_Matrix = Xpetra::IO::Read(Kn_file, node_map); + + if ((M0inv_file == "") && (M0_file != "")) { + // nodal mass matrix + RCP M0_Matrix = Xpetra::IO::Read(M0_file, node_map); + // build lumped mass matrix inverse (M0inv_Matrix) + RCP diag = Utilities::GetLumpedMatrixDiagonal(M0_Matrix); + RCP M0inv_MatrixWrap = rcp(new CrsMatrixWrap(node_map, node_map, 0)); + RCP M0inv_CrsMatrix = M0inv_MatrixWrap->getCrsMatrix(); + Teuchos::ArrayRCP rowPtr; + Teuchos::ArrayRCP colInd; + Teuchos::ArrayRCP values; + Teuchos::ArrayRCP diags = diag->getData(0); + size_t nodeNumElements = node_map->getNodeNumElements(); + M0inv_CrsMatrix->allocateAllValues(nodeNumElements, rowPtr, colInd, values); + SC ONE = Teuchos::ScalarTraits::one(); + for (size_t i = 0; i < nodeNumElements; i++) { + rowPtr[i] = i; colInd[i] = i; values[i] = ONE / diags[i]; + } + rowPtr[nodeNumElements] = nodeNumElements; + M0inv_CrsMatrix->setAllValues(rowPtr, colInd, values); + M0inv_CrsMatrix->expertStaticFillComplete(node_map, node_map); + M0inv_Matrix = Teuchos::rcp_dynamic_cast(M0inv_MatrixWrap); + } else if (M0inv_file != "") { + M0inv_Matrix = Xpetra::IO::Read(M0inv_file, node_map); + } + + // coordinates + RCP::magnitudeType, LO, GO, NO> > coords = Xpetra::IO::magnitudeType, LO, GO, NO>::ReadMultiVector(coords_file, node_map); + + RCP nullspace = Teuchos::null; + if (nullspace_file != "") + nullspace = Xpetra::IO::ReadMultiVector(nullspace_file, edge_map); + + RCP material = Teuchos::null; + if (material_file != "") { + material = Xpetra::IO::ReadMultiVector(material_file, node_map); + } + + // set parameters + Teuchos::ParameterList params, belosParams; + Teuchos::updateParametersFromXmlFileAndBroadcast(xml, Teuchos::Ptr(¶ms), *comm); + Teuchos::updateParametersFromXmlFileAndBroadcast(belos_xml,Teuchos::Ptr(&belosParams),*comm); + belosParams.sublist(belosSolverType).set("Maximum Iterations", maxIts); + belosParams.sublist(belosSolverType).set("Convergence Tolerance",tol); + + // setup LHS, RHS + RCP X, X0, B; + if (rhs_file == "") { + B = MultiVectorFactory::Build(edge_map,1); + RCP vec = MultiVectorFactory::Build(edge_map,1); + vec -> putScalar(Teuchos::ScalarTraits::one()); + SM_Matrix->apply(*vec,*B); + } else + B = Xpetra::IO::ReadMultiVector(rhs_file, edge_map); + + X = MultiVectorFactory::Build(edge_map,1); + X -> putScalar(Teuchos::ScalarTraits::zero()); + if (x0_file != "") + X0 = Xpetra::IO::ReadMultiVector(x0_file, edge_map); + + comm->barrier(); + tm = Teuchos::null; + + std::map inputs; + inputs["SM"] = &SM_Matrix; + inputs["D0"] = &D0_Matrix; + inputs["M1"] = &M1_Matrix; + inputs["Ms"] = &Ms_Matrix; + inputs["M0inv"] = &M0inv_Matrix; + inputs["Kn"] = &Kn_Matrix; + + inputs["coordinates"] = &coords; + inputs["nullspace"] = &nullspace; + inputs["material"] = &material; + + inputs["B"] = &B; + inputs["X"] = &X; + inputs["X0"] = &X0; + + inputs["params"] = ¶ms; + inputs["solverName"] = &solverName; + inputs["belosSolverType"] = &belosSolverType; + inputs["precType"] = &precType; + inputs["belosParams"] = &belosParams; + inputs["numResolves"] = &numResolves; + + inputs["lib"] = &lib; + inputs["comm"] = &comm; + inputs["out"] = &out; + + bool success = SetupSolveWrappers::SetupSolve(inputs); + + globalTimeMonitor = Teuchos::null; + + if (printTimings) { + if (use_stacked_timer) { + stacked_timer->stop("Maxwell Driver"); + Teuchos::StackedTimer::OutputOptions options; + options.output_fraction = options.output_histogram = options.output_minmax = true; + stacked_timer->report(*out, comm, options); + } else { RCP reportParams = rcp(new Teuchos::ParameterList); if (timingsFormat == "yaml") { reportParams->set("Report format", "YAML"); // "Table" or "YAML" @@ -870,31 +768,14 @@ int MainWrappers::main_(Teuchos::Command std::ios_base::fmtflags ff(out->flags()); if (timingsFormat == "table-fixed") *out << std::fixed; else * out << std::scientific; - if (use_stacked_timer) { - stacked_timer->stop("Maxwell Driver"); - Teuchos::StackedTimer::OutputOptions options; - options.output_fraction = options.output_histogram = options.output_minmax = true; - stacked_timer->report(*out, comm, options); - } else - TimeMonitor::report(comm.ptr(), *out, filter, reportParams); - *out << std::setiosflags(ff); + TimeMonitor::report(comm.ptr(), *out, filter, reportParams); + *out << std::setiosflags(ff); } - - TimeMonitor::clearCounters(); -#endif // #ifdef HAVE_MUELU_BELOS } - TEUCHOS_STANDARD_CATCH_STATEMENTS(verbose, std::cerr, success); - - return ( success ? EXIT_SUCCESS : EXIT_FAILURE ); -#else - return EXIT_SUCCESS; -#endif -} // main + TimeMonitor::clearCounters(); -template -int main_(Teuchos::CommandLineProcessor &clp, Xpetra::UnderlyingLib lib, int argc, char *argv[]) { - return MainWrappers::main_(clp, lib, argc, argv); + return ( success ? EXIT_SUCCESS : EXIT_FAILURE ); } //- -- --------------------------------------------------------