heat_dg.hpp

Go to the documentation of this file.

/********************************************************************************
  Copyright (C) 2024 by the DUBeat authors.
 
  This file is part of DUBeat.
 
  DUBeat is free software; you can 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 3 of the License, or
  (at your option) any later version.
 
  DUBeat is distributed in the hope that it will be useful,
  but WITHOUT ANY WARRANTY; without even the implied warranty of
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  Lesser General Public License for more details.
 
  You should have received a copy of the GNU Lesser General Public License
  along with DUBeat.  If not, see <http://www.gnu.org/licenses/>.
********************************************************************************/
 
#ifndef HEAT_DG_HPP_
#define HEAT_DG_HPP_
 
#include <math.h>
 
#include <memory>
#include <vector>
 
#include "../source/DUBValues.hpp"
#include "../source/DUB_FEM_handler.hpp"
#include "../source/assemble_DG.hpp"
#include "../source/face_handler_DG.hpp"
#include "../source/model_DG_t.hpp"
#include "../source/volume_handler_DG.hpp"
#include "source/core_model.hpp"
#include "source/geometry/mesh_handler.hpp"
#include "source/init.hpp"
#include "source/io/data_writer.hpp"
#include "source/numerics/bc_handler.hpp"
#include "source/numerics/linear_solver_handler.hpp"
#include "source/numerics/preconditioner_handler.hpp"
#include "source/numerics/tools.hpp"
 
namespace DUBeat::models
{
  namespace heat_DG
  {
    class ExactSolution : public lifex::utils::FunctionDirichlet
    {
    public:
      ExactSolution()
        : lifex::utils::FunctionDirichlet()
      {}
 
      virtual double
      value(const dealii::Point<lifex::dim> &p,
            const unsigned int /*component*/ = 0) const override
      {
        if (lifex::dim == 2)
          return std::sin(2 * M_PI * p[0]) * std::sin(2 * M_PI * p[1]) *
                 std::exp(-5 * this->get_time());
        else
          return std::sin(2 * M_PI * p[0] + M_PI / 4) *
                 std::sin(2 * M_PI * p[1] + M_PI / 4) *
                 std::sin(2 * M_PI * p[2] + M_PI / 4) *
                 std::exp(-5 * this->get_time());
      }
    };
 
    class RightHandSide : public lifex::Function<lifex::dim>
    {
    public:
      RightHandSide()
        : lifex::Function<lifex::dim>()
      {}
 
      virtual double
      value(const dealii::Point<lifex::dim> &p,
            const unsigned int /*component*/ = 0) const override
      {
        if (lifex::dim == 2)
          return std::sin(2 * M_PI * p[0]) * std::sin(2 * M_PI * p[1]) *
                 std::exp(-5 * this->get_time()) * (8 * M_PI * M_PI - 5);
        else
          return std::sin(2 * M_PI * p[0] + M_PI / 4) *
                 std::sin(2 * M_PI * p[1] + M_PI / 4) *
                 std::sin(2 * M_PI * p[2] + M_PI / 4) *
                 std::exp(-5 * this->get_time()) * (12 * M_PI * M_PI - 5);
      }
    };
 
    class BCNeumann : public lifex::Function<lifex::dim>
    {
    public:
      BCNeumann()
        : lifex::Function<lifex::dim>()
      {}
 
