linear_problem.h

Go to the documentation of this file.

// ---------------------------------------------------------------------
//
// Copyright (C) 2022 by Luca Heltai
//
// This file is part of the FSI-suite platform, based on the deal.II library.
//
// The FSI-suite platform 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 3.0 of the License,
// or (at your option) any later version. The full text of the license can be
// found in the file LICENSE at the top level of the FSI-suite platform
// distribution.
//
// ---------------------------------------------------------------------
 
// Make sure we don't redefine things
#ifndef pdes_linear_problem_h
#define pdes_linear_problem_h
 
#include <deal.II/base/conditional_ostream.h>
#include <deal.II/base/convergence_table.h>
#include <deal.II/base/function.h>
#include <deal.II/base/function_parser.h>
#include <deal.II/base/parameter_acceptor.h>
#include <deal.II/base/quadrature_lib.h>
#include <deal.II/base/timer.h>
 
#include <deal.II/distributed/grid_refinement.h>
#include <deal.II/distributed/shared_tria.h>
#include <deal.II/distributed/tria.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/fe_values.h>
#include <deal.II/fe/fe_values_extractors.h>
#include <deal.II/fe/mapping_q_generic.h>
 
#include <deal.II/grid/grid_generator.h>
#include <deal.II/grid/tria.h>
 
#include <deal.II/lac/affine_constraints.h>
#include <deal.II/lac/dynamic_sparsity_pattern.h>
#include <deal.II/lac/full_matrix.h>
#include <deal.II/lac/generic_linear_algebra.h>
#include <deal.II/lac/precondition.h>
#include <deal.II/lac/solver_cg.h>
#include <deal.II/lac/sparse_matrix.h>
#include <deal.II/lac/vector.h>
 
#include <deal.II/meshworker/copy_data.h>
#include <deal.II/meshworker/scratch_data.h>
 
#include <deal.II/numerics/data_out.h>
#include <deal.II/numerics/matrix_tools.h>
#include <deal.II/numerics/vector_tools.h>
 
#include <deal.II/sundials/arkode.h>
 
#include <boost/signals2.hpp>
 
#include <fstream>
#include <iostream>
 
#include "lac.h"
#include "parsed_lac/amg.h"
#include "parsed_lac/inverse_operator.h"
#include "parsed_tools/boundary_conditions.h"
#include "parsed_tools/constants.h"
#include "parsed_tools/convergence_table.h"
#include "parsed_tools/data_out.h"
#include "parsed_tools/finite_element.h"
#include "parsed_tools/function.h"
#include "parsed_tools/grid_generator.h"
#include "parsed_tools/grid_refinement.h"
 
namespace PDEs
{
  using namespace dealii;
 
  enum class EvolutionType
  {
    steady_state = 1 << 0, //< Steady state problem
    quasi_static = 1 << 1, //< Quasi static problem
    transient    = 1 << 2, //< Transient problem
  };
 
  template <int dim, int spacedim = dim, class LacType = LAC::LAdealii>
  class LinearProblem : public ParameterAcceptor
  {
  public:
    LinearProblem(const std::string &component_names = "u",
                  const std::string &problem_name    = "");
 
    virtual ~LinearProblem() = default;
 
    boost::signals2::signal<void()> check_consistency_call_back;
 
    virtual void
    run();
 
    void
    run_steady_state();
 
    void
    run_quasi_static();
 
    void
    run_transient();
 
    using ARKode = typename SUNDIALS::ARKode<typename LacType::BlockVector>;
 
    virtual void
    setup_transient(ARKode &arkode);
 
    static constexpr bool lac_is_dealii =
      std::is_same<LAC::LAdealii, LacType>::value;
 
    static constexpr bool lac_is_petsc =
      std::is_same<LAC::LAPETSc, LacType>::value;
 
    static constexpr bool lac_is_trilinos =
      std::is_same<LAC::LATrilinos, LacType>::value;
 
    using Triangulation = typename std::conditional<
      dim == 1,
      dealii::parallel::shared::Triangulation<dim, spacedim>,
      dealii::parallel::distributed::Triangulation<dim, spacedim>>::type;
 
    using CopyData = MeshWorker::CopyData<1, 1, 1>;
 
    using ScratchData = MeshWorker::ScratchData<dim, spacedim>;
 
