interface.h

Go to the documentation of this file.

/*
  Copyright (C) 2017 - 2022 by the authors of the ASPECT code.

  This file is part of ASPECT.

  ASPECT is free software; you can redistribute it and/or modify
  it under the terms of the GNU General Public License as published by
  the Free Software Foundation; either version 2, or (at your option)
  any later version.

  ASPECT 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 General Public License for more details.

  You should have received a copy of the GNU General Public License
  along with ASPECT; see the file LICENSE.  If not see
  <http://www.gnu.org/licenses/>.
*/

#ifndef _aspect_simulator_assemblers_interface_h
#define _aspect_simulator_assemblers_interface_h

#include <aspect/global.h>
#include <aspect/heating_model/interface.h>
#include <aspect/material_model/interface.h>

#include <deal.II/fe/fe_values.h>

namespace aspect
{
  using namespace dealii;

  template <int dim>
  class Simulator;

  struct AdvectionField;

  namespace internal
  {
    namespace Assembly
    {
      namespace Scratch
      {
        template <int dim>
        struct ScratchBase
        {
          ScratchBase()
            :
            cell(),
            face_number(numbers::invalid_unsigned_int)
          {}

          ScratchBase(const ScratchBase &scratch)
            :
            cell(scratch.cell),
            face_number(scratch.face_number)
          {}

          virtual ~ScratchBase ()  = default;

          typename DoFHandler<dim>::active_cell_iterator cell;

          unsigned face_number;
        };

        template <int dim>
        struct StokesPreconditioner: public ScratchBase<dim>
        {
          StokesPreconditioner (const FiniteElement<dim> &finite_element,
                                const Quadrature<dim>    &quadrature,
                                const Mapping<dim>       &mapping,
                                const UpdateFlags         update_flags,
                                const unsigned int        n_compositional_fields,
                                const unsigned int        stokes_dofs_per_cell,
                                const bool                add_compaction_pressure,
                                const bool                rebuild_matrix);
          StokesPreconditioner (const StokesPreconditioner &scratch);

          ~StokesPreconditioner () override;

          FEValues<dim> finite_element_values;

          void reinit (const typename DoFHandler<dim>::active_cell_iterator &cell_ref);

          std::vector<types::global_dof_index> local_dof_indices;
          std::vector<unsigned int>            dof_component_indices;
          std::vector<SymmetricTensor<2,dim>> grads_phi_u;
          std::vector<double>                  div_phi_u;
          std::vector<double>                  phi_p;
          std::vector<double>                  phi_p_c;
          std::vector<Tensor<1,dim>>          grad_phi_p;

          MaterialModel::MaterialModelInputs<dim> material_model_inputs;
          MaterialModel::MaterialModelOutputs<dim> material_model_outputs;

          const bool rebuild_stokes_matrix;
        };



        template <int dim>
        struct StokesSystem : public StokesPreconditioner<dim>
        {
          StokesSystem (const FiniteElement<dim> &finite_element,
                        const Mapping<dim>       &mapping,
                        const Quadrature<dim>    &quadrature,
                        const Quadrature<dim-1>  &face_quadrature,
                        const UpdateFlags         update_flags,
                        const UpdateFlags         face_update_flags,
                        const unsigned int        n_compositional_fields,
                        const unsigned int        stokes_dofs_per_cell,
                        const bool                add_compaction_pressure,
                        const bool                use_reference_density_profile,
                        const bool                rebuild_stokes_matrix,
                        const bool                rebuild_newton_stokes_matrix);

          StokesSystem (const StokesSystem<dim> &scratch);

          FEFaceValues<dim> face_finite_element_values;

          using StokesPreconditioner<dim>::reinit;

          void reinit (const typename DoFHandler<dim>::active_cell_iterator &cell_ref,
                       const unsigned face_number_ref);

          std::vector<Tensor<1,dim>>          phi_u;
          std::vector<Tensor<1,dim>>          velocity_values;
          std::vector<double>                  velocity_divergence;
          std::vector<Tensor<1,dim>>          temperature_gradients;

          MaterialModel::MaterialModelInputs<dim> face_material_model_inputs;
          MaterialModel::MaterialModelOutputs<dim> face_material_model_outputs;

          std::vector<double> reference_densities;
          std::vector<double> reference_densities_depth_derivative;