      virtual double
      value(const dealii::Point<lifex::dim> &p,
            const unsigned int /*component*/ = 0) const override
      {
        if (lifex::dim == 2)
          return -2 * M_PI * std::sin(2 * M_PI * p[0]) *
                   std::cos(2 * M_PI * p[1]) * std::exp(-5 * this->get_time()) *
                   (std::abs(p[1]) < 1e-10) +
                 2 * M_PI * std::cos(2 * M_PI * p[0]) *
                   std::sin(2 * M_PI * p[1]) * std::exp(-5 * this->get_time()) *
                   (std::abs(p[0] - 1) < 1e-10) +
                 2 * M_PI * std::sin(2 * M_PI * p[0]) *
                   std::cos(2 * M_PI * p[1]) * std::exp(-5 * this->get_time()) *
                   (std::abs(p[1] - 1) < 1e-10) -
                 2 * M_PI * std::cos(2 * M_PI * p[0]) *
                   std::sin(2 * M_PI * p[1]) * std::exp(-5 * this->get_time()) *
                   (std::abs(p[0]) < 1e-10);
        else
          return 2 * M_PI * std::cos(2 * M_PI * p[0] + M_PI / 4) *
                   std::sin(2 * M_PI * p[1] + M_PI / 4) *
                   std::sin(2 * M_PI * p[2] + M_PI / 4) *
                   std::exp(-5 * this->get_time()) *
                   (std::abs(p[0] - 1) < 1e-10) -
                 2 * M_PI * std::cos(2 * M_PI * p[0] + M_PI / 4) *
                   std::sin(2 * M_PI * p[1] + M_PI / 4) *
                   std::sin(2 * M_PI * p[2] + M_PI / 4) *
                   std::exp(-5 * this->get_time()) * (std::abs(p[0]) < 1e-10) +
                 2 * M_PI * std::sin(2 * M_PI * p[0] + M_PI / 4) *
                   std::cos(2 * M_PI * p[1] + M_PI / 4) *
                   std::sin(2 * M_PI * p[2] + M_PI / 4) *
                   std::exp(-5 * this->get_time()) *
                   (std::abs(p[1] - 1) < 1e-10) -
                 2 * M_PI * std::sin(2 * M_PI * p[0] + M_PI / 4) *
                   std::cos(2 * M_PI * p[1] + M_PI / 4) *
                   std::sin(2 * M_PI * p[2] + M_PI / 4) *
                   std::exp(-5 * this->get_time()) * (std::abs(p[1]) < 1e-10) +
                 2 * M_PI * std::sin(2 * M_PI * p[0] + M_PI / 4) *
                   std::sin(2 * M_PI * p[1] + M_PI / 4) *
                   std::cos(2 * M_PI * p[2] + M_PI / 4) *
                   std::exp(-5 * this->get_time()) *
                   (std::abs(p[2] - 1) < 1e-10) -
                 2 * M_PI * std::sin(2 * M_PI * p[0] + M_PI / 4) *
                   std::sin(2 * M_PI * p[1] + M_PI / 4) *
                   std::cos(2 * M_PI * p[2] + M_PI / 4) *
                   std::exp(-5 * this->get_time()) * (std::abs(p[2]) < 1e-10);
      }
    };
 
    class GradExactSolution : public lifex::Function<lifex::dim>
    {
    public:
      GradExactSolution()
        : lifex::Function<lifex::dim>()
      {}
 
      virtual double
      value(const dealii::Point<lifex::dim> &p,
            const unsigned int               component = 0) const override
      {
        if (lifex::dim == 2)
          {
            if (component == 0) // x
              return 2 * M_PI * std::cos(2 * M_PI * p[0]) *
                     std::sin(2 * M_PI * p[1]) *
                     std::exp(-5 * this->get_time());
            else // y
              return 2 * M_PI * std::sin(2 * M_PI * p[0]) *
                     std::cos(2 * M_PI * p[1]) *
                     std::exp(-5 * this->get_time());
          }
        else // dim=3
          {
            if (component == 0) // x
              return 2 * M_PI * std::cos(2 * M_PI * p[0] + M_PI / 4) *
                     std::sin(2 * M_PI * p[1] + M_PI / 4) *
                     std::sin(2 * M_PI * p[2] + M_PI / 4) *
                     std::exp(-5 * this->get_time());
            if (component == 1) // y
              return 2 * M_PI * std::sin(2 * M_PI * p[0] + M_PI / 4) *
                     std::cos(2 * M_PI * p[1] + M_PI / 4) *
                     std::sin(2 * M_PI * p[2] + M_PI / 4) *
                     std::exp(-5 * this->get_time());
            else // z
              return 2 * M_PI * std::sin(2 * M_PI * p[0] + M_PI / 4) *
                     std::sin(2 * M_PI * p[1] + M_PI / 4) *
                     std::cos(2 * M_PI * p[2] + M_PI / 4) *
                     std::exp(-5 * this->get_time());
          }
      }
    };
  } // namespace heat_DG
 
  template <class basis>
  class HeatDG : public ModelDG_t<basis>
  {
  public:
    HeatDG<basis>()
      : ModelDG_t<basis>("Heat")
    {
      this->u_ex      = std::make_shared<heat_DG::ExactSolution>();
      this->grad_u_ex = std::make_shared<heat_DG::GradExactSolution>();
      this->f_ex      = std::make_shared<heat_DG::RightHandSide>();
      this->g_n       = std::make_shared<heat_DG::BCNeumann>();
    }
 
  private:
    void
    assemble_system() override;
  };
 
