structured_data.h

Go to the documentation of this file.

/*
  Copyright (C) 2014 - 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_structured_data_h
#define _aspect_structured_data_h

#include <aspect/global.h>
#include <aspect/simulator_access.h>

#include <array>

namespace aspect
{
  namespace Utilities
  {
    using namespace ::Utilities;

    template <int dim>
    class StructuredDataLookup
    {
      public:
        StructuredDataLookup(const unsigned int n_components,
                             const double scale_factor);

        explicit StructuredDataLookup(const double scale_factor);

        void reinit(const std::vector<std::string> &column_names,
                    std::vector<std::vector<double>> &&coordinate_values,
                    std::vector<Table<dim,double>> &&data_table,
                    const MPI_Comm mpi_communicator = MPI_COMM_SELF,
                    const unsigned int root_process = numbers::invalid_unsigned_int);

        void
        load_ascii(const std::string &filename,
                   const MPI_Comm communicator);

        void
        load_netcdf(const std::string &filename, const std::vector<std::string> &data_column_names = {});


        void
        load_file(const std::string &filename,
                  const MPI_Comm communicator);

        double
        get_data(const Point<dim> &position,
                 const unsigned int component,
                 const bool crash_if_not_in_range = false) const;

        Tensor<1,dim>
        get_gradients(const Point<dim> &position,
                      const unsigned int component);

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

        bool
        has_equidistant_coordinates() const;

        const std::vector<double> &
        get_interpolation_point_coordinates(const unsigned int dimension) const;

        unsigned int
        get_column_index_from_name(const std::string &column_name) const;

        std::string
        get_column_name_from_index(const unsigned int column_index) const;

        double get_maximum_component_value(const unsigned int component) const;

        unsigned int get_number_of_coordinates(const unsigned int dimension) const;

      private:
        unsigned int n_components;

        std::vector<std::string> data_component_names;

        std::vector<std::unique_ptr<Function<dim>>> data;

        std::array<std::vector<double>,dim> coordinate_values;

        std::vector<double> maximum_component_value;

        TableIndices<dim> table_points;

        const double scale_factor;

        bool coordinate_values_are_equidistant;

        TableIndices<dim>
        compute_table_indices(const TableIndices<dim> &sizes, const std::size_t idx) const;

    };

    template <int dim>
    class AsciiDataBase
    {
      public:
        AsciiDataBase();

        static
        void
        declare_parameters (ParameterHandler  &prm,
                            const std::string &default_directory,
                            const std::string &default_filename,
                            const std::string &subsection_name = "Ascii data model");

        void
        parse_parameters (ParameterHandler &prm,
                          const std::string &subsection_name = "Ascii data model");

        std::string data_directory;

        std::string data_file_name;

        double scale_factor;
    };

    template <int dim>
    class AsciiDataBoundary : public Utilities::AsciiDataBase<dim>, public SimulatorAccess<dim>
    {
      public:
        AsciiDataBoundary();

        virtual
        void
        initialize (const std::set<types::boundary_id> &boundary_ids,
                    const unsigned int components);

        void
        update();

        double
        get_data_component (const types::boundary_id             boundary_indicator,
                            const Point<dim>                    &position,
                            const unsigned int                   component) const;

        double
        get_maximum_component_value (const types::boundary_id boundary_indicator,
                                     const unsigned int       component) const;

        Tensor<1,dim-1>
        vector_gradient(const types::boundary_id boundary_indicator,
                        const Point<dim>        &p,
                        const unsigned int       component) const;

        static
        void
        declare_parameters (ParameterHandler  &prm,
                            const std::string &default_directory,
                            const std::string &default_filename,
                            const std::string &subsection_name = "Ascii data model",
                            const bool declare_time_dependent_parameters = true);

        void
        parse_parameters (ParameterHandler &prm,
                          const std::string &subsection_name = "Ascii data model",
                          const bool parse_time_dependent_parameters = true);

      protected:
        int current_file_number;

        int first_data_file_number;

        bool decreasing_file_order;

        double data_file_time_step;

        double time_weight;

        bool time_dependent;

        std::map<types::boundary_id,
            std::unique_ptr<aspect::Utilities::StructuredDataLookup<dim-1>>> lookups;

        std::map<types::boundary_id,
            std::unique_ptr<aspect::Utilities::StructuredDataLookup<dim-1>>> old_lookups;

        void
        update_data (const types::boundary_id boundary_id,
                     const bool reload_both_files);

        void
        end_time_dependence ();

        std::string
        create_filename (const int filenumber,
                         const types::boundary_id boundary_id) const;
    };



    template <int dim>
    class AsciiDataInitial : public Utilities::AsciiDataBase<dim>, public SimulatorAccess<dim>
    {
      public:
        AsciiDataInitial();

        virtual
        void
        initialize (const unsigned int components);


        double
        get_data_component (const Point<dim> &position,
                            const unsigned int component) const;

        static
        void
        declare_parameters (ParameterHandler  &prm,
                            const std::string &default_directory,
                            const std::string &default_filename,
                            const std::string &subsection_name = "Ascii data model");

        void
        parse_parameters (ParameterHandler &prm,
                          const std::string &subsection_name = "Ascii data model");

      protected:
        std::unique_ptr<aspect::Utilities::StructuredDataLookup<dim>> lookup;

        std::unique_ptr<aspect::Utilities::StructuredDataLookup<3>> slice_lookup;

        bool slice_data;

        Tensor<2,3> rotation_matrix;
    };


    template <int dim>
    class AsciiDataLayered : public Utilities::AsciiDataBase<dim>, public SimulatorAccess<dim>
    {
      public:
        AsciiDataLayered();

        virtual
        void
        initialize (const unsigned int components);


        double
        get_data_component (const Point<dim> &position,
                            const unsigned int component) const;


        static
        void
        declare_parameters (ParameterHandler  &prm,
                            const std::string &default_directory,
                            const std::string &default_filename,
                            const std::string &subsection_name = "Ascii data model");

        void
        parse_parameters (ParameterHandler &prm,
                          const std::string &subsection_name = "Ascii data model");

      protected:
        std::vector<std::unique_ptr<aspect::Utilities::StructuredDataLookup<dim-1>>> lookups;

      private:

        std::string data_directory;

        std::vector<std::string> data_file_names;

        unsigned int number_of_layer_boundaries;

        std::string interpolation_scheme;


    };


    template <int dim>
    class AsciiDataProfile : public Utilities::AsciiDataBase<dim>
    {
      public:
        AsciiDataProfile();

        virtual
        void
        initialize (const MPI_Comm communicator);


        double
        get_data_component (const Point<1>                      &position,
                            const unsigned int                   component) const;

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

        const std::vector<double> &
        get_interpolation_point_coordinates() const;

        unsigned int
        get_column_index_from_name(const std::string &column_name) const;

        unsigned int
        maybe_get_column_index_from_name(const std::string &column_name) const;

        std::string
        get_column_name_from_index(const unsigned int column_index) const;
      protected:
        std::unique_ptr<aspect::Utilities::StructuredDataLookup<1>> lookup;
    };
  }
}

#endif