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