project_Hcurl_to_L2.hpp

/******************************************************************************
 * Copyright (C) Siarhei Uzunbajakau, 2023.
 *
 * This program is free software. You can use, modify, and redistribute it under
 * the terms of the GNU Lesser General Public License as published by the Free
 * Software Foundation, either version 3 or (at your option) any later version.
 * This program is distributed without any warranty.
 *
 * Refer to COPYING.LESSER for more details.
 ******************************************************************************/
 
#ifndef ProjectHcurlToL2_H__
#define ProjectHcurlToL2_H__
 
#define BOOST_ALLOW_DEPRECATED_HEADERS
 
#include <deal.II/base/timer.h>
#include <deal.II/base/work_stream.h>
 
#include <deal.II/grid/tria.h>
 
#include <deal.II/dofs/dof_handler.h>
#include <deal.II/dofs/dof_renumbering.h>
#include <deal.II/dofs/dof_tools.h>
 
#include <deal.II/lac/full_matrix.h>
#include <deal.II/lac/sparse_direct.h>
 
#include <deal.II/lac/precondition.h>
#include <deal.II/lac/solver_cg.h>
#include <deal.II/lac/solver_control.h>
 
#include <deal.II/fe/fe_dgq.h>
 
#include <deal.II/fe/fe_values.h>
#include <deal.II/fe/mapping_q1.h>
 
#include <deal.II/numerics/data_out.h>
#include <deal.II/numerics/matrix_tools.h>
#include <deal.II/numerics/vector_tools.h>
 
#include <fstream>
#include <iomanip>
#include <ios>
#include <iostream>
#include <string>
 
#include "constants.hpp"
#include "static_vector_input.hpp"
 
#define SE scratch_data.se
 
#define TMR(__name) TimerOutput::Scope timer_section(timer, __name)
 
using namespace dealii;
 
namespace StaticVectorSolver {
template<int stage = 1>
class ProjectHcurlToL2
{
public:
  ProjectHcurlToL2() = delete;
 
  ProjectHcurlToL2(unsigned int p,
                   unsigned int mapping_degree,
                   const Triangulation<2>& triangulation_Hcurl,
                   const DoFHandler<2>& dof_handler_Hcurl,
                   const Vector<double>& solution_Hcurl,
                   std::string fname = "L2",
                   const Function<2>* exact_solution = nullptr,
                   bool print_time_tables = false,
                   bool project_exact_solution = false,
                   bool log_cg_convergence = false,
                   bool write_higher_order_cells = false);
 
  double get_L2_norm() { return L2_norm; };
 
  double get_Linfty_norm() { return Linfty_norm; }
 
  unsigned int get_n_cells() const
  {
    return static_cast<unsigned int>(triangulation_Hcurl.n_active_cells());
  }
 
  unsigned int get_n_dofs() const
  {
    return static_cast<unsigned int>(dof_handler_L2.n_dofs());
  }
 
  void clear()
  {
    system_matrix.clear();
    system_rhs.reinit(0);
  }
 
  const Triangulation<2>& get_tria() const { return triangulation_Hcurl; }
 
  const DoFHandler<2>& get_dof_handler() const { return dof_handler_L2; }
 
  const Vector<double>& get_solution() const { return solution_L2; }
 
  void save_matrix_and_rhs_to_csv(std::string fname) const;
 
private:
  void setup();
  void assemble();
  void solve();
  void save() const;
  void compute_error_norms();
  void project_exact_solution_fcn();
 
  const std::string fname;
 
  const DoFHandler<2>& dof_handler_Hcurl;
  const Vector<double>& solution_Hcurl;
 
  const Triangulation<2>& triangulation_Hcurl;
  const FE_DGQ<2> fe_L2;
  DoFHandler<2> dof_handler_L2;
 
  SparsityPattern sparsity_pattern;
  SparseMatrix<double> system_matrix;
 
  Vector<double> solution_L2;
  Vector<double> system_rhs;
 
  Vector<double> projected_exact_solution;
 
