Skip to content

Remove NewtonSolver from hyperelasticity demo #3929

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
wants to merge 8 commits into
base: main
Choose a base branch
from
Open

Remove NewtonSolver from hyperelasticity demo #3929

Changes from all commits
Commits
Show all changes
8 commits
Select commit Hold shift + click to select a range
29d9376
Use userspace solver
chrisrichardson Sep 17, 2025
4845826
More trimming
chrisrichardson Sep 17, 2025
a1a1a51
Incorporate solver into main.cpp
chrisrichardson Sep 17, 2025
bf2fefc
Reduce class variables
chrisrichardson Sep 17, 2025
1038836
Simplify
chrisrichardson Sep 17, 2025
877f583
Simplify
chrisrichardson Sep 17, 2025
83cb6c3
Remove external files
chrisrichardson Sep 17, 2025
55ede5c
Merge branch 'main' into chris/hyperelasticity-demo
chrisrichardson Sep 19, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
199 changes: 138 additions & 61 deletions cpp/demo/hyperelasticity/main.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -28,7 +28,6 @@
#include <dolfinx/la/petsc.h>
#include <dolfinx/mesh/Mesh.h>
#include <dolfinx/mesh/cell_types.h>
#include <dolfinx/nls/NewtonSolver.h>
#include <numbers>
#include <petscmat.h>
#include <petscsys.h>
Expand All @@ -47,9 +46,10 @@ class HyperElasticProblem
HyperElasticProblem(fem::Form<T>& L, fem::Form<T>& J,
const std::vector<fem::DirichletBC<T>>& bcs)
: _l(L), _j(J), _bcs(bcs.begin(), bcs.end()),
_b(L.function_spaces()[0]->dofmap()->index_map,
L.function_spaces()[0]->dofmap()->index_map_bs()),
_matA(la::petsc::Matrix(fem::petsc::create_matrix(J, "aij"), false))
_b_vec(L.function_spaces()[0]->dofmap()->index_map,
L.function_spaces()[0]->dofmap()->index_map_bs()),
_matJ(la::petsc::Matrix(fem::petsc::create_matrix(J, "aij"), false)),
_solver(L.function_spaces()[0]->dofmap()->index_map->comm())
{
auto map = L.function_spaces()[0]->dofmap()->index_map;
const int bs = L.function_spaces()[0]->dofmap()->index_map_bs();
Expand All @@ -59,82 +59,163 @@ class HyperElasticProblem
std::int64_t size_global = bs * map->size_global();
VecCreateGhostBlockWithArray(map->comm(), bs, size_local, size_global,
ghosts.size(), ghosts.data(),
_b.array().data(), &_b_petsc);
_b_vec.array().data(), &_b);

// Create linear solver. Default to LU.
_solver.set_options_prefix("nls_solve_");
la::petsc::options::set("nls_solve_ksp_type", "preonly");
la::petsc::options::set("nls_solve_pc_type", "lu");
_solver.set_from_options();
}

/// Destructor
virtual ~HyperElasticProblem()
{
if (_b_petsc)
VecDestroy(&_b_petsc);
assert(_b);
VecDestroy(&_b);
}

/// @brief Form
/// @return
auto form()
/// @brief Newton Solver
/// @param x Solution vector
/// @return Iteration count and flag indicating convergence
std::pair<int, bool> solve(Vec x)
{
return [](Vec x)
int iteration = 0;
PetscReal residual0 = 0;

auto converged
= [&iteration, &residual0, this](const Vec r) -> std::pair<double, bool>
{
VecGhostUpdateBegin(x, INSERT_VALUES, SCATTER_FORWARD);
VecGhostUpdateEnd(x, INSERT_VALUES, SCATTER_FORWARD);
PetscReal residual = 0;
VecNorm(r, NORM_2, &residual);
Copy link
Contributor

@michalhabera michalhabera Sep 26, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Should we use NORM_INFINITY for residual norms instead (as a naive fix)? I always struggle that NORM_2 grows with numbers of elements so for a finer mesh you'll violate atol very soon.

Btw. are C++ demos run in CI? Could not find the job.

Copy link
Contributor

@schnellerhase schnellerhase Sep 27, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Btw. are C++ demos run in CI? Could not find the job.

https://github.com/FEniCS/dolfinx/blob/main/.github/workflows/ccpp.yml#L230 - CI step has a confusing name.


// Relative residual
const double relative_residual = residual / residual0;

// Output iteration number and residual
spdlog::info("Newton iteration {}"
": r (abs) = {} (tol = {}), r (rel) = {} (tol = {})",
iteration, residual, atol, relative_residual, rtol);

// Return true if convergence criterion is met
if (relative_residual < rtol or residual < atol)
return {residual, true};
else
return {residual, false};
};

assert(_b);
F(x);

auto [residual, newton_converged] = converged(_b);

_solver.set_operators(_matJ.mat(), _matJ.mat());

Vec dx;
MatCreateVecs(_matJ.mat(), &dx, nullptr);

int max_it = 50;
int krylov_iterations = 0;

// Start iterations
while (!newton_converged and iteration < max_it)
{
// Compute Jacobian
assert(_matJ.mat());
J(x, _matJ.mat());

// Perform linear solve and update total number of Krylov iterations
krylov_iterations += _solver.solve(dx, _b);

// Update solution
double relaxation_parameter = 1.0;
VecAXPY(x, -relaxation_parameter, dx);

// Increment iteration count
++iteration;

// Compute F
F(x);

// Initialize residual0
if (iteration == 1)
VecNorm(dx, NORM_2, &residual0);

// Test for convergence
std::tie(residual, newton_converged) = converged(_b);
}

if (newton_converged)
{
spdlog::info("Newton solver finished in {} iterations and {} linear "
"solver iterations.",
iteration, krylov_iterations);
}
else
throw std::runtime_error("Newton solver did not converge.");

VecDestroy(&dx);

return {iteration, newton_converged};
}

