utilities.h

Go to the documentation of this file.

/*
  Copyright (C) 2011 - 2024 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_utilities_h
#define _aspect_material_model_utilities_h

#include <aspect/global.h>
#include <deal.II/base/point.h>
#include <deal.II/base/symmetric_tensor.h>
#include <deal.II/fe/component_mask.h>
#include <deal.II/base/signaling_nan.h>
#include <deal.II/base/parameter_handler.h>

#include <mpi.h>

#include <limits>


namespace aspect
{
  template <int dim> class SimulatorAccess;
  namespace Utilities
  {
    using namespace ::Utilities;

    template <int dim>
    class StructuredDataLookup;
  }
  namespace MaterialModel
  {
    template <int dim> class MaterialModelOutputs;
    template <int dim> struct EquationOfStateOutputs;

    namespace MaterialUtilities
    {
      namespace Lookup
      {
        class MaterialLookup
        {
          public:

            double
            specific_heat(const double temperature,
                          const double pressure) const;

            double
            density(const double temperature,
                    const double pressure) const;

            double
            thermal_expansivity(const double temperature,
                                const double pressure) const;

            double
            seismic_Vp(const double temperature,
                       const double pressure) const;

            double
            seismic_Vs(const double temperature,
                       const double pressure) const;

            double
            enthalpy(const double temperature,
                     const double pressure) const;

            double
            dHdT (const double temperature,
                  const double pressure) const;

            double
            dHdp (const double temperature,
                  const double pressure) const;

            std::array<std::pair<double, unsigned int>,2>
            enthalpy_derivatives(const std::vector<double> &temperatures,
                                 const std::vector<double> &pressures,
                                 const unsigned int n_substeps = 1) const;

            double
            dRhodp (const double temperature,
                    const double pressure) const;

            unsigned int
            dominant_phase (const double temperature,
                            const double pressure) const;

            bool
            has_dominant_phase() const;

            std::vector<std::string>
            phase_volume_column_names() const;

            double
            phase_volume_fraction(const int phase_id,
                                  const double temperature,
                                  const double pressure) const;

            std::array<double,2>
            get_pT_steps() const;


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

          protected:
            double
            value (const double temperature,
                   const double pressure,
                   const Table<2, double> &values,
                   const bool interpol) const;

            unsigned int
            value (const double temperature,
                   const double pressure,
                   const Table<2, unsigned int> &values) const;

            double get_nT(const double temperature) const;

            double get_np(const double pressure) const;

            ::Table<2,double> density_values;
            ::Table<2,double> thermal_expansivity_values;
            ::Table<2,double> specific_heat_values;
            ::Table<2,double> vp_values;
            ::Table<2,double> vs_values;
            ::Table<2,double> enthalpy_values;
            ::Table<2,unsigned int> dominant_phase_indices;

            std::vector<std::string> phase_column_names;
            std::vector<::Table<2,double>> phase_volume_fractions;

            double delta_press;
            double min_press;
            double max_press;
            double delta_temp;
            double min_temp;
            double max_temp;
            unsigned int n_temperature;
            unsigned int n_pressure;
            unsigned int n_phases;
            unsigned int n_columns;
            bool interpolation;
            bool has_dominant_phase_column;
            std::vector<std::string> dominant_phase_names;
        };

        class HeFESToReader : public MaterialLookup
        {
          public:
            HeFESToReader(const std::string &material_filename,
                          const std::string &derivatives_filename,
                          const bool interpol,
                          const MPI_Comm comm);
        };

        class PerplexReader : public MaterialLookup
        {
          public:
            PerplexReader(const std::string &filename,
                          const bool interpol,
                          const MPI_Comm comm);
        };

        class EntropyReader
        {
          public:

            void
            initialize(const MPI_Comm comm,
                       const std::string &data_directory,
                       const std::string &material_file_name);

            double
            specific_heat(const double entropy,
                          const double pressure) const;

            double
            density(const double entropy,
                    const double pressure) const;

            double
            thermal_expansivity(const double entropy,
                                const double pressure) const;

            double
            temperature(const double entropy,
                        const double pressure) const;

            double
            seismic_vp(const double entropy,
                       const double pressure) const;

            double
            seismic_vs(const double entropy,
                       const double pressure) const;

            Tensor<1, 2>
            density_gradient(const double entropy,
                             const double pressure) const;

          private:
            std::unique_ptr<Utilities::StructuredDataLookup<2>> material_lookup;
        };
      }

      std::vector<double>
      compute_only_composition_fractions(const std::vector<double> &compositional_fields,
                                         const std::vector<unsigned int> &indices_to_use,
                                         const double minimum_fraction = -std::numeric_limits<double>::max());

      std::vector<double>
      compute_composition_fractions(const std::vector<double> &compositional_fields,
                                    const ComponentMask &field_mask = ComponentMask(),
                                    const double minimum_fraction = -std::numeric_limits<double>::max());

      std::vector<double>
      compute_volumes_from_masses(const std::vector<double> &masses,
                                  const std::vector<double> &densities,
                                  const bool return_as_fraction);

      enum CompositionalAveragingOperation
      {
        harmonic,
        arithmetic,
        geometric,
        maximum_composition
      };



      CompositionalAveragingOperation
      parse_compositional_averaging_operation (const std::string &parameter_name,
                                               const ParameterHandler &prm);



      double average_value (const std::vector<double> &volume_fractions,
                            const std::vector<double> &parameter_values,
                            const CompositionalAveragingOperation &average_type);



      template <int dim>
      void
      fill_averaged_equation_of_state_outputs(const EquationOfStateOutputs<dim> &eos_outputs,
                                              const std::vector<double> &mass_fractions,
                                              const std::vector<double> &volume_fractions,
                                              const unsigned int i,
                                              MaterialModelOutputs<dim> &out);



      namespace PhaseUtilities
      {
        enum PhaseAveragingOperation
        {
          arithmetic,
          logarithmic
        };
      }



      void
      reaction_progress_modify_values (const std::vector<double> &reaction_progress_values,
                                       const std::vector<unsigned int> &reaction_progress_mapping,
                                       const std::vector<unsigned int> &n_phase_transitions_per_composition,
                                       std::vector<double> &parameter_values,
                                       const unsigned int composition_index,
                                       const PhaseUtilities::PhaseAveragingOperation operation = PhaseUtilities::arithmetic);



      double phase_average_value (const std::vector<double> &phase_function_values,
                                  const std::vector<unsigned int> &n_phase_transitions_per_composition,
                                  const std::vector<double> &parameter_values,
                                  const unsigned int composition_index,
                                  const PhaseUtilities::PhaseAveragingOperation operation = PhaseUtilities::arithmetic);
      std::pair<bool, double>
      get_reaction_progress_for_phase_transition(const std::vector<double> &reaction_progress_values,
                                                 const std::vector<unsigned> &reaction_progress_mapping,
                                                 const unsigned int phase_transition_index);

      template <int dim>
      struct PhaseFunctionInputs
      {
        PhaseFunctionInputs(const double temperature,
                            const double pressure,
                            const double depth,
                            const double pressure_depth_derivative,
                            const unsigned int phase_transition_index);

        double temperature;
        double pressure;
        double depth;
        double pressure_depth_derivative;

        unsigned int phase_transition_index;
      };

      template <int dim>
      class PhaseFunctionDiscrete: public ::aspect::SimulatorAccess<dim>
      {
        public:

          void initialize();

          double compute_value (const PhaseFunctionInputs<dim> &in) const;

          double compute_derivative () const;

          unsigned int n_phase_transitions () const;

          unsigned int n_phases () const;

          unsigned int n_phases_over_all_chemical_compositions () const;

          const std::vector<unsigned int> &
          n_phase_transitions_for_each_chemical_composition () const;

          const std::vector<unsigned int> &
          n_phases_for_each_chemical_composition () const;

          const std::vector<unsigned int> &
          n_phase_transitions_for_each_composition () const;

          const std::vector<unsigned int> &
          n_phases_for_each_composition () const;

          static
          void
          declare_parameters (ParameterHandler &prm);

          void
          parse_parameters (ParameterHandler &prm);


        private:
          std::string data_directory;

          std::vector<std::string> material_file_names;

          std::vector<double> minimum_temperature;

          std::vector<double> maximum_temperature;

          std::vector<double> interval_temperature;

          std::vector<double> minimum_pressure;

          std::vector<double> maximum_pressure;

          std::vector<double> interval_pressure;

          std::vector<std::unique_ptr<Utilities::StructuredDataLookup<2>>> material_lookup;

          std::vector<unsigned int> transition_indicators;

          std::unique_ptr<std::vector<unsigned int>> n_phase_transitions_per_composition;

          std::vector<unsigned int> n_phases_per_composition;

          std::vector<unsigned int> n_phase_transitions_per_chemical_composition;

          std::vector<unsigned int> n_phases_per_chemical_composition;

          unsigned int n_phases_total;

          unsigned int n_phases_total_chemical_compositions;
      };

      template <int dim>
      class PhaseFunction: public ::aspect::SimulatorAccess<dim>
      {
        public:
          double compute_value (const PhaseFunctionInputs<dim> &in) const;

          double compute_derivative (const PhaseFunctionInputs<dim> &in) const;

          unsigned int n_phase_transitions () const;

          unsigned int n_phases () const;

          unsigned int n_phases_over_all_chemical_compositions () const;

          double get_transition_slope (const unsigned int phase_transition_index) const;

          double get_transition_depth (const unsigned int phase_transition_index) const;

          const std::vector<unsigned int> &
          n_phase_transitions_for_each_chemical_composition () const;

          const std::vector<unsigned int> &
          n_phases_for_each_chemical_composition () const;

          const std::vector<unsigned int> &
          n_phase_transitions_for_each_composition () const;

          const std::vector<unsigned int> &
          n_phases_for_each_composition () const;

          static
          void
          declare_parameters (ParameterHandler &prm);

          void
          parse_parameters (ParameterHandler &prm);


        private:
          std::vector<double> transition_depths;
          std::vector<double> transition_pressures;
          std::vector<double> transition_temperatures;
          std::vector<double> transition_widths;
          std::vector<double> transition_pressure_widths;
          std::vector<double> transition_slopes;
          std::vector<double> transition_temperature_upper_limits;
          std::vector<double> transition_temperature_lower_limits;

          bool use_depth_instead_of_pressure;

          std::unique_ptr<std::vector<unsigned int>> n_phase_transitions_per_composition;

          std::vector<unsigned int> n_phases_per_composition;

          std::vector<unsigned int> n_phase_transitions_per_chemical_composition;

          std::vector<unsigned int> n_phases_per_chemical_composition;

          unsigned int n_phases_total;

          unsigned int n_phases_total_chemical_compositions;
      };
    }
  }
}


#endif