  template <class basis>
  void
  HeatDG<basis>::assemble_system()
  {
    this->matrix = 0;
    this->rhs    = 0;
 
    // The method is needed to define how the system matrix and rhs term are
    // defined for the heat problem with Neumann boundary conditions. The full
    // matrix is composed by different sub-matrices that are called with simple
    // capital letters. We refer here to the DG_Assemble methods for their
    // definition.
 
    // See DG_Assemble::local_V().
    dealii::FullMatrix<double> V(this->dofs_per_cell, this->dofs_per_cell);
    // See DG_Assemble::local_M().
    dealii::FullMatrix<double> M(this->dofs_per_cell, this->dofs_per_cell);
    // See DG_Assemble::local_SC().
    dealii::FullMatrix<double> SC(this->dofs_per_cell, this->dofs_per_cell);
    // See DG_Assemble::local_IC().
    dealii::FullMatrix<double> IC(this->dofs_per_cell, this->dofs_per_cell);
    // Transpose of IC.
    dealii::FullMatrix<double> IC_t(this->dofs_per_cell, this->dofs_per_cell);
    // See DG_Assemble::local_IN().
    dealii::FullMatrix<double> IN(this->dofs_per_cell, this->dofs_per_cell);
    // Transpose of IN.
    dealii::FullMatrix<double> IN_t(this->dofs_per_cell, this->dofs_per_cell);
    // See DG_Assemble::local_SN().
    dealii::FullMatrix<double> SN(this->dofs_per_cell, this->dofs_per_cell);
 
    // See DG_Assemble::local_rhs().
    dealii::Vector<double> cell_rhs(this->dofs_per_cell);
    // See DG_Assemble::local_rhs_edge_neumann().
    dealii::Vector<double> cell_rhs_edge(this->dofs_per_cell);
    // See DG_Assemble::local_u0_M_rhs().
    dealii::Vector<double> u0_rhs(this->dofs_per_cell);
 
    std::vector<lifex::types::global_dof_index> dof_indices(
      this->dofs_per_cell);
    dealii::IndexSet owned_dofs = this->dof_handler.locally_owned_dofs();
    const double    &alpha_bdf  = this->bdf_handler.get_alpha();
 
    for (const auto &cell : this->dof_handler.active_cell_iterators())
      {
        if (cell->is_locally_owned())
          {
            this->assemble->reinit(cell);
            dof_indices = this->dof_handler.get_dof_indices(cell);
 
            V = this->assemble->local_V();
            M = this->assemble->local_M();
            M /= this->prm_time_step;
            M *= alpha_bdf;
 
            cell_rhs = this->assemble->local_rhs(this->f_ex);
            u0_rhs =
              this->assemble->local_u0_M_rhs(this->solution_bdf, dof_indices);
            u0_rhs /= this->prm_time_step;
 
            this->matrix.add(dof_indices, V);
            this->matrix.add(dof_indices, M);
            this->rhs.add(dof_indices, cell_rhs);
            this->rhs.add(dof_indices, u0_rhs);
 
            for (const auto &edge : cell->face_indices())
              {
                this->assemble->reinit(cell, edge);
                std::vector<lifex::types::global_dof_index> dof_indices_neigh(
                  this->dofs_per_cell);
 
                if (!cell->at_boundary(edge))
                  {
                    SC = this->assemble->local_SC(this->prm_stability_coeff);
                    std::tie(IC, IC_t) =
                      this->assemble->local_IC(this->prm_penalty_coeff);
                    this->matrix.add(dof_indices, SC);
                    this->matrix.add(dof_indices, IC);
                    this->matrix.add(dof_indices, IC_t);
 
                    const auto neighcell = cell->neighbor(edge);
                    dof_indices_neigh =
                      this->dof_handler.get_dof_indices(neighcell);
 
                    std::tie(IN, IN_t) =
                      this->assemble->local_IN(this->prm_penalty_coeff);
                    SN = this->assemble->local_SN(this->prm_stability_coeff);
 
                    this->matrix.add(dof_indices, dof_indices_neigh, IN);
                    this->matrix.add(dof_indices_neigh, dof_indices, IN_t);
                    this->matrix.add(dof_indices, dof_indices_neigh, SN);
                  }
                else
                  {
                    cell_rhs_edge =
                      this->assemble->local_rhs_edge_neumann(this->g_n);
                    this->rhs.add(dof_indices, cell_rhs_edge);
                  }
              }
          }
      }
 
    this->matrix.compress(dealii::VectorOperation::add);
    this->rhs.compress(dealii::VectorOperation::add);
  }
} // namespace DUBeat::models
 
#endif /* HEAT_DG_HPP_*/