diff --git a/feddlib/problems/Solver/NonLinearSchwarzSolver/NonLinearSchwarzOperator_decl.hpp b/feddlib/problems/Solver/NonLinearSchwarzSolver/NonLinearSchwarzOperator_decl.hpp index cdc89f25..214f2f21 100644 --- a/feddlib/problems/Solver/NonLinearSchwarzSolver/NonLinearSchwarzOperator_decl.hpp +++ b/feddlib/problems/Solver/NonLinearSchwarzSolver/NonLinearSchwarzOperator_decl.hpp @@ -156,7 +156,7 @@ class NonLinearSchwarzOperator : public SchwarzOperator, 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_ diff --git a/feddlib/problems/Solver/NonLinearSchwarzSolver/NonLinearSchwarzOperator_def.hpp b/feddlib/problems/Solver/NonLinearSchwarzSolver/NonLinearSchwarzOperator_def.hpp index f86b24ad..ac280674 100644 --- a/feddlib/problems/Solver/NonLinearSchwarzSolver/NonLinearSchwarzOperator_def.hpp +++ b/feddlib/problems/Solver/NonLinearSchwarzSolver/NonLinearSchwarzOperator_def.hpp @@ -54,7 +54,7 @@ NonLinearSchwarzOperator::NonLinearSchwarzOperator(CommPtr seria y_{Teuchos::rcp(new FEDD::BlockMultiVector(1))}, localJacobianGhosts_{Teuchos::rcp(new FEDD::BlockMatrix(1))}, mapOverlappingGhostsLocal_{}, mapVecFieldOverlappingGhostsLocal_{}, relNewtonTol_{}, absNewtonTol_{}, maxNumIts_{}, - combinationMode_{CombinationMode::Full}, + combinationMode_{}, multiplicity_{Teuchos::rcp(new FEDD::BlockMultiVector(1))}, mapRepeatedMpiTmp_{}, mapUniqueMpiTmp_{}, mapVecFieldRepeatedMpiTmp_{}, mapVecFieldUniqueMpiTmp_{}, pointsRepTmp_{}, pointsUniTmp_{}, bcFlagRepTmp_{}, bcFlagUniTmp_{}, elementsCTmp_{}, @@ -106,7 +106,7 @@ template int NonLinearSchwarzOperatorgetParameterList()->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 int NonLinearSchwarzOperatorMpiComm_->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::Build( diff --git a/feddlib/problems/Solver/NonLinearSchwarzSolver/SimpleOverlappingOperator_def.hpp b/feddlib/problems/Solver/NonLinearSchwarzSolver/SimpleOverlappingOperator_def.hpp index 2f04bd4f..df994178 100644 --- a/feddlib/problems/Solver/NonLinearSchwarzSolver/SimpleOverlappingOperator_def.hpp +++ b/feddlib/problems/Solver/NonLinearSchwarzSolver/SimpleOverlappingOperator_def.hpp @@ -33,11 +33,18 @@ SimpleOverlappingOperator::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::CombinationType::Averaging; - } else if (!this->ParameterList_->get("Combine Mode", "Restricted").compare("Full")) { + } else if (combineModeTemp == "Full") { this->Combine_ = OverlappingOperator::CombinationType::Full; - } else if (!this->ParameterList_->get("Combine Mode", "Restricted").compare("Restricted")) { + } else if (combineModeTemp == "Restricted") { + this->Combine_ = OverlappingOperator::CombinationType::Restricted; + } else { + if (this->MpiComm_->getRank() == 0) { + std::cerr << "\nInvalid Recombination Mode in SimpleOverlappingOperator: \"" << combineModeTemp + << "\". Defaulting to \"Restricted\"" << std::endl; + } this->Combine_ = OverlappingOperator::CombinationType::Restricted; } } @@ -120,7 +127,7 @@ void SimpleOverlappingOperator::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::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 void SimpleOverlappingOperator::apply(const XMultiVector &x, XMultiVector &y, bool usePreconditionerOnly, ETransp mode, SC alpha, @@ -383,10 +391,10 @@ void SimpleOverlappingOperator::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(); diff --git a/feddlib/problems/Solver/NonLinearSolver_def.hpp b/feddlib/problems/Solver/NonLinearSolver_def.hpp index fe94a711..3d8b0f73 100644 --- a/feddlib/problems/Solver/NonLinearSolver_def.hpp +++ b/feddlib/problems/Solver/NonLinearSolver_def.hpp @@ -630,7 +630,6 @@ void NonLinearSolver::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> rhsCombineOperator; Teuchos::RCP> simpleCombineOperator; auto variantString = problem.getParameterList()->get("Nonlin Schwarz Variant", "Additive"); @@ -700,6 +699,19 @@ void NonLinearSolver::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::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); diff --git a/feddlib/problems/examples/nonLinElasWithNonLinSchwarz/main.cpp b/feddlib/problems/examples/nonLinElasWithNonLinSchwarz/main.cpp index d95b6a0a..f1ecb488 100644 --- a/feddlib/problems/examples/nonLinElasWithNonLinSchwarz/main.cpp +++ b/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>(new MultiVector(domain->getMapUnique())); + rankVec->putScalar(comm->getRank()); + if (parameterListAll->sublist("General").get("ParaViewExport", false)) { Teuchos::RCP> exParaDisp(new ExporterParaView()); Teuchos::RCP> exportSolutionU = elasticity.getSolution()->getBlock(0); + Teuchos::RCP> 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); }