diff --git a/feddlib/core/Utils/FEDDUtils.hpp b/feddlib/core/Utils/FEDDUtils.hpp index 4f461390..f2b8a8b9 100644 --- a/feddlib/core/Utils/FEDDUtils.hpp +++ b/feddlib/core/Utils/FEDDUtils.hpp @@ -291,8 +291,13 @@ inline void print(const double s, const Teuchos::RCP> & } } -template void logVec(const std::vector vec, const Teuchos::RCP> &comm, int rank = 0) { +template +void logVec(const std::vector &vec, const Teuchos::RCP> &comm, int rank = 0) { if (comm->getRank() == rank) { + if (vec.empty()) { + std::cout << "You are trying to print an empty vector!" << std::endl; + return; + } std::cout << "[ "; for (auto it = vec.begin(); it != vec.end() - 1; it++) { std::cout << *it << ", "; diff --git a/feddlib/problems/Solver/NonLinearSchwarzSolver/CoarseNonLinearSchwarzOperator_decl.hpp b/feddlib/problems/Solver/NonLinearSchwarzSolver/CoarseNonLinearSchwarzOperator_decl.hpp index 2e4501d7..074fce4d 100644 --- a/feddlib/problems/Solver/NonLinearSchwarzSolver/CoarseNonLinearSchwarzOperator_decl.hpp +++ b/feddlib/problems/Solver/NonLinearSchwarzSolver/CoarseNonLinearSchwarzOperator_decl.hpp @@ -116,8 +116,6 @@ class CoarseNonLinearSchwarzOperator : public IPOUHarmonicCoarseOperator::one(), SC beta = ScalarTraits::zero()) const override; - BlockMatrixPtrFEDD getCoarseJacobian() const; - void exportCoarseBasis(); std::vector getRunStats() const; @@ -133,8 +131,6 @@ class CoarseNonLinearSchwarzOperator : public IPOUHarmonicCoarseOperator +#include #include #include #include @@ -20,6 +21,7 @@ #include #include #include +#include #include #include #include @@ -48,10 +50,13 @@ CoarseNonLinearSchwarzOperator::CoarseNonLinearSchwarzOperator(N ParameterListPtr parameterList) : IPOUHarmonicCoarseOperator(problem->system_->getBlock(0, 0)->getXpetraMatrix(), parameterList), problem_{problem}, x_{Teuchos::rcp(new FEDD::BlockMultiVector(1))}, - y_{Teuchos::rcp(new FEDD::BlockMultiVector(1))}, - coarseJacobian_{Teuchos::rcp(new FEDD::BlockMatrix(1))}, relNewtonTol_{}, absNewtonTol_{}, - maxNumIts_{}, solutionTmp_{Teuchos::rcp(new FEDD::BlockMultiVector(1))}, - systemTmp_{Teuchos::rcp(new FEDD::BlockMatrix(1))}, totalIters_{0} { + y_{Teuchos::rcp(new FEDD::BlockMultiVector(1))}, relNewtonTol_{}, absNewtonTol_{}, maxNumIts_{}, + solutionTmp_{Teuchos::rcp(new FEDD::BlockMultiVector(1))}, + systemTmp_{Teuchos::rcp(new FEDD::BlockMatrix(1))}, + rhsTmp_{Teuchos::rcp(new FEDD::BlockMultiVector(1))}, + sourceTermTmp_{Teuchos::rcp(new FEDD::BlockMultiVector(1))}, + previousSolutionTmp_{Teuchos::rcp(new FEDD::BlockMultiVector(1))}, + residualVecTmp_{Teuchos::rcp(new FEDD::BlockMultiVector(1))}, totalIters_{0} { // Ensure that the mesh object has been initialized and a dual graph generated auto domainPtr_vec = problem_->getDomainVector(); @@ -71,7 +76,15 @@ template int CoarseNonLinearSchwarzOper auto domainVec = problem_->getDomainVector(); auto mesh = domainVec.at(0)->getMesh(); - auto repeatedMap = domainVec.at(0)->getMapRepeated()->getXpetraMap(); + ConstXMapPtr repeatedNodesMap; + ConstXMapPtr repeatedDofsMap; + if (problem_->getDofsPerNode(0) > 1) { + repeatedNodesMap = domainVec.at(0)->getMapRepeated()->getXpetraMap(); + repeatedDofsMap = domainVec.at(0)->getMapVecFieldRepeated()->getXpetraMap(); + } else { + repeatedNodesMap = domainVec.at(0)->getMapRepeated()->getXpetraMap(); + repeatedDofsMap = repeatedNodesMap; + } // Extract info from parameterList relNewtonTol_ = @@ -82,14 +95,16 @@ template int CoarseNonLinearSchwarzOper auto dimension = problem_->getParameterList()->sublist("Parameter").get("Dimension", 2); auto dofsPerNode = problem_->getDofsPerNode(0); totalIters_ = 0; + // Initialize the underlying IPOUHarmonicCoarseOperator object - // TODO: kho dofsMaps need to be built properly for problems with more than one dof per node + // dofs of a node are lumped together + int dofOrdering = 0; auto dofsMaps = Teuchos::ArrayRCP(1); - dofsMaps[0] = domainVec.at(0)->getMapRepeated()->getXpetraMap(); + BuildDofMaps(repeatedDofsMap, dofsPerNode, dofOrdering, repeatedNodesMap, dofsMaps); // Build nodeList auto nodeListFEDD = mesh->getPointsRepeated(); - auto nodeList = Xpetra::MultiVectorFactory::Build(repeatedMap, nodeListFEDD->at(0).size()); + auto nodeList = Xpetra::MultiVectorFactory::Build(repeatedNodesMap, nodeListFEDD->at(0).size()); for (auto i = 0; i < nodeListFEDD->size(); i++) { // pointsRepeated are distributed for (auto j = 0; j < nodeListFEDD->at(0).size(); j++) { @@ -106,27 +121,79 @@ template int CoarseNonLinearSchwarzOper } else { FROSCH_ASSERT(false, "Null Space Type unknown."); } - auto nullSpaceBasis = BuildNullSpace(dimension, nullSpaceType, repeatedMap, dofsPerNode, dofsMaps, + // nullSpaceBasis is a multiVector with one vector for each function in the nullspace e.g. translations and + // rotations for elasticity. Each vector contains the nullspace function on the repeated map. + auto nullSpaceBasis = BuildNullSpace(dimension, nullSpaceType, repeatedDofsMap, dofsPerNode, dofsMaps, implicit_cast(nodeList)); - // Build the vector of Dirichlet node indices + // See FROSch::FindOneEntryOnlyRowsGlobal() for reference auto dirichletBoundaryDofsVec = Teuchos::rcp(new std::vector(0)); int block = 0; int loc = 0; - for (auto i = 0; i < mesh->bcFlagRep_->size(); i++) { + auto out = Teuchos::VerboseObjectBase::getDefaultOStream(); + /* FEDD::logGreen("repeatedNodesmap", this->MpiComm_); */ + /* repeatedNodesMap->describe(*out, VERB_EXTREME); */ + /* FEDD::logGreen("bcFlagRep_", this->MpiComm_); */ + /* FEDD::logVec(*mesh->bcFlagRep_, this->MpiComm_); */ + /* std::cout << std::flush; */ + /* this->MpiComm_->barrier(); */ + /* this->MpiComm_->barrier(); */ + /* this->MpiComm_->barrier(); */ + /* FEDD::logVec(*mesh->bcFlagRep_, this->MpiComm_, 1); */ + /* std::cout << std::flush; */ + /* this->MpiComm_->barrier(); */ + /* this->MpiComm_->barrier(); */ + /* this->MpiComm_->barrier(); */ + /* FEDD::logVec(*mesh->bcFlagRep_, this->MpiComm_, 2); */ + /* std::cout << std::flush; */ + /* this->MpiComm_->barrier(); */ + /* this->MpiComm_->barrier(); */ + /* this->MpiComm_->barrier(); */ + /* FEDD::logVec(*mesh->bcFlagRep_, this->MpiComm_, 3); */ + /* std::cout << std::flush; */ + /* this->MpiComm_->barrier(); */ + /* this->MpiComm_->barrier(); */ + /* this->MpiComm_->barrier(); */ + for (auto i = 0; i < repeatedNodesMap->getLocalNumElements(); i++) { auto flag = mesh->bcFlagRep_->at(i); + // The vector vecFlag_ contains the flags that have been set with addBC(). + // loc: local index in vecFlag_ of the flag if found. Is set by the function. + // block: block id in which to search. Needs to be provided to the function. if (problem_->getBCFactory()->findFlag(flag, block, loc)) { - dirichletBoundaryDofsVec->push_back(repeatedMap->getGlobalElement(i)); + for (auto j = 0; j < dofsPerNode; j++) { + dirichletBoundaryDofsVec->push_back(repeatedDofsMap->getGlobalElement(dofsPerNode * i + j)); + } } } + /* FEDD::logGreen("dirichletBoundaryDofsVec", this->MpiComm_); */ + /* FEDD::logVec(*dirichletBoundaryDofsVec, this->MpiComm_); */ + /* std::cout << std::flush; */ + /* this->MpiComm_->barrier(); */ + /* this->MpiComm_->barrier(); */ + /* this->MpiComm_->barrier(); */ + /* FEDD::logVec(*dirichletBoundaryDofsVec, this->MpiComm_, 1); */ + /* std::cout << std::flush; */ + /* this->MpiComm_->barrier(); */ + /* this->MpiComm_->barrier(); */ + /* this->MpiComm_->barrier(); */ + /* FEDD::logVec(*dirichletBoundaryDofsVec, this->MpiComm_, 2); */ + /* std::cout << std::flush; */ + /* this->MpiComm_->barrier(); */ + /* this->MpiComm_->barrier(); */ + /* this->MpiComm_->barrier(); */ + /* FEDD::logVec(*dirichletBoundaryDofsVec, this->MpiComm_, 3); */ + /* std::cout << std::flush; */ + /* this->MpiComm_->barrier(); */ + /* this->MpiComm_->barrier(); */ + /* this->MpiComm_->barrier(); */ // Convert std::vector to ArrayRCP auto dirichletBoundaryDofs = arcp(dirichletBoundaryDofsVec); // This builds the coarse spaces, assembles the coarse solve map and does symbolic factorization of the coarse // problem IPOUHarmonicCoarseOperator::initialize( - dimension, dofsPerNode, repeatedMap, dofsMaps, nullSpaceBasis, implicit_cast(nodeList), - dirichletBoundaryDofs); + dimension, dofsPerNode, repeatedNodesMap, dofsMaps, nullSpaceBasis, + implicit_cast(nodeList), dirichletBoundaryDofs); return 0; } @@ -139,9 +206,19 @@ void CoarseNonLinearSchwarzOperator::apply(const BlockMultiVecto "input map does not correspond to domain map of nonlinear operator"); x_ = x; + MapConstPtrFEDD uniqueDofsMap; + if (problem_->getDofsPerNode(0) > 1) { + uniqueDofsMap = problem_->getDomainVector().at(0)->getMapVecFieldUnique(); + } else { + uniqueDofsMap = problem_->getDomainVector().at(0)->getMapUnique(); + } // Save problem state - solutionTmp_ = problem_->getSolution(); systemTmp_ = problem_->system_; + solutionTmp_ = problem_->getSolution(); + rhsTmp_ = problem_->rhs_; + sourceTermTmp_ = problem_->sourceTerm_; + previousSolutionTmp_ = problem_->previousSolution_; + residualVecTmp_ = problem_->residualVec_; // Reset problem problem_->initializeProblem(); @@ -156,10 +233,14 @@ void CoarseNonLinearSchwarzOperator::apply(const BlockMultiVecto Xpetra::MultiVectorFactory::Build(this->GatheringMaps_[this->GatheringMaps_.size() - 1], 1); auto coarseDeltaG0 = Xpetra::MultiVectorFactory::Build(this->GatheringMaps_[this->GatheringMaps_.size() - 1], 1); - auto tempMV = Teuchos::rcp(new FEDD::MultiVector(problem_->getDomain(0)->getMapUnique(), 1)); + auto tempMV = Teuchos::rcp(new FEDD::MultiVector(uniqueDofsMap, 1)); auto deltaG0 = Teuchos::rcp(new FEDD::BlockMultiVector(1)); deltaG0->addBlock(tempMV, 0); + // Need to initialize the rhs_ and set boundary values in rhs_ + problem_->assemble(); + problem_->setBoundaries(); + // Set solution_ to be u+P_0*g_0. g_0 is zero on the boundary, simulating P_0 locally // x_ corresponds to u and problem_->solution_ corresponds to g_0 // this = alpha*xTmp + beta*this @@ -222,7 +303,7 @@ void CoarseNonLinearSchwarzOperator::apply(const BlockMultiVecto // Set solution_ to g_i problem_->solution_->update(-ST::one(), *x_, ST::one()); totalIters_ += nlIts; - FEDD::logGreen("==> Terminated coarse Newton", this->MpiComm_); + FEDD::logGreen("Terminated coarse Newton", this->MpiComm_); if (problem_->getParameterList()->sublist("Parameter").get("Cancel MaxNonLinIts", false)) { TEUCHOS_TEST_FOR_EXCEPTION(nlIts == maxNumIts_, std::runtime_error, @@ -230,14 +311,15 @@ void CoarseNonLinearSchwarzOperator::apply(const BlockMultiVecto "step. Still we cancel here."); } - // The currently assembled Jacobian is from the previous Newton iteration. The error this causes is negligable. - // Worth it since a reassemble is avoided. - coarseJacobian_->addBlock(problem_->system_->getBlock(0, 0), 0, 0); y->getBlockNonConst(0)->update(alpha, problem_->getSolution()->getBlock(0), beta); // Restore problem state problem_->initializeProblem(); - problem_->solution_ = solutionTmp_; problem_->system_ = systemTmp_; + problem_->solution_ = solutionTmp_; + problem_->rhs_ = rhsTmp_; + problem_->sourceTerm_ = sourceTermTmp_; + problem_->previousSolution_ = previousSolutionTmp_; + problem_->residualVec_ = residualVecTmp_; } template @@ -264,16 +346,11 @@ void CoarseNonLinearSchwarzOperator::apply(const XMultiVector &x "This apply() overload should not be used with nonlinear operators"); } -template -Teuchos::RCP> -CoarseNonLinearSchwarzOperator::getCoarseJacobian() const { - return coarseJacobian_; -} - template void CoarseNonLinearSchwarzOperator::exportCoarseBasis() { auto comm = problem_->getComm(); + FEDD::logGreen("Exporting coarse basis functions", comm); FEDD::ParameterListPtr_Type pList = problem_->getParameterList(); FEDD::ParameterListPtr_Type pLCoarse = sublist(pList, "Coarse Nonlinear Schwarz"); @@ -281,6 +358,8 @@ void CoarseNonLinearSchwarzOperator::exportCoarseBasis() { TEUCHOS_TEST_FOR_EXCEPTION(!pLCoarse->isParameter("RCP(Phi)"), std::runtime_error, "No parameter to extract Phi pointer."); + // This will have numNodes*dofs rows and numCoarseBasisFunctions columns. Each column is a coarse basis function + // that has a value for each dof at each node. Teuchos::RCP> phiXpetra = this->Phi_; // Convert to a FEDD matrix object @@ -343,14 +422,15 @@ void CoarseNonLinearSchwarzOperator::exportCoarseBasis() { Teuchos::RCP> exportPhi = phiBlockMV->getBlock(i)->getVector(j); std::string varName = fileName + std::to_string(j); - if (pList->get("DofsPerNode" + std::to_string(i + 1), 1) > 1) + if (pLCoarse->get("DofsPerNode" + std::to_string(i + 1), 1) > 1) { exporter->addVariable(exportPhi, varName, "Vector", dofsPerNode, domain->getMapUnique()); - else + } else { exporter->addVariable(exportPhi, varName, "Scalar", 1, domain->getMapUnique()); + } } Teuchos::RCP> exportSumPhi = phiSumBlockMV->getBlock(i); std::string varName = fileName + "Sum"; - if (pList->get("DofsPerNode" + std::to_string(i + 1), 1) > 1) + if (pLCoarse->get("DofsPerNode" + std::to_string(i + 1), 1) > 1) exporter->addVariable(exportSumPhi, varName, "Vector", dofsPerNode, domain->getMapUnique()); else exporter->addVariable(exportSumPhi, varName, "Scalar", 1, domain->getMapUnique()); diff --git a/feddlib/problems/Solver/NonLinearSchwarzSolver/NonLinearSchwarzOperator_decl.hpp b/feddlib/problems/Solver/NonLinearSchwarzSolver/NonLinearSchwarzOperator_decl.hpp index ef7ce09b..cdc89f25 100644 --- a/feddlib/problems/Solver/NonLinearSchwarzSolver/NonLinearSchwarzOperator_decl.hpp +++ b/feddlib/problems/Solver/NonLinearSchwarzSolver/NonLinearSchwarzOperator_decl.hpp @@ -166,7 +166,6 @@ class NonLinearSchwarzOperator : public SchwarzOperator, public BlockMultiVectorPtrFEDD solutionTmp_; BlockMultiVectorPtrFEDD rhsTmp_; BlockMultiVectorPtrFEDD sourceTermTmp_; - std::vector rhsFuncVecTmp_; BlockMultiVectorPtrFEDD previousSolutionTmp_; BlockMultiVectorPtrFEDD residualVecTmp_; // FE assembly factory for global and local assembly diff --git a/feddlib/problems/Solver/NonLinearSchwarzSolver/NonLinearSchwarzOperator_def.hpp b/feddlib/problems/Solver/NonLinearSchwarzSolver/NonLinearSchwarzOperator_def.hpp index ef13653e..ffba01c3 100644 --- a/feddlib/problems/Solver/NonLinearSchwarzSolver/NonLinearSchwarzOperator_def.hpp +++ b/feddlib/problems/Solver/NonLinearSchwarzSolver/NonLinearSchwarzOperator_def.hpp @@ -61,7 +61,7 @@ NonLinearSchwarzOperator::NonLinearSchwarzOperator(CommPtr seria systemTmp_{Teuchos::rcp(new FEDD::BlockMatrix(1))}, solutionTmp_{Teuchos::rcp(new FEDD::BlockMultiVector(1))}, rhsTmp_{Teuchos::rcp(new FEDD::BlockMultiVector(1))}, - sourceTermTmp_{Teuchos::rcp(new FEDD::BlockMultiVector(1))}, rhsFuncVecTmp_{0}, + sourceTermTmp_{Teuchos::rcp(new FEDD::BlockMultiVector(1))}, previousSolutionTmp_{Teuchos::rcp(new FEDD::BlockMultiVector(1))}, residualVecTmp_{Teuchos::rcp(new FEDD::BlockMultiVector(1))}, feFactoryTmp_{Teuchos::rcp(new FEDD::FE())}, @@ -181,10 +181,10 @@ void NonLinearSchwarzOperator::apply(const BlockMultiVectorPtrFE solutionTmp_ = problem_->getSolution(); rhsTmp_ = problem_->rhs_; sourceTermTmp_ = problem_->sourceTerm_; - rhsFuncVecTmp_ = problem_->rhsFuncVec_; feFactoryTmp_ = problem_->feFactory_; previousSolutionTmp_ = problem_->previousSolution_; residualVecTmp_ = problem_->residualVec_; + // Do not need to store rhsFuncVec_ because reseting problem does not erase it // ================= Replace shared objects =============================== // Done to "trick" FEDD::Problem assembly routines to assemble locally on the overlapping subdomain @@ -214,8 +214,6 @@ void NonLinearSchwarzOperator::apply(const BlockMultiVectorPtrFE // Problems block vectors and matrices need to be reinitialized problem_->initializeProblem(); - // the rhsFuncVec_ does not care about being local or global - problem_->rhsFuncVec_ = rhsFuncVecTmp_; if (problem_->getDofsPerNode(0) > 1) { x_->getBlockNonConst(0)->replaceMap(FEDDMapVecFieldOverlappingGhostsLocal); @@ -261,14 +259,14 @@ void NonLinearSchwarzOperator::apply(const BlockMultiVectorPtrFE double absResidual = 1.; int nlIts = 0; - // Need to initialize the rhs_ and set boundary values in rhs_ + // Need to initialize the rhs_ and set boundary values in rhs_ problem_->assemble(); problem_->setBoundaries(); // Set solution_ to be u+P_i*g_i. g_i is zero on the boundary, simulating P_i locally // this = alpha*xTmp + beta*this problem_->solution_->update(ST::one(), *x_, ST::one()); - + // Need to update solution_ within each iteration to assemble at u+P_i*g_i but update only g_i // This is necessary since u is nonzero on the artificial (interface) zero Dirichlet boundary // It would be more efficient to only store u on the boundary and update this value in each iteration @@ -307,7 +305,7 @@ void NonLinearSchwarzOperator::apply(const BlockMultiVectorPtrFE // Set solution_ to g_i problem_->solution_->update(-ST::one(), *x_, ST::one()); - FEDD::logGreen("==> Terminated inner Newton", this->MpiComm_); + FEDD::logGreen("Terminated inner Newton", this->MpiComm_); totalIters_ += nlIts; if (problem_->getParameterList()->sublist("Parameter").get("Cancel MaxNonLinIts", false)) { @@ -363,7 +361,6 @@ void NonLinearSchwarzOperator::apply(const BlockMultiVectorPtrFE this->replaceMapAndExportProblem(); this->MpiComm_->barrier(); // Restore system state - // TODO: kho maybe these reinits are not needed? problem_->initializeProblem(); problem_->system_ = systemTmp_; problem_->solution_ = solutionTmp_; @@ -450,7 +447,7 @@ template string NonLinearSchwarzOperato template void NonLinearSchwarzOperator::replaceMapAndExportProblem() { - // TODO: KHo make this work for block and vector-valued systems + // TODO: kho make this work for block and vector-valued systems auto domainPtr_vec = problem_->getDomainVector(); MapConstPtrFEDD mapUnique; MapConstPtrFEDD mapOverlappingGhosts; diff --git a/feddlib/problems/Solver/NonLinearSchwarzSolver/SimpleCoarseOperator_decl.hpp b/feddlib/problems/Solver/NonLinearSchwarzSolver/SimpleCoarseOperator_decl.hpp index bae97b2b..ec6adbd0 100644 --- a/feddlib/problems/Solver/NonLinearSchwarzSolver/SimpleCoarseOperator_decl.hpp +++ b/feddlib/problems/Solver/NonLinearSchwarzSolver/SimpleCoarseOperator_decl.hpp @@ -21,8 +21,10 @@ /*! Declaration of simple coarse operator This class is just a wrapper around an existing NonlinearCoarseOperator providing it with an alternative apply() method - that evaluates the tanget. This design choice facilitates reusing the underlying CoarseOperator for evaluation of the - coarse problem and its tangent. + that evaluates the tanget. The alternative apply() method is provided by the underlying CoarseOperator object. The + initialize() method uses the default assignment operator of CoarseOperator to copy the member variables (which are + mostly pointers) of the passed CoarseOperator object. By passing a NonlinearCoarseOperator object, its coarse Jacobian + and coarse space can be used in the apply() method of CoarseOperator through the SimpleCoarseOperator object. @brief Implements the coarse tangent $D\mathcal{F}(u) = P_0(R_0DF(u_0)P_0)^{-1}R_0DF(u_0)$ from the nonlinear Schwarz approach @@ -105,7 +107,7 @@ class SimpleCoarseOperator : public CoarseOperator { ~SimpleCoarseOperator() = default; int initialize() override; - int initialize(Teuchos::RCP> inputOp); + int initialize(const Teuchos::RCP> inputOp); int compute() override; diff --git a/feddlib/problems/Solver/NonLinearSchwarzSolver/SimpleCoarseOperator_def.hpp b/feddlib/problems/Solver/NonLinearSchwarzSolver/SimpleCoarseOperator_def.hpp index a05f9ed9..2a8d6ee3 100644 --- a/feddlib/problems/Solver/NonLinearSchwarzSolver/SimpleCoarseOperator_def.hpp +++ b/feddlib/problems/Solver/NonLinearSchwarzSolver/SimpleCoarseOperator_def.hpp @@ -2,6 +2,7 @@ #define SIMPLECOARSEOPERATPR_DEF_HPP #include "SimpleCoarseOperator_decl.hpp" +#include "feddlib/core/Utils/FEDDUtils.hpp" #include #include #include @@ -36,10 +37,10 @@ template int SimpleCoarseOperator -int SimpleCoarseOperator::initialize(Teuchos::RCP> inputOp) { +int SimpleCoarseOperator::initialize(const Teuchos::RCP> inputOp) { // Shallow copy into underlying CoarseOperator object // For this to work FROSch has to be modified (considering master branch at commit - // 5d204ffb426c04e52f76bb49c3ba59627406ac00) by making LevelID_ in FROSch_SchwarzOperator_decl.hpp non-const. This + // 986c0264d5bc81983e78227ed9375dd1105bde10) by making LevelID_ in FROSch_SchwarzOperator_decl.hpp non-const. This // ensures that default assignment operators are generated by the compiler (and can be used here). CoarseOperator::operator=(*inputOp); return 0; diff --git a/feddlib/problems/Solver/NonLinearSolver_def.hpp b/feddlib/problems/Solver/NonLinearSolver_def.hpp index 9d34a92e..7f251675 100644 --- a/feddlib/problems/Solver/NonLinearSolver_def.hpp +++ b/feddlib/problems/Solver/NonLinearSolver_def.hpp @@ -774,10 +774,12 @@ void NonLinearSolver::solveNonLinearSchwarz(NonLinearProblem_Typ print("================= Nonlinear Schwarz terminated =========================", mpiComm); print("\n\nOuter Newton:", mpiComm, 25); print("Rel. residual:", mpiComm, 15); + print("Abs. residual:", mpiComm, 15); print("Iters:", mpiComm, 15); print("GMRES iters:", mpiComm, 15); print("\n", mpiComm, 25); print(relativeResidual, mpiComm, 15, 2); + print(residual, mpiComm, 15, 2); print(outerNonLinIts, mpiComm, 15, 2); print(gmresIts, mpiComm, 15, 2); print("\n\nInner Newton:", mpiComm, 25); diff --git a/feddlib/problems/examples/nonLinElasWithNonLinSchwarz/CMakeLists.txt b/feddlib/problems/examples/nonLinElasWithNonLinSchwarz/CMakeLists.txt index 7bd72550..0a8a56b2 100644 --- a/feddlib/problems/examples/nonLinElasWithNonLinSchwarz/CMakeLists.txt +++ b/feddlib/problems/examples/nonLinElasWithNonLinSchwarz/CMakeLists.txt @@ -24,11 +24,13 @@ TRIBITS_COPY_FILES_TO_BINARY_DIR(mesh_nonLinElasWithNonLinSchwarz SOURCE_FILES square_9.mesh SOURCE_FILES square_900.mesh SOURCE_FILES square_2500.mesh + SOURCE_FILES beam2D.mesh DEST_DIR ${CMAKE_CURRENT_BINARY_DIR} DEST_FILES square.mesh DEST_FILES square_9.mesh DEST_FILES square_900.mesh DEST_FILES square_2500.mesh + DEST_FILES beam2D.mesh EXEDEPS nonLinElasWithNonLinSchwarz ) diff --git a/feddlib/problems/examples/nonLinElasWithNonLinSchwarz/main.cpp b/feddlib/problems/examples/nonLinElasWithNonLinSchwarz/main.cpp index 3bf10f80..d95b6a0a 100644 --- a/feddlib/problems/examples/nonLinElasWithNonLinSchwarz/main.cpp +++ b/feddlib/problems/examples/nonLinElasWithNonLinSchwarz/main.cpp @@ -102,10 +102,6 @@ int main(int argc, char *argv[]) { myCLP.setOption("length", &length, "length of domain."); string xmlProblemFile = "parametersProblem.xml"; myCLP.setOption("problemfile", &xmlProblemFile, ".xml file with Inputparameters."); - string xmlPrecFile = "parametersPrec.xml"; - myCLP.setOption("precfile", &xmlPrecFile, ".xml file with Inputparameters."); - string xmlSolverFile = "parametersSolver.xml"; - myCLP.setOption("solverfile", &xmlSolverFile, ".xml file with Inputparameters."); string xmlSchwarzSolverFile = "parametersSolverNonLinSchwarz.xml"; myCLP.setOption("schwarzsolverfile", &xmlSchwarzSolverFile, ".xml file with Inputparameters."); @@ -118,8 +114,6 @@ int main(int argc, char *argv[]) { } ParameterListPtr_Type parameterListProblem = Teuchos::getParametersFromXmlFile(xmlProblemFile); - ParameterListPtr_Type parameterListPrec = Teuchos::getParametersFromXmlFile(xmlPrecFile); - ParameterListPtr_Type parameterListSolver = Teuchos::getParametersFromXmlFile(xmlSolverFile); ParameterListPtr_Type parameterListSchwarzSolver = Teuchos::getParametersFromXmlFile(xmlSchwarzSolverFile); int dim = parameterListProblem->sublist("Parameter").get("Dimension", 3); @@ -132,8 +126,6 @@ int main(int argc, char *argv[]) { int n; ParameterListPtr_Type parameterListAll(new Teuchos::ParameterList(*parameterListProblem)); - parameterListAll->setParameters(*parameterListPrec); - parameterListAll->setParameters(*parameterListSolver); parameterListAll->setParameters(*parameterListSchwarzSolver); int minNumberSubdomains = (int)length; @@ -200,22 +192,31 @@ int main(int argc, char *argv[]) { Teuchos::RCP> bcFactory(new BCBuilder()); if (meshType == "structured") { - if (dim == 2) + if (dim == 2) { + /* bcFactory->addBC(zeroDirichlet2D, 1, 0, domain, "Dirichlet", dim); */ bcFactory->addBC(zeroDirichlet2D, 2, 0, domain, "Dirichlet", dim); - else if (dim == 3) + /* bcFactory->addBC(zeroDirichlet2D, 3, 0, domain, "Dirichlet", dim); */ + bcFactory->addBC(zeroDirichlet2D, 4, 0, domain, "Dirichlet", dim); + } else if (dim == 3) { bcFactory->addBC(zeroDirichlet3D, 2, 0, domain, "Dirichlet", dim); + } } else if (meshType == "unstructured") { - if (dim == 2) + if (dim == 2) { // For unstructured meshes and surface force set zeroBC at flag 1 since the surface force is applied at flag // 3. For volume forces any boundary flag can be taken. This will affect which side of the geometry is // stationary - bcFactory->addBC(zeroDirichlet2D, 1, 0, domain, "Dirichlet", dim); - else if (dim == 3) + /* bcFactory->addBC(zeroDirichlet2D, 1, 0, domain, "Dirichlet", dim); */ + bcFactory->addBC(zeroDirichlet2D, 2, 0, domain, "Dirichlet", dim); + /* bcFactory->addBC(zeroDirichlet2D, 3, 0, domain, "Dirichlet", dim); */ + bcFactory->addBC(zeroDirichlet2D, 4, 0, domain, "Dirichlet", dim); + + } else if (dim == 3) { bcFactory->addBC(zeroDirichlet3D, 1, 0, domain, "Dirichlet", dim); + } } // The current global solution must be set as the Dirichlet BC on the ghost nodes for nonlinear Schwarz solver to // correctly solve on the subdomains - std::vector funcParams {static_cast(dim)}; + std::vector funcParams{static_cast(dim)}; bcFactory->addBC(Helper::currentSolutionDirichlet, -99, 0, domain, "Dirichlet", dim, funcParams); NonLinElasticity elasticity(domain, FEType, parameterListAll); @@ -232,7 +233,8 @@ int main(int argc, char *argv[]) { // corresponding to the specified boundary /* elasticity.addRhsFunction(rhsY); */ // Use this in conjuction with Source Type = Volume. It applies a constant force at every node - // Applying an elongating force results in a stable solution for much larger deformations than a compressing force + // Applying an elongating force results in a stable solution for much larger deformations than a compressing + // force elasticity.addRhsFunction(rhs2D); else if (dim == 3) elasticity.addRhsFunction(rhsX); diff --git a/feddlib/problems/examples/nonLinElasWithNonLinSchwarz/parametersProblem.xml b/feddlib/problems/examples/nonLinElasWithNonLinSchwarz/parametersProblem.xml index 34469757..cf5aa621 100644 --- a/feddlib/problems/examples/nonLinElasWithNonLinSchwarz/parametersProblem.xml +++ b/feddlib/problems/examples/nonLinElasWithNonLinSchwarz/parametersProblem.xml @@ -5,18 +5,18 @@ - + - + - + - + @@ -41,7 +41,10 @@ - + + + + diff --git a/feddlib/problems/examples/nonLinElasWithNonLinSchwarz/parametersSolverNonLinSchwarz.xml b/feddlib/problems/examples/nonLinElasWithNonLinSchwarz/parametersSolverNonLinSchwarz.xml index 3118435a..1a28dd7c 100644 --- a/feddlib/problems/examples/nonLinElasWithNonLinSchwarz/parametersSolverNonLinSchwarz.xml +++ b/feddlib/problems/examples/nonLinElasWithNonLinSchwarz/parametersSolverNonLinSchwarz.xml @@ -1,9 +1,9 @@ - - - + + + @@ -23,7 +23,7 @@ - + @@ -41,15 +41,18 @@ - - - + + + + + + diff --git a/feddlib/problems/examples/nonLinElasticity/parametersProblem.xml b/feddlib/problems/examples/nonLinElasticity/parametersProblem.xml index 5772eec4..ae76f117 100644 --- a/feddlib/problems/examples/nonLinElasticity/parametersProblem.xml +++ b/feddlib/problems/examples/nonLinElasticity/parametersProblem.xml @@ -8,15 +8,15 @@ - + - - + + - + @@ -41,13 +41,14 @@ - + + - + diff --git a/feddlib/problems/specific/NonLinElasticity_decl.hpp b/feddlib/problems/specific/NonLinElasticity_decl.hpp index 4d61e736..a498ebf1 100644 --- a/feddlib/problems/specific/NonLinElasticity_decl.hpp +++ b/feddlib/problems/specific/NonLinElasticity_decl.hpp @@ -58,38 +58,38 @@ class NonLinElasticity : public NonLinearProblem { //@} ~NonLinElasticity(); - virtual void info(); + void info() override; - virtual void assemble( std::string type = "" ) const; + void assemble( std::string type = "" ) const override; void reAssemble(std::string type) const; - virtual void reAssemble( BlockMultiVectorPtr_Type previousSolution ) const{}; + void reAssemble( BlockMultiVectorPtr_Type previousSolution ) const override{}; virtual void reAssemble( MatrixPtr_Type& massmatrix, std::string type ) const; - virtual void reAssembleExtrapolation(BlockMultiVectorPtrArray_Type previousSolutions); + void reAssembleExtrapolation(BlockMultiVectorPtrArray_Type previousSolutions) override; - virtual void calculateNonLinResidualVec(std::string type, double time=0.) const; + void calculateNonLinResidualVec(std::string type, double time=0.) const override; - virtual void getValuesOfInterest( vec_dbl_Type& values ){}; + void getValuesOfInterest( vec_dbl_Type& values ) override {}; - virtual void computeValuesOfInterestAndExport() {}; + void computeValuesOfInterestAndExport() override {}; // virtual void assembleExternal( std::string type ){}; - Teuchos::RCP< Thyra::LinearOpBase > create_W_op() const; + Teuchos::RCP< Thyra::LinearOpBase > create_W_op() const override; - Teuchos::RCP > create_W_prec() const; + Teuchos::RCP > create_W_prec() const override; void reInitSpecificProblemVectors(const MapConstPtr_Type newMap) override; private: - virtual void evalModelImpl( + void evalModelImpl( const ::Thyra::ModelEvaluatorBase::InArgs &inArgs, const ::Thyra::ModelEvaluatorBase::OutArgs &outArgs - ) const; + ) const override; mutable MultiVectorPtr_Type u_rep_; double E_; diff --git a/feddlib/problems/specific/NonLinElasticity_def.hpp b/feddlib/problems/specific/NonLinElasticity_def.hpp index cd29aaaf..e527f190 100644 --- a/feddlib/problems/specific/NonLinElasticity_def.hpp +++ b/feddlib/problems/specific/NonLinElasticity_def.hpp @@ -99,7 +99,6 @@ void NonLinElasticity::reAssemble(std::string type) const { MultiVectorPtr_Type f = Teuchos::rcp( new MultiVector_Type( this->getDomain(0)->getMapVecFieldRepeated(), 1 ) ); MatrixPtr_Type W = Teuchos::rcp(new Matrix_Type( this->getDomain(0)->getMapVecFieldUnique(), this->getDomain(0)->getDimension() * this->getDomain(0)->getApproxEntriesPerRow() ) ); - //TODO: kho is it necessary to assemble Jacobian and rhs at the same time or is this done for efficiency? (not iterate over all elements twice) this->feFactory_->assemblyElasticityJacobianAndStressAceFEM(this->dim_, this->getDomain(0)->getFEType(), W, f, u_rep_, this->parameterList_, C_); MultiVectorPtr_Type fUnique = Teuchos::rcp( new MultiVector_Type( this->getDomain(0)->getMapVecFieldUnique(), 1 ) ); diff --git a/feddlib/problems/tests/nonLinSchwarz/main.cpp b/feddlib/problems/tests/nonLinSchwarz/main.cpp index 33010101..ca5bf165 100644 --- a/feddlib/problems/tests/nonLinSchwarz/main.cpp +++ b/feddlib/problems/tests/nonLinSchwarz/main.cpp @@ -122,7 +122,7 @@ int main(int argc, char *argv[]) { string meshType = parameterListProblem->sublist("Parameter").get("Mesh Type", "structured"); string FEType = parameterListProblem->sublist("Parameter").get("Discretization", "P1"); auto overlap = parameterListSolver->get("Overlap", 1); - int numProcsCoarseSolve = 0; // TODO: kho need to implement this + int numProcsCoarseSolve = 0; int size = comm->getSize() - numProcsCoarseSolve; // ########################