    using BlockVectorType = typename LacType::BlockVector;
 
    using VectorType = typename BlockVectorType::BlockType;
 
    using BlockMatrixType = typename LacType::BlockSparseMatrix;
 
    virtual void
    assemble_system_one_cell(
      const typename DoFHandler<dim, spacedim>::active_cell_iterator &cell,
      ScratchData                                                    &scratch,
      CopyData                                                       &copy);
 
    virtual void
    copy_one_cell(const CopyData &copy);
 
    virtual void
    solve();
 
    virtual void
    estimate(Vector<float> &error_per_cell) const;
 
    void
    mark(const Vector<float> &error_per_cell);
 
    void
    refine();
 
    virtual void
    setup_system();
 
    virtual void
    custom_estimator(dealii::Vector<float> &error_per_cell) const;
 
    boost::signals2::signal<void()> add_constraints_call_back;
 
    boost::signals2::signal<void()> setup_system_call_back;
 
    boost::signals2::signal<void()> output_results_call_back;
 
    virtual void
    assemble_system();
 
    boost::signals2::signal<void()> assemble_system_call_back;
 
    virtual void
    output_results(const unsigned cycle) const;
 
    virtual void
    print_system_info() const;
 
    boost::signals2::signal<void(ParsedTools::DataOut<dim, spacedim> &)>
      add_data_vector;
 
    boost::signals2::signal<void(const double       &time,
                                 const double       &time_step,
                                 const unsigned int &time_step_number)>
      advance_time_call_back;
 
 
    const std::string component_names;
 
    const unsigned int n_components;
 
    const std::string problem_name;
 
    const std::string section_name;
 
    MPI_Comm mpi_communicator;
 
    const unsigned int mpi_rank;
 
    const unsigned int mpi_size;
 
    int number_of_threads = 1;
 
    unsigned int verbosity_level = 4;
 
    dealii::ConditionalOStream pcout;
 
    mutable TimerOutput timer;
 
    EvolutionType evolution_type;
 
    ParsedTools::GridGenerator<dim, spacedim> grid_generator;
 
    ParsedTools::GridRefinement grid_refinement;
 
    Triangulation triangulation;
 
 
    ParsedTools::FiniteElement<dim, spacedim> finite_element;
 
    std::unique_ptr<Mapping<dim, spacedim>> mapping;
 
    Quadrature<dim> cell_quadrature;
 
    Quadrature<dim - 1> face_quadrature;
 
 
    DoFHandler<dim, spacedim> dof_handler;
 
    AffineConstraints<double> constraints;
 
    std::vector<types::global_dof_index> dofs_per_block;
 
    std::vector<IndexSet> locally_owned_dofs;
 
    std::vector<IndexSet> locally_relevant_dofs;
 
    typename LacType::BlockSparsityPattern sparsity;
 
    typename LacType::BlockSparseMatrix matrix;
 
    typename LacType::BlockSparseMatrix mass_matrix;
 
    typename LacType::BlockVector locally_relevant_solution;
 
    typename LacType::BlockVector solution;
 
    typename LacType::BlockVector rhs;
 
    Vector<float> error_per_cell;
 
    ParsedLAC::InverseOperator inverse_operator;
 
    typename LacType::AMG preconditioner;
 
    ParsedLAC::InverseOperator mass_inverse_operator;
 
    typename LacType::AMG mass_preconditioner;
 
    ParsedTools::Function<spacedim> forcing_term;
 
    ParsedTools::Function<spacedim> exact_solution;
 
    ParsedTools::Function<spacedim> initial_value;
 
    ParsedTools::BoundaryConditions<spacedim> boundary_conditions;
 
 
    mutable ParsedTools::ConvergenceTable error_table;
 
    mutable ParsedTools::DataOut<dim, spacedim> data_out;
 
    double start_time = 0.0;
 
    double end_time = 1.0;
 
    double desired_start_step_size = .0625;
 
    unsigned int output_frequency = 1;
 
    ParsedTools::Proxy<
      typename SUNDIALS::ARKode<typename LacType::BlockVector>::AdditionalData>
      ark_ode_data;
 
    boost::signals2::signal<void(ARKode &)> setup_arkode_call_back;
  };
} // namespace PDEs
#endif