/// Compute F at current point x
auto F()
void F(const Vec x)
{
return [&](const Vec x, Vec)
{
// Assemble b and update ghosts
std::span b(_b.array());
std::ranges::fill(b, 0);
fem::assemble_vector(b, _l);
VecGhostUpdateBegin(_b_petsc, ADD_VALUES, SCATTER_REVERSE);
VecGhostUpdateEnd(_b_petsc, ADD_VALUES, SCATTER_REVERSE);

// Set bcs
Vec x_local;
VecGhostGetLocalForm(x, &x_local);
PetscInt n = 0;
VecGetSize(x_local, &n);
const T* _x = nullptr;
VecGetArrayRead(x_local, &_x);
std::ranges::for_each(_bcs, [b, x = std::span(_x, n)](auto& bc)
{ bc.get().set(b, x, -1); });
VecRestoreArrayRead(x_local, &_x);
};
VecGhostUpdateBegin(x, INSERT_VALUES, SCATTER_FORWARD);
VecGhostUpdateEnd(x, INSERT_VALUES, SCATTER_FORWARD);

// Assemble b and update ghosts
std::span b(_b_vec.array());
std::ranges::fill(b, 0);
fem::assemble_vector(b, _l);
VecGhostUpdateBegin(_b, ADD_VALUES, SCATTER_REVERSE);
VecGhostUpdateEnd(_b, ADD_VALUES, SCATTER_REVERSE);

// Set bcs
Vec x_local;
VecGhostGetLocalForm(x, &x_local);
PetscInt n = 0;
VecGetSize(x_local, &n);
const T* _x = nullptr;
VecGetArrayRead(x_local, &_x);
std::ranges::for_each(_bcs, [b, x = std::span(_x, n)](auto& bc)
{ bc.get().set(b, x, -1); });
VecRestoreArrayRead(x_local, &_x);
}

/// Compute J = F' at current point x
auto J()
void J(const Vec, Mat A)
{
return [&](const Vec, Mat A)
{
MatZeroEntries(A);
fem::assemble_matrix(la::petsc::Matrix::set_block_fn(A, ADD_VALUES), _j,
_bcs);
MatAssemblyBegin(A, MAT_FLUSH_ASSEMBLY);
MatAssemblyEnd(A, MAT_FLUSH_ASSEMBLY);
fem::set_diagonal(la::petsc::Matrix::set_fn(A, INSERT_VALUES),
*_j.function_spaces()[0], _bcs);
MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
};
MatZeroEntries(A);
fem::assemble_matrix(la::petsc::Matrix::set_block_fn(A, ADD_VALUES), _j,
_bcs);
MatAssemblyBegin(A, MAT_FLUSH_ASSEMBLY);
MatAssemblyEnd(A, MAT_FLUSH_ASSEMBLY);
fem::set_diagonal(la::petsc::Matrix::set_fn(A, INSERT_VALUES),
*_j.function_spaces()[0], _bcs);
MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
}

/// RHS vector
Vec vector() { return _b_petsc; }
/// @brief Relative convergence tolerance.
double rtol = 1e-9;

/// Jacobian matrix
Mat matrix() { return _matA.mat(); }
/// @brief Absolute convergence tolerance.
double atol = 1e-10;

private:
fem::Form<T>& _l;
fem::Form<T>& _j;
std::vector<std::reference_wrapper<const fem::DirichletBC<T>>> _bcs;
la::Vector<T> _b;
Vec _b_petsc = nullptr;
la::petsc::Matrix _matA;
la::Vector<T> _b_vec;
Vec _b = nullptr;

// Jacobian matrix
la::petsc::Matrix _matJ;

// Linear solver
dolfinx::la::petsc::KrylovSolver _solver;
};

int main(int argc, char* argv[])
Expand Down Expand Up @@ -247,15 +328,11 @@ int main(int argc, char* argv[])
fem::DirichletBC<T>(u_rotation, bdofs_right)};

HyperElasticProblem problem(L, a, bcs);
nls::petsc::NewtonSolver newton_solver(mesh->comm());
newton_solver.setF(problem.F(), problem.vector());
newton_solver.setJ(problem.J(), problem.matrix());
newton_solver.set_form(problem.form());
newton_solver.rtol = 10 * std::numeric_limits<T>::epsilon();
newton_solver.atol = 10 * std::numeric_limits<T>::epsilon();
problem.rtol = 10 * std::numeric_limits<T>::epsilon();
problem.atol = 10 * std::numeric_limits<T>::epsilon();

la::petsc::Vector _u(la::petsc::create_vector_wrap(*u->x()), false);
auto [niter, success] = newton_solver.solve(_u.vec());
auto [niter, success] = problem.solve(_u.vec());
std::cout << "Number of Newton iterations: " << niter << std::endl;

// Compute Cauchy stress. Construct appropriate Basix element for
Expand Down
Loading