When applying the laplace operator to my solution vector, what is the
difference between forming an explicit matrix containing the laplace
operator by using
template <int dim>
void LaplaceProblem<dim>::assemble_mass_matrix (){
TimerOutput::Scope t(computing_timer, "Mass matrix assembly");
mass_matrix = 0.;
system_matrix = 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_mass_matrix(dofs_per_cell,
dofs_per_cell),
cell_matrix(dofs_per_cell,
dofs_per_cell);
std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell
);
for (const auto &cell : dof_handler.active_cell_iterators())
if (cell->is_locally_owned())
{
cell_mass_matrix = 0.;
cell_matrix = 0.;
fe_values.reinit(cell);
for (unsigned int q = 0; q < n_q_points; ++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) += -1.
* fe_values.shape_grad (i, q)
* fe_values.shape_grad (j, q)
* fe_values.JxW (q);
cell_mass_matrix(i, j) += (fe_values.shape_value
(i, q)
* fe_values.
shape_value(j, q))
* fe_values.JxW(q);
}
}
}
cell->get_dof_indices(local_dof_indices);
constraints.distribute_local_to_global(cell_mass_matrix,
local_dof_indices,
mass_matrix);
constraints.distribute_local_to_global(cell_matrix,
local_dof_indices,
system_matrix);
}
mass_matrix.compress(VectorOperation::add);
system_matrix.compress(VectorOperation::add);
}
and assembling the result vector directly by getting the gradients from the
solution vector via get_function_gradients() and multiplying it afterwards
with shape_grad(), by using
template <int dim>
void LaplaceProblem<dim>::assemble_rhs(const vector_t &y, const double
local_time, vector_t &dst){
TimerOutput::Scope t(computing_timer, "RHS assembly");
dst.reinit(locally_owned_dofs, mpi_communicator);
dst = 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();
Vector<double> cell_rhs(dofs_per_cell);
std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell
);
std::vector<Tensor<1, dim>> solution_gradients(n_q_points);
for (const auto &cell : dof_handler.active_cell_iterators())
if (cell->is_locally_owned())
{
cell_rhs = 0.;
fe_values.reinit(cell);
fe_values.get_function_gradients(y, solution_gradients);
for (unsigned int q = 0; q < n_q_points; ++q)
{
for (unsigned int i = 0; i < dofs_per_cell; ++i)
{
cell_rhs(i) = (-1.
* solution_gradients[q]
* fe_values.shape_grad(i, q))
* fe_values.JxW(q);
}
}
cell->get_dof_indices(local_dof_indices);
constraints.distribute_local_to_global(cell_rhs,
local_dof_indices,
dst);
}
dst.compress(VectorOperation::add);
}
I wrote a short test program which should solve the diffusion equation
using the time-stepping class, and implemented both methods. When
calculating the matrix and applying it to my solution vector, I get a
different result compared to reading the gradients from the solution with
get_function_gradients() and multiplying it with the gradients returned by
shape_grad(). The results obtained from the matrix multiplication are
correct compared to the expected solution, the results obtained from the
direct approach are not. Why is that?
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 dealii+unsubscr...@googlegroups.com.
To view this discussion on the web visit
https://groups.google.com/d/msgid/dealii/f1c59802-44ba-4398-85db-3a2a27f316ca%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>
#include <deal.II/base/time_stepping.h>
#include <deal.II/base/utilities.h>
#include <deal.II/base/conditional_ostream.h>
#include <deal.II/base/index_set.h>
#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/lac/sparse_direct.h>
#include <deal.II/lac/la_parallel_vector.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/matrix_free/fe_evaluation.h>
#include <deal.II/matrix_free/matrix_free.h>
#include <deal.II/matrix_free/operators.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.001;
constexpr unsigned int fe_degree = 2;
constexpr unsigned int n_q_points_1d = fe_degree + 2;
using Number = double;
using vector_t = dealii::LinearAlgebra::distributed::Vector<Number>;
constexpr bool use_matrix = true;
namespace Step40
{
using namespace dealii;
void print_vector(const vector_t &input_vector, const std::string filename){
std::ofstream data_out(filename);
for(size_t i = 0; i < input_vector.size(); ++i)
data_out << input_vector(i) << '\n';
data_out.close();
}
template <class MatrixType>
class InverseMatrix : public Subscriptor
{
public:
InverseMatrix(const MatrixType & m,
const IndexSet &locally_owned_dofs);
void vmult(vector_t &dst,
const vector_t &src) const;
private:
const SmartPointer<const MatrixType> matrix;
const IndexSet locally_owned_dofs;
};
template <class MatrixType>
InverseMatrix<MatrixType>::InverseMatrix(
const MatrixType & m,
const IndexSet &locally_owned_dofs)
: matrix(&m),
locally_owned_dofs(locally_owned_dofs)
{}
template <class MatrixType>
void InverseMatrix<MatrixType>::vmult(
vector_t & dst,
const vector_t &src) const
{
SolverControl solver_control(src.size(), 1e-6 * src.l2_norm());
SolverCG<vector_t> cg(solver_control);
dst.reinit(src);
dst = 0;
cg.solve(*matrix, dst, src, PreconditionIdentity());
}
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 Solution : public Function<dim>
{
public:
Solution(const double time) : Function<dim>(1), solution_local_time(time) {
}
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;
private:
const double solution_local_time;
};
template <int dim>
double Solution<dim>::value(const Point<dim> &p, const unsigned int) const{
const double x = p[0];
const double y = p[1];
return exp(-2 * M_PI * M_PI * solution_local_time) * sin(M_PI * x) * sin(M_PI * y);
}
template <int dim>
Tensor<1, dim> Solution<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] = exp(-2 * M_PI * M_PI * solution_local_time) * M_PI * cos(M_PI * x) * sin(M_PI * y);
return_value[1] = exp(-2 * M_PI * M_PI * solution_local_time) * M_PI * sin(M_PI * x) * cos(M_PI * y);
return return_value;
}
template <int dim>
class LaplaceProblem
{
public:
LaplaceProblem();
void run_with_MF();
private:
void setup_system();
void assemble_mass_matrix();
void assemble_rhs(const vector_t &y, const double local_time, vector_t &dst);
void output_results(const unsigned int cycle, const double time) const;
vector_t evaluate_diffusion(const double time,
const vector_t &y);
unsigned int embedded_explicit_method(const TimeStepping::runge_kutta_method method,
const unsigned int n_time_steps,
const double initial_time,
const double final_time);
MPI_Comm mpi_communicator;
parallel::distributed::Triangulation<dim> triangulation;
FE_Q<dim> fe;
MappingQGeneric<dim> mapping;
DoFHandler<dim> dof_handler;
IndexSet locally_owned_dofs;
IndexSet locally_relevant_dofs;
AffineConstraints<double> constraints;
vector_t locally_relevant_solution;
SparsityPattern sparsity_pattern;
SparseMatrix<double> mass_matrix, system_matrix;
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(fe_degree)
, mapping(fe.degree + 1)
, 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);
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);
sparsity_pattern.copy_from(dsp);
mass_matrix.reinit(sparsity_pattern);
system_matrix.reinit(sparsity_pattern);
}
template <int dim>
void LaplaceProblem<dim>::assemble_mass_matrix (){
TimerOutput::Scope t(computing_timer, "Mass matrix assembly");
mass_matrix = 0.;
system_matrix = 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_mass_matrix(dofs_per_cell,
dofs_per_cell),
cell_matrix(dofs_per_cell,
dofs_per_cell);
std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
for (const auto &cell : dof_handler.active_cell_iterators())
if (cell->is_locally_owned())
{
cell_mass_matrix = 0.;
cell_matrix = 0.;
fe_values.reinit(cell);
for (unsigned int q = 0; q < n_q_points; ++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) += -1.
* fe_values.shape_grad (i, q)
* fe_values.shape_grad (j, q)
* fe_values.JxW (q);
cell_mass_matrix(i, j) += (fe_values.shape_value(i, q)
* fe_values.shape_value(j, q))
* fe_values.JxW(q);
}
}
}
cell->get_dof_indices(local_dof_indices);
constraints.distribute_local_to_global(cell_mass_matrix,
local_dof_indices,
mass_matrix);
constraints.distribute_local_to_global(cell_matrix,
local_dof_indices,
system_matrix);
}
mass_matrix.compress(VectorOperation::add);
system_matrix.compress(VectorOperation::add);
}
template <int dim>
void LaplaceProblem<dim>::assemble_rhs(const vector_t &y, const double local_time, vector_t &dst){
TimerOutput::Scope t(computing_timer, "RHS assembly");
dst.reinit(locally_owned_dofs, mpi_communicator);
dst = 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();
Vector<double> cell_rhs(dofs_per_cell);
std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
std::vector<Tensor<1, dim>> solution_gradients(n_q_points);
for (const auto &cell : dof_handler.active_cell_iterators())
if (cell->is_locally_owned())
{
cell_rhs = 0.;
fe_values.reinit(cell);
fe_values.get_function_gradients(y, solution_gradients);
for (unsigned int q = 0; q < n_q_points; ++q)
{
for (unsigned int i = 0; i < dofs_per_cell; ++i)
{
cell_rhs(i) = (-1.
* solution_gradients[q]
* fe_values.shape_grad(i, q))
* fe_values.JxW(q);
}
}
cell->get_dof_indices(local_dof_indices);
constraints.distribute_local_to_global(cell_rhs,
local_dof_indices,
dst);
}
dst.compress(VectorOperation::add);
}
template <int dim>
vector_t LaplaceProblem<dim>::evaluate_diffusion(const double time,
const vector_t &y){
vector_t tmp,
value;
tmp.reinit(dof_handler.locally_owned_dofs (),
mpi_communicator);
value.reinit(dof_handler.locally_owned_dofs (),
mpi_communicator);
tmp = 0.;
value = 0.;
InverseMatrix local_inverse_matrix(mass_matrix, locally_owned_dofs);
if(use_matrix)
system_matrix.vmult(tmp, y);
else
assemble_rhs(y, time, tmp);
local_inverse_matrix.vmult(value, tmp);
return value;
}
template <int dim>
unsigned int LaplaceProblem<dim>::embedded_explicit_method(const TimeStepping::runge_kutta_method method,
const unsigned int n_time_steps,
const double initial_time,
const double final_time){
TimerOutput::Scope t(computing_timer, "DOPRI step");
double local_time_step = (final_time - initial_time) / static_cast<double>(n_time_steps);
double time = initial_time;
const double coarsen_param = 1.2;
const double refine_param = 0.8;
const double min_delta = 1e-8;
const double max_delta = 10 * local_time_step;
const double refine_tol = 1e-1;
const double coarsen_tol = 1e-5;
TimeStepping::EmbeddedExplicitRungeKutta<vector_t>
embedded_explicit_runge_kutta(method,
coarsen_param,
refine_param,
min_delta,
max_delta,
refine_tol,
coarsen_tol);
TimeStepping::ExplicitRungeKutta <vector_t> explicit_runge_kutta(TimeStepping::FORWARD_EULER);
//output_results(0, method);
unsigned int n_steps = 0;
vector_t local_solution(locally_owned_dofs,
mpi_communicator);
local_solution = locally_relevant_solution;
while (time < final_time)
{
if (time + local_time_step > final_time)
local_time_step = final_time - time;
time = embedded_explicit_runge_kutta.evolve_one_time_step(
[this](const double time, const vector_t &y) {
return this->evaluate_diffusion(time, y);
},
time,
local_time_step,
local_solution);
constraints.distribute(local_solution);
local_time_step = embedded_explicit_runge_kutta.get_status().delta_t_guess;
++n_steps;
//pcout << "Time: " << time << '\n';
}
// for(size_t i = 0; i < n_time_steps; ++i){
// time = explicit_runge_kutta.evolve_one_time_step(
// [this](const double time, const vector_t &y){
// return this->evaluate_diffusion (time, y);
// },
// time,
// time_step,
// local_solution);
// constraints.distribute(local_solution);
// }
locally_relevant_solution = local_solution;
return n_steps;
}
template <int dim>
void LaplaceProblem<dim>::output_results(const unsigned int cycle, const double time) const
{
DataOut<dim> data_out;
data_out.attach_dof_handler(dof_handler);
vector_t exact_solution;
exact_solution.reinit(dof_handler.locally_owned_dofs (),
mpi_communicator);
VectorTools::interpolate(dof_handler,
Solution<dim>(time),
exact_solution);
data_out.add_data_vector(exact_solution, "Exact_solution");
locally_relevant_solution.update_ghost_values();
data_out.add_data_vector(locally_relevant_solution, "u");
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_with_MF(){
{
const unsigned int n_vect_number =
VectorizedArray<Number>::n_array_elements;
const unsigned int n_vect_bits = 8 * sizeof(Number) * n_vect_number;
pcout << "Running with "
<< Utilities::MPI::n_mpi_processes(MPI_COMM_WORLD)
<< " MPI processes" << std::endl;
pcout << "Vectorization over " << n_vect_number << " "
<< (std::is_same<Number, double>::value ? "doubles" : "floats")
<< " = " << n_vect_bits << " bits ("
<< Utilities::System::get_current_vectorization_level()
<< "), VECTORIZATION_LEVEL=" << DEAL_II_COMPILER_VECTORIZATION_LEVEL
<< std::endl;
}
const unsigned int n_cycles = 8;
double local_time = 0.;
double local_end_time = 0.;
GridGenerator::hyper_cube(triangulation);
triangulation.refine_global(5);
setup_system();
assemble_mass_matrix ();
vector_t local_solution(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;
output_results (0, 0.);
// assemble_rhs (local_solution, 0., locally_relevant_solution);
// output_results (1, 0.);
for (unsigned int cycle = 0; cycle < n_cycles; ++cycle)
{
pcout << "Cycle " << cycle << ':' << std::endl;
local_time = cycle * time_step;
local_end_time = (cycle + 1) * time_step;
//triangulation.refine_global(1);
Timer timer;
unsigned int timestep_number = 1;
timestep_number = embedded_explicit_method(TimeStepping::DOPRI,
1,
local_time,
local_end_time);
pcout << "Number of time steps: " << timestep_number << '\n';
timestep_number = 1;
if (Utilities::MPI::n_mpi_processes(mpi_communicator) <= 32)
{
output_results(cycle + 1, local_end_time);
}
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;
Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, dealii::numbers::invalid_unsigned_int);
LaplaceProblem<2> laplace_problem_2d;
laplace_problem_2d.run_with_MF();
}
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;
}
