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> points_per_direction;
 
        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::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