interface.h

Go to the documentation of this file.

/*
  Copyright (C) 2011 - 2023 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_material_model_interface_h
#define _aspect_material_model_interface_h

#include <aspect/global.h>
#include <aspect/plugins.h>
#include <aspect/material_model/utilities.h>

#include <deal.II/base/point.h>
#include <deal.II/base/quadrature.h>
#include <deal.II/base/symmetric_tensor.h>
#include <deal.II/base/parameter_handler.h>
#include <deal.II/dofs/dof_handler.h>
#include <deal.II/dofs/dof_accessor.h>
#include <deal.II/fe/mapping.h>
#include <deal.II/fe/fe_values.h>
#include <deal.II/fe/component_mask.h>
#include <deal.II/numerics/data_postprocessor.h>
#include <deal.II/base/signaling_nan.h>

namespace aspect
{
  template <int dim>
  struct Introspection;


  namespace MaterialModel
  {
    using namespace dealii;

    namespace NonlinearDependence
    {
      enum Dependence
      {
        uninitialized        = 0,

        none                 = 1,
        temperature          = 2,
        pressure             = 4,
        strain_rate          = 8,
        compositional_fields = 16,

        any_variable         = temperature | pressure | strain_rate | compositional_fields
      };


      inline Dependence operator | (const Dependence d1,
                                    const Dependence d2)
      {
        return Dependence(static_cast<int>(d1) | static_cast<int>(d2));
      }

      inline Dependence operator |= (Dependence &d1,
                                     const Dependence d2)
      {
        d1 = (d1 | d2);
        return d1;
      }

      struct ModelDependence
      {
        Dependence viscosity;

        Dependence density;

        Dependence compressibility;

        Dependence specific_heat;

        Dependence thermal_conductivity;

        ModelDependence ();
      };

      bool
      identifies_single_variable(const Dependence dependence);
    }

    namespace MaterialProperties
    {
      enum Property
      {
        uninitialized                  = 0,

        none                           = 1,
        viscosity                      = 2,
        density                        = 4,
        thermal_expansion_coefficient  = 8,
        specific_heat                  = 16,
        thermal_conductivity           = 32,
        compressibility                = 64,
        entropy_derivative_pressure    = 128,
        entropy_derivative_temperature = 256,
        reaction_terms                 = 512,
        reaction_rates                 = 1024,
        additional_outputs             = 2048,

        equation_of_state_properties   = density |
                                         thermal_expansion_coefficient |
                                         specific_heat |
                                         compressibility |
                                         entropy_derivative_pressure |
                                         entropy_derivative_temperature,
        all_properties                 = equation_of_state_properties |
                                         viscosity |
                                         thermal_conductivity |
                                         reaction_terms |
                                         reaction_rates |
                                         additional_outputs
      };

      inline Property operator | (const Property d1,
                                  const Property d2)
      {
        return Property(static_cast<int>(d1) | static_cast<int>(d2));
      }

      inline Property operator |= (Property &d1,
                                   const Property d2)
      {
        return (d1 | d2);
      }
    }

    template <int dim> class AdditionalMaterialInputs;

    template <int dim>
    struct MaterialModelInputs
    {
      MaterialModelInputs(const unsigned int n_points,
                          const unsigned int n_comp);

      MaterialModelInputs(const DataPostprocessorInputs::Vector<dim> &input_data,
                          const Introspection<dim> &introspection,
                          const bool compute_strain_rate = true);


      MaterialModelInputs(const FEValuesBase<dim,dim> &fe_values,
                          const typename DoFHandler<dim>::active_cell_iterator &cell,
                          const Introspection<dim> &introspection,
                          const LinearAlgebra::BlockVector &solution_vector,
                          const bool compute_strain_rate = true);

      MaterialModelInputs (const MaterialModelInputs &source);

      MaterialModelInputs (MaterialModelInputs &&)  noexcept = default;

      MaterialModelInputs &operator= (const MaterialModelInputs &source) = delete;

      MaterialModelInputs &operator= (MaterialModelInputs &&) = default;

      void reinit(const FEValuesBase<dim,dim>                          &fe_values,
                  const typename DoFHandler<dim>::active_cell_iterator &cell,
                  const Introspection<dim>                             &introspection,
                  const LinearAlgebra::BlockVector                     &solution_vector,
                  const bool                                            compute_strain_rate = true);

      unsigned int n_evaluation_points() const;

      bool requests_property(const MaterialProperties::Property &property) const;

      std::vector<Point<dim>> position;

      std::vector<double> temperature;

      std::vector<double> pressure;

      std::vector<Tensor<1,dim>> pressure_gradient;

      std::vector<Tensor<1,dim>> velocity;

      std::vector<std::vector<double>> composition;

      std::vector<SymmetricTensor<2,dim>> strain_rate;

      typename DoFHandler<dim>::active_cell_iterator current_cell;

      MaterialProperties::Property requested_properties;

      std::vector<std::unique_ptr<AdditionalMaterialInputs<dim>>> additional_inputs;

      template <class AdditionalInputType>
      AdditionalInputType *get_additional_input();

      template <class AdditionalInputType>
      const AdditionalInputType *get_additional_input() const;
    };


    template <int dim>     class AdditionalMaterialOutputs;


    template <int dim>
    struct MaterialModelOutputs
    {
      MaterialModelOutputs (const unsigned int n_points,
                            const unsigned int n_comp);

      MaterialModelOutputs (const MaterialModelOutputs &source);

      MaterialModelOutputs (MaterialModelOutputs &&)  noexcept = default;

      MaterialModelOutputs &operator= (const MaterialModelOutputs &source) = delete;

      MaterialModelOutputs &operator= (MaterialModelOutputs &&)  noexcept = default;

      unsigned int n_evaluation_points() const;

      std::vector<double> viscosities;

      std::vector<double> densities;

      std::vector<double> thermal_expansion_coefficients;

      std::vector<double> specific_heat;

      std::vector<double> thermal_conductivities;

      std::vector<double> compressibilities;

      std::vector<double> entropy_derivative_pressure;

      std::vector<double> entropy_derivative_temperature;

      std::vector<std::vector<double>> reaction_terms;

      std::vector<std::unique_ptr<AdditionalMaterialOutputs<dim>>> additional_outputs;

      template <class AdditionalOutputType>
      AdditionalOutputType *get_additional_output();

      template <class AdditionalOutputType>
      const AdditionalOutputType *get_additional_output() const;

      void move_additional_outputs_from(MaterialModelOutputs<dim> &other);
    };



    namespace MaterialAveraging
    {
      enum AveragingOperation
      {
        none,
        arithmetic_average,
        harmonic_average,
        geometric_average,
        pick_largest,
        project_to_Q1,
        log_average,
        harmonic_average_only_viscosity,
        geometric_average_only_viscosity,
        project_to_Q1_only_viscosity
      };


      std::string get_averaging_operation_names ();

      AveragingOperation parse_averaging_operation_name (const std::string &s);

      template <int dim>
      void average (const AveragingOperation operation,
                    const typename DoFHandler<dim>::active_cell_iterator &cell,
                    const Quadrature<dim>         &quadrature_formula,
                    const Mapping<dim>            &mapping,
                    MaterialModelOutputs<dim>     &values_out);

      void average_property (const AveragingOperation  operation,
                             const FullMatrix<double>      &projection_matrix,
                             const FullMatrix<double>      &expansion_matrix,
                             std::vector<double>           &values_out);
    }


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

        virtual void
        fill (const LinearAlgebra::BlockVector &solution,
              const FEValuesBase<dim>          &fe_values,
              const Introspection<dim>         &introspection) = 0;
    };


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

        virtual void average (const MaterialAveraging::AveragingOperation /*operation*/,
                              const FullMatrix<double>  &/*projection_matrix*/,
                              const FullMatrix<double>  &/*expansion_matrix*/)
        {}
    };


    template <int dim>
    class NamedAdditionalMaterialOutputs : public AdditionalMaterialOutputs<dim>
    {
      public:
        NamedAdditionalMaterialOutputs(const std::vector<std::string> &output_names);

        NamedAdditionalMaterialOutputs(const std::vector<std::string> &output_names,
                                       const unsigned int n_points);

        ~NamedAdditionalMaterialOutputs() override;

        const std::vector<std::string> &get_names() const;

        virtual std::vector<double> get_nth_output(const unsigned int idx) const;

        void average (const MaterialAveraging::AveragingOperation /*operation*/,
                      const FullMatrix<double>  &/*projection_matrix*/,
                      const FullMatrix<double>  &/*expansion_matrix*/) override
        {}


        std::vector<std::vector<double>> output_values;

      private:
        const std::vector<std::string> names;
    };


    template <int dim>
    class SeismicAdditionalOutputs : public NamedAdditionalMaterialOutputs<dim>
    {
      public:
        SeismicAdditionalOutputs(const unsigned int n_points);

        std::vector<double> get_nth_output(const unsigned int idx) const override;

        std::vector<double> vs;

        std::vector<double> vp;
    };


    template <int dim>
    class ReactionRateOutputs : public NamedAdditionalMaterialOutputs<dim>
    {
      public:
        ReactionRateOutputs (const unsigned int n_points,
                             const unsigned int n_comp);

        std::vector<double> get_nth_output(const unsigned int idx) const override;

        std::vector<std::vector<double>> reaction_rates;
    };



    template <int dim>
    class PhaseOutputs : public NamedAdditionalMaterialOutputs<dim>
    {
      public:
        PhaseOutputs(const unsigned int n_points);
    };



    template <int dim>
    class PrescribedFieldOutputs : public NamedAdditionalMaterialOutputs<dim>
    {
      public:
        PrescribedFieldOutputs (const unsigned int n_points,
                                const unsigned int n_comp);

        std::vector<double> get_nth_output(const unsigned int idx) const override;

        std::vector<std::vector<double>> prescribed_field_outputs;
    };

    template <int dim>
    class PrescribedTemperatureOutputs : public NamedAdditionalMaterialOutputs<dim>
    {
      public:
        PrescribedTemperatureOutputs (const unsigned int n_points);

        std::vector<double> get_nth_output(const unsigned int idx) const override;

        std::vector<double> prescribed_temperature_outputs;
    };

    template <int dim>
    class AdditionalMaterialOutputsStokesRHS: public AdditionalMaterialOutputs<dim>
    {
      public:
        AdditionalMaterialOutputsStokesRHS(const unsigned int n_points)
          : rhs_u(n_points), rhs_p(n_points), rhs_melt_pc(n_points)
        {}

        ~AdditionalMaterialOutputsStokesRHS() override
          = default;

        void average (const MaterialAveraging::AveragingOperation /*operation*/,
                      const FullMatrix<double>  &/*projection_matrix*/,
                      const FullMatrix<double>  &/*expansion_matrix*/) override
        {
          // TODO: not implemented
        }

        std::vector<Tensor<1,dim>> rhs_u;

        std::vector<double> rhs_p;

        std::vector<double> rhs_melt_pc;
    };



    template <int dim>
    class PrescribedPlasticDilation : public NamedAdditionalMaterialOutputs<dim>
    {
      public:
        explicit PrescribedPlasticDilation (const unsigned int n_points);

        std::vector<double> get_nth_output(const unsigned int idx) const override;

        std::vector<double> dilation;
    };



    template <int dim>
    class ElasticOutputs: public AdditionalMaterialOutputs<dim>
    {
      public:
        ElasticOutputs(const unsigned int n_points)
          : elastic_force(n_points, numbers::signaling_nan<SymmetricTensor<2,dim>>())
          , viscoelastic_strain_rate(n_points, numbers::signaling_nan<SymmetricTensor<2,dim>>())
        {}

        ~ElasticOutputs() override
          = default;

        void average (const MaterialAveraging::AveragingOperation operation,
                      const FullMatrix<double>  &/*projection_matrix*/,
                      const FullMatrix<double>  &/*expansion_matrix*/) override
        {
          AssertThrow(operation == MaterialAveraging::AveragingOperation::none,ExcNotImplemented());
          return;
        }

        std::vector<SymmetricTensor<2,dim>> elastic_force;

        std::vector<SymmetricTensor<2,dim>> viscoelastic_strain_rate;
    };



    template <int dim>
    class EnthalpyOutputs : public AdditionalMaterialOutputs<dim>
    {
      public:
        EnthalpyOutputs(const unsigned int n_points)
          : enthalpies_of_fusion(n_points, numbers::signaling_nan<double>())
        {}

        virtual ~EnthalpyOutputs()
          = default;

        void average (const MaterialAveraging::AveragingOperation operation,
                      const FullMatrix<double>  &/*projection_matrix*/,
                      const FullMatrix<double>  &/*expansion_matrix*/) override
        {
          AssertThrow(operation == MaterialAveraging::AveragingOperation::none,ExcNotImplemented());
          return;
        }

        std::vector<double> enthalpies_of_fusion;
    };



    template <int dim>
    class Interface
    {
      public:
        using MaterialModelInputs = MaterialModel::MaterialModelInputs<dim>;
        using MaterialModelOutputs = MaterialModel::MaterialModelOutputs<dim>;

        virtual ~Interface() = default;

        virtual
        void
        initialize ();

        virtual void update ();

        const NonlinearDependence::ModelDependence &
        get_model_dependence () const;

        virtual bool is_compressible () const = 0;
        virtual
        void evaluate (const MaterialModel::MaterialModelInputs<dim> &in,
                       MaterialModel::MaterialModelOutputs<dim> &out) const = 0;
        static
        void
        declare_parameters (ParameterHandler &prm);

        virtual
        void
        parse_parameters (ParameterHandler &prm);
        virtual
        void
        create_additional_named_outputs (MaterialModelOutputs &outputs) const;


        virtual
        void
        fill_additional_material_model_inputs(MaterialModel::MaterialModelInputs<dim> &input,
                                              const LinearAlgebra::BlockVector        &solution,
                                              const FEValuesBase<dim>                 &fe_values,
                                              const Introspection<dim>                &introspection) const;

      protected:
        NonlinearDependence::ModelDependence model_dependence;
    };

    template <int dim>
    void
    register_material_model (const std::string &name,
                             const std::string &description,
                             void (*declare_parameters_function) (ParameterHandler &),
                             std::unique_ptr<Interface<dim>> (*factory_function) ());

    template <int dim>
    std::unique_ptr<Interface<dim>>
    create_material_model (const std::string &model_name);


    template <int dim>
    std::unique_ptr<Interface<dim>>
    create_material_model (ParameterHandler &prm);


    template <int dim>
    std::string
    get_valid_model_names_pattern ();


    template <int dim>
    void
    declare_parameters (ParameterHandler &prm);



    template <int dim>
    void
    write_plugin_graph (std::ostream &output_stream);



