Hi all,
I've encountered a residual jump at each restart step within the inner
iteration of a restarted GMRES while using geometrical multigrid
precondition and matrix free operator.
To simplify the test, I modify the step-66
<https://dealii.org/current/doxygen/deal.II/step_66.html> to solve the
Gelfand problem using both *GMG + Matrix-free* method and *AMG + sparse
matrix* respectively with a restarted GMRES(4) (*max_n_tmp_vectors=6*).
For the GMG case, there're clearly two jumps within one Newton step solving:
DEAL:GMRES::Check 0 60.1874
DEAL:GMRES::Starting value 60.1874
DEAL:GMRES::Check 1 0.0551893
DEAL:GMRES::Check 2 0.000396485
DEAL:GMRES::Check 3 5.95803e-06
DEAL:GMRES::Check 4 1.09095e-07
DEAL:GMRES::Check 4 1.00888e-05 <-- JUMP HERE
DEAL:GMRES::Check 5 2.19635e-07
DEAL:GMRES::Check 6 3.17790e-09
DEAL:GMRES::Check 7 6.61810e-11
DEAL:GMRES::Check 8 1.82112e-12
DEAL:GMRES::Check 8 6.59447e-12 <-- JUMP HERE
DEAL:GMRES::Check 9 1.15358e-13
DEAL:GMRES::Convergence step 9 value 1.15358e-13
And for the AMG case, the residual goes all the way down normally:
DEAL:GMRES::Check 0 25.2492
DEAL:GMRES::Starting value 25.2492
DEAL:GMRES::Check 1 1.02964
DEAL:GMRES::Check 2 0.464888
DEAL:GMRES::Check 3 0.263190
DEAL:GMRES::Check 4 0.0703056
DEAL:GMRES::Check 4 0.0703056
DEAL:GMRES::Check 5 0.0580013
DEAL:GMRES::Check 6 0.0363219
DEAL:GMRES::Check 7 0.0330487
DEAL:GMRES::Check 8 0.0222448
DEAL:GMRES::Check 8 0.0222448
DEAL:GMRES::Check 9 0.0200449
DEAL:GMRES::Check 10 0.0173175
DEAL:GMRES::Check 11 0.0106881
DEAL:GMRES::Check 12 0.00388264
DEAL:GMRES::Check 12 0.00388264
DEAL:GMRES::Check 13 0.00140426
DEAL:GMRES::Check 14 0.00116016
DEAL:GMRES::Check 15 0.000288553
DEAL:GMRES::Check 16 0.000219209
DEAL:GMRES::Check 16 0.000219209
DEAL:GMRES::Check 17 0.000169839
DEAL:GMRES::Check 18 6.81943e-05
DEAL:GMRES::Check 19 6.79830e-05
DEAL:GMRES::Check 20 1.91041e-05
...
Complete test code & log file is attached below, the first half of the log
file is *GMG* and the latter is *AMG.*
I'm not an expert in the linear solvers but GMRES(m) residual should
remains unchanged during each restart since all its value are untouched,
just throw away the Arnoldi basis? And thus that seem to me like a
misbehave related to *PreconditionMG*.
Does anyone have idea on how this happens or why it should just jump if I
get it wrong. Thanks!
Regards,
Chengjiang Yin
P.S. I've already opened an issue
<https://github.com/dealii/dealii/issues/13297> about this on the github
repo and no reply has been received till now. Sorry for the duplicate and I
will keep update on both side
--
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 [email protected].
To view this discussion on the web visit
https://groups.google.com/d/msgid/dealii/fc0137c1-a480-4a67-b7e5-ae78fdf798f6n%40googlegroups.com.
================================================================================
START DATE: 2022/1/27, TIME: 6:06:48
--------------------------------------------------------------------------------
Running with 32 MPI processes
Vectorization over 2 doubles = 128 bits (SSE2), VECTORIZATION_LEVEL=1
Finite element space: FE_Q<3>(2)
================================================================================
JobId fe0cba80120c Thu Jan 27 06:06:48 2022
--------------------------------------------------------------------------------
Cycle 0
--------------------------------------------------------------------------------
Set up system...
Triangulation: 28672 cells
DoFHandler: 232609 DoFs
Solve using Newton's method...
DEAL:GMRES::Check 0 60.1874
DEAL:GMRES::Starting value 60.1874
DEAL:GMRES::Check 1 0.0551893
DEAL:GMRES::Check 2 0.000396485
DEAL:GMRES::Check 3 5.95803e-06
DEAL:GMRES::Check 4 1.09095e-07
DEAL:GMRES::Check 4 1.00888e-05
DEAL:GMRES::Check 5 2.19635e-07
DEAL:GMRES::Check 6 3.17790e-09
DEAL:GMRES::Check 7 6.61810e-11
DEAL:GMRES::Check 8 1.82112e-12
DEAL:GMRES::Check 8 6.59447e-12
DEAL:GMRES::Check 9 1.15358e-13
DEAL:GMRES::Convergence step 9 value 1.15358e-13
Nstep 1, errf = 4.85787e-05, errx = 60.2426, it = 9
DEAL:GMRES::Check 0 0.517403
DEAL:GMRES::Starting value 0.517403
DEAL:GMRES::Check 1 0.000713779
DEAL:GMRES::Check 2 4.69099e-06
DEAL:GMRES::Check 3 7.05096e-08
DEAL:GMRES::Check 4 1.62490e-09
DEAL:GMRES::Check 4 8.61865e-08
DEAL:GMRES::Check 5 1.87820e-09
DEAL:GMRES::Check 6 2.64609e-11
DEAL:GMRES::Check 7 5.05381e-13
DEAL:GMRES::Convergence step 7 value 5.05381e-13
Nstep 2, errf = 3.28133e-09, errx = 0.518129, it = 7
DEAL:GMRES::Check 0 4.01821e-05
DEAL:GMRES::Starting value 4.01821e-05
DEAL:GMRES::Check 1 5.29173e-08
DEAL:GMRES::Check 2 3.46129e-10
DEAL:GMRES::Check 3 5.27488e-12
DEAL:GMRES::Check 4 1.38412e-13
DEAL:GMRES::Convergence step 4 value 1.38412e-13
Nstep 3, errf = 6.35346e-13, errx = 4.02406e-05, it = 4
Convergence step 3 value 6.35346e-13 (used wall time: 0.201269 s)
Time for setup+solve (CPU/Wall) 0.462676/0.479048 s
Output results...
H1 seminorm: 0.588667
+---------------------------------------------+------------+------------+
| Total wallclock time elapsed since start | 0.951s | |
| | | |
| Section | no. calls | wall time | % of total |
+---------------------------------+-----------+------------+------------+
| assemble right hand side | 3 | 0.00257s | 0.27% |
| compute residual | 3 | 0.00313s | 0.33% |
| compute update | 3 | 0.194s | 20.% |
| make grid | 1 | 0.282s | 30.% |
| setup system | 1 | 0.277s | 29.% |
| solve | 1 | 0.201s | 21.% |
+---------------------------------+-----------+------------+------------+
================================================================================
START DATE: 2022/1/27, TIME: 6:06:49
--------------------------------------------------------------------------------
Running with 32 MPI processes
Vectorization over 2 doubles = 128 bits (SSE2), VECTORIZATION_LEVEL=1
Finite element space: FE_Q<3>(2)
================================================================================
JobId fe0cba80120c Thu Jan 27 06:06:48 2022
--------------------------------------------------------------------------------
Cycle 0
--------------------------------------------------------------------------------
Set up system...
Triangulation: 28672 cells
DoFHandler: 232609 DoFs
Solve using Newton's method...
DEAL:GMRES::Check 0 25.2492
DEAL:GMRES::Starting value 25.2492
DEAL:GMRES::Check 1 1.02964
DEAL:GMRES::Check 2 0.464888
DEAL:GMRES::Check 3 0.263190
DEAL:GMRES::Check 4 0.0703056
DEAL:GMRES::Check 4 0.0703056
DEAL:GMRES::Check 5 0.0580013
DEAL:GMRES::Check 6 0.0363219
DEAL:GMRES::Check 7 0.0330487
DEAL:GMRES::Check 8 0.0222448
DEAL:GMRES::Check 8 0.0222448
DEAL:GMRES::Check 9 0.0200449
DEAL:GMRES::Check 10 0.0173175
DEAL:GMRES::Check 11 0.0106881
DEAL:GMRES::Check 12 0.00388264
DEAL:GMRES::Check 12 0.00388264
DEAL:GMRES::Check 13 0.00140426
DEAL:GMRES::Check 14 0.00116016
DEAL:GMRES::Check 15 0.000288553
DEAL:GMRES::Check 16 0.000219209
DEAL:GMRES::Check 16 0.000219209
DEAL:GMRES::Check 17 0.000169839
DEAL:GMRES::Check 18 6.81943e-05
DEAL:GMRES::Check 19 6.79830e-05
DEAL:GMRES::Check 20 1.91041e-05
DEAL:GMRES::Check 20 1.91041e-05
DEAL:GMRES::Check 21 1.69033e-05
DEAL:GMRES::Check 22 6.16672e-06
DEAL:GMRES::Check 23 5.80567e-06
DEAL:GMRES::Check 24 5.00262e-06
DEAL:GMRES::Check 24 5.00262e-06
DEAL:GMRES::Check 25 4.73400e-06
DEAL:GMRES::Check 26 4.73083e-06
DEAL:GMRES::Check 27 1.15946e-06
DEAL:GMRES::Check 28 5.78294e-07
DEAL:GMRES::Check 28 5.78294e-07
DEAL:GMRES::Check 29 4.76543e-07
DEAL:GMRES::Check 30 2.31646e-07
DEAL:GMRES::Check 31 2.30765e-07
DEAL:GMRES::Check 32 1.46000e-07
DEAL:GMRES::Check 32 1.46000e-07
DEAL:GMRES::Check 33 1.40243e-07
DEAL:GMRES::Check 34 7.90235e-08
DEAL:GMRES::Check 35 4.37569e-08
DEAL:GMRES::Check 36 2.08236e-08
DEAL:GMRES::Check 36 2.08236e-08
DEAL:GMRES::Check 37 1.83929e-08
DEAL:GMRES::Check 38 1.03665e-08
DEAL:GMRES::Check 39 6.46205e-09
DEAL:GMRES::Check 40 3.44757e-09
DEAL:GMRES::Check 40 3.44757e-09
DEAL:GMRES::Check 41 1.71066e-09
DEAL:GMRES::Check 42 1.71042e-09
DEAL:GMRES::Check 43 9.55959e-10
DEAL:GMRES::Check 44 3.76580e-10
DEAL:GMRES::Check 44 3.76580e-10
DEAL:GMRES::Check 45 2.27942e-10
DEAL:GMRES::Check 46 1.22379e-10
DEAL:GMRES::Check 47 1.20322e-10
DEAL:GMRES::Check 48 4.47285e-11
DEAL:GMRES::Check 48 4.47284e-11
DEAL:GMRES::Check 49 4.15619e-11
DEAL:GMRES::Check 50 1.68414e-11
DEAL:GMRES::Check 51 7.70878e-12
DEAL:GMRES::Check 52 3.08794e-12
DEAL:GMRES::Check 52 3.08799e-12
DEAL:GMRES::Check 53 2.42904e-12
DEAL:GMRES::Check 54 1.62947e-12
DEAL:GMRES::Check 55 6.95513e-13
DEAL:GMRES::Convergence step 55 value 6.95513e-13
Nstep 1, errf = 4.85787e-05, errx = 60.2426, it = 55
DEAL:GMRES::Check 0 0.200010
DEAL:GMRES::Starting value 0.200010
DEAL:GMRES::Check 1 0.00788639
DEAL:GMRES::Check 2 0.00506430
DEAL:GMRES::Check 3 0.00246018
DEAL:GMRES::Check 4 0.000551852
DEAL:GMRES::Check 4 0.000551852
DEAL:GMRES::Check 5 0.000451143
DEAL:GMRES::Check 6 0.000358859
DEAL:GMRES::Check 7 0.000353628
DEAL:GMRES::Check 8 0.000259076
DEAL:GMRES::Check 8 0.000259076
DEAL:GMRES::Check 9 0.000234279
DEAL:GMRES::Check 10 0.000141028
DEAL:GMRES::Check 11 0.000112595
DEAL:GMRES::Check 12 5.66726e-05
DEAL:GMRES::Check 12 5.66726e-05
DEAL:GMRES::Check 13 4.12704e-05
DEAL:GMRES::Check 14 2.34796e-05
DEAL:GMRES::Check 15 2.20837e-05
DEAL:GMRES::Check 16 1.96805e-05
DEAL:GMRES::Check 16 1.96805e-05
DEAL:GMRES::Check 17 1.94617e-05
DEAL:GMRES::Check 18 7.22846e-06
DEAL:GMRES::Check 19 4.62679e-06
DEAL:GMRES::Check 20 2.94292e-06
DEAL:GMRES::Check 20 2.94292e-06
DEAL:GMRES::Check 21 1.97684e-06
DEAL:GMRES::Check 22 1.97221e-06
DEAL:GMRES::Check 23 1.73070e-06
DEAL:GMRES::Check 24 8.22979e-08
DEAL:GMRES::Check 24 8.22979e-08
DEAL:GMRES::Check 25 3.26905e-08
DEAL:GMRES::Check 26 2.83243e-08
DEAL:GMRES::Check 27 1.14731e-08
DEAL:GMRES::Check 28 8.48691e-09
DEAL:GMRES::Check 28 8.48691e-09
DEAL:GMRES::Check 29 6.66122e-09
DEAL:GMRES::Check 30 6.21653e-09
DEAL:GMRES::Check 31 1.50705e-09
DEAL:GMRES::Check 32 6.14413e-10
DEAL:GMRES::Check 32 6.14413e-10
DEAL:GMRES::Check 33 5.63526e-10
DEAL:GMRES::Check 34 2.07749e-10
DEAL:GMRES::Check 35 1.78125e-10
DEAL:GMRES::Check 36 1.26029e-10
DEAL:GMRES::Check 36 1.26029e-10
DEAL:GMRES::Check 37 1.23350e-10
DEAL:GMRES::Check 38 4.04264e-11
DEAL:GMRES::Check 39 1.56014e-11
DEAL:GMRES::Check 40 7.15859e-12
DEAL:GMRES::Check 40 7.15859e-12
DEAL:GMRES::Check 41 6.03967e-12
DEAL:GMRES::Check 42 5.47824e-12
DEAL:GMRES::Check 43 3.81547e-12
DEAL:GMRES::Check 44 1.86419e-12
DEAL:GMRES::Check 44 1.86419e-12
DEAL:GMRES::Check 45 1.64310e-12
DEAL:GMRES::Check 46 4.26039e-13
DEAL:GMRES::Convergence step 46 value 4.26039e-13
Nstep 2, errf = 3.28133e-09, errx = 0.518129, it = 46
DEAL:GMRES::Check 0 1.51009e-05
DEAL:GMRES::Starting value 1.51009e-05
DEAL:GMRES::Check 1 6.06440e-07
DEAL:GMRES::Check 2 4.32493e-07
DEAL:GMRES::Check 3 1.90244e-07
DEAL:GMRES::Check 4 4.24134e-08
DEAL:GMRES::Check 4 4.24134e-08
DEAL:GMRES::Check 5 3.35372e-08
DEAL:GMRES::Check 6 2.95403e-08
DEAL:GMRES::Check 7 2.92521e-08
DEAL:GMRES::Check 8 2.18082e-08
DEAL:GMRES::Check 8 2.18082e-08
DEAL:GMRES::Check 9 1.95304e-08
DEAL:GMRES::Check 10 1.02972e-08
DEAL:GMRES::Check 11 8.60819e-09
DEAL:GMRES::Check 12 4.60425e-09
DEAL:GMRES::Check 12 4.60425e-09
DEAL:GMRES::Check 13 4.08036e-09
DEAL:GMRES::Check 14 1.93057e-09
DEAL:GMRES::Check 15 1.65261e-09
DEAL:GMRES::Check 16 1.64997e-09
DEAL:GMRES::Check 16 1.64997e-09
DEAL:GMRES::Check 17 1.64919e-09
DEAL:GMRES::Check 18 5.74881e-10
DEAL:GMRES::Check 19 3.84769e-10
DEAL:GMRES::Check 20 1.80949e-10
DEAL:GMRES::Check 20 1.80949e-10
DEAL:GMRES::Check 21 1.44165e-10
DEAL:GMRES::Check 22 8.90440e-11
DEAL:GMRES::Check 23 2.89740e-11
DEAL:GMRES::Check 24 8.66621e-12
DEAL:GMRES::Check 24 8.66621e-12
DEAL:GMRES::Check 25 6.22335e-12
DEAL:GMRES::Check 26 6.11289e-12
DEAL:GMRES::Check 27 2.01745e-12
DEAL:GMRES::Check 28 1.62864e-12
DEAL:GMRES::Check 28 1.62864e-12
DEAL:GMRES::Check 29 1.61566e-12
DEAL:GMRES::Check 30 6.46141e-13
DEAL:GMRES::Convergence step 30 value 6.46141e-13
Nstep 3, errf = 7.95997e-14, errx = 4.02406e-05, it = 30
Convergence step 3 value 7.95997e-14 (used wall time: 3.29803 s)
Time for setup+solve (CPU/Wall) 3.52526/3.56872 s
Output results...
H1 seminorm: 0.588667
+---------------------------------------------+------------+------------+
| Total wallclock time elapsed since start | 3.92s | |
| | | |
| Section | no. calls | wall time | % of total |
+---------------------------------+-----------+------------+------------+
| assemble right hand side | 3 | 0.00178s | 0.00% |
| compute residual | 3 | 0.00199s | 0.00% |
| compute update | 3 | 3.29s | 84.% |
| make grid | 1 | 0.219s | 5.6% |
| setup system | 1 | 0.270s | 6.9% |
| solve | 1 | 3.30s | 84.% |
+---------------------------------+-----------+------------+------------+
#include <deal.II/base/function.h>
#include <deal.II/base/quadrature_lib.h>
#include <deal.II/base/timer.h>
#include <deal.II/base/vectorization.h>
#include <deal.II/dofs/dof_accessor.h>
#include <deal.II/dofs/dof_handler.h>
#include <deal.II/dofs/dof_tools.h>
#include <deal.II/fe/fe_q.h>
#include <deal.II/fe/mapping_q_generic.h>
#include <deal.II/grid/grid_generator.h>
#include <deal.II/grid/grid_out.h>
#include <deal.II/grid/manifold_lib.h>
#include <deal.II/grid/tria.h>
#include <deal.II/grid/tria_accessor.h>
#include <deal.II/grid/tria_iterator.h>
#include <deal.II/lac/affine_constraints.h>
#include <deal.II/lac/precondition.h>
// MODIFIED: use gmres rather than cg
#include <deal.II/lac/solver_gmres.h>
#include <deal.II/lac/trilinos_precondition.h>
#include <deal.II/lac/vector.h>
#include <deal.II/matrix_free/fe_evaluation.h>
#include <deal.II/matrix_free/matrix_free.h>
#include <deal.II/matrix_free/operators.h>
#include <deal.II/matrix_free/tools.h>
#include <deal.II/multigrid/mg_coarse.h>
#include <deal.II/multigrid/mg_constrained_dofs.h>
#include <deal.II/multigrid/mg_matrix.h>
#include <deal.II/multigrid/mg_smoother.h>
#include <deal.II/multigrid/mg_tools.h>
#include <deal.II/multigrid/mg_transfer_matrix_free.h>
#include <deal.II/multigrid/multigrid.h>
#include <deal.II/numerics/data_out.h>
#include <deal.II/numerics/vector_tools.h>
// MODIFIED: add boost null sink
#include <boost/iostreams/device/null.hpp>
#include <fstream>
#include <iostream>
namespace Step66 {
using namespace dealii;
// MODIFIED: add null sink object
boost::iostreams::stream<boost::iostreams::null_sink> null_log{
boost::iostreams::null_sink{}};
template <int dim, int fe_degree, typename number>
class JacobianOperator : public MatrixFreeOperators::Base<
dim, LinearAlgebra::distributed::Vector<number>> {
public:
using value_type = number;
using FECellIntegrator =
FEEvaluation<dim, fe_degree, fe_degree + 1, 1, number>;
JacobianOperator();
virtual void clear() override;
void evaluate_newton_step(
const LinearAlgebra::distributed::Vector<number> &newton_step);
virtual void compute_diagonal() override;
// MODIFIED: compute real matrix within given matrix and constraint
template <typename MatrixType>
void compute_matrix(MatrixType &matrix,
const AffineConstraints<number> &constraint =
AffineConstraints<number>()) const {
matrix = 0.;
MatrixFreeTools::compute_matrix(
*this->data, constraint, matrix,
&JacobianOperator<dim, fe_degree, number>::local_compute_diagonal,
this);
}
private:
virtual void apply_add(
LinearAlgebra::distributed::Vector<number> &dst,
const LinearAlgebra::distributed::Vector<number> &src) const override;
void local_apply(
const MatrixFree<dim, number> &data,
LinearAlgebra::distributed::Vector<number> &dst,
const LinearAlgebra::distributed::Vector<number> &src,
const std::pair<unsigned int, unsigned int> &cell_range) const;
void local_compute_diagonal(FECellIntegrator &integrator) const;
Table<2, VectorizedArray<number>> nonlinear_values;
};
template <int dim, int fe_degree, typename number>
JacobianOperator<dim, fe_degree, number>::JacobianOperator()
: MatrixFreeOperators::Base<dim,
LinearAlgebra::distributed::Vector<number>>() {}
template <int dim, int fe_degree, typename number>
void JacobianOperator<dim, fe_degree, number>::clear() {
nonlinear_values.reinit(0, 0);
MatrixFreeOperators::Base<
dim, LinearAlgebra::distributed::Vector<number>>::clear();
}
template <int dim, int fe_degree, typename number>
void JacobianOperator<dim, fe_degree, number>::evaluate_newton_step(
const LinearAlgebra::distributed::Vector<number> &newton_step) {
const unsigned int n_cells = this->data->n_cell_batches();
FECellIntegrator phi(*this->data);
nonlinear_values.reinit(n_cells, phi.n_q_points);
for (unsigned int cell = 0; cell < n_cells; ++cell) {
phi.reinit(cell);
phi.read_dof_values_plain(newton_step);
phi.evaluate(EvaluationFlags::values);
for (unsigned int q = 0; q < phi.n_q_points; ++q) {
nonlinear_values(cell, q) = std::exp(phi.get_value(q));
}
}
}
template <int dim, int fe_degree, typename number>
void JacobianOperator<dim, fe_degree, number>::local_apply(
const MatrixFree<dim, number> &data,
LinearAlgebra::distributed::Vector<number> &dst,
const LinearAlgebra::distributed::Vector<number> &src,
const std::pair<unsigned int, unsigned int> &cell_range) const {
FECellIntegrator phi(data);
for (unsigned int cell = cell_range.first; cell < cell_range.second; ++cell) {
AssertDimension(nonlinear_values.size(0),
phi.get_matrix_free().n_cell_batches());
AssertDimension(nonlinear_values.size(1), phi.n_q_points);
phi.reinit(cell);
phi.gather_evaluate(src,
EvaluationFlags::values | EvaluationFlags::gradients);
for (unsigned int q = 0; q < phi.n_q_points; ++q) {
phi.submit_value(-nonlinear_values(cell, q) * phi.get_value(q), q);
phi.submit_gradient(phi.get_gradient(q), q);
}
phi.integrate_scatter(EvaluationFlags::values | EvaluationFlags::gradients,
dst);
}
}
template <int dim, int fe_degree, typename number>
void JacobianOperator<dim, fe_degree, number>::apply_add(
LinearAlgebra::distributed::Vector<number> &dst,
const LinearAlgebra::distributed::Vector<number> &src) const {
this->data->cell_loop(&JacobianOperator::local_apply, this, dst, src);
}
template <int dim, int fe_degree, typename number>
void JacobianOperator<dim, fe_degree, number>::local_compute_diagonal(
FECellIntegrator &phi) const {
AssertDimension(nonlinear_values.size(0),
phi.get_matrix_free().n_cell_batches());
AssertDimension(nonlinear_values.size(1), phi.n_q_points);
const unsigned int cell = phi.get_current_cell_index();
phi.evaluate(EvaluationFlags::values | EvaluationFlags::gradients);
for (unsigned int q = 0; q < phi.n_q_points; ++q) {
phi.submit_value(-nonlinear_values(cell, q) * phi.get_value(q), q);
phi.submit_gradient(phi.get_gradient(q), q);
}
phi.integrate(EvaluationFlags::values | EvaluationFlags::gradients);
}
template <int dim, int fe_degree, typename number>
void JacobianOperator<dim, fe_degree, number>::compute_diagonal() {
this->inverse_diagonal_entries.reset(
new DiagonalMatrix<LinearAlgebra::distributed::Vector<number>>());
LinearAlgebra::distributed::Vector<number> &inverse_diagonal =
this->inverse_diagonal_entries->get_vector();
this->data->initialize_dof_vector(inverse_diagonal);
MatrixFreeTools::compute_diagonal(*this->data, inverse_diagonal,
&JacobianOperator::local_compute_diagonal,
this);
for (auto &diagonal_element : inverse_diagonal) {
diagonal_element =
(std::abs(diagonal_element) > 1.0e-10) ? (1.0 / diagonal_element) : 1.0;
}
}
// MODIFIED: add solver method to choose to use gmg or amg
enum class SolverMethod { use_gmg, use_amg };
template <int dim, int fe_degree>
class GelfandProblem {
public:
GelfandProblem();
void run();
SolverMethod solver_method;
private:
void make_grid();
void setup_system();
void evaluate_residual(
LinearAlgebra::distributed::Vector<double> &dst,
const LinearAlgebra::distributed::Vector<double> &src) const;
void local_evaluate_residual(
const MatrixFree<dim, double> &data,
LinearAlgebra::distributed::Vector<double> &dst,
const LinearAlgebra::distributed::Vector<double> &src,
const std::pair<unsigned int, unsigned int> &cell_range) const;
void assemble_rhs();
double compute_residual(const double alpha);
void compute_update();
void solve();
double compute_solution_norm() const;
void output_results(const unsigned int cycle) const;
// MODIFIED: add function to assemble real matrix by using compute_matrix
void assemble_matrix_real() {
if (matrix_real.m() != dof_handler.n_locally_owned_dofs() ||
matrix_real.n() != dof_handler.n_locally_owned_dofs()) {
TrilinosWrappers::SparsityPattern sp(dof_handler.locally_owned_dofs(),
MPI_COMM_WORLD);
DoFTools::make_sparsity_pattern(dof_handler, sp, constraints);
sp.compress();
matrix_real.reinit(sp);
}
system_matrix.compute_matrix(matrix_real, constraints);
}
TrilinosWrappers::SparseMatrix matrix_real;
parallel::distributed::Triangulation<dim> triangulation;
const MappingQGeneric<dim> mapping;
FE_Q<dim> fe;
DoFHandler<dim> dof_handler;
AffineConstraints<double> constraints;
using SystemMatrixType = JacobianOperator<dim, fe_degree, double>;
SystemMatrixType system_matrix;
MGConstrainedDoFs mg_constrained_dofs;
using LevelMatrixType = JacobianOperator<dim, fe_degree, float>;
MGLevelObject<LevelMatrixType> mg_matrices;
MGLevelObject<LinearAlgebra::distributed::Vector<float>> mg_solution;
MGTransferMatrixFree<dim, float> mg_transfer;
LinearAlgebra::distributed::Vector<double> solution;
LinearAlgebra::distributed::Vector<double> newton_update;
LinearAlgebra::distributed::Vector<double> system_rhs;
unsigned int linear_iterations;
ConditionalOStream pcout;
TimerOutput computing_timer;
};
template <int dim, int fe_degree>
GelfandProblem<dim, fe_degree>::GelfandProblem()
: triangulation(MPI_COMM_WORLD,
Triangulation<dim>::limit_level_difference_at_vertices,
parallel::distributed::Triangulation<
dim>::construct_multigrid_hierarchy),
mapping(fe_degree),
fe(fe_degree),
dof_handler(triangulation),
pcout(std::cout, Utilities::MPI::this_mpi_process(MPI_COMM_WORLD) == 0),
computing_timer(MPI_COMM_WORLD, pcout, TimerOutput::never,
TimerOutput::wall_times) {}
template <int dim, int fe_degree>
void GelfandProblem<dim, fe_degree>::make_grid() {
TimerOutput::Scope t(computing_timer, "make grid");
SphericalManifold<dim> boundary_manifold;
TransfiniteInterpolationManifold<dim> inner_manifold;
GridGenerator::hyper_ball(triangulation);
triangulation.set_all_manifold_ids(1);
triangulation.set_all_manifold_ids_on_boundary(0);
triangulation.set_manifold(0, boundary_manifold);
inner_manifold.initialize(triangulation);
triangulation.set_manifold(1, inner_manifold);
// MODIFIED: fixed refine level
triangulation.refine_global(4);
}
template <int dim, int fe_degree>
void GelfandProblem<dim, fe_degree>::setup_system() {
TimerOutput::Scope t(computing_timer, "setup system");
system_matrix.clear();
mg_matrices.clear_elements();
dof_handler.distribute_dofs(fe);
dof_handler.distribute_mg_dofs();
IndexSet locally_relevant_dofs;
DoFTools::extract_locally_relevant_dofs(dof_handler, locally_relevant_dofs);
constraints.clear();
constraints.reinit(locally_relevant_dofs);
DoFTools::make_hanging_node_constraints(dof_handler, constraints);
VectorTools::interpolate_boundary_values(
dof_handler, 0, Functions::ZeroFunction<dim>(), constraints);
constraints.close();
{
typename MatrixFree<dim, double>::AdditionalData additional_data;
additional_data.tasks_parallel_scheme =
MatrixFree<dim, double>::AdditionalData::partition_color;
additional_data.mapping_update_flags =
(update_values | update_gradients | update_JxW_values |
update_quadrature_points);
auto system_mf_storage = std::make_shared<MatrixFree<dim, double>>();
system_mf_storage->reinit(mapping, dof_handler, constraints,
QGauss<1>(fe.degree + 1), additional_data);
system_matrix.initialize(system_mf_storage);
}
system_matrix.initialize_dof_vector(solution);
system_matrix.initialize_dof_vector(newton_update);
system_matrix.initialize_dof_vector(system_rhs);
const unsigned int nlevels = triangulation.n_global_levels();
mg_matrices.resize(0, nlevels - 1);
mg_solution.resize(0, nlevels - 1);
std::set<types::boundary_id> dirichlet_boundary;
dirichlet_boundary.insert(0);
mg_constrained_dofs.initialize(dof_handler);
mg_constrained_dofs.make_zero_boundary_constraints(dof_handler,
dirichlet_boundary);
mg_transfer.initialize_constraints(mg_constrained_dofs);
mg_transfer.build(dof_handler);
for (unsigned int level = 0; level < nlevels; ++level) {
IndexSet relevant_dofs;
DoFTools::extract_locally_relevant_level_dofs(dof_handler, level,
relevant_dofs);
AffineConstraints<double> level_constraints;
level_constraints.reinit(relevant_dofs);
level_constraints.add_lines(
mg_constrained_dofs.get_boundary_indices(level));
level_constraints.close();
typename MatrixFree<dim, float>::AdditionalData additional_data;
additional_data.tasks_parallel_scheme =
MatrixFree<dim, float>::AdditionalData::partition_color;
additional_data.mapping_update_flags =
(update_values | update_gradients | update_JxW_values |
update_quadrature_points);
additional_data.mg_level = level;
auto mg_mf_storage_level = std::make_shared<MatrixFree<dim, float>>();
mg_mf_storage_level->reinit(mapping, dof_handler, level_constraints,
QGauss<1>(fe.degree + 1), additional_data);
mg_matrices[level].initialize(mg_mf_storage_level, mg_constrained_dofs,
level);
mg_matrices[level].initialize_dof_vector(mg_solution[level]);
}
}
template <int dim, int fe_degree>
void GelfandProblem<dim, fe_degree>::evaluate_residual(
LinearAlgebra::distributed::Vector<double> &dst,
const LinearAlgebra::distributed::Vector<double> &src) const {
auto matrix_free = system_matrix.get_matrix_free();
matrix_free->cell_loop(&GelfandProblem::local_evaluate_residual, this, dst,
src, true);
}
template <int dim, int fe_degree>
void GelfandProblem<dim, fe_degree>::local_evaluate_residual(
const MatrixFree<dim, double> &data,
LinearAlgebra::distributed::Vector<double> &dst,
const LinearAlgebra::distributed::Vector<double> &src,
const std::pair<unsigned int, unsigned int> &cell_range) const {
FEEvaluation<dim, fe_degree, fe_degree + 1, 1, double> phi(data);
for (unsigned int cell = cell_range.first; cell < cell_range.second; ++cell) {
phi.reinit(cell);
phi.read_dof_values_plain(src);
phi.evaluate(EvaluationFlags::values | EvaluationFlags::gradients);
for (unsigned int q = 0; q < phi.n_q_points; ++q) {
phi.submit_value(-std::exp(phi.get_value(q)), q);
phi.submit_gradient(phi.get_gradient(q), q);
}
phi.integrate_scatter(EvaluationFlags::values | EvaluationFlags::gradients,
dst);
}
}
template <int dim, int fe_degree>
void GelfandProblem<dim, fe_degree>::assemble_rhs() {
TimerOutput::Scope t(computing_timer, "assemble right hand side");
evaluate_residual(system_rhs, solution);
system_rhs *= -1.0;
}
template <int dim, int fe_degree>
double GelfandProblem<dim, fe_degree>::compute_residual(const double alpha) {
TimerOutput::Scope t(computing_timer, "compute residual");
LinearAlgebra::distributed::Vector<double> residual;
LinearAlgebra::distributed::Vector<double> evaluation_point;
system_matrix.initialize_dof_vector(residual);
system_matrix.initialize_dof_vector(evaluation_point);
evaluation_point = solution;
if (alpha > 1e-12) {
evaluation_point.add(alpha, newton_update);
}
evaluate_residual(residual, evaluation_point);
return residual.l2_norm();
}
template <int dim, int fe_degree>
void GelfandProblem<dim, fe_degree>::compute_update() {
TimerOutput::Scope t(computing_timer, "compute update");
solution.update_ghost_values();
system_matrix.evaluate_newton_step(solution);
newton_update = 0.0;
// MODIFIED: enable solver check log & set solver to GMRES(4)
SolverControl solver_control(100, 1.e-12, true);
SolverGMRES<LinearAlgebra::distributed::Vector<double>>::AdditionalData
solver_data;
solver_data.max_n_tmp_vectors = 6;
SolverGMRES<LinearAlgebra::distributed::Vector<double>> gmres(solver_control,
solver_data);
// MODIFIED: add switch for solver method
switch (solver_method) {
case SolverMethod::use_gmg: {
mg_transfer.interpolate_to_mg(dof_handler, mg_solution, solution);
using SmootherType =
PreconditionChebyshev<LevelMatrixType,
LinearAlgebra::distributed::Vector<float>>;
mg::SmootherRelaxation<SmootherType,
LinearAlgebra::distributed::Vector<float>>
mg_smoother;
MGLevelObject<typename SmootherType::AdditionalData> smoother_data;
smoother_data.resize(0, triangulation.n_global_levels() - 1);
for (unsigned int level = 0; level < triangulation.n_global_levels();
++level) {
if (level > 0) {
smoother_data[level].smoothing_range = 15.;
smoother_data[level].degree = 4;
smoother_data[level].eig_cg_n_iterations = 10;
} else {
smoother_data[0].smoothing_range = 1e-3;
smoother_data[0].degree = numbers::invalid_unsigned_int;
smoother_data[0].eig_cg_n_iterations = mg_matrices[0].m();
}
mg_matrices[level].evaluate_newton_step(mg_solution[level]);
mg_matrices[level].compute_diagonal();
smoother_data[level].preconditioner =
mg_matrices[level].get_matrix_diagonal_inverse();
}
mg_smoother.initialize(mg_matrices, smoother_data);
MGCoarseGridApplySmoother<LinearAlgebra::distributed::Vector<float>>
mg_coarse;
mg_coarse.initialize(mg_smoother);
mg::Matrix<LinearAlgebra::distributed::Vector<float>> mg_matrix(
mg_matrices);
MGLevelObject<MatrixFreeOperators::MGInterfaceOperator<LevelMatrixType>>
mg_interface_matrices;
mg_interface_matrices.resize(0, triangulation.n_global_levels() - 1);
for (unsigned int level = 0; level < triangulation.n_global_levels();
++level) {
mg_interface_matrices[level].initialize(mg_matrices[level]);
}
mg::Matrix<LinearAlgebra::distributed::Vector<float>> mg_interface(
mg_interface_matrices);
Multigrid<LinearAlgebra::distributed::Vector<float>> mg(
mg_matrix, mg_coarse, mg_transfer, mg_smoother, mg_smoother);
mg.set_edge_matrices(mg_interface, mg_interface);
PreconditionMG<dim, LinearAlgebra::distributed::Vector<float>,
MGTransferMatrixFree<dim, float>>
preconditioner(dof_handler, mg, mg_transfer);
gmres.solve(system_matrix, newton_update, system_rhs, preconditioner);
break;
}
case SolverMethod::use_amg: {
assemble_matrix_real();
TrilinosWrappers::PreconditionAMG preconditioner;
preconditioner.initialize(matrix_real);
gmres.solve(system_matrix, newton_update, system_rhs, preconditioner);
break;
}
}
constraints.distribute(newton_update);
linear_iterations = solver_control.last_step();
solution.zero_out_ghost_values();
}
template <int dim, int fe_degree>
void GelfandProblem<dim, fe_degree>::solve() {
TimerOutput::Scope t(computing_timer, "solve");
const unsigned int itmax = 10;
const double TOLf = 1e-12;
const double TOLx = 1e-10;
Timer solver_timer;
solver_timer.start();
for (unsigned int newton_step = 1; newton_step <= itmax; ++newton_step) {
assemble_rhs();
compute_update();
const double ERRx = newton_update.l2_norm();
const double ERRf = compute_residual(1.0);
solution.add(1.0, newton_update);
pcout << " Nstep " << newton_step << ", errf = " << ERRf
<< ", errx = " << ERRx << ", it = " << linear_iterations << std::endl;
if (ERRf < TOLf || ERRx < TOLx) {
solver_timer.stop();
pcout << "Convergence step " << newton_step << " value " << ERRf
<< " (used wall time: " << solver_timer.wall_time() << " s)"
<< std::endl;
break;
} else if (newton_step == itmax) {
solver_timer.stop();
pcout << "WARNING: No convergence of Newton's method after "
<< newton_step << " steps." << std::endl;
break;
}
}
}
template <int dim, int fe_degree>
double GelfandProblem<dim, fe_degree>::compute_solution_norm() const {
solution.update_ghost_values();
Vector<float> norm_per_cell(triangulation.n_active_cells());
VectorTools::integrate_difference(
mapping, dof_handler, solution, Functions::ZeroFunction<dim>(),
norm_per_cell, QGauss<dim>(fe.degree + 2), VectorTools::H1_seminorm);
solution.zero_out_ghost_values();
return VectorTools::compute_global_error(triangulation, norm_per_cell,
VectorTools::H1_seminorm);
}
template <int dim, int fe_degree>
void GelfandProblem<dim, fe_degree>::output_results(
const unsigned int cycle) const {
if (triangulation.n_global_active_cells() > 1e6) return;
solution.update_ghost_values();
DataOut<dim> data_out;
data_out.attach_dof_handler(dof_handler);
data_out.add_data_vector(solution, "solution");
Vector<float> subdomain(triangulation.n_active_cells());
for (unsigned int i = 0; i < subdomain.size(); ++i) {
subdomain(i) = triangulation.locally_owned_subdomain();
}
data_out.add_data_vector(subdomain, "subdomain");
data_out.build_patches(mapping, fe.degree, DataOut<dim>::curved_inner_cells);
DataOutBase::VtkFlags flags;
flags.compression_level = DataOutBase::VtkFlags::best_speed;
data_out.set_flags(flags);
data_out.write_vtu_with_pvtu_record(
"./", "solution_" + std::to_string(dim) + "d", cycle, MPI_COMM_WORLD, 3);
solution.zero_out_ghost_values();
}
template <int dim, int fe_degree>
void GelfandProblem<dim, fe_degree>::run() {
{
const unsigned int n_ranks =
Utilities::MPI::n_mpi_processes(MPI_COMM_WORLD);
const unsigned int n_vect_doubles = VectorizedArray<double>::size();
const unsigned int n_vect_bits = 8 * sizeof(double) * n_vect_doubles;
std::string DAT_header = "START DATE: " + Utilities::System::get_date() +
", TIME: " + Utilities::System::get_time();
std::string MPI_header = "Running with " + std::to_string(n_ranks) +
" MPI process" + (n_ranks > 1 ? "es" : "");
std::string VEC_header =
"Vectorization over " + std::to_string(n_vect_doubles) +
" doubles = " + std::to_string(n_vect_bits) + " bits (" +
Utilities::System::get_current_vectorization_level() +
"), VECTORIZATION_LEVEL=" +
std::to_string(DEAL_II_COMPILER_VECTORIZATION_LEVEL);
std::string SOL_header = "Finite element space: " + fe.get_name();
pcout << std::string(80, '=') << std::endl;
pcout << DAT_header << std::endl;
pcout << std::string(80, '-') << std::endl;
pcout << MPI_header << std::endl;
pcout << VEC_header << std::endl;
pcout << SOL_header << std::endl;
pcout << std::string(80, '=') << std::endl;
}
// MODIFIED: attach deallog output
{
if (Utilities::MPI::this_mpi_process(MPI_COMM_WORLD) == 0) {
deallog.attach(std::cout);
deallog.depth_file(3);
} else {
deallog.attach(null_log);
}
}
for (unsigned int cycle = 0; cycle < 1; ++cycle) {
pcout << std::string(80, '-') << std::endl;
pcout << "Cycle " << cycle << std::endl;
pcout << std::string(80, '-') << std::endl;
if (cycle == 0) {
make_grid();
} else {
triangulation.refine_global(1);
}
Timer timer;
pcout << "Set up system..." << std::endl;
setup_system();
pcout << " Triangulation: " << triangulation.n_global_active_cells()
<< " cells" << std::endl;
pcout << " DoFHandler: " << dof_handler.n_dofs() << " DoFs"
<< std::endl;
pcout << std::endl;
pcout << "Solve using Newton's method..." << std::endl;
solve();
pcout << std::endl;
timer.stop();
pcout << "Time for setup+solve (CPU/Wall) " << timer.cpu_time() << "/"
<< timer.wall_time() << " s" << std::endl;
pcout << std::endl;
pcout << "Output results..." << std::endl;
const double norm = compute_solution_norm();
// MODIFIED: disable output
// output_results(cycle);
pcout << " H1 seminorm: " << norm << std::endl;
pcout << std::endl;
computing_timer.print_summary();
computing_timer.reset();
}
}
} // namespace Step66
int main(int argc, char *argv[]) {
try {
using namespace Step66;
Utilities::MPI::MPI_InitFinalize mpi_init(argc, argv, 1);
// MODIFIED: solve with gmg & amg respectively
{
GelfandProblem<3, 2> gelfand_problem;
gelfand_problem.solver_method = SolverMethod::use_gmg;
gelfand_problem.run();
}
{
GelfandProblem<3, 2> gelfand_problem;
gelfand_problem.solver_method = SolverMethod::use_amg;
gelfand_problem.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;
}
