newton.h

Go to the documentation of this file.

/*
  Copyright (C) 2016 - 2022 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_newton_h
#define _aspect_newton_h

#include <aspect/simulator/assemblers/interface.h>
#include <aspect/simulator_access.h>
#include <aspect/global.h>
#include <aspect/material_model/interface.h>

#include <deal.II/base/signaling_nan.h>

namespace aspect
{
  using namespace dealii;

  namespace MaterialModel
  {
    template <int dim>
    class MaterialModelDerivatives : public AdditionalMaterialOutputs<dim>
    {
      public:
        MaterialModelDerivatives (const unsigned int n_points);

        std::vector<double> viscosity_derivative_wrt_pressure;
        std::vector<SymmetricTensor<2,dim>> viscosity_derivative_wrt_strain_rate;

    };
  }


  namespace Newton
  {
    struct Parameters
    {

      enum Stabilization
      {
        none = 0,
        symmetric = 1,
        PD = 2,
        SPD = symmetric | PD
      };


      friend
      Stabilization
      operator| (const Stabilization a,
                 const Stabilization b)
      {
        return static_cast<Stabilization>(
                 static_cast<int>(a) | static_cast<int>(b));
      }


      friend
      Stabilization
      operator& (const Stabilization a,
                 const Stabilization b)
      {
        return static_cast<Stabilization>(
                 static_cast<int>(a) & static_cast<int>(b));
      }


      static void declare_parameters (ParameterHandler &prm);

      void parse_parameters (ParameterHandler &prm);

      double              newton_derivative_scaling_factor;

      Stabilization       preconditioner_stabilization;
      Stabilization       velocity_block_stabilization;

      bool                use_Newton_failsafe;

      double              nonlinear_switch_tolerance;

      bool                use_Eisenstat_Walker_method_for_Picard_iterations;
      unsigned int        max_pre_newton_nonlinear_iterations;
      unsigned int        max_newton_line_search_iterations;
      bool                use_newton_residual_scaling_method;
      double              maximum_linear_stokes_solver_tolerance;
      double              SPD_safety_factor;
    };


    std::string
    to_string(const Newton::Parameters::Stabilization preconditioner_stabilization);
  }



  template <int dim>
  class NewtonHandler : public SimulatorAccess<dim>
  {
    public:
      void set_assemblers (Assemblers::Manager<dim> &assemblers) const;

      static void create_material_model_outputs(MaterialModel::MaterialModelOutputs<dim> &output);

      Newton::Parameters parameters;
  };


  namespace Assemblers
  {
    template <int dim>
    class NewtonInterface : public aspect::Assemblers::Interface<dim>,
      public SimulatorAccess<dim>
    {
      public:
        void
        create_additional_material_model_outputs(MaterialModel::MaterialModelOutputs<dim> &outputs) const override;
    };

    template <int dim>
    class NewtonStokesPreconditioner : public NewtonInterface<dim>
    {
      public:
        void
        execute (internal::Assembly::Scratch::ScratchBase<dim>  &scratch_base,
                 internal::Assembly::CopyData::CopyDataBase<dim> &data_base) const override;
    };

    template <int dim>
    class NewtonStokesIncompressibleTerms : public NewtonInterface<dim>
    {
      public:
        void
        execute (internal::Assembly::Scratch::ScratchBase<dim>  &scratch_base,
                 internal::Assembly::CopyData::CopyDataBase<dim> &data_base) const override;

        void create_additional_material_model_outputs(MaterialModel::MaterialModelOutputs<dim> &outputs) const override;
    };

    template <int dim>
    class NewtonStokesCompressibleStrainRateViscosityTerm : public Assemblers::Interface<dim>,
      public SimulatorAccess<dim>
    {
      public:
        void
        execute(internal::Assembly::Scratch::ScratchBase<dim>  &scratch_base,
                internal::Assembly::CopyData::CopyDataBase<dim> &data_base) const override;
    };

    template <int dim>
    class NewtonStokesReferenceDensityCompressibilityTerm : public Assemblers::Interface<dim>,
      public SimulatorAccess<dim>
    {
      public:
        void
        execute(internal::Assembly::Scratch::ScratchBase<dim>  &scratch_base,
                internal::Assembly::CopyData::CopyDataBase<dim> &data_base) const override;
    };

    template <int dim>
    class NewtonStokesImplicitReferenceDensityCompressibilityTerm : public Assemblers::Interface<dim>,
      public SimulatorAccess<dim>
    {
      public:
        void
        execute(internal::Assembly::Scratch::ScratchBase<dim>  &scratch_base,
                internal::Assembly::CopyData::CopyDataBase<dim> &data_base) const override;
    };

    template <int dim>
    class NewtonStokesIsentropicCompressionTerm : public Assemblers::Interface<dim>,
      public SimulatorAccess<dim>
    {
      public:
        void
        execute(internal::Assembly::Scratch::ScratchBase<dim>  &scratch_base,
                internal::Assembly::CopyData::CopyDataBase<dim> &data_base) const override;
    };

    template <int dim>
    class NewtonStokesProjectedDensityFieldTerm : public Assemblers::Interface<dim>,
      public SimulatorAccess<dim>
    {
      public:
        void
        execute(internal::Assembly::Scratch::ScratchBase<dim>   &scratch,
                internal::Assembly::CopyData::CopyDataBase<dim> &data) const override;
    };
  }
}

#endif