// --------------------- template function definitions ----------------------------------

    template <int dim>
    template <class AdditionalInputType>
    AdditionalInputType *MaterialModelInputs<dim>::get_additional_input()
    {
      for (unsigned int i=0; i<additional_inputs.size(); ++i)
        {
          AdditionalInputType *result = dynamic_cast<AdditionalInputType *> (additional_inputs[i].get());
          if (result)
            return result;
        }
      return nullptr;
    }


    template <int dim>
    template <class AdditionalInputType>
    const AdditionalInputType *MaterialModelInputs<dim>::get_additional_input() const
    {
      for (unsigned int i=0; i<additional_inputs.size(); ++i)
        {
          const AdditionalInputType *result = dynamic_cast<const AdditionalInputType *> (additional_inputs[i].get());
          if (result)
            return result;
        }
      return nullptr;
    }


    template <int dim>
    template <class AdditionalOutputType>
    AdditionalOutputType *MaterialModelOutputs<dim>::get_additional_output()
    {
      for (unsigned int i=0; i<additional_outputs.size(); ++i)
        {
          AdditionalOutputType *result = dynamic_cast<AdditionalOutputType *> (additional_outputs[i].get());
          if (result)
            return result;
        }
      return nullptr;
    }


    template <int dim>
    template <class AdditionalOutputType>
    const AdditionalOutputType *MaterialModelOutputs<dim>::get_additional_output() const
    {
      for (unsigned int i=0; i<additional_outputs.size(); ++i)
        {
          const AdditionalOutputType *result = dynamic_cast<const AdditionalOutputType *> (additional_outputs[i].get());
          if (result)
            return result;
        }
      return nullptr;
    }


    template <int dim>
    void MaterialModelOutputs<dim>::move_additional_outputs_from(MaterialModelOutputs<dim> &other)
    {
      Assert(this->additional_outputs.empty(), ExcMessage("Destination of move needs to be empty!"));
      this->additional_outputs = std::move(other.additional_outputs);
    }


#define ASPECT_REGISTER_MATERIAL_MODEL(classname,name,description) \
  template class classname<2>; \
  template class classname<3>; \
  namespace ASPECT_REGISTER_MATERIAL_MODEL_ ## classname \
  { \
    aspect::internal::Plugins::RegisterHelper<aspect::MaterialModel::Interface<2>,classname<2>> \
    dummy_ ## classname ## _2d (&aspect::MaterialModel::register_material_model<2>, \
                                name, description); \
    aspect::internal::Plugins::RegisterHelper<aspect::MaterialModel::Interface<3>,classname<3>> \
    dummy_ ## classname ## _3d (&aspect::MaterialModel::register_material_model<3>, \
                                name, description); \
  }
  }
}


#endif