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

#include <aspect/global.h>
#include <aspect/initial_temperature/plate_cooling.h>
#include <aspect/adiabatic_conditions/interface.h>
#include <aspect/boundary_temperature/interface.h>
#include <aspect/geometry_model/box.h>
#include <aspect/geometry_model/two_merged_boxes.h>
#include <aspect/geometry_model/spherical_shell.h>
#include <aspect/geometry_model/chunk.h>
#include <deal.II/base/parameter_handler.h>

#include <cmath>

namespace aspect
{
  namespace InitialTemperature
  {
    template <int dim>
    PlateCooling<dim>::AsciiData ()
    {}


    template <int dim>
    void
    PlateCooling<dim>::initialize ()
    {
      const std::set<types::boundary_id> boundary_ids = this->get_fixed_temperature_boundary_indicators();
      Utilities::AsciiDataBoundary<dim>::initialize(boundary_ids, 1);
    }


    template <int dim>
    void
    PlateCooling<dim>::update ()
    {
      Interface<dim>::update ();
      Utilities::AsciiDataBoundary<dim>::update();
    }


    template <int dim>
    double
    PlateCooling<dim>::initial_temperature (const types::boundary_id boundary_indicator,
                                     const Point<dim> &position) const
    {
      double age_top = Utilities::AsciiDataBoundary<dim>::get_data_component(boundary_indicator,
                                                                             position,
                                                                             0);
    //  return age_top;
    //}

    //template <int dim>
    //double
    //PlateCooling<dim>::
    //initial_temperature (const Point<dim> &position) const
    //{ 
      // convert input ages to seconds
      //const double age_top =    (this->convert_output_to_years() ? age_top_boundary_layer * year_in_seconds
      //                           : age_top_boundary_layer);
      
      // then, get the temperature at the top and bottom boundary of the model
      // if no boundary temperature is prescribed simply use the adiabatic.
      // This implementation assumes that the top and bottom boundaries have
      // prescribed temperatures and minimal_temperature() returns the value
      // at the surface and maximal_temperature() the value at the bottom.
      const double T_surface = this->get_boundary_temperature_manager().minimal_temperature(this->get_fixed_temperature_boundary_indicators());
      const double T_bottom = this->get_boundary_temperature_manager().maximal_temperature(this->get_fixed_temperature_boundary_indicators());
      
      // get a representative profile of the compositional fields as an input
      // for the material model
      const double depth = this->get_geometry_model().depth(position);

      // look up material properties
      MaterialModel::MaterialModelInputs<dim> in(1, this->n_compositional_fields());
      MaterialModel::MaterialModelOutputs<dim> out(1, this->n_compositional_fields());
      in.position[0]=position;
      in.temperature[0]=this->get_adiabatic_conditions().temperature(position);
      in.temperature[0]=this->get_adiabatic_conditions().temperature(position);
      in.pressure[0]=this->get_adiabatic_conditions().pressure(position);
      in.velocity[0]= Tensor<1,dim> ();
      for (unsigned int c=0; c<this->n_compositional_fields(); ++c)
        in.composition[0][c] = function->value(Point<1>(depth),c);
      in.strain_rate.resize(0); // adiabat has strain=0.
      this->get_material_model().evaluate(in, out);

      const double kappa = ( (this->get_parameters().formulation_temperature_equation ==
                              Parameters<dim>::Formulation::TemperatureEquation::reference_density_profile)
                             ?
                             out.thermal_conductivities[0] /
                             (this->get_adiabatic_conditions().density(in.position[0]) * out.specific_heat[0])
                             :
                             out.thermal_conductivities[0] / (out.densities[0] * out.specific_heat[0])
                           );

      // analytical solution for the thermal boundary layer from half-space cooling model
      if (this->get_geometry_model().depth(position) > zlo)
        {
          const double temperature_profile = T_bottom;
          return temperature_profile;
        }
      else
        {
          const double pi = 3.1415;
          const double exp_fac = -kappa * std::pow(pi, 2) * age_top / std::pow(zlo, 2);
          const double sin_fac = pi / zlo;
          unsigned int n = 1;
          double sum_terms = 0;
          while (n <= m)
            {
              sum_terms += 1/(double)n * std::exp(std::pow((double)n, 2) * exp_fac) * std::sin((double)n * (this->get_geometry_model().depth(position) * sin_fac));
              n += 1;
            }
          const double temperature_profile = T_surface + (T_bottom - T_surface) * ((this->get_geometry_model().depth(position)) / zlo + 2/pi * sum_terms);
          return temperature_profile;
        }
    }
    template <int dim>
    void
    PlateCooling<dim>::declare_parameters (ParameterHandler &prm)
    {
      prm.enter_subsection ("Initial temperature model");
      {
        Utilities::AsciiDataBoundary<dim>::declare_parameters(prm,
                                                              "$ASPECT_SOURCE_DIR/data/boundary-temperature/ascii-data/test/",
                                                              "box_2d_%s.%d.txt");
        prm.enter_subsection("Plate cooling");
        {
          prm.declare_entry ("Age top boundary layer", "0.",
                             Patterns::Double (0.),
                             "The age of the upper thermal boundary layer, used for the calculation "
                             "of the half-space cooling model temperature. Units: years if the "
                             "'Use years in output instead of seconds' parameter is set; "
                             "seconds otherwise.");
          prm.declare_entry ("Depth to plate base", "125.e3",
                             Patterns::Anything (),
                             "The depth to the base of the plate. Units: m");
          prm.declare_entry ("Number of terms in sum", "5", 
                             Patterns::Integer (1),
                             "The number of terms to include in the summation");
          //Utilities::AsciiDataBoundary<dim>::declare_parameters(prm,
          //                                                      "$ASPECT_SOURCE_DIR/data/boundary-temperature/ascii-data/test/",
          //                                                      "box_2d_%s.%d.txt");
                             
          prm.enter_subsection("Function");
          {
            Functions::ParsedFunction<1>::declare_parameters (prm, 1);
          }
          prm.leave_subsection();
        }
        prm.leave_subsection ();
      }
      prm.leave_subsection ();
    }