  AffineConstraints<double> constraints;
 
  const Function<2>* exact_solution;
 
  const unsigned int mapping_degree;
  const bool project_exact_solution;
  const bool log_cg_convergence;
  const bool write_higher_order_cells;
 
  Vector<double> L2_per_cell;
  double L2_norm;
 
  Vector<double> Linfty_per_cell;
  double Linfty_norm;
 
  // ----------------------------------------------------------------------------
  // These structures and functions are related to the Work Stream algorithm.
  // See article "WorkStream – A Design Pattern for Multicore-Enabled Finite
  // Element Computations." by BRUNO TURCKSIN, MARTIN KRONBICHLER,
  // WOLFGANG BANGERTH for more details.
  // ----------------------------------------------------------------------------
 
  using IteratorTuple =
    std::tuple<typename DoFHandler<2>::active_cell_iterator,
               typename DoFHandler<2>::active_cell_iterator>;
 
  using IteratorPair = SynchronousIterators<IteratorTuple>;
 
  struct AssemblyScratchData
  {
    AssemblyScratchData(const FiniteElement<2>& fe,
                        const DoFHandler<2>& dof_handr_Hcurl,
                        const Vector<double>& dofs_Hcurl,
                        unsigned int mapping_degree);
 
    AssemblyScratchData(const AssemblyScratchData& scratch_data);
 
    MappingQ<2> mapping;
    Constants::QuadratureTableVector<2> qt;
    FEValues<2> fe_values_L2;
    FEValues<2> fe_values_Hcurl;
 
    const unsigned int dofs_per_cell;
    const unsigned int n_q_points;
 
    std::vector<std::vector<Tensor<1, 2>>> vector_gradients;
    Tensor<1, 1> curl_vec_in_Hcurl;
 
    const FEValuesExtractors::Scalar se;
 
    const DoFHandler<2>& dof_hand_Hcurl;
    const Vector<double>& dofs_Hcurl;
  };
 
  struct AssemblyCopyData
  {
    FullMatrix<double> cell_matrix;
    Vector<double> cell_rhs;
    std::vector<types::global_dof_index> local_dof_indices;
  };
 
  void system_matrix_local(const IteratorPair& IP,
                           AssemblyScratchData& scratch_data,
                           AssemblyCopyData& copy_data);
 
  void copy_local_to_global(const AssemblyCopyData& copy_data);
 
