melt.h

Go to the documentation of this file.

/*
  Copyright (C) 2016 - 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_melt_h
#define _aspect_melt_h

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

namespace aspect
{
  namespace MaterialModel
  {
    template <int dim>
    class MeltInputs : public AdditionalMaterialInputs<dim>
    {
      public:
        MeltInputs (const unsigned int n_points);

        std::vector<double> compaction_pressures;

        std::vector<Tensor<1,dim>> fluid_velocities;

        void fill (const LinearAlgebra::BlockVector &solution,
                   const FEValuesBase<dim>          &fe_values,
                   const Introspection<dim>         &introspection) override;
    };

    template <int dim>
    class MeltOutputs : public AdditionalMaterialOutputs<dim>
    {
      public:
        MeltOutputs (const unsigned int n_points,
                     const unsigned int n_comp);

        std::vector<double> compaction_viscosities;

        std::vector<double> fluid_viscosities;

        std::vector<double> permeabilities;

        std::vector<double> fluid_densities;

        std::vector<Tensor<1,dim>> fluid_density_gradients;

        void average (const MaterialAveraging::AveragingOperation operation,
                      const FullMatrix<double>  &projection_matrix,
                      const FullMatrix<double>  &expansion_matrix) override;
    };

    template <int dim>
    class MeltFractionModel
    {
      public:
        virtual ~MeltFractionModel () = default;

        virtual void melt_fractions (const MaterialModel::MaterialModelInputs<dim> &in,
                                     std::vector<double> &melt_fractions,
                                     const MaterialModel::MaterialModelOutputs<dim> *out = nullptr) const = 0;

        template <typename ModelType>
        static
        bool is_melt_fraction_model (const ModelType &model_object);

        template <typename ModelType>
        static
        const MeltFractionModel<dim> &
        as_melt_fraction_model (const ModelType &model_object);
    };



    template <int dim>
    template <typename ModelType>
    inline
    bool
    MeltFractionModel<dim>::is_melt_fraction_model (const ModelType &model_object)
    {
      return (dynamic_cast<const MeltFractionModel<dim>*>(&model_object)
              != nullptr);
    }


    template <int dim>
    template <typename ModelType>
    inline
    const MeltFractionModel<dim> &
    MeltFractionModel<dim>::as_melt_fraction_model (const ModelType &model_object)
    {
      Assert (is_melt_fraction_model(model_object) == true,
              ExcMessage ("This function can only be called for model objects "
                          "whose types are derived from MeltFractionModel."));

      return dynamic_cast<const MeltFractionModel<dim>&>(model_object);
    }



    template <int dim>
    class MeltInterface: public MaterialModel::Interface<dim>
    {
      public:
        virtual double reference_darcy_coefficient () const = 0;

        double p_c_scale (const MaterialModel::MaterialModelInputs<dim> &inputs,
                          const MaterialModel::MaterialModelOutputs<dim> &outputs,
                          const MeltHandler<dim> &melt_handler,
                          const bool consider_is_melt_cell) const;
    };


  }



  namespace Assemblers
  {
    template <int dim>
    class MeltInterface : 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 MeltStokesPreconditioner : public MeltInterface<dim>
    {
      public:
        void
        execute(internal::Assembly::Scratch::ScratchBase<dim>   &scratch_base,
                internal::Assembly::CopyData::CopyDataBase<dim> &data_base) const override;
    };

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


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

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

        std::vector<double>
        compute_residual(internal::Assembly::Scratch::ScratchBase<dim> &scratch_base) const override;
    };

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

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


  namespace Melt
  {
    template <int dim>
    struct Parameters
    {
      static void declare_parameters (ParameterHandler &prm);

      void parse_parameters (ParameterHandler &prm);

      double melt_scaling_factor_threshold;

      bool heat_advection_by_melt;

      bool use_discontinuous_p_c;

      bool average_melt_velocity;
    };
  }


  template <int dim>
  class MeltHandler: public SimulatorAccess<dim>
  {
    public:
      MeltHandler(ParameterHandler &prm);

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

      void edit_finite_element_variables(const Parameters<dim> &parameters,
                                         std::vector<VariableDeclaration<dim>> &variables);

      void set_assemblers (Assemblers::Manager<dim> &assemblers) const;


      void initialize() const;

      void initialize_simulator (const Simulator<dim> &simulator_object) override;

      void compute_melt_variables(LinearAlgebra::BlockSparseMatrix &system_matrix,
                                  LinearAlgebra::BlockVector &solution,
                                  LinearAlgebra::BlockVector &system_rhs) const;

      bool is_porosity (const AdvectionField &advection_field) const;

      void apply_free_surface_stabilization_with_melt (const double free_surface_theta,
                                                       const typename DoFHandler<dim>::active_cell_iterator &cell,
                                                       internal::Assembly::Scratch::StokesSystem<dim>       &scratch,
                                                       internal::Assembly::CopyData::StokesSystem<dim>      &data) const;

      void add_current_constraints(AffineConstraints<double> &constraints);

      bool is_melt_cell(const typename DoFHandler<dim>::active_cell_iterator &cell) const;

      double limited_darcy_coefficient(const double K_D,
                                       const bool is_melt_cell) const;

      const BoundaryFluidPressure::Interface<dim> &
      get_boundary_fluid_pressure () const;

      Melt::Parameters<dim> melt_parameters;

    private:
      const std::unique_ptr<aspect::BoundaryFluidPressure::Interface<dim>> boundary_fluid_pressure;

      std::vector<bool> is_melt_cell_vector;

      AffineConstraints<double> current_constraints;

  };

}

#endif