          const bool rebuild_newton_stokes_matrix;
        };



        template <int dim>
        struct AdvectionSystem: public ScratchBase<dim>
        {
          AdvectionSystem (const FiniteElement<dim> &finite_element,
                           const FiniteElement<dim> &advection_element,
                           const Mapping<dim>       &mapping,
                           const Quadrature<dim>    &quadrature,
                           const Quadrature<dim-1>  &face_quadrature,
                           const UpdateFlags         update_flags,
                           const UpdateFlags         face_update_flags,
                           const unsigned int        n_compositional_fields,
                           const typename Simulator<dim>::AdvectionField     &field);
          AdvectionSystem (const AdvectionSystem &scratch);

          FEValues<dim> finite_element_values;

          void reinit (const typename DoFHandler<dim>::active_cell_iterator &cell_ref);

          std::unique_ptr<FEFaceValues<dim>>    face_finite_element_values;
          std::unique_ptr<FEFaceValues<dim>>    neighbor_face_finite_element_values;
          std::unique_ptr<FESubfaceValues<dim>> subface_finite_element_values;

          std::vector<types::global_dof_index>   local_dof_indices;

          std::vector<double>         phi_field;
          std::vector<Tensor<1,dim>> grad_phi_field;
          std::vector<double>         laplacian_phi_field;
          std::vector<double>         face_phi_field;
          std::vector<Tensor<1,dim>> face_grad_phi_field;
          std::vector<double>         neighbor_face_phi_field;
          std::vector<Tensor<1,dim>> neighbor_face_grad_phi_field;

          std::vector<Tensor<1,dim>> old_velocity_values;
          std::vector<Tensor<1,dim>> old_old_velocity_values;

          std::vector<double>         old_pressure;
          std::vector<double>         old_old_pressure;
          std::vector<Tensor<1,dim>> old_pressure_gradients;
          std::vector<Tensor<1,dim>> old_old_pressure_gradients;

          std::vector<SymmetricTensor<2,dim>> old_strain_rates;
          std::vector<SymmetricTensor<2,dim>> old_old_strain_rates;

          std::vector<double>         old_temperature_values;
          std::vector<double>         old_old_temperature_values;

          std::vector<double>         old_field_values;
          std::vector<double>         old_old_field_values;
          std::vector<Tensor<1,dim>> old_field_grads;
          std::vector<Tensor<1,dim>> old_old_field_grads;
          std::vector<double>         old_field_laplacians;
          std::vector<double>         old_old_field_laplacians;

          std::vector<std::vector<double>> old_composition_values;
          std::vector<std::vector<double>> old_old_composition_values;

          std::vector<double>         current_temperature_values;
          std::vector<Tensor<1,dim>> current_velocity_values;
          std::vector<Tensor<1,dim>> face_current_velocity_values;
          std::vector<Tensor<1,dim>> mesh_velocity_values;
          std::vector<Tensor<1,dim>> face_mesh_velocity_values;

          std::vector<SymmetricTensor<2,dim>> current_strain_rates;
          std::vector<std::vector<double>> current_composition_values;
          std::vector<double>         current_velocity_divergences;

          MaterialModel::MaterialModelInputs<dim> material_model_inputs;
          MaterialModel::MaterialModelOutputs<dim> material_model_outputs;

          MaterialModel::MaterialModelInputs<dim> face_material_model_inputs;
          MaterialModel::MaterialModelOutputs<dim> face_material_model_outputs;

          MaterialModel::MaterialModelInputs<dim> neighbor_face_material_model_inputs;
          MaterialModel::MaterialModelOutputs<dim> neighbor_face_material_model_outputs;

          HeatingModel::HeatingModelOutputs heating_model_outputs;
          HeatingModel::HeatingModelOutputs face_heating_model_outputs;
          HeatingModel::HeatingModelOutputs neighbor_face_heating_model_outputs;

          const typename Simulator<dim>::AdvectionField *advection_field;

          double artificial_viscosity;
        };
      }



      namespace CopyData
      {
        template <int dim>
        struct CopyDataBase
        {
          virtual ~CopyDataBase () = default;
        };

        template <int dim>
        struct StokesPreconditioner: public CopyDataBase<dim>
        {
          StokesPreconditioner (const unsigned int stokes_dofs_per_cell);

          StokesPreconditioner (const StokesPreconditioner &data);