  //-----------------------------------------------------------------------------
  //-----------------------------------------------------------------------------
  //-----------------------------------------------------------------------------
};
 
template<int stage>
ProjectHcurlToL2<stage>::ProjectHcurlToL2(
  unsigned int p,
  unsigned int mapping_degree,
  const Triangulation<2>& triangulation_Hcurl,
  const DoFHandler<2>& dof_handler_Hcurl,
  const Vector<double>& solution_Hcurl,
  std::string fname,
  const Function<2>* exact_solution,
  bool print_time_tables,
  bool project_exact_solution,
  bool log_cg_convergence,
  bool write_higher_order_cells)
  : fname(fname)
  , dof_handler_Hcurl(dof_handler_Hcurl)
  , solution_Hcurl(solution_Hcurl)
  , triangulation_Hcurl(triangulation_Hcurl)
  , fe_L2(p)
  , exact_solution(exact_solution)
  , mapping_degree(mapping_degree)
  , project_exact_solution(project_exact_solution)
  , log_cg_convergence(log_cg_convergence)
  , write_higher_order_cells(write_higher_order_cells)
{
  TimerOutput::OutputFrequency tf =
    (print_time_tables) ? TimerOutput::summary : TimerOutput::never;
 
  TimerOutput timer(std::cout, tf, TimerOutput::cpu_and_wall_times_grouped);
 
  {
    TMR("Setup");
    setup();
  }
  {
    TMR("Assemble");
    assemble();
  }
  {
    TMR("Solve");
    solve();
  }
 
  if (exact_solution) {
    {
      TMR("Compute error norms");
      compute_error_norms();
    }
 
    if (project_exact_solution) {
      {
        TMR("Project exact solution");
        project_exact_solution_fcn();
      }
    }
  }
 
  {
    TMR("Save");
    save();
  }
}
 
template<int stage>
void
ProjectHcurlToL2<stage>::save_matrix_and_rhs_to_csv(std::string fname) const
{
  std::ofstream ofs_matrix(fname + "_matrix.csv");
  std::ofstream ofs_rhs(fname + "_rhs.csv");
 
  for (unsigned int i = 0; i < system_matrix.m(); ++i) {
    ofs_rhs << system_rhs(i);
    if (i < (system_matrix.m() - 1))
      ofs_rhs << "\n";
 
    for (unsigned int j = 0; j < system_matrix.n(); ++j) {
      ofs_matrix << std::scientific << std::setprecision(16)
                 << system_matrix.el(i, j);
 
      if (j < (system_matrix.m() - 1))
        ofs_matrix << ", ";
    }
    if (i < (system_matrix.m() - 1))
      ofs_matrix << "\n";
  }
 
  ofs_rhs.close();
  ofs_matrix.close();
}
 
template<int stage>
void
ProjectHcurlToL2<stage>::setup()
{
  constraints.close();
 
  dof_handler_L2.reinit(triangulation_Hcurl);
  dof_handler_L2.distribute_dofs(fe_L2);
  constraints.close();
 
  DynamicSparsityPattern dsp(dof_handler_L2.n_dofs(), dof_handler_L2.n_dofs());
  DoFTools::make_sparsity_pattern(dof_handler_L2, dsp, constraints, false);
 
  sparsity_pattern.copy_from(dsp);
  system_matrix.reinit(sparsity_pattern);
  solution_L2.reinit(dof_handler_L2.n_dofs());
  system_rhs.reinit(dof_handler_L2.n_dofs());
 
  if (project_exact_solution)
    projected_exact_solution.reinit(dof_handler_L2.n_dofs());
 
  if (exact_solution) {
    L2_per_cell.reinit(triangulation_Hcurl.n_active_cells());
    Linfty_per_cell.reinit(triangulation_Hcurl.n_active_cells());
  }
}
 
template<int stage>
void
ProjectHcurlToL2<stage>::assemble()
{
  WorkStream::run(
    IteratorPair(IteratorTuple(dof_handler_L2.begin_active(),
                               dof_handler_Hcurl.begin_active())),
    IteratorPair(IteratorTuple(dof_handler_L2.end(), dof_handler_Hcurl.end())),
    *this,
    &ProjectHcurlToL2<stage>::system_matrix_local,
    &ProjectHcurlToL2<stage>::copy_local_to_global,
    AssemblyScratchData(
      fe_L2, dof_handler_Hcurl, solution_Hcurl, mapping_degree),
    AssemblyCopyData());
}
 
template<int stage>
ProjectHcurlToL2<stage>::AssemblyScratchData::AssemblyScratchData(
  const FiniteElement<2>& fe,
  const DoFHandler<2>& dof_hand_Hcurl,
  const Vector<double>& dofs_Hcurl,
  unsigned int mapping_degree)
  : mapping(mapping_degree)
  , qt(fe.degree)
  , fe_values_L2(mapping,
                 fe,
                 QGauss<2>(qt.sim()),
                 update_values | update_quadrature_points | update_JxW_values)
  , fe_values_Hcurl(mapping,
                    dof_hand_Hcurl.get_fe(),
                    QGauss<2>(qt.sim()),
                    update_gradients)
  , dofs_per_cell(fe_values_L2.dofs_per_cell)
  , n_q_points(fe_values_L2.get_quadrature().size())
  , vector_gradients(n_q_points, std::vector<Tensor<1, 2>>(2))
  , se(0)
  , dof_hand_Hcurl(dof_hand_Hcurl)
  , dofs_Hcurl(dofs_Hcurl)
{
}
 
template<int stage>
ProjectHcurlToL2<stage>::AssemblyScratchData::AssemblyScratchData(
  const AssemblyScratchData& scratch_data)
  : mapping(scratch_data.mapping.get_degree())
  , qt(scratch_data.qt)
  , fe_values_L2(mapping,
                 scratch_data.fe_values_L2.get_fe(),
                 scratch_data.fe_values_L2.get_quadrature(),
                 update_values | update_quadrature_points | update_JxW_values)
  , fe_values_Hcurl(mapping,
                    scratch_data.fe_values_Hcurl.get_fe(),
                    scratch_data.fe_values_Hcurl.get_quadrature(),
                    update_gradients)
  , dofs_per_cell(fe_values_L2.dofs_per_cell)
  , n_q_points(fe_values_L2.get_quadrature().size())
  , vector_gradients(n_q_points, std::vector<Tensor<1, 2>>(2))
  , se(0)
  , dof_hand_Hcurl(scratch_data.dof_hand_Hcurl)
  , dofs_Hcurl(scratch_data.dofs_Hcurl)
{
}
 
template<int stage>
void
ProjectHcurlToL2<stage>::system_matrix_local(const IteratorPair& IP,
                                             AssemblyScratchData& scratch_data,
                                             AssemblyCopyData& copy_data)
{
  // See the color box
  // Recipe for projections from H(curl) to L_2 nr. 14
 
  copy_data.cell_matrix.reinit(scratch_data.dofs_per_cell,
                               scratch_data.dofs_per_cell);
 
  copy_data.cell_rhs.reinit(scratch_data.dofs_per_cell);
 
  copy_data.local_dof_indices.resize(scratch_data.dofs_per_cell);
 
  scratch_data.fe_values_L2.reinit(std::get<0>(*IP));
  scratch_data.fe_values_Hcurl.reinit(std::get<1>(*IP));
 
  scratch_data.fe_values_Hcurl.get_function_gradients(
    scratch_data.dofs_Hcurl, scratch_data.vector_gradients);
 
  for (unsigned int q_index = 0; q_index < scratch_data.n_q_points; ++q_index) {
    for (unsigned int i = 0; i < scratch_data.dofs_per_cell; ++i) {
      for (unsigned int j = 0; j < scratch_data.dofs_per_cell; ++j) {
        copy_data.cell_matrix(i, j) +=
          scratch_data.fe_values_L2[SE].value(i, q_index) *
          scratch_data.fe_values_L2[SE].value(j, q_index) *
          scratch_data.fe_values_L2.JxW(q_index);
      }
      scratch_data.curl_vec_in_Hcurl[0] =
        scratch_data.vector_gradients[q_index][1][0] -
        scratch_data.vector_gradients[q_index][0][1];
 
      copy_data.cell_rhs(i) += scratch_data.curl_vec_in_Hcurl[0] *
                               scratch_data.fe_values_L2[SE].value(i, q_index) *
                               scratch_data.fe_values_L2.JxW(q_index);
    }
  }
 
  std::get<0>(*IP)->get_dof_indices(copy_data.local_dof_indices);
}
 
template<int stage>
void
ProjectHcurlToL2<stage>::copy_local_to_global(const AssemblyCopyData& copy_data)
{
  constraints.distribute_local_to_global(copy_data.cell_matrix,
                                         copy_data.cell_rhs,
                                         copy_data.local_dof_indices,
                                         system_matrix,
                                         system_rhs);
}
 
template<int stage>
void
ProjectHcurlToL2<stage>::solve()
{
  SolverControl control(
    1000 * system_rhs.size(), 1e-12 * system_rhs.l2_norm(), false, false);
 
  if (log_cg_convergence)
    control.enable_history_data();
 
  GrowingVectorMemory<Vector<double>> memory;
  SolverCG<Vector<double>> cg(control, memory);
 
  PreconditionSSOR<SparseMatrix<double>> preconditioner;
  preconditioner.initialize(system_matrix, 1.2);
 
  cg.solve(system_matrix, solution_L2, system_rhs, preconditioner);
 
  if (log_cg_convergence) {
    const std::vector<double> history_data = control.get_history_data();
 
    std::ofstream ofs(fname + "_cg_convergence.csv");
 
    unsigned int i = 1;
    for (auto item : history_data) {
      ofs << i << ", " << item << "\n";
      i++;
    }
    ofs.close();
  }
}
 
template<int stage>
void
ProjectHcurlToL2<stage>::save() const
{
  std::vector<std::string> solution_names(1, "ScalarField");
  std::vector<DataComponentInterpretation::DataComponentInterpretation>
    data_component_interpretation(
      1, DataComponentInterpretation::component_is_scalar);
 
  DataOut<2> data_out;
  data_out.attach_dof_handler(dof_handler_L2);
  data_out.add_data_vector(solution_L2,
                           solution_names,
                           DataOut<2>::type_dof_data,
                           data_component_interpretation);
 
  if (project_exact_solution && exact_solution) {
    solution_names[0] = "ScalarFieldExact";
 
    data_out.add_data_vector(projected_exact_solution,
                             solution_names,
                             DataOut<2>::type_dof_data,
                             data_component_interpretation);
  }
 
  if (exact_solution) {
    solution_names[0] = "L2norm";
 
    data_out.add_data_vector(L2_per_cell,
                             solution_names,
                             DataOut<2>::type_cell_data,
                             data_component_interpretation);
 
    solution_names[0] = "LinftyNorm";
 
    data_out.add_data_vector(Linfty_per_cell,
                             solution_names,
                             DataOut<2>::type_cell_data,
                             data_component_interpretation);
  }
 
  std::ofstream ofs;
 
  if (write_higher_order_cells) {
    DataOutBase::VtkFlags flags;
    flags.write_higher_order_cells = true;
    data_out.set_flags(flags);
 
    const MappingQ<2> mapping(mapping_degree);
 
    data_out.build_patches(mapping,
                           fe_L2.degree + 2,
                           DataOut<2>::CurvedCellRegion::curved_inner_cells);
 
    ofs.open(fname + ".vtu");
    data_out.write_vtu(ofs);
 
  } else {
 
    data_out.build_patches();
 
    ofs.open(fname + ".vtk");
    data_out.write_vtk(ofs);
  }
 
  ofs.close();
}
 
template<int stage>
void
ProjectHcurlToL2<stage>::compute_error_norms()
{
  Weight<2, stage> weight;
  const Function<2, double>* mask = &weight;
 
  Constants::QuadratureTableScalar<2> qt(dof_handler_L2.get_fe().degree + 1);
 
  VectorTools::integrate_difference(MappingQ<2>(mapping_degree),
                                    dof_handler_L2,
                                    solution_L2,
                                    *exact_solution,
                                    L2_per_cell,
                                    QGauss<2>(qt.enorm()),
                                    VectorTools::L2_norm,
                                    mask);
 
  L2_norm = VectorTools::compute_global_error(
    triangulation_Hcurl, L2_per_cell, VectorTools::L2_norm);
 
  VectorTools::integrate_difference(MappingQ<2>(mapping_degree),
                                    dof_handler_L2,
                                    solution_L2,
                                    *exact_solution,
                                    Linfty_per_cell,
                                    QGauss<2>(1),
                                    VectorTools::Linfty_norm,
                                    mask // & B_mask
  );
 
  Linfty_norm = Linfty_per_cell.linfty_norm();
}
 
template<int stage>
void
ProjectHcurlToL2<stage>::project_exact_solution_fcn()
{
  Constants::QuadratureTableScalar<2> qt(fe_L2.degree + 1);
 
  AffineConstraints<double> constraints_empty;
  constraints_empty.close();
 
  VectorTools::project(MappingQ<2>(mapping_degree),
                       dof_handler_L2,
                       constraints_empty,
                       QGauss<2>(qt.sim()),
                       *exact_solution,
                       projected_exact_solution);
}
 
} // namespace StaticVectorSolver
 
#endif