dynamic_core.h

Go to the documentation of this file.

/*
  Copyright (C) 2011 - 2021 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>

namespace aspect
{
  namespace BoundaryTemperature
  {

    namespace internal
    {
      struct 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;
      };
    }


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

        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;

        void update() override;

        const internal::CoreData &
        get_core_data() const;

        bool
        is_OES_used() const;


      private:

        internal::CoreData core_data;

        double inner_temperature;

        double outer_temperature;

        double init_dT_dt;

        double init_dR_dt;

        double init_dX_dt;

        bool is_first_call;

        double Rc;

        double CpRho;

        double X_init;

        double Delta;

        double g;

        double P_CMB;

        double P_Core;

        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;

        int max_steps;

        double dTa;

        std::string name_OES;
        struct str_data_OES
        {
          double t;
          double w;
        };
        std::vector<struct str_data_OES> data_OES;
        void read_data_OES();
        double get_OES(double t) const;



        bool solve_time_step(double &X, double &T, double &R);

        double get_dT(double r) const;

        double get_Tc(double r) const;

        double get_Ts(double r) const;

        double get_solidus(double X,double p) const;

        double get_initial_Ri(double T);

        double get_X(double r) const;

        double get_Mass(double r) const;

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

        double get_Rho(double r) const;

        double get_g(double r) const;

        double get_T(double Tc, double r) const;

        double get_Pressure(double r) const;

        double get_gravity_potential(double r) const;

        void get_specific_heating(double Tc, double &Qs,double &Es);

        void get_radio_heating(double Tc, double &Qr, double &Er);

        void get_gravity_heating(double Tc, double r,double X,double &Qg,double &Eg);

        void get_adiabatic_heating(double Tc, double &Ek, double &Qk);

        void get_latent_heating(double Tc, double r, double &El, double &Ql);

        void get_heat_solution(double Tc, double r, double X, double &Eh);

        double get_radioheating_rate() const;

        void update_core_data();

    };
  }
}


#endif