From c478b2ffe9d08af19dfa4e4af3e170e01faf26b2 Mon Sep 17 00:00:00 2001 From: Roy Stogner Date: Mon, 12 Dec 2022 11:04:47 -0600 Subject: [PATCH 1/5] enforce_constraints_exactly() overload w/o System& This can be called from within a solver even if it doesn't have access to the System. --- include/base/dof_map.h | 16 ++++++ include/solvers/linear_solver.h | 21 ++++++++ include/solvers/petsc_linear_solver.h | 36 ++++++++++++- src/base/dof_map_constraints.C | 49 +++++++++-------- src/solvers/linear_solver.C | 16 ++++++ src/solvers/petsc_linear_solver.C | 76 ++++++++++++++++++++++----- 6 files changed, 178 insertions(+), 36 deletions(-) diff --git a/include/base/dof_map.h b/include/base/dof_map.h index 6c16c14629c..03a3bf31b7f 100644 --- a/include/base/dof_map.h +++ b/include/base/dof_map.h @@ -1251,6 +1251,20 @@ class DofMap : public ReferenceCountedObject, */ void constrain_nothing (std::vector & dofs) const; + /** + * Constrains the numeric vector \p v, which represents a solution defined on + * the mesh. This may need to be used after a linear solve, if your linear + * solver's solutions do not satisfy your DoF constraints to a tight enough + * tolerance. + * + * If \p homogeneous == true, heterogeneous constraints are enforced + * as if they were homogeneous. This might be appropriate for e.g. a + * vector representing a difference between two + * heterogeneously-constrained solutions. + */ + void enforce_constraints_exactly (NumericVector & v, + bool homogeneous = false) const; + /** * Constrains the numeric vector \p v, which represents a solution defined on * the mesh. This may need to be used after a linear solve, if your linear @@ -1263,6 +1277,8 @@ class DofMap : public ReferenceCountedObject, * as if they were homogeneous. This might be appropriate for e.g. a * vector representing a difference between two * heterogeneously-constrained solutions. + * + * Will be deprecated eventually; use the non-System& overload. */ void enforce_constraints_exactly (const System & system, NumericVector * v = nullptr, diff --git a/include/solvers/linear_solver.h b/include/solvers/linear_solver.h index cb3ec645869..2045fdce960 100644 --- a/include/solvers/linear_solver.h +++ b/include/solvers/linear_solver.h @@ -40,6 +40,7 @@ template class SparseMatrix; template class NumericVector; template class ShellMatrix; template class Preconditioner; +class DofMap; class System; class SolverConfiguration; enum SolverPackage : int; @@ -152,6 +153,26 @@ class LinearSolver : public ReferenceCountedObject>, virtual void restrict_solve_to (const std::vector * const dofs, const SubsetSolveMode subset_solve_mode=SUBSET_ZERO); + /** + * After calling this method, all successive solves will be + * restricted to unconstrained dofs. This mode can be disabled by + * calling this method with \p dof_map being a \p nullptr. + * + * This does a pure submatrix \ subvector solve, with no correction + * of the rhs to account for ignored columns (so the application of + * constraints to the matrix should be done symmetrically + * beforehand) and no post-solution application of ignored rows (so + * homogeneous or heterogeneous primal or heterogeneous adjoint + * constraints should be applied manually afterward) + */ + virtual void restrict_solve_to_unconstrained (const DofMap * dof_map); + + /** + * Returns the current restriction (or lack thereof) state with + * regard to constrained/unconstrained dofs. + */ + virtual const DofMap * restrictions_to_unconstrained(); + /** * This function calls the solver \p _solver_type preconditioned * with the \p _preconditioner_type preconditioner. diff --git a/include/solvers/petsc_linear_solver.h b/include/solvers/petsc_linear_solver.h index d14e237a68f..f0cdf663a6d 100644 --- a/include/solvers/petsc_linear_solver.h +++ b/include/solvers/petsc_linear_solver.h @@ -147,6 +147,26 @@ class PetscLinearSolver : public LinearSolver virtual void restrict_solve_to (const std::vector * const dofs, const SubsetSolveMode subset_solve_mode=SUBSET_ZERO) override; + /** + * After calling this method, all successive solves will be + * restricted to unconstrained dofs. This mode can be disabled by + * calling this method with \p dof_map being a \p nullptr. + * + * This does a pure submatrix \ subvector solve, with no correction + * of the rhs to account for ignored columns (so the application of + * constraints to the matrix should be done symmetrically + * beforehand) and no post-solution application of ignored rows (so + * homogeneous or heterogeneous primal or heterogeneous adjoint + * constraints should be applied manually afterward) + */ + virtual void restrict_solve_to_unconstrained (const DofMap * dof_map) override; + + /** + * Returns the current restriction (or lack thereof) state with + * regard to constrained/unconstrained dofs. + */ + virtual const DofMap * restrictions_to_unconstrained() override; + /** * Call the Petsc solver. It calls the method below, using the * same matrix for the system and preconditioner matrices. @@ -344,9 +364,21 @@ class PetscLinearSolver : public LinearSolver WrappedPetsc _restrict_solve_to_is_complement; /** - * \returns The local size of \p _restrict_solve_to_is. + * DoFMap to use if we're restricting solves to unconstrained DoFs. + * (\p nullptr means solve on constrained DoFs too) + */ + const DofMap * _restrict_to_unconstrained_dofmap; + + /** + * PETSc index set containing unconstrained dofs on which to solve + * (\p nullptr means solve on all dofs). + */ + WrappedPetsc _unconstrained_dofs_is; + + /** + * \returns The local size of \p is. */ - PetscInt restrict_solve_to_is_local_size() const; + PetscInt local_is_size(IS & is) const; /** * Creates \p _restrict_solve_to_is_complement to contain all diff --git a/src/base/dof_map_constraints.C b/src/base/dof_map_constraints.C index be96c3b07f5..3c8ba4985df 100644 --- a/src/base/dof_map_constraints.C +++ b/src/base/dof_map_constraints.C @@ -2895,8 +2895,7 @@ void DofMap::constrain_nothing (std::vector & dofs) const -void DofMap::enforce_constraints_exactly (const System & system, - NumericVector * v, +void DofMap::enforce_constraints_exactly (NumericVector & v, bool homogeneous) const { parallel_object_only(); @@ -2906,13 +2905,10 @@ void DofMap::enforce_constraints_exactly (const System & system, LOG_SCOPE("enforce_constraints_exactly()","DofMap"); - if (!v) - v = system.solution.get(); - NumericVector * v_local = nullptr; // will be initialized below NumericVector * v_global = nullptr; // will be initialized below std::unique_ptr> v_built; - if (v->type() == SERIAL) + if (v.type() == SERIAL) { v_built = NumericVector::build(this->comm()); v_built->init(this->n_dofs(), this->n_local_dofs(), true, PARALLEL); @@ -2920,39 +2916,38 @@ void DofMap::enforce_constraints_exactly (const System & system, for (dof_id_type i=v_built->first_local_index(); ilast_local_index(); i++) - v_built->set(i, (*v)(i)); + v_built->set(i, v(i)); v_built->close(); v_global = v_built.get(); - v_local = v; + v_local = &v; libmesh_assert (v_local->closed()); } - else if (v->type() == PARALLEL) + else if (v.type() == PARALLEL) { v_built = NumericVector::build(this->comm()); - v_built->init (v->size(), v->local_size(), + v_built->init (v.size(), v.local_size(), this->get_send_list(), true, GHOSTED); - v->localize(*v_built, this->get_send_list()); + v.localize(*v_built, this->get_send_list()); v_built->close(); v_local = v_built.get(); - v_global = v; + v_global = &v; } - else if (v->type() == GHOSTED) + else if (v.type() == GHOSTED) { - v_local = v; - v_global = v; + v_local = &v; + v_global = &v; } - else // unknown v->type() - libmesh_error_msg("ERROR: Unsupported NumericVector type == " << Utility::enum_to_string(v->type())); + else // unknown v.type() + libmesh_error_msg("ERROR: Unsupported NumericVector type == " << Utility::enum_to_string(v.type())); // We should never hit these asserts because we should error-out in // else clause above. Just to be sure we don't try to use v_local // and v_global uninitialized... libmesh_assert(v_local); libmesh_assert(v_global); - libmesh_assert_equal_to (this, &(system.get_dof_map())); for (const auto & [constrained_dof, constraint_row] : _dof_constraints) { @@ -2975,14 +2970,26 @@ void DofMap::enforce_constraints_exactly (const System & system, // If the old vector was serial, we probably need to send our values // to other processors - if (v->type() == SERIAL) + if (v.type() == SERIAL) { #ifndef NDEBUG v_global->close(); #endif - v_global->localize (*v); + v_global->localize (v); } - v->close(); + v.close(); +} + +void DofMap::enforce_constraints_exactly (const System & system, + NumericVector * v, + bool homogeneous) const +{ + libmesh_assert_equal_to (this, &(system.get_dof_map())); + + if (!v) + v = system.solution.get(); + + this->enforce_constraints_exactly(*v, homogeneous); } void DofMap::enforce_constraints_on_residual (const NonlinearImplicitSystem & system, diff --git a/src/solvers/linear_solver.C b/src/solvers/linear_solver.C index 7f54e304c1b..aa49c17c82b 100644 --- a/src/solvers/linear_solver.C +++ b/src/solvers/linear_solver.C @@ -143,6 +143,22 @@ LinearSolver::restrict_solve_to(const std::vector * const dofs, } +template +void +LinearSolver::restrict_solve_to_unconstrained(const DofMap * dof_map) +{ + if (dof_map != nullptr) + libmesh_not_implemented(); +} + + +template +const DofMap * LinearSolver::restrictions_to_unconstrained() +{ + return nullptr; // No support unless this is overridden +} + + template std::pair LinearSolver::adjoint_solve (SparseMatrix & mat, NumericVector & sol, diff --git a/src/solvers/petsc_linear_solver.C b/src/solvers/petsc_linear_solver.C index 8056f7aad55..431a9b5c1fa 100644 --- a/src/solvers/petsc_linear_solver.C +++ b/src/solvers/petsc_linear_solver.C @@ -95,6 +95,8 @@ PetscLinearSolver::PetscLinearSolver(const libMesh::Parallel::Communicator & LinearSolver(comm_in), _restrict_solve_to_is(nullptr), _restrict_solve_to_is_complement(nullptr), + _restrict_to_unconstrained_dofmap(nullptr), + _unconstrained_dofs_is(nullptr), _subset_solve_mode(SUBSET_ZERO) { if (this->n_processors() == 1) @@ -118,6 +120,10 @@ void PetscLinearSolver::clear () if (_restrict_solve_to_is_complement) _restrict_solve_to_is_complement.reset_to_zero(); + _restrict_to_unconstrained_dofmap = nullptr; + if (_unconstrained_dofs_is) + _unconstrained_dofs_is.reset_to_zero(); + // Previously we only called KSPDestroy(), we did not reset _ksp // to nullptr, so that behavior is maintained here. _ksp.destroy(); @@ -349,6 +355,20 @@ PetscLinearSolver::restrict_solve_to (const std::vector * const } +template +void +PetscLinearSolver::restrict_solve_to_unconstrained (const DofMap * dof_map) +{ + _restrict_to_unconstrained_dofmap = dof_map; +} + + +template +const DofMap * PetscLinearSolver::restrictions_to_unconstrained() +{ + return _restrict_to_unconstrained_dofmap; +} + template std::pair @@ -521,11 +541,40 @@ PetscLinearSolver::solve_base (SparseMatrix * matrix, WrappedPetsc subsolution; WrappedPetsc scatter; + // Set up restrictions based on DoF constraints, if asked. + if (_restrict_to_unconstrained_dofmap) + { + std::vector unconstrained_dofs; + for (dof_id_type i = _restrict_to_unconstrained_dofmap->first_dof(), + end_i = _restrict_to_unconstrained_dofmap->end_dof(); + i != end_i; ++i) + if (!_restrict_to_unconstrained_dofmap->is_constrained_dof(i)) + unconstrained_dofs.push_back(cast_int(i)); + + ierr = ISCreateGeneral(this->comm().get(), + cast_int(unconstrained_dofs.size()), + unconstrained_dofs.data(), PETSC_COPY_VALUES, + _unconstrained_dofs_is.get()); + } + // Restrict rhs and solution vectors and set operators. The input // matrix works as the preconditioning matrix. - if (_restrict_solve_to_is) + if (_restrict_solve_to_is || _unconstrained_dofs_is) { - PetscInt is_local_size = this->restrict_solve_to_is_local_size(); + // We'll try to handle the combination eventually, but not + // today. + IS restrict_to_is = [this]() { + if (_restrict_solve_to_is) + { + if (_unconstrained_dofs_is) + libmesh_not_implemented_msg("Can't doubly restrict a solve"); + return *_restrict_solve_to_is.get(); + } + else + return *_unconstrained_dofs_is.get(); + }(); + + PetscInt is_local_size = this->local_is_size(restrict_to_is); ierr = VecCreate(this->comm().get(), subrhs.get()); LIBMESH_CHKERR(ierr); @@ -541,15 +590,15 @@ PetscLinearSolver::solve_base (SparseMatrix * matrix, ierr = VecSetFromOptions(subsolution); LIBMESH_CHKERR(ierr); - ierr = VecScatterCreate(rhs->vec(), _restrict_solve_to_is, subrhs, nullptr, scatter.get()); + ierr = VecScatterCreate(rhs->vec(), restrict_to_is, subrhs, nullptr, scatter.get()); LIBMESH_CHKERR(ierr); VecScatterBeginEnd(this->comm(), scatter, rhs->vec(), subrhs, INSERT_VALUES, SCATTER_FORWARD); VecScatterBeginEnd(this->comm(), scatter, solution->vec(), subsolution, INSERT_VALUES, SCATTER_FORWARD); ierr = LibMeshCreateSubMatrix(mat, - _restrict_solve_to_is, - _restrict_solve_to_is, + restrict_to_is, + restrict_to_is, MAT_INITIAL_MATRIX, submat.get()); LIBMESH_CHKERR(ierr); @@ -557,8 +606,8 @@ PetscLinearSolver::solve_base (SparseMatrix * matrix, if (precond) { ierr = LibMeshCreateSubMatrix(const_cast *>(precond)->mat(), - _restrict_solve_to_is, - _restrict_solve_to_is, + restrict_to_is, + restrict_to_is, MAT_INITIAL_MATRIX, subprecond.get()); LIBMESH_CHKERR(ierr); @@ -568,7 +617,7 @@ PetscLinearSolver::solve_base (SparseMatrix * matrix, // system, we will now change the right hand side to compensate // for this. Note that this is not necessary if \p SUBSET_ZERO // has been selected. - if (_subset_solve_mode!=SUBSET_ZERO) + else if (_subset_solve_mode!=SUBSET_ZERO) { this->create_complement_is(rhs_in); PetscInt is_complement_local_size = @@ -662,7 +711,7 @@ PetscLinearSolver::solve_base (SparseMatrix * matrix, } // Solve the linear system - if (_restrict_solve_to_is) + if (_restrict_solve_to_is || _unconstrained_dofs_is) { ierr = solve_func (_ksp, subrhs, subsolution); LIBMESH_CHKERR(ierr); @@ -702,7 +751,10 @@ PetscLinearSolver::solve_base (SparseMatrix * matrix, default: libmesh_error_msg("Invalid subset solve mode = " << _subset_solve_mode); } + } + if (_restrict_solve_to_is || _unconstrained_dofs_is) + { VecScatterBeginEnd(this->comm(), scatter, subsolution, solution->vec(), INSERT_VALUES, SCATTER_REVERSE); if (precond && this->_preconditioner) @@ -1026,12 +1078,10 @@ PetscLinearSolver::create_complement_is (const NumericVector & vec_in) template PetscInt -PetscLinearSolver::restrict_solve_to_is_local_size() const +PetscLinearSolver::local_is_size(IS & is) const { - libmesh_assert(_restrict_solve_to_is); - PetscInt s; - int ierr = ISGetLocalSize(_restrict_solve_to_is, &s); + int ierr = ISGetLocalSize(is, &s); LIBMESH_CHKERR(ierr); return s; From bad5ff986d95c0e5bcd6123ada2279e016e2f345 Mon Sep 17 00:00:00 2001 From: Roy Stogner Date: Wed, 14 Dec 2022 13:05:54 -0600 Subject: [PATCH 2/5] Allow DiffSolver to restrict to unconstrained Still need a PetscDiffSolver implementation --- include/solvers/diff_solver.h | 9 +++++++++ include/solvers/newton_solver.h | 9 +++++++++ include/solvers/petsc_diff_solver.h | 11 +++++++++++ src/solvers/newton_solver.C | 22 ++++++++++++++++++++++ src/solvers/time_solver.C | 6 ++++++ 5 files changed, 57 insertions(+) diff --git a/include/solvers/diff_solver.h b/include/solvers/diff_solver.h index a8e3cb9abe5..54e6f8e0ad2 100644 --- a/include/solvers/diff_solver.h +++ b/include/solvers/diff_solver.h @@ -139,6 +139,15 @@ class DiffSolver : public ReferenceCountedObject, */ sys_type & system () { return _system; } + /** + * After calling this method with \p true, all successive solves + * will be restricted to unconstrained dofs, with constraints + * applied separately as necessary. + */ + virtual void restrict_solves_to_unconstrained(bool restricting) = 0; + + virtual bool get_restrict_solves_to_unconstrained() = 0; + /** * Each linear solver step should exit after \p max_linear_iterations * is exceeded. diff --git a/include/solvers/newton_solver.h b/include/solvers/newton_solver.h index 92ea32f8dd9..f620d23676b 100644 --- a/include/solvers/newton_solver.h +++ b/include/solvers/newton_solver.h @@ -88,6 +88,15 @@ class NewtonSolver : public DiffSolver return *_linear_solver; } + /** + * After calling this method with \p true, all successive solves + * will be restricted to unconstrained dofs, with constraints + * applied separately as necessary. + */ + virtual void restrict_solves_to_unconstrained(bool restricting) override; + + virtual bool get_restrict_solves_to_unconstrained() override; + /** * If this is set to true, the solver is forced to test the residual * after each Newton step, and to reduce the length of its steps diff --git a/include/solvers/petsc_diff_solver.h b/include/solvers/petsc_diff_solver.h index 7d8ed99dfbe..33f9300370a 100644 --- a/include/solvers/petsc_diff_solver.h +++ b/include/solvers/petsc_diff_solver.h @@ -86,6 +86,17 @@ class PetscDiffSolver : public DiffSolver */ void clear (); + /** + * After calling this method with \p true, all successive solves + * will be restricted to unconstrained dofs, with constraints + * applied separately as necessary. + */ + virtual void restrict_solves_to_unconstrained(bool restricting) override + { if (restricting) libmesh_not_implemented(); } + + virtual bool get_restrict_solves_to_unconstrained() override + { return false; } + /** * This method performs a solve. What occurs in * this method will depend on the PETSc SNES settings. diff --git a/src/solvers/newton_solver.C b/src/solvers/newton_solver.C index 9e3d9af6209..32d064e43af 100644 --- a/src/solvers/newton_solver.C +++ b/src/solvers/newton_solver.C @@ -270,12 +270,34 @@ void NewtonSolver::reinit() { Parent::reinit(); + const DofMap * d = _linear_solver->restrictions_to_unconstrained(); + + // Clear what might be invalidated by system changes _linear_solver->clear(); + // Then "un-clear" a few things + _linear_solver->restrict_solve_to_unconstrained(d); _linear_solver->init_names(_system); } +void NewtonSolver::restrict_solves_to_unconstrained(bool restricting) +{ + if (restricting) + { + const DofMap * d = &this->system().get_dof_map(); + this->get_linear_solver().restrict_solve_to_unconstrained(d); + } + else + this->get_linear_solver().restrict_solve_to_unconstrained(nullptr); +} + + +bool NewtonSolver::get_restrict_solves_to_unconstrained() +{ + return this->get_linear_solver().restrictions_to_unconstrained(); +} + unsigned int NewtonSolver::solve() { diff --git a/src/solvers/time_solver.C b/src/solvers/time_solver.C index 238311c6c84..22d94414e6f 100644 --- a/src/solvers/time_solver.C +++ b/src/solvers/time_solver.C @@ -58,7 +58,13 @@ void TimeSolver::reinit () this->diff_solver()->reinit(); libmesh_assert(this->linear_solver().get()); + const DofMap * d = this->linear_solver()->restrictions_to_unconstrained(); + + // Clear what might be invalidated by system changes this->linear_solver()->clear(); + + // Then "un-clear" a few things + this->linear_solver()->restrict_solve_to_unconstrained(d); if (libMesh::on_command_line("--solver-system-names")) this->linear_solver()->init((_system.name()+"_").c_str()); else From ac32cf0e5d85730154b35fbc7d4c922e3271fd68 Mon Sep 17 00:00:00 2001 From: Roy Stogner Date: Wed, 14 Dec 2022 13:06:49 -0600 Subject: [PATCH 3/5] ImplicitSystem API to restrict to unconstrained This helps encapsulate some of the tricky distinctions between system subclasses and their associated solvers --- include/systems/diff_system.h | 7 +++++++ include/systems/implicit_system.h | 7 +++++++ src/systems/diff_system.C | 7 +++++++ src/systems/implicit_system.C | 20 ++++++++++++++++---- 4 files changed, 37 insertions(+), 4 deletions(-) diff --git a/include/systems/diff_system.h b/include/systems/diff_system.h index e688e8ac19a..9356f40e278 100644 --- a/include/systems/diff_system.h +++ b/include/systems/diff_system.h @@ -304,6 +304,13 @@ class DifferentiableSystem : public ImplicitSystem, */ virtual void side_postprocess (DiffContext &) {} + /** + * After calling this method with \p true, all successive solves + * will be restricted to unconstrained dofs, with constraints + * applied separately as necessary. + */ + virtual void restrict_solves_to_unconstrained(bool restricting) override; + /** * For a given second order (in time) variable var, this method will return * the index to the corresponding "dot" variable. For FirstOrderUnsteadySolver diff --git a/include/systems/implicit_system.h b/include/systems/implicit_system.h index edc21e82970..504eb7cc670 100644 --- a/include/systems/implicit_system.h +++ b/include/systems/implicit_system.h @@ -308,6 +308,13 @@ class ImplicitSystem : public ExplicitSystem */ SparseMatrix & get_system_matrix(); + /** + * After calling this method with \p true, all successive solves + * will be restricted to unconstrained dofs, with constraints + * applied separately as necessary. + */ + virtual void restrict_solves_to_unconstrained(bool restricting); + /** * The system matrix. Implicit systems are characterized by * the need to solve the linear system Ax=b. This is the diff --git a/src/systems/diff_system.C b/src/systems/diff_system.C index 30566cce607..cdc79d95b49 100644 --- a/src/systems/diff_system.C +++ b/src/systems/diff_system.C @@ -176,6 +176,13 @@ std::pair DifferentiableSystem::get_linear_solve_parameters( } +void DifferentiableSystem::restrict_solves_to_unconstrained(bool restricting) { + // We need to restrict both the linear solver that's used for + // adjoint and sensitivity solves *and* the forward solver. + ImplicitSystem::restrict_solves_to_unconstrained(restricting); + this->time_solver->diff_solver()->restrict_solves_to_unconstrained(restricting); +} + void DifferentiableSystem::add_second_order_dot_vars() { diff --git a/src/systems/implicit_system.C b/src/systems/implicit_system.C index 4a73147de9e..1d8a42d34b1 100644 --- a/src/systems/implicit_system.C +++ b/src/systems/implicit_system.C @@ -108,6 +108,14 @@ void ImplicitSystem::disable_cache () { } +void ImplicitSystem::restrict_solves_to_unconstrained(bool restricting) { + if (restricting) + this->get_linear_solver()->restrict_solve_to_unconstrained(&this->get_dof_map()); + else + this->get_linear_solver()->restrict_solve_to_unconstrained(nullptr); +} + + std::pair ImplicitSystem::sensitivity_solve (const ParameterVector & parameters) @@ -152,7 +160,8 @@ ImplicitSystem::sensitivity_solve (const ParameterVector & parameters) totalrval.second += rval.second; } - // The linear solver may not have fit our constraints exactly + // The linear solver may not have fit our constraints exactly, or we + // might have disabled solves of constrained dofs #ifdef LIBMESH_ENABLE_CONSTRAINTS for (auto p : make_range(parameters.size())) this->get_dof_map().enforce_constraints_exactly @@ -203,7 +212,8 @@ ImplicitSystem::adjoint_solve (const QoISet & qoi_indices) totalrval.second += rval.second; } - // The linear solver may not have fit our constraints exactly + // The linear solver may not have fit our constraints exactly, or we + // might have disabled solves of constrained dofs #ifdef LIBMESH_ENABLE_CONSTRAINTS for (auto i : make_range(this->n_qois())) if (qoi_indices.has_index(i)) @@ -348,7 +358,8 @@ ImplicitSystem::weighted_sensitivity_adjoint_solve (const ParameterVector & para totalrval.second += rval.second; } - // The linear solver may not have fit our constraints exactly + // The linear solver may not have fit our constraints exactly, or we + // might have disabled solves of constrained dofs #ifdef LIBMESH_ENABLE_CONSTRAINTS for (auto i : make_range(this->n_qois())) if (qoi_indices.has_index(i)) @@ -432,7 +443,8 @@ ImplicitSystem::weighted_sensitivity_solve (const ParameterVector & parameters_i double(solver_params.second), solver_params.first); - // The linear solver may not have fit our constraints exactly + // The linear solver may not have fit our constraints exactly, or we + // might have disabled solves of constrained dofs #ifdef LIBMESH_ENABLE_CONSTRAINTS this->get_dof_map().enforce_constraints_exactly (*this, &this->get_weighted_sensitivity_solution(), From 4d8dedb4e10faf427df82c634b70f496e74d8b42 Mon Sep 17 00:00:00 2001 From: Roy Stogner Date: Wed, 14 Dec 2022 13:20:19 -0600 Subject: [PATCH 4/5] A little test coverage for restricted solves --- examples/adjoints/adjoints_ex3/adjoints_ex3.C | 19 +++++++++++++++++++ examples/adjoints/adjoints_ex3/run.sh | 2 +- 2 files changed, 20 insertions(+), 1 deletion(-) diff --git a/examples/adjoints/adjoints_ex3/adjoints_ex3.C b/examples/adjoints/adjoints_ex3/adjoints_ex3.C index e3c7813ab40..5ed745f7624 100644 --- a/examples/adjoints/adjoints_ex3/adjoints_ex3.C +++ b/examples/adjoints/adjoints_ex3/adjoints_ex3.C @@ -790,6 +790,17 @@ int main (int argc, char ** argv) equation_systems.init (); + if (infile.search("--restrict_solves_to_unconstrained")) + { + libMesh::out << "Restricting solves to unconstrained DoFs" << std::endl; + system.restrict_solves_to_unconstrained(true); + libmesh_assert_equal_to(system.get_linear_solver()->restrictions_to_unconstrained(), + &system.get_dof_map()); + libmesh_assert(system.get_time_solver().diff_solver()->get_restrict_solves_to_unconstrained()); + } + else + libMesh::out << "Not restricting solves to unconstrained DoFs" << std::endl; + // Print information about the mesh and system to the screen. mesh.print_info(); equation_systems.print_info(); @@ -830,9 +841,17 @@ int main (int argc, char ** argv) << std::endl << std::endl; + libmesh_assert_equal_to(system.get_linear_solver()->restrictions_to_unconstrained(), + &system.get_dof_map()); + libmesh_assert(system.get_time_solver().diff_solver()->get_restrict_solves_to_unconstrained()); + // Solve the forward system system.solve(); + libmesh_assert_equal_to(system.get_linear_solver()->restrictions_to_unconstrained(), + &system.get_dof_map()); + libmesh_assert(system.get_time_solver().diff_solver()->get_restrict_solves_to_unconstrained()); + // Write the output write_output(equation_systems, 0, a_step, "primal", param); diff --git a/examples/adjoints/adjoints_ex3/run.sh b/examples/adjoints/adjoints_ex3/run.sh index f3534b0ccd6..59a0378811e 100755 --- a/examples/adjoints/adjoints_ex3/run.sh +++ b/examples/adjoints/adjoints_ex3/run.sh @@ -8,7 +8,7 @@ example_name=adjoints_ex3 example_dir=examples/adjoints/$example_name -run_example "$example_name" $options +run_example "$example_name" --restrict_solves_to_unconstrained $options # This example needs better solvers on larger meshes + processor counts #benchmark_example 1 "$example_name" "max_adaptivesteps=18 -pc_type bjacobi -sub_pc_factor_levels 4" From 61c30df99caf49ad11187acc3e70236efdb578c3 Mon Sep 17 00:00:00 2001 From: Roy Stogner Date: Wed, 14 Dec 2022 17:33:20 -0600 Subject: [PATCH 5/5] PetscDiffSolver support for restricted solves --- include/solvers/petsc_diff_solver.h | 29 ++++- src/solvers/petsc_diff_solver.C | 171 +++++++++++++++++++++++----- 2 files changed, 168 insertions(+), 32 deletions(-) diff --git a/include/solvers/petsc_diff_solver.h b/include/solvers/petsc_diff_solver.h index 33f9300370a..622d395d9cf 100644 --- a/include/solvers/petsc_diff_solver.h +++ b/include/solvers/petsc_diff_solver.h @@ -91,11 +91,9 @@ class PetscDiffSolver : public DiffSolver * will be restricted to unconstrained dofs, with constraints * applied separately as necessary. */ - virtual void restrict_solves_to_unconstrained(bool restricting) override - { if (restricting) libmesh_not_implemented(); } + virtual void restrict_solves_to_unconstrained(bool restricting) override; - virtual bool get_restrict_solves_to_unconstrained() override - { return false; } + virtual bool get_restrict_solves_to_unconstrained() override; /** * This method performs a solve. What occurs in @@ -104,6 +102,24 @@ class PetscDiffSolver : public DiffSolver */ virtual unsigned int solve () override; + /* + * Scatter connecting full vectors with unconstrained subvectors + */ + WrappedPetsc scatter; + + /* + * Submatrix of Jacobian with unconstrained-unconstrained coupling + */ + WrappedPetsc submat; + + bool submat_created; + + /** + * PETSc index set containing unconstrained dofs on which to solve + * (\p nullptr means solve on all dofs). + */ + WrappedPetsc _unconstrained_dofs_is; + protected: /** @@ -120,6 +136,11 @@ class PetscDiffSolver : public DiffSolver #endif #endif + /** + * Are we restricting solves to unconstrained DoFs? + */ + bool _restrict_to_unconstrained; + private: /** diff --git a/src/solvers/petsc_diff_solver.C b/src/solvers/petsc_diff_solver.C index fa330db8fa7..14298cf252c 100644 --- a/src/solvers/petsc_diff_solver.C +++ b/src/solvers/petsc_diff_solver.C @@ -40,10 +40,10 @@ extern "C" // Function to hand to PETSc's SNES, // which monitors convergence at X PetscErrorCode - __libmesh_petsc_diff_solver_monitor (SNES snes, - PetscInt its, - PetscReal fnorm, - void * ctx) + _libmesh_petsc_diff_solver_monitor (SNES snes, + PetscInt its, + PetscReal fnorm, + void * ctx) { PetscDiffSolver & solver = *(static_cast (ctx)); @@ -84,7 +84,7 @@ extern "C" // Functions to hand to PETSc's SNES, // which compute the residual or jacobian at X PetscErrorCode - __libmesh_petsc_diff_solver_residual (SNES, Vec x, Vec r, void * ctx) + _libmesh_petsc_diff_solver_residual (SNES, Vec x, Vec r, void * ctx) { libmesh_assert(x); libmesh_assert(r); @@ -106,9 +106,20 @@ extern "C" // DiffSystem assembles from the solution and into the rhs, so swap // those with our input vectors before assembling. They'll probably // already be references to the same vectors, but PETSc might do - // something tricky. - X_input.swap(X_system); - R_input.swap(R_system); + // something tricky. ... or we might do something tricky. If + // we're solving only on constrained DoFs, we'll need to do some + // scatters here. + if (solver.get_restrict_solves_to_unconstrained()) + { + VecScatterBeginEnd(solver.comm(), solver.scatter, + X_input.vec(), X_system.vec(), + INSERT_VALUES, SCATTER_REVERSE); + } + else + { + X_input.swap(X_system); + R_input.swap(R_system); + } // We may need to localize a parallel solution sys.update(); @@ -121,8 +132,17 @@ extern "C" R_system.close(); // Swap back - X_input.swap(X_system); - R_input.swap(R_system); + if (solver.get_restrict_solves_to_unconstrained()) + { + VecScatterBeginEnd(solver.comm(), solver.scatter, + R_system.vec(), R_input.vec(), + INSERT_VALUES, SCATTER_FORWARD); + } + else + { + X_input.swap(X_system); + R_input.swap(R_system); + } // No errors, we hope return 0; @@ -130,11 +150,11 @@ extern "C" PetscErrorCode - __libmesh_petsc_diff_solver_jacobian (SNES, - Vec x, - Mat libmesh_dbg_var(j), - Mat pc, - void * ctx) + _libmesh_petsc_diff_solver_jacobian (SNES, + Vec x, + Mat libmesh_dbg_var(j), + Mat pc, + void * ctx) { libmesh_assert(x); libmesh_assert(j); @@ -159,9 +179,20 @@ extern "C" // DiffSystem assembles from the solution and into the jacobian, so // swap those with our input vectors before assembling. They'll // probably already be references to the same vectors, but PETSc - // might do something tricky. - X_input.swap(X_system); - J_input.swap(J_system); + // might do something tricky. ... or we might do something + // tricky. If we're solving only on constrained DoFs, we'll need + // to do some scatters here. + if (solver.get_restrict_solves_to_unconstrained()) + { + VecScatterBeginEnd(solver.comm(), solver.scatter, + X_input.vec(), X_system.vec(), + INSERT_VALUES, SCATTER_REVERSE); + } + else + { + X_input.swap(X_system); + J_input.swap(J_system); + } // We may need to localize a parallel solution sys.update(); @@ -174,8 +205,24 @@ extern "C" J_system.close(); // Swap back - X_input.swap(X_system); - J_input.swap(J_system); + if (solver.get_restrict_solves_to_unconstrained()) + { + PetscInt ierr = + LibMeshCreateSubMatrix(J_system.mat(), + solver._unconstrained_dofs_is, + solver._unconstrained_dofs_is, + solver.submat_created ? + MAT_REUSE_MATRIX : + MAT_INITIAL_MATRIX, + solver.submat.get()); + LIBMESH_CHKERR(ierr); + solver.submat_created = true; + } + else + { + X_input.swap(X_system); + J_input.swap(J_system); + } // No errors, we hope return 0; @@ -185,7 +232,8 @@ extern "C" PetscDiffSolver::PetscDiffSolver (sys_type & s) - : Parent(s) + : Parent(s), submat_created(false), + _restrict_to_unconstrained(false) { } @@ -280,6 +328,17 @@ DiffSolver::SolveResult convert_solve_result(SNESConvergedReason r) } +void PetscDiffSolver::restrict_solves_to_unconstrained(bool restricting) +{ + _restrict_to_unconstrained = restricting; +} + + +bool PetscDiffSolver::get_restrict_solves_to_unconstrained() +{ + return _restrict_to_unconstrained; +} + unsigned int PetscDiffSolver::solve() { @@ -292,20 +351,76 @@ unsigned int PetscDiffSolver::solve() PetscVector & r = *(cast_ptr *>(_system.rhs)); + WrappedPetsc subrhs, subsolution; + PetscErrorCode ierr = 0; - ierr = SNESSetFunction (_snes, r.vec(), - __libmesh_petsc_diff_solver_residual, this); - LIBMESH_CHKERR(ierr); + if (_restrict_to_unconstrained) + { + std::vector unconstrained_dofs; + const DofMap & dofmap = this->system().get_dof_map(); + for (dof_id_type i = dofmap.first_dof(), + end_i = dofmap.end_dof(); + i != end_i; ++i) + if (!dofmap.is_constrained_dof(i)) + unconstrained_dofs.push_back(cast_int(i)); + + const PetscInt is_local_size = + cast_int(unconstrained_dofs.size()); + + ierr = ISCreateGeneral(this->comm().get(), + cast_int(unconstrained_dofs.size()), + unconstrained_dofs.data(), PETSC_COPY_VALUES, + _unconstrained_dofs_is.get()); + + ierr = VecCreate(this->comm().get(), subrhs.get()); + LIBMESH_CHKERR(ierr); + ierr = VecSetSizes(subrhs, is_local_size, PETSC_DECIDE); + LIBMESH_CHKERR(ierr); + ierr = VecSetFromOptions(subrhs); + LIBMESH_CHKERR(ierr); - ierr = SNESSetJacobian (_snes, jac.mat(), jac.mat(), - __libmesh_petsc_diff_solver_jacobian, this); + ierr = VecCreate(this->comm().get(), subsolution.get()); + LIBMESH_CHKERR(ierr); + ierr = VecSetSizes(subsolution, is_local_size, PETSC_DECIDE); + LIBMESH_CHKERR(ierr); + ierr = VecSetFromOptions(subsolution); + LIBMESH_CHKERR(ierr); + + ierr = VecScatterCreate(r.vec(), _unconstrained_dofs_is, + subrhs, nullptr, scatter.get()); + LIBMESH_CHKERR(ierr); + + VecScatterBeginEnd(this->comm(), scatter, r.vec(), subrhs, INSERT_VALUES, SCATTER_FORWARD); + VecScatterBeginEnd(this->comm(), scatter, x.vec(), subsolution, INSERT_VALUES, SCATTER_FORWARD); + + LIBMESH_CHKERR(ierr); + + ierr = SNESSetFunction (_snes, subrhs, + _libmesh_petsc_diff_solver_residual, this); + LIBMESH_CHKERR(ierr); + + ierr = SNESSetJacobian (_snes, submat, submat, + _libmesh_petsc_diff_solver_jacobian, this); + } + else + { + ierr = SNESSetFunction (_snes, r.vec(), + _libmesh_petsc_diff_solver_residual, this); + LIBMESH_CHKERR(ierr); + + ierr = SNESSetJacobian (_snes, jac.mat(), jac.mat(), + _libmesh_petsc_diff_solver_jacobian, this); + } LIBMESH_CHKERR(ierr); ierr = SNESSetFromOptions(_snes); LIBMESH_CHKERR(ierr); - ierr = SNESSolve (_snes, PETSC_NULL, x.vec()); + if (_restrict_to_unconstrained) + ierr = SNESSolve (_snes, PETSC_NULL, x.vec()); + else + ierr = SNESSolve (_snes, PETSC_NULL, x.vec()); LIBMESH_CHKERR(ierr); #ifdef LIBMESH_ENABLE_CONSTRAINTS @@ -334,7 +449,7 @@ void PetscDiffSolver::setup_petsc_data() ierr = SNESCreate(this->comm().get(), _snes.get()); LIBMESH_CHKERR(ierr); - ierr = SNESMonitorSet (_snes, __libmesh_petsc_diff_solver_monitor, + ierr = SNESMonitorSet (_snes, _libmesh_petsc_diff_solver_monitor, this, PETSC_NULL); LIBMESH_CHKERR(ierr);