dynamic_core.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 doc/COPYING.  If not see
  <http://www.gnu.org/licenses/>.
*/



#ifndef _aspect_boundary_temperature_dynamic_core_h
#define _aspect_boundary_temperature_dynamic_core_h

#include <aspect/boundary_temperature/interface.h>
#include <aspect/simulator_access.h>

#include <tuple>

namespace aspect
{
  namespace BoundaryTemperature
  {

    namespace internal
    {
      struct CoreData
      {
        CoreData ();

        double Qs,Qr,Qg,Qk,Ql;

        double Es,Er,Eg,Ek,El,Eh;

        double Ri,Ti,Xi;

        double Q,H;

        double dt;

        double dR_dt,dT_dt,dX_dt;

        double Q_OES;

        bool is_initialized;
      };

      using SolveTimeStepResult = std::tuple<double, double, double>;
    }


    template <int dim>
    class DynamicCore : public Interface<dim>, public aspect::SimulatorAccess<dim>
    {
      public:
        DynamicCore();

        void
        update() override;

        const internal::CoreData &
        get_core_data() const;

        bool
        is_OES_used() const;

        double
        boundary_temperature (const types::boundary_id boundary_indicator,
                              const Point<dim> &location) const override;

        double
        minimal_temperature (const std::set<types::boundary_id> &fixed_boundary_ids) const override;

        double
        maximal_temperature (const std::set<types::boundary_id> &fixed_boundary_ids) const override;

        static
        void
        declare_parameters (ParameterHandler &prm);

        void
        parse_parameters (ParameterHandler &prm) override;

        template <class Archive>
        void
        serialize (Archive &ar, const unsigned int version);

        void
        save (std::map<std::string, std::string> &status_strings) const override;

        void
        load (const std::map<std::string, std::string> &status_strings) override;

      private:

        internal::CoreData core_data;

        double inner_temperature;

        double outer_temperature;

        types::boundary_id inner_boundary_id;
        types::boundary_id outer_boundary_id;

        double init_dT_dt;

        double init_dR_dt;

        double init_dX_dt;

        bool is_first_call;

        double Rc;

        double X_init;

        double Delta;

        double P_CMB;

        double Tm0;
        double Tm1;
        double Tm2;
        double Theta;
        bool composition_dependency;

        bool use_bw11;

        //Variables for formulation of \cite NPB+04
        double K0;

        double Alpha;

        double Rho_0;

        double Rho_cen;

        double Lh;

        double Rh;

        double Beta_c;

        double k_c;

        double Cp;

        unsigned int n_radioheating_elements;

        std::vector<double> heating_rate;

        std::vector<double> half_life;

        std::vector<double> initial_concentration;

        double L;
        double D;

        double Mc;

        unsigned int max_steps;

        double dTa;

        std::string name_OES;
        struct str_data_OES
        {
          double t;
          double w;

          template <class Archive>
          void serialize (Archive &ar, const unsigned int)
          {
            ar &t
            & w;
          }
        };
        std::vector<struct str_data_OES> data_OES;
        void read_data_OES();
        double compute_OES(double t) const;

        internal::SolveTimeStepResult solve_time_step() const;

        double compute_dT(const double r) const;

        double compute_Tc(const double r) const;

        double compute_Ts(const double r) const;

        double compute_solidus(const double X, const double pressure) const;

        double compute_initial_Ri(const double T) const;

        double compute_X(const double r) const;

        double compute_mass(const double r) const;

        double fun_Sn(const double B, const double R, const unsigned int n) const;

        double compute_rho(const double r) const;

        double compute_T(const double Tc, const double r) const;

        double compute_pressure(const double r) const;

        double compute_gravity_potential(const double r) const;

        std::pair<double,double>
        compute_specific_heating(const double Tc) const;

        std::pair<double,double>
        compute_radio_heating(const double Tc) const;

        std::pair<double,double>
        compute_gravity_heating(const double Tc, const double r, const double X) const;

        std::pair<double,double>
        compute_adiabatic_heating(const double Tc) const;

        std::pair<double,double>
        compute_latent_heating(const double Tc, const double r) const;

        double
        compute_heat_solution(const double Tc, const double r, const double X) const;

        double compute_radioheating_rate() const;

        void update_core_data();

    };
  }
}


#endif