    template <int dim>
    void
    PlateCooling<dim>::parse_parameters (ParameterHandler &prm)
    {
      // initialize the function parser. unfortunately, we can't get it
      // via SimulatorAccess from the simulator itself because at the
      // current point the SimulatorAccess hasn't been initialized
      // yet. so get it from the parameter file directly.
      prm.enter_subsection ("Compositional fields");
      const unsigned int n_compositional_fields = prm.get_integer ("Number of fields");
      prm.leave_subsection ();

      prm.enter_subsection ("Initial temperature model");
      {
        Utilities::AsciiDataBoundary<dim>::parse_parameters(prm);
        prm.enter_subsection("Plate cooling");
        {
          age_top_boundary_layer = prm.get_double ("Age top boundary layer");
          zlo = prm.get_double ("Depth to plate base");
          m = prm.get_integer ("Number of terms in sum");
          //Utilities::AsciiDataBoundary<dim>::parse_parameters(prm);
          if (n_compositional_fields > 0)
            {
              prm.enter_subsection("Function");
              try
                {
                  function
                    = std_cxx14::make_unique<Functions::ParsedFunction<1>>(n_compositional_fields);
                  function->parse_parameters (prm);
                }
              catch (...)
                {
                  std::cerr << "ERROR: FunctionParser failed to parse\n"
                            << "\t'Initial temperature model.Adiabatic.Function'\n"
                            << "with expression\n"
                            << "\t'" << prm.get("Function expression") << "'"
                            << "More information about the cause of the parse error \n"
                            << "is shown below.\n";
                  throw;
                }

              prm.leave_subsection();
            }
        }
        prm.leave_subsection ();
      }
      prm.leave_subsection ();
    }

  }
}

// explicit instantiations
namespace aspect
{
  namespace InitialTemperature
  {
    ASPECT_REGISTER_INITIAL_TEMPERATURE_MODEL(PlateCooling,
                                              "plate cooling",
                                              "Temperature is prescribed as a plate cooling model "
                                              "profile with upper and lower thermal boundary layers, "
                                              "whose ages are given as input parameters.")
  }
}