I tried to use the AD framework in my own codes (mainly for solving
nonlinear equations), and therefore wrote a small test program for solving
the linear diffusion equation (after I can calculate the jacobian by hand
for that equation). Thus, my residual for the AD-jacobian is -0.5 * (\nabla
u_1 + \nabla u_0)\cdot\nabla\varphi - \frac{1}{dt}(u_1 - u_0)\cdot\varphi,
with u_0 my current solution, and u_1 my solution for the next step. When
creating the jacobian by hand, I use the same approach as shown in example
33, and assume that u_1 = u_0 in my first step, resulting in
-0.5\cdot(\nabla\varphi_i\cdot\nabla\varphi_j) -
\frac{1}{dt}\cdot\varphi_i\cdot\varphi_j, and my right hand side is \nabla
u_0\cdot\nabla\varphi.
Solving both systems should give me a value for \delta u, i.e. the update
for my current solution, such that u_1 = \delta u + u_0 (after the system
is linear).
Though, when running the code in the attachment, I get the correct result
(i.e. \delta u) for my hand-assembled system, but the solution I get for
the AD-generated jacobian is u_1, not \delta u. Did I misinterpret
something, or is there another bug in the code?
Thanks!
--
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/e556dfdd-d5e8-453c-896a-0fb2f3cc73b2%40googlegroups.com.
/* ---------------------------------------------------------------------
*
* Copyright (C) 2009 - 2019 by the deal.II authors
*
* This file is part of the deal.II library.
*
* The deal.II library is free software; you can use it, redistribute
* it, and/or modify it under the terms of the GNU Lesser General
* Public License as published by the Free Software Foundation; either
* version 2.1 of the License, or (at your option) any later version.
* The full text of the license can be found in the file LICENSE.md at
* the top level directory of deal.II.
*
* ---------------------------------------------------------------------
*
* Author: Wolfgang Bangerth, Texas A&M University, 2009, 2010
* Timo Heister, University of Goettingen, 2009, 2010
*/
#include <deal.II/base/quadrature_lib.h>
#include <deal.II/base/function.h>
#include <deal.II/base/timer.h>
#include <deal.II/lac/generic_linear_algebra.h>
namespace LA
{
#undef DEAL_II_WITH_PETSC
#if defined(DEAL_II_WITH_PETSC) && !defined(DEAL_II_PETSC_WITH_COMPLEX) && \
!(defined(DEAL_II_WITH_TRILINOS) && defined(FORCE_USE_OF_TRILINOS))
using namespace dealii::LinearAlgebraPETSc;
# define USE_PETSC_LA
#elif defined(DEAL_II_WITH_TRILINOS)
using namespace dealii::LinearAlgebraTrilinos;
#else
# error DEAL_II_WITH_PETSC or DEAL_II_WITH_TRILINOS required
#endif
} // namespace LA
#include <deal.II/lac/vector.h>
#include <deal.II/lac/full_matrix.h>
#include <deal.II/lac/solver_cg.h>
#include <deal.II/lac/affine_constraints.h>
#include <deal.II/lac/dynamic_sparsity_pattern.h>
#include <deal.II/lac/sparsity_tools.h>
#include <deal.II/lac/linear_operator.h>
#include <deal.II/lac/parpack_solver.h>
#include <deal.II/lac/solver_gmres.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/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/data_out.h>
#include <deal.II/numerics/error_estimator.h>
#include <deal.II/base/utilities.h>
#include <deal.II/base/conditional_ostream.h>
#include <deal.II/base/index_set.h>
#include <deal.II/differentiation/ad.h>
#include <deal.II/distributed/tria.h>
#include <deal.II/distributed/grid_refinement.h>
#include <fstream>
#include <iostream>
constexpr double time_step = 0.0001;
namespace Step40
{
using namespace dealii;
template <int dim>
class InitialCondition : public Function<dim>
{
public:
InitialCondition() : Function<dim>(1){
}
virtual double value(const Point<dim> &p, const unsigned int component) const override;
virtual Tensor<1, dim> gradient(const Point<dim> &p, const unsigned int component) const override;
};
template <int dim>
double InitialCondition<dim>::value(const Point<dim> &p, const unsigned int) const{
const double x = p[0];
const double y = p[1];
return sin(M_PI * x) * sin(M_PI * y);
}
template <int dim>
Tensor<1, dim> InitialCondition<dim>::gradient(const Point<dim> &p, const unsigned int) const{
const double x = p[0];
const double y = p[1];
Tensor<1, dim> return_value;
return_value[0] = M_PI * cos(M_PI * x) * sin(M_PI * y);
return_value[1] = M_PI * sin(M_PI * x) * cos(M_PI * y);
return return_value;
}
template <int dim>
class LaplaceProblem
{
public:
LaplaceProblem();
void run();
private:
void setup_system();
void assemble_jacobian_by_hand_dynamic();
void assemble_jacobian_by_AD_dynamic();
void solve();
void solve_with_newton_hand();
void solve_with_newton_AD();
void refine_grid();
void output_results(const unsigned int cycle) const;
MPI_Comm mpi_communicator;
parallel::distributed::Triangulation<dim> triangulation;
FE_Q<dim> fe;
DoFHandler<dim> dof_handler;
IndexSet locally_owned_dofs;
IndexSet locally_relevant_dofs;
AffineConstraints<double> constraints;
LA::MPI::SparseMatrix system_matrix, jacobian_matrix_AD, jacobian_matrix_hand, mass_matrix;
LA::MPI::Vector locally_relevant_solution,
locally_relevant_solution_AD,
locally_relevant_solution_hand;
LA::MPI::Vector system_rhs, residual_rhs_AD, residual_rhs_hand;
ConditionalOStream pcout;
TimerOutput computing_timer;
};
template <int dim>
LaplaceProblem<dim>::LaplaceProblem()
: mpi_communicator(MPI_COMM_WORLD)
, triangulation(mpi_communicator,
typename Triangulation<dim>::MeshSmoothing(
Triangulation<dim>::smoothing_on_refinement |
Triangulation<dim>::smoothing_on_coarsening))
, fe(2)
, dof_handler(triangulation)
, pcout(std::cout,
(Utilities::MPI::this_mpi_process(mpi_communicator) == 0))
, computing_timer(mpi_communicator,
pcout,
TimerOutput::summary,
TimerOutput::wall_times)
{}
template <int dim>
void LaplaceProblem<dim>::setup_system()
{
TimerOutput::Scope t(computing_timer, "setup");
dof_handler.distribute_dofs(fe);
locally_owned_dofs = dof_handler.locally_owned_dofs();
DoFTools::extract_locally_relevant_dofs(dof_handler, locally_relevant_dofs);
locally_relevant_solution.reinit(locally_owned_dofs,
locally_relevant_dofs,
mpi_communicator);
locally_relevant_solution_AD.reinit(locally_owned_dofs,
locally_relevant_dofs,
mpi_communicator);
locally_relevant_solution_hand.reinit(locally_owned_dofs,
locally_relevant_dofs,
mpi_communicator);
system_rhs.reinit(locally_owned_dofs, mpi_communicator);
residual_rhs_hand.reinit(locally_owned_dofs, mpi_communicator);
residual_rhs_AD.reinit(locally_owned_dofs, mpi_communicator);
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();
DynamicSparsityPattern dsp(locally_relevant_dofs);
DoFTools::make_sparsity_pattern(dof_handler, dsp, constraints, false);
SparsityTools::distribute_sparsity_pattern(
dsp,
dof_handler.compute_n_locally_owned_dofs_per_processor(),
mpi_communicator,
locally_relevant_dofs);
system_matrix.reinit(locally_owned_dofs,
locally_owned_dofs,
dsp,
mpi_communicator);
jacobian_matrix_AD.reinit(locally_owned_dofs,
locally_owned_dofs,
dsp,
mpi_communicator);
jacobian_matrix_hand.reinit(locally_owned_dofs,
locally_owned_dofs,
dsp,
mpi_communicator);
mass_matrix.reinit(locally_owned_dofs,
locally_owned_dofs,
dsp,
mpi_communicator);
}
template <int dim>
void LaplaceProblem<dim>::assemble_jacobian_by_AD_dynamic()
{
TimerOutput::Scope t(computing_timer, "jacobian assembly, by AD");
jacobian_matrix_hand = 0.;
jacobian_matrix_AD = 0.;
residual_rhs_hand = 0.;
residual_rhs_AD = 0.;
const QGauss<dim> quadrature_formula(fe.degree + 1);
FEValues<dim> fe_values(fe,
quadrature_formula,
update_values | 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<types::global_dof_index> local_dof_indices(dofs_per_cell);
namespace AD = dealii::Differentiation::AD;
typedef AD::VectorFunction<dim, AD::NumberTypes::sacado_dfad, double> ADHelper;
typedef typename ADHelper::ad_type ADNumberType;
typedef typename ADHelper::scalar_type ScalarNumberType;
std::vector<Tensor<1, dim>> old_solution_gradients(n_q_points);
std::vector<double> old_solution_values(n_q_points);
const size_t tape_no = 1;
for (const auto &cell : dof_handler.active_cell_iterators())
if (cell->is_locally_owned())
{
cell_matrix = 0.;
cell_rhs = 0.;
fe_values.reinit(cell);
cell->get_dof_indices(local_dof_indices);
fe_values.get_function_gradients(locally_relevant_solution_AD, old_solution_gradients);
fe_values.get_function_values(locally_relevant_solution_AD, old_solution_values);
const size_t n_independent_variables = local_dof_indices.size();
const size_t n_dependent_variables = dofs_per_cell;
ADHelper ad_helper(n_independent_variables, n_dependent_variables);
ad_helper.set_tape_buffer_sizes();
std::vector<double> dof_values(n_independent_variables);
for(size_t i = 0; i < n_independent_variables; ++i)
dof_values[i] = 0;
const bool is_recording = ad_helper.start_recording_operations(tape_no,
true,
true);
if(is_recording){
ad_helper.register_independent_variables(dof_values);
const std::vector<ADNumberType> dof_values_ad = ad_helper.get_sensitive_variables();
std::vector<ADNumberType> residual_ad(n_dependent_variables, ADNumberType(0.0));
//Getting data
std::vector<Tensor<1, dim, ADNumberType>> gradients(n_q_points);
std::vector<ADNumberType> values(n_q_points);
//Filling data with values
for(size_t q = 0; q < n_q_points; ++q){
for(size_t k = 0; k < dofs_per_cell; ++k){
gradients[q] += dof_values_ad[k] * fe_values.shape_grad(k, q);
values[q] += dof_values_ad[k] * fe_values.shape_value(k, q);
}
}
for(size_t q = 0; q < n_q_points; ++q){
for(size_t i = 0; i < dofs_per_cell; ++i){
residual_ad[i] += (-0.5
* (gradients[q]
+ old_solution_gradients[q])
* fe_values.shape_grad(i, q)
- 1./time_step
* (values[q]
- old_solution_values[q]
)
* fe_values.shape_value(i, q)
)
* fe_values.JxW(q);
}
}
ad_helper.register_dependent_variables(residual_ad);
ad_helper.stop_recording_operations(false);
}
else{
Assert(is_recording == true, ExcInternalError());
pcout << "Using recorded tape\n";
ad_helper.activate_recorded_tape(tape_no);
ad_helper.set_independent_variables(dof_values);
}
ad_helper.compute_values(cell_rhs);
cell_rhs *= -1.0;
ad_helper.compute_jacobian(cell_matrix);
constraints.distribute_local_to_global(cell_matrix,
cell_rhs,
local_dof_indices,
jacobian_matrix_AD,
residual_rhs_AD);
}
jacobian_matrix_AD.compress(VectorOperation::add);
residual_rhs_AD.compress(VectorOperation::add);
}
template <int dim>
void LaplaceProblem<dim>::assemble_jacobian_by_hand_dynamic()
{
TimerOutput::Scope t(computing_timer, "jacobian assembly, by hand");
jacobian_matrix_hand = 0.;
jacobian_matrix_AD = 0.;
residual_rhs_hand = 0.;
residual_rhs_AD = 0.;
const QGauss<dim> quadrature_formula(fe.degree + 1);
FEValues<dim> fe_values(fe,
quadrature_formula,
update_values | 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<types::global_dof_index> local_dof_indices(dofs_per_cell);
std::vector<Tensor<1, dim>> old_solution_gradients(n_q_points);
std::vector<double> old_solution_values(n_q_points);
for (const auto &cell : dof_handler.active_cell_iterators())
if (cell->is_locally_owned())
{
cell_matrix = 0.;
cell_rhs = 0.;
fe_values.reinit(cell);
fe_values.get_function_gradients(locally_relevant_solution_hand, old_solution_gradients);
fe_values.get_function_values(locally_relevant_solution_hand, old_solution_values);
for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
{
for (unsigned int i = 0; i < dofs_per_cell; ++i)
{
for (unsigned int j = 0; j < dofs_per_cell; ++j)
cell_matrix(i, j) += (-0.5
* fe_values.shape_grad(i, q_point)
* fe_values.shape_grad(j, q_point)
- 1./time_step
* fe_values.shape_value(i, q_point)
* fe_values.shape_value(j, q_point))
* fe_values.JxW(q_point);
cell_rhs(i) += -1.
* -1.
* (old_solution_gradients[q_point]
* fe_values.shape_grad(i, q_point))
* fe_values.JxW(q_point);
}
}
cell->get_dof_indices(local_dof_indices);
constraints.distribute_local_to_global(cell_matrix,
cell_rhs,
local_dof_indices,
jacobian_matrix_hand,
residual_rhs_hand);
}
jacobian_matrix_hand.compress(VectorOperation::add);
residual_rhs_hand.compress(VectorOperation::add);
}
template <int dim>
void LaplaceProblem<dim>::solve_with_newton_hand()
{
TimerOutput::Scope t(computing_timer, "solve Hand-Newton");
LA::MPI::Vector completely_distributed_solution(locally_owned_dofs,
mpi_communicator);
SolverControl solver_control(dof_handler.n_dofs(), 1e-12);
#ifdef USE_PETSC_LA
LA::SolverGMRES solver(solver_control, mpi_communicator);
#else
LA::SolverGMRES solver(solver_control);
#endif
LA::MPI::PreconditionAMG preconditioner;
LA::MPI::PreconditionAMG::AdditionalData data;
#ifdef USE_PETSC_LA
data.symmetric_operator = true;
#else
/* Trilinos defaults are good */
#endif
preconditioner.initialize(jacobian_matrix_hand, data);
solver.solve(jacobian_matrix_hand,
completely_distributed_solution,
residual_rhs_hand,
preconditioner);
pcout << "Handmade jacobian solved in " << solver_control.last_step() << " iterations."
<< std::endl;
constraints.distribute(completely_distributed_solution);
locally_relevant_solution_hand += completely_distributed_solution;//Only correct for linear systems
}
template <int dim>
void LaplaceProblem<dim>::solve_with_newton_AD()
{
TimerOutput::Scope t(computing_timer, "solve AD-Newton");
LA::MPI::Vector completely_distributed_solution(locally_owned_dofs,
mpi_communicator);
SolverControl solver_control(dof_handler.n_dofs(), 1e-12);
#ifdef USE_PETSC_LA
LA::SolverGMRES solver(solver_control, mpi_communicator);
#else
LA::SolverGMRES solver(solver_control);
#endif
LA::MPI::PreconditionAMG preconditioner;
LA::MPI::PreconditionAMG::AdditionalData data;
#ifdef USE_PETSC_LA
data.symmetric_operator = true;
#else
/* Trilinos defaults are good */
#endif
preconditioner.initialize(jacobian_matrix_AD, data);
solver.solve(jacobian_matrix_AD,
completely_distributed_solution,
residual_rhs_AD,
preconditioner);
pcout << "AD-jacobian solved in " << solver_control.last_step() << " iterations."
<< std::endl;
constraints.distribute(completely_distributed_solution);
locally_relevant_solution_AD /*+=*/= completely_distributed_solution;//Only correct for linear systems
}
template <int dim>
void LaplaceProblem<dim>::refine_grid()
{
return;
}
template <int dim>
void LaplaceProblem<dim>::output_results(const unsigned int cycle) const
{
DataOut<dim> data_out;
data_out.attach_dof_handler(dof_handler);
data_out.add_data_vector(locally_relevant_solution, "u");
data_out.add_data_vector(locally_relevant_solution_AD, "u_AD");
data_out.add_data_vector(locally_relevant_solution_hand, "u_hand");
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();
data_out.write_vtu_with_pvtu_record("", "solution", cycle);
}
template <int dim>
void LaplaceProblem<dim>::run()
{
pcout << "Running with "
#ifdef USE_PETSC_LA
<< "PETSc"
#else
<< "Trilinos"
#endif
<< " on " << Utilities::MPI::n_mpi_processes(mpi_communicator)
<< " MPI rank(s)..." << std::endl;
const unsigned int n_cycles = 8;
for (unsigned int cycle = 0; cycle < n_cycles; ++cycle)
{
pcout << "Cycle " << cycle << ':' << std::endl;
if (cycle == 0)
{
GridGenerator::hyper_cube(triangulation);
triangulation.refine_global(5);
}
else
refine_grid();
//triangulation.refine_global(1);
if(cycle == 0){
setup_system();
LA::MPI::Vector local_solution;
local_solution.reinit(dof_handler.locally_owned_dofs(), mpi_communicator);
VectorTools::project(dof_handler,
constraints,
QGauss<dim>(fe.degree + 1),
InitialCondition<dim>(),
local_solution);
locally_relevant_solution = local_solution;
locally_relevant_solution_AD = local_solution;
locally_relevant_solution_hand = local_solution;
}
if (Utilities::MPI::n_mpi_processes(mpi_communicator) <= 32 && cycle == 0)
{
TimerOutput::Scope t(computing_timer, "output");
output_results(cycle);
}
pcout << " Number of active cells: "
<< triangulation.n_global_active_cells() << std::endl
<< " Number of degrees of freedom: " << dof_handler.n_dofs()
<< std::endl;
assemble_jacobian_by_AD_dynamic();
solve_with_newton_AD();
assemble_jacobian_by_hand_dynamic();
solve_with_newton_hand();
if (Utilities::MPI::n_mpi_processes(mpi_communicator) <= 32)
{
TimerOutput::Scope t(computing_timer, "output");
output_results(cycle + 1);
std::ofstream hand_out(std::string("rhs_hand_") + std::to_string(cycle) + std::string(".txt"));
std::ofstream AD_out(std::string("rhs_AD_") + std::to_string(cycle) + std::string(".txt"));
for(size_t i = 0; i < residual_rhs_hand.size(); ++i){
hand_out << residual_rhs_hand[i] << '\n';
}
for(size_t i = 0; i < residual_rhs_AD.size(); ++i){
AD_out << residual_rhs_AD[i] << '\n';
}
}
computing_timer.print_summary();
computing_timer.reset();
pcout << std::endl;
}
}
} // namespace Step40
int main(int argc, char *argv[])
{
std::cout << "Hello World\n";
try
{
using namespace dealii;
using namespace Step40;
#ifdef DEAL_II_WITH_PETSC
Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1);
#else
Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, dealii::numbers::invalid_unsigned_int);
#endif
LaplaceProblem<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;
}
