Updated recombination options

feddlib/problems/Solver/NonLinearSchwarzSolver/NonLinearSchwarzOperator_decl.hpp

@@ -156,7 +156,7 @@ class NonLinearSchwarzOperator : public SchwarzOperator<SC, LO, GO, NO>, public
     // Vectors for saving repeated and unique points
     FEDD::vec2D_dbl_ptr_Type pointsRepTmp_;
     FEDD::vec2D_dbl_ptr_Type pointsUniTmp_;
-    // Vectors for saving the boundary conditios
+    // Vectors for saving the boundary conditions
     FEDD::vec_int_ptr_Type bcFlagRepTmp_;
     FEDD::vec_int_ptr_Type bcFlagUniTmp_;
     // Vector of elements for saving elementsC_

feddlib/problems/Solver/NonLinearSchwarzSolver/NonLinearSchwarzOperator_def.hpp

@@ -54,7 +54,7 @@ NonLinearSchwarzOperator<SC, LO, GO, NO>::NonLinearSchwarzOperator(CommPtr seria
       y_{Teuchos::rcp(new FEDD::BlockMultiVector<SC, LO, GO, NO>(1))},
       localJacobianGhosts_{Teuchos::rcp(new FEDD::BlockMatrix<SC, LO, GO, NO>(1))}, mapOverlappingGhostsLocal_{},
       mapVecFieldOverlappingGhostsLocal_{}, relNewtonTol_{}, absNewtonTol_{}, maxNumIts_{},
-      combinationMode_{CombinationMode::Full},
+      combinationMode_{},
       multiplicity_{Teuchos::rcp(new FEDD::BlockMultiVector<SC, LO, GO, NO>(1))}, mapRepeatedMpiTmp_{},
       mapUniqueMpiTmp_{}, mapVecFieldRepeatedMpiTmp_{}, mapVecFieldUniqueMpiTmp_{}, pointsRepTmp_{}, pointsUniTmp_{},
       bcFlagRepTmp_{}, bcFlagUniTmp_{}, elementsCTmp_{},
@@ -106,7 +106,7 @@ template <class SC, class LO, class GO, class NO> int NonLinearSchwarzOperator<S
         problem_->getParameterList()->sublist("Inner Newton Nonlinear Schwarz").get("Absolute Tolerance", 1.0e-6);
     maxNumIts_ = problem_->getParameterList()->sublist("Inner Newton Nonlinear Schwarz").get("Max Iterations", 10);
     totalIters_ = 0;
-    auto combineModeTemp = problem_->getParameterList()->get("Combine Mode", "Full");
+    auto combineModeTemp = problem_->getParameterList()->get("Combine Mode", "Restricted");
 
     if (combineModeTemp == "Averaging") {
         combinationMode_ = CombinationMode::Averaging;
@@ -116,9 +116,9 @@ template <class SC, class LO, class GO, class NO> int NonLinearSchwarzOperator<S
         combinationMode_ = CombinationMode::Restricted;
     } else {
         if (this->MpiComm_->getRank() == 0) {
-            std::cout << "\nInvalid Recombination Mode. Defaulting to Addition" << std::endl;
+            std::cerr << "\nInvalid Recombination Mode in NonLinearSchwarzOperator: \"" << combineModeTemp << "\". Defaulting to \"Restricted\"" << std::endl;
         }
-        combinationMode_ = CombinationMode::Full;
+        combinationMode_ = CombinationMode::Restricted;
     }
 
     mapOverlappingGhostsLocal_ = Xpetra::MapFactory<LO, GO, NO>::Build(

feddlib/problems/Solver/NonLinearSchwarzSolver/SimpleOverlappingOperator_def.hpp

@@ -33,11 +33,18 @@ SimpleOverlappingOperator<SC, LO, GO, NO>::SimpleOverlappingOperator(NonLinearPr
       uniqueMap_(), importerUniqueToGhosts_(), x_Ghosts_(), y_unique_(), y_Ghosts_(),
       bcFlagOverlappingGhosts_(), problem_(problem) {
     // Override the combine mode of the FROSch operator base object from the nonlinear Schwarz configuration
-    if (!this->ParameterList_->get("Combine Mode", "Restricted").compare("Averaging")) {
+    auto combineModeTemp = problem_->getParameterList()->get("Combine Mode", "Restricted");
+    if (combineModeTemp == "Averaging") {
         this->Combine_ = OverlappingOperator<SC, LO, GO, NO>::CombinationType::Averaging;
-    } else if (!this->ParameterList_->get("Combine Mode", "Restricted").compare("Full")) {
+    } else if (combineModeTemp == "Full") {
         this->Combine_ = OverlappingOperator<SC, LO, GO, NO>::CombinationType::Full;
-    } else if (!this->ParameterList_->get("Combine Mode", "Restricted").compare("Restricted")) {
+    } else if (combineModeTemp == "Restricted") {
+        this->Combine_ = OverlappingOperator<SC, LO, GO, NO>::CombinationType::Restricted;
+    } else {
+        if (this->MpiComm_->getRank() == 0) {
+            std::cerr << "\nInvalid Recombination Mode in SimpleOverlappingOperator: \"" << combineModeTemp
+                      << "\". Defaulting to \"Restricted\"" << std::endl;
+        }
         this->Combine_ = OverlappingOperator<SC, LO, GO, NO>::CombinationType::Restricted;
     }
 }
@@ -120,7 +127,7 @@ void SimpleOverlappingOperator<SC, LO, GO, NO>::apply(const XMultiVector &x, XMu
 
     //  Apply DF(u_i)
     this->OverlappingMatrix_->apply(*x_Ghosts_, *x_Ghosts_, mode, ST::one(), ST::zero());
-   // Set solution on ghost points to zero so that column entries in (R_iDF(u_i)P_i)^-1 corresponding to ghost nodes do
+    // Set solution on ghost points to zero so that column entries in (R_iDF(u_i)P_i)^-1 corresponding to ghost nodes do
     // not affect the solution
     for (int i = 0; i < bcFlagOverlappingGhosts_->size(); i++) {
         if (bcFlagOverlappingGhosts_->at(i) == -99) {
@@ -180,7 +187,8 @@ void SimpleOverlappingOperator<SC, LO, GO, NO>::apply(const XMultiVector &x, XMu
     y.update(alpha, *y_unique_, beta);
 }
 
-// The standard way to apply this operator is with the local Jacobians --> usePreonditionerOnly is not required for this.
+// The standard way to apply this operator is with the local Jacobians --> usePreonditionerOnly is not required for
+// this.
 template <class SC, class LO, class GO, class NO>
 void SimpleOverlappingOperator<SC, LO, GO, NO>::apply(const XMultiVector &x, XMultiVector &y,
                                                       bool usePreconditionerOnly, ETransp mode, SC alpha,
@@ -383,10 +391,10 @@ void SimpleOverlappingOperator<SC, LO, GO, NO>::describe(FancyOStream &out, cons
                                 << rowvals[j] << ") ";
                         }
                     } // globally or locally indexed
-                } // vl == VERB_EXTREME
+                }     // vl == VERB_EXTREME
                 out << endl;
             } // for each row r on this process
-        } // if (myRank == curRank)
+        }     // if (myRank == curRank)
 
         // Give output time to complete
         comm->barrier();

feddlib/problems/Solver/NonLinearSolver_def.hpp

@@ -630,7 +630,6 @@ void NonLinearSolver<SC, LO, GO, NO>::solveNonLinearSchwarz(NonLinearProblem_Typ
 
     // When using auto, this is an rcp to a const NonLinearSumOperator. We need non-const here to ensure that the
     // non-const (nonlinear) apply() overload is called
-  // TODO: kho finish implementing this for diffferent variants. Use an enum here? 
     Teuchos::RCP<FROSch::NonLinearCombineOperator<SC, LO, GO, NO>> rhsCombineOperator;
     Teuchos::RCP<FROSch::CombineOperator<SC, LO, GO, NO>> simpleCombineOperator;
     auto variantString = problem.getParameterList()->get("Nonlin Schwarz Variant", "Additive");
@@ -700,6 +699,19 @@ void NonLinearSolver<SC, LO, GO, NO>::solveNonLinearSchwarz(NonLinearProblem_Typ
     int outerNonLinIts = 0;
     int maxOuterNonLinIts = problem.getParameterList()->sublist("Parameter").get("MaxNonLinIts", 10);
 
+    // Print solver settings
+    logGreen("Nonlinear Schwarz solver settings", mpiComm);
+    print("\tUse ASPEN: ", mpiComm);
+    if (mpiComm->getRank() == 0) {
+        std::cout << boolalpha << useASPEN;
+    }
+    print("\n\tSolver variant: " + variantString, mpiComm);
+    print("\n\tCombine mode: " + problem.getParameterList()->get("Combine Mode", "Restricted"), mpiComm);
+    print("\n\tOverlap: " + std::to_string(problem.getParameterList()->get("Overlap", 1)), mpiComm);
+    print("\n\tNum. levels: " + std::to_string(numLevels), mpiComm);
+    print("\n\tRel. tol: " + std::to_string(outerTol), mpiComm);
+    print("\n\tMax outer Newton iters.: " + std::to_string(maxOuterNonLinIts) + "\n", mpiComm);
+
     // Compute the residual
     problem.calculateNonLinResidualVec("reverse");
     auto residual0 = problem.calculateResidualNorm();
@@ -726,15 +738,11 @@ void NonLinearSolver<SC, LO, GO, NO>::solveNonLinearSchwarz(NonLinearProblem_Typ
             logGreen("Building ASPEN tangent", mpiComm);
             localJacobian = nonLinearSchwarzOp->getLocalJacobianGhosts()->getBlock(0, 0)->getXpetraMatrix();
         } else {
-
-            logGreen("Building ASPIN tangent with FROSch", mpiComm);
+            logGreen("Building ASPIN tangent", mpiComm);
             problem.assemble("Newton");
             problem.setBoundariesSystem();
             // Compute D\mathcal{F}(u) using FROSch and DF(u)
             auto jacobian = problem.getSystem()->getBlock(0, 0)->getXpetraMatrix();
-            logGreen("Global num entries: " +
-                         std::to_string(problem.getSystem()->getBlock(0, 0)->getXpetraMatrix()->getGlobalNumEntries()),
-                     mpiComm);
             jacobianGhosts->setAllToScalar(ST::zero());
             jacobianGhosts->resumeFill();
             jacobianGhosts->doImport(*jacobian, *uniqueToOverlappingGhostsImporter, Xpetra::ADD);

feddlib/problems/examples/nonLinElasWithNonLinSchwarz/main.cpp

@@ -1,6 +1,7 @@
 #include "feddlib/core/FE/Domain.hpp"
 #include "feddlib/core/FEDDCore.hpp"
 #include "feddlib/core/General/ExporterParaView.hpp"
+#include "feddlib/core/LinearAlgebra/MultiVector_decl.hpp"
 #include "feddlib/core/Mesh/MeshPartitioner.hpp"
 #include "feddlib/problems/Solver/NonLinearSolver.hpp"
 #include "feddlib/problems/specific/NonLinElasticity.hpp"
@@ -76,7 +77,6 @@ void zeroDirichlet3D(double *x, double *res, double t, const double *parameters)
     return;
 }
 
-
 typedef unsigned UN;
 typedef default_sc SC;
 typedef default_lo LO;
@@ -256,16 +256,21 @@ int main(int argc, char *argv[]) {
     elasticitySolver.solve(elasticity);
     FEDD_TIMER_STOP(SolveTimer);
 
+    auto rankVec = Teuchos::RCP<MultiVector<SC, LO, GO, NO>>(new MultiVector<SC, LO, GO, NO>(domain->getMapUnique()));
+    rankVec->putScalar(comm->getRank());
+
     if (parameterListAll->sublist("General").get("ParaViewExport", false)) {
         Teuchos::RCP<ExporterParaView<SC, LO, GO, NO>> exParaDisp(new ExporterParaView<SC, LO, GO, NO>());
 
         Teuchos::RCP<const MultiVector<SC, LO, GO, NO>> exportSolutionU = elasticity.getSolution()->getBlock(0);
+        Teuchos::RCP<const MultiVector<SC, LO, GO, NO>> rank = rankVec;
 
         exParaDisp->setup("displacement", domain->getMesh(), domain->getFEType());
 
         UN dofsPerNode = dim;
 
         exParaDisp->addVariable(exportSolutionU, "u", "Vector", dofsPerNode, domain->getMapUnique());
+    exParaDisp->addVariable(rank, "rank", "Scalar", 1, domain->getMapUnique());
         exParaDisp->save(0.0);
         Teuchos::TimeMonitor::report(std::cout);
     }