          ~StokesPreconditioner () override = default;
          StokesPreconditioner<dim> &operator= (const StokesPreconditioner<dim> &data) = default;

          FullMatrix<double> local_matrix;
          std::vector<types::global_dof_index> local_dof_indices;

          void extract_stokes_dof_indices(const std::vector<types::global_dof_index> &all_dof_indices,
                                          const Introspection<dim>                   &introspection,
                                          const FiniteElement<dim>                   &finite_element);
        };

        template <int dim>
        struct StokesSystem : public StokesPreconditioner<dim>
        {
          StokesSystem (const unsigned int        stokes_dofs_per_cell,
                        const bool                do_pressure_rhs_compatibility_modification);
          StokesSystem (const StokesSystem<dim> &data);

          ~StokesSystem () override = default;
          StokesSystem<dim> &operator= (const StokesSystem<dim> &data) = default;

          Vector<double> local_rhs;
          Vector<double> local_pressure_shape_function_integrals;
        };

        template <int dim>
        struct AdvectionSystem: public CopyDataBase<dim>
        {
          AdvectionSystem (const FiniteElement<dim> &finite_element,
                           const bool                field_is_discontinuous);

          FullMatrix<double>          local_matrix;

          std::vector<FullMatrix<double>>         local_matrices_int_ext;
          std::vector<FullMatrix<double>>         local_matrices_ext_int;
          std::vector<FullMatrix<double>>         local_matrices_ext_ext;

          Vector<double>              local_rhs;

          std::vector<bool>               assembled_matrices;

          std::vector<types::global_dof_index>   local_dof_indices;

          std::vector<std::vector<types::global_dof_index>>   neighbor_dof_indices;
        };
      }
    }
  }

  namespace Assemblers
  {

    unsigned int
    n_interface_matrices (const ReferenceCell &reference_cell);

    unsigned int
    nth_interface_matrix (const ReferenceCell &reference_cell,
                          const unsigned int face);

    unsigned int
    nth_interface_matrix (const ReferenceCell &reference_cell,
                          const unsigned int face,
                          const unsigned int sub_face);

    template <int dim>
    class Interface
    {
      public:
        virtual ~Interface () = default;

        virtual
        void
        execute(internal::Assembly::Scratch::ScratchBase<dim> &scratch,
                internal::Assembly::CopyData::CopyDataBase<dim> &data) const = 0;

        virtual
        void
        create_additional_material_model_outputs(MaterialModel::MaterialModelOutputs<dim> &) const;

        virtual
        std::vector<double>
        compute_residual(internal::Assembly::Scratch::ScratchBase<dim> &) const;
    };



    template <int dim>
    class AdvectionStabilizationInterface
    {
      public:
        virtual ~AdvectionStabilizationInterface ();

        virtual
        std::vector<double>
        advection_prefactors(internal::Assembly::Scratch::ScratchBase<dim> &scratch_base) const;

        virtual
        std::vector<double>
        diffusion_prefactors(internal::Assembly::Scratch::ScratchBase<dim> &scratch_base) const;
    };



    template <int dim>
    class Manager
    {
      public:

        void reset ();

        std::vector<std::unique_ptr<Assemblers::Interface<dim>>> stokes_preconditioner;

        std::vector<std::unique_ptr<Assemblers::Interface<dim>>> stokes_system;

        std::vector<std::unique_ptr<Assemblers::Interface<dim>>> stokes_system_on_boundary_face;

        std::vector<std::vector<std::unique_ptr<Assemblers::Interface<dim>>>> advection_system;

        std::vector<std::vector<std::unique_ptr<Assemblers::Interface<dim>>>> advection_system_on_boundary_face;

        std::vector<std::vector<std::unique_ptr<Assemblers::Interface<dim>>>> advection_system_on_interior_face;

        struct Properties
        {
          Properties ();

          bool need_face_material_model_data;

          bool need_face_finite_element_evaluation;

          bool need_viscosity;

          UpdateFlags needed_update_flags;
        };

        Properties stokes_preconditioner_assembler_properties;
        Properties stokes_system_assembler_properties;
        Properties stokes_system_assembler_on_boundary_face_properties;
        std::vector<Properties> advection_system_assembler_properties;
        std::vector<Properties> advection_system_assembler_on_face_properties;
    };
  }
}


#endif