Dear all,
I am currently trying to implement the PETSc SNES solver
(https://petsc.org/release/docs/manual/snes/) into my code.
As a starting point, I am trying to implement it into the step-15 from the
tutorials.
Attached can be found a MWE of the current version of the code which
follows an older implementation
(https://github.com/gujans/dealii-snes/blob/master/step-15.cc).
The code performs the following steps:
1 - create the mesh;
2 - initialize the PETScWrappers::MPI matrix/vectors;
3 - set the boundary values;
4 - initialize the SNES solver;
5 - assemble the FormFunction and FormJacobian using the assemble_rhs and
assemble_system as functions to set the values to the system matrix and
rhs;
6 - solve the problem using the SNES solver.
At the moment the code runs but it produces a solution which is constant
zero, moreover the SNES iteration count returns zero as well.
Am I missing something while applying the boundary conditions?
How can I pass the correct assembled matrices to the SNES Solver?
To be noted that if I pass to the solver the PETScWrappers::MPI::Vector
present_solution, which is initialized during the setup_system(), instead
of a Vec test_x which i created just to store the solution from the solver,
it returns the following error:
[0]PETSC ERROR: ----- Error Message --------
[0]PETSC ERROR: Object is in wrong state
[0]PETSC ERROR: Not for unassembled vector
Lastly, I also found an old post in the group
(https://groups.google.com/g/dealii/c/pg9bj_z8dbk/m/OTtfdCY0AAAJ) which
hinted another possible approach to the one i am using, however i am not
sure how to extract PETSc Vec object from PETScWrappers::MPI::Vector
without passing from a PETScWrapper::VectorBase.
Thanks in advance,
Lorenzo
--
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see
https://groups.google.com/d/forum/dealii?hl=en
---
You received this message because you are subscribed to the Google Groups
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email
to dealii+unsubscr...@googlegroups.com.
To view this discussion on the web visit
https://groups.google.com/d/msgid/dealii/6a940fd0-b6af-4e1e-9746-e72100b586d5n%40googlegroups.com.
// Base step 15 include files
#include <deal.II/base/quadrature_lib.h>
#include <deal.II/base/function.h>
#include <deal.II/base/logstream.h>
#include <deal.II/base/utilities.h>
#include <deal.II/lac/vector.h>
#include <deal.II/lac/full_matrix.h>
#include <deal.II/lac/sparse_matrix.h>
#include <deal.II/lac/dynamic_sparsity_pattern.h>
#include <deal.II/lac/solver_cg.h>
#include <deal.II/lac/precondition.h>
#include <deal.II/lac/affine_constraints.h>
#include <deal.II/grid/tria.h>
#include <deal.II/grid/grid_generator.h>
#include <deal.II/grid/tria_accessor.h>
#include <deal.II/grid/tria_iterator.h>
#include <deal.II/grid/manifold_lib.h>
#include <deal.II/grid/grid_refinement.h>
#include <deal.II/dofs/dof_handler.h>
#include <deal.II/dofs/dof_accessor.h>
#include <deal.II/dofs/dof_tools.h>
#include <deal.II/fe/fe_values.h>
#include <deal.II/fe/fe_q.h>
#include <deal.II/numerics/vector_tools.h>
#include <deal.II/numerics/matrix_tools.h>
#include <deal.II/numerics/data_out.h>
#include <deal.II/numerics/error_estimator.h>
#include <fstream>
#include <iostream>
// Added include files for the SNES version
#include <deal.II/base/conditional_ostream.h>
#include <deal.II/base/mpi.h>
#include <deal.II/lac/petsc_vector.h>
#include <deal.II/lac/petsc_sparse_matrix.h>
#include <deal.II/lac/petsc_solver.h>
#include <deal.II/lac/petsc_precondition.h>
#include <deal.II/grid/grid_tools.h>
#include <deal.II/grid/grid_out.h>
#include <deal.II/dofs/dof_renumbering.h>
#include <petscsnes.h>
#include <deal.II/numerics/solution_transfer.h>
#include <deal.II/distributed/solution_transfer.h>
namespace Step15
{
// Setup Classe
using namespace dealii;
template <int dim>
class MinimalSurfaceProblem
{
public:
MinimalSurfaceProblem();
~MinimalSurfaceProblem();
void run();
private:
void setup_system(const bool initial_step);
void set_boundary_values();
void assemble_system(PETScWrappers::MPI::Vector &present_solution);
void assemble_rhs(PETScWrappers::MPI::Vector &present_solution,
PETScWrappers::MPI::Vector &system_rhs);
void refine_mesh();
static PetscErrorCode FormFunction(SNES snes, Vec x, Vec f, void* ctx);
static PetscErrorCode FormJacobian(SNES snes, Vec x, Mat jac, Mat B, void*
ctx);
MPI_Comm mpi_communicator;
const unsigned int n_mpi_processes;
const unsigned int this_mpi_process;
IndexSet locally_owned_dofs;
ConditionalOStream pcout;
Triangulation<dim> triangulation;
DoFHandler<dim> dof_handler;
FE_Q<dim> fe;
AffineConstraints<double> hanging_node_constraints;
SparsityPattern sparsity_pattern;
PETScWrappers::MPI::SparseMatrix system_matrix;
PETScWrappers::MPI::Vector present_solution;
PETScWrappers::MPI::Vector system_rhs;
Vec test_x;
};
// Costructor
template <int dim>
MinimalSurfaceProblem<dim>::MinimalSurfaceProblem()
: mpi_communicator(MPI_COMM_WORLD)
, n_mpi_processes(Utilities::MPI::n_mpi_processes(mpi_communicator))
, this_mpi_process(Utilities::MPI::this_mpi_process(mpi_communicator))
, pcout(std::cout, (this_mpi_process == 0))
, dof_handler(triangulation)
, fe(2)
{}
// Destructor
template <int dim>
MinimalSurfaceProblem<dim>::~MinimalSurfaceProblem()
{
dof_handler.clear();
}
// Boundary values
template <int dim>
class BoundaryValues : public Function<dim>
{
public:
virtual double value(const Point<dim>& p,
const unsigned int component = 0) const override;
};
template <int dim>
double BoundaryValues<dim>::value(const Point<dim>& p,
const unsigned int /*component*/) const
{
return std::sin(2 * numbers::PI * (p[0] + p[1]));
}
// Function used to create the FormFunction Input for the SNESSetFunction
template <int dim>
PetscErrorCode MinimalSurfaceProblem<dim>::FormFunction(SNES snes, Vec x, Vec
f, void* ctx)
{
auto p_ctx = reinterpret_cast<MinimalSurfaceProblem<dim>*>(ctx);
PETScWrappers::VectorBase xx(x);
PETScWrappers::VectorBase ff(f);
PETScWrappers::MPI::Vector x_wrap(p_ctx->mpi_communicator, xx,
xx.local_size());
PETScWrappers::MPI::Vector f_wrap(p_ctx->mpi_communicator, ff,
ff.local_size());
p_ctx->assemble_rhs(x_wrap, f_wrap);
return 0;
}
// Function used to create the FormFunction Input for the SNESSetJacobian
template <int dim>
PetscErrorCode MinimalSurfaceProblem<dim>::FormJacobian(SNES snes, Vec x, Mat
jac, Mat B, void* ctx)
{
auto p_ctx = reinterpret_cast<MinimalSurfaceProblem<dim>*>(ctx);
PETScWrappers::VectorBase xx(x);
PETScWrappers::MPI::Vector x_wrap(p_ctx->mpi_communicator, xx,
xx.local_size());
p_ctx->assemble_system(x_wrap);
/*
Assemble matrix
*/
PetscErrorCode ierr;
ierr = MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
ierr = MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
if (jac != B) {
ierr = MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
ierr = MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
}
return 0;
}
// Setup System
template <int dim>
void MinimalSurfaceProblem<dim>::setup_system(const bool initial_step)
{
GridTools::partition_triangulation(n_mpi_processes, triangulation);
if (initial_step)
{
dof_handler.distribute_dofs(fe);
DoFRenumbering::subdomain_wise(dof_handler);
const std::vector<IndexSet> locally_owned_dofs_per_proc =
DoFTools::locally_owned_dofs_per_subdomain(dof_handler);
locally_owned_dofs = locally_owned_dofs_per_proc[this_mpi_process];
present_solution.reinit(locally_owned_dofs, mpi_communicator);
}
const std::vector<IndexSet> locally_owned_dofs_per_proc =
DoFTools::locally_owned_dofs_per_subdomain(dof_handler);
locally_owned_dofs = locally_owned_dofs_per_proc[this_mpi_process];
hanging_node_constraints.clear();
DoFTools::make_hanging_node_constraints(dof_handler,
hanging_node_constraints);
VectorTools::interpolate_boundary_values(dof_handler,
0,
ZeroFunction<dim>(),
hanging_node_constraints);
hanging_node_constraints.close();
system_rhs.reinit(locally_owned_dofs, mpi_communicator);
DynamicSparsityPattern dsp(dof_handler.n_dofs());
DoFTools::make_sparsity_pattern(dof_handler, dsp, hanging_node_constraints,
false);
hanging_node_constraints.condense(dsp);
sparsity_pattern.copy_from(dsp);
system_matrix.reinit(locally_owned_dofs,
locally_owned_dofs,
sparsity_pattern,
mpi_communicator);
}
// Assemble bilinear form
template <int dim>
void MinimalSurfaceProblem<dim>::assemble_system(PETScWrappers::MPI::Vector
&present_solution)
{
const QGauss<dim> quadrature_formula(fe.degree + 1);
system_matrix = 0;
FEValues<dim> fe_values(fe,
quadrature_formula,
update_gradients |
update_quadrature_points |
update_JxW_values);
const unsigned int dofs_per_cell = fe.dofs_per_cell;
const unsigned int n_q_points = quadrature_formula.size();
FullMatrix<double> cell_matrix(dofs_per_cell, dofs_per_cell);
std::vector<Tensor<1, dim>> old_solution_gradients(n_q_points);
std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
for (const auto& cell : dof_handler.active_cell_iterators())
if (cell->subdomain_id() == this_mpi_process)
{
cell_matrix = 0;
fe_values.reinit(cell);
fe_values.get_function_gradients(present_solution,
old_solution_gradients);
for (unsigned int q = 0; q < n_q_points; ++q)
{
const double coeff = 1.0 / std::sqrt(1 +
old_solution_gradients[q] * old_solution_gradients[q]);
for (unsigned int i = 0; i < dofs_per_cell; ++i)
{
for (unsigned int j = 0; j < dofs_per_cell; ++j)
cell_matrix(i, j) +=
(((fe_values.shape_grad(i, q) // ((\nabla \phi_i
* coeff // * a_n
* fe_values.shape_grad(j, q)) // * \nabla
\phi_j)
- // -
(fe_values.shape_grad(i, q) // (\nabla \phi_i
* coeff * coeff * coeff // * a_n^3
* (fe_values.shape_grad(j, q) // * (\nabla
\phi_j
* old_solution_gradients[q]) // *
\nabla u_n)
* old_solution_gradients[q])) // * \nabla
u_n)))
* fe_values.JxW(q)); // * dx
}
}
cell->get_dof_indices(local_dof_indices);
hanging_node_constraints.distribute_local_to_global(cell_matrix,
local_dof_indices,
system_matrix);
}
system_matrix.compress(VectorOperation::add);
std::map<types::global_dof_index, double> boundary_values;
VectorTools::interpolate_boundary_values(dof_handler,
0,
Functions::ZeroFunction<dim>(),
boundary_values);
MatrixTools::apply_boundary_values(boundary_values,
system_matrix,
present_solution,
system_rhs);
}
// Assemble RHS
template <int dim>
void MinimalSurfaceProblem<dim>::assemble_rhs(PETScWrappers::MPI::Vector
&present_solution, PETScWrappers::MPI::Vector &system_rhs)
{
const QGauss<dim> quadrature_formula(fe.degree + 1);
system_rhs = 0;
FEValues<dim> fe_values(fe,
quadrature_formula,
update_gradients | update_quadrature_points |
update_JxW_values);
const unsigned int dofs_per_cell = fe.dofs_per_cell;
const unsigned int n_q_points = quadrature_formula.size();
FullMatrix<double> cell_matrix(dofs_per_cell, dofs_per_cell);
Vector<double> cell_rhs(dofs_per_cell);
std::vector<Tensor<1, dim>> old_solution_gradients(n_q_points);
std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
for (const auto& cell : dof_handler.active_cell_iterators())
if (cell->subdomain_id() == this_mpi_process)
{
cell_rhs = 0;
fe_values.reinit(cell);
fe_values.get_function_gradients(present_solution,
old_solution_gradients);
for (unsigned int q = 0; q < n_q_points; ++q)
{
const double coeff = 1.0 / std::sqrt(1 +
old_solution_gradients[q] * old_solution_gradients[q]);
for (unsigned int i = 0; i < dofs_per_cell; ++i)
{
cell_rhs(i) -= (fe_values.shape_grad(i, q) // \nabla \phi_i
* coeff // * a_n
* old_solution_gradients[q] // * u_n
* fe_values.JxW(q)); // * dx
}
}
cell->get_dof_indices(local_dof_indices);
hanging_node_constraints.distribute_local_to_global(cell_rhs,
local_dof_indices,
system_rhs);
}
system_rhs.compress(VectorOperation::add);
}
// Setup boundary values
template <int dim>
void MinimalSurfaceProblem<dim>::set_boundary_values()
{
std::map<types::global_dof_index, double> boundary_values;
VectorTools::interpolate_boundary_values(dof_handler,
0,
BoundaryValues<dim>(),
boundary_values);
for (auto &boundary_value : boundary_values)
present_solution(boundary_value.first) = boundary_value.second;
}
// Run
template <int dim>
void MinimalSurfaceProblem<dim>::run()
{
GridGenerator::hyper_ball(triangulation);
triangulation.refine_global(2);
setup_system(true);
set_boundary_values();
// SNES solver
SNES snes;
PetscErrorCode ierr;
ierr = SNESCreate(mpi_communicator, &snes);
ierr = VecCreate(mpi_communicator,&test_x);
ierr = VecSetSizes(test_x,PETSC_DECIDE,present_solution.size());
ierr = VecSetFromOptions(test_x);
ierr = SNESSetFunction(snes, system_rhs, FormFunction, this);
ierr = SNESSetJacobian(snes, system_matrix, system_matrix, FormJacobian,
this);
ierr = SNESSetFromOptions(snes);
ierr = SNESSolve(snes, NULL, test_x);
PETScWrappers::VectorBase test_xx(test_x);
PETScWrappers::MPI::Vector output_solution(mpi_communicator, test_xx,
test_xx.local_size());
PetscInt its;
SNESGetIterationNumber(snes, &its);
PetscPrintf(MPI_COMM_WORLD, "Number of SNES iterations = %D\n", its);
ierr = SNESDestroy(&snes);
ierr = VecDestroy(&test_x);
DataOut<dim> data_out;
data_out.attach_dof_handler(dof_handler);
data_out.add_data_vector(output_solution, "solution");
data_out.build_patches();
const std::string filename = "solution.vtk";
std::ofstream output(filename.c_str());
data_out.write_vtk(output);
}
} // namespace Step15
// @sect4{The main function}
int main(int argc, char** argv)
{
try
{
using namespace dealii;
using namespace Step15;
Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1);
deallog.depth_console(0);
MinimalSurfaceProblem<2> laplace_problem_2d;
laplace_problem_2d.run();
}
catch (std::exception &exc)
{
std::cerr << std::endl
<< std::endl
<< "----------------------------------------------------"
<< std::endl;
std::cerr << "Exception on processing: " << std::endl
<< exc.what() << std::endl
<< "Aborting!" << std::endl
<< "----------------------------------------------------"
<< std::endl;
return 1;
}
catch (...)
{
std::cerr << std::endl
<< std::endl
<< "----------------------------------------------------"
<< std::endl;
std::cerr << "Unknown exception!" << std::endl
<< "Aborting!" << std::endl
<< "----------------------------------------------------"
<< std::endl;
return 1;
}
return 0;
}
