simulator.h

Go to the documentation of this file.

/*
  Copyright (C) 2011 - 2026 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_simulator_h
#define _aspect_simulator_h

#include <fstream>

#include <deal.II/base/timer.h>
#include <deal.II/base/parameter_handler.h>
#include <deal.II/base/conditional_ostream.h>
#include <deal.II/base/symmetric_tensor.h>

DEAL_II_DISABLE_EXTRA_DIAGNOSTICS

#include <deal.II/lac/affine_constraints.h>

#include <deal.II/distributed/tria.h>

#include <deal.II/dofs/dof_handler.h>
#include <deal.II/dofs/dof_tools.h>

#include <deal.II/fe/fe_system.h>
#include <deal.II/fe/mapping.h>
#include <deal.II/base/tensor_function.h>

DEAL_II_ENABLE_EXTRA_DIAGNOSTICS

#include <aspect/global.h>
#include <aspect/linear_algebra_types.h>
#include <aspect/simulator_access.h>
#include <aspect/lateral_averaging.h>
#include <aspect/simulator_signals.h>
#include <aspect/material_model/interface.h>
#include <aspect/heating_model/interface.h>
#include <aspect/prescribed_dilation/interface.h>
#include <aspect/geometry_model/initial_topography_model/interface.h>
#include <aspect/geometry_model/interface.h>
#include <aspect/gravity_model/interface.h>
#include <aspect/boundary_temperature/interface.h>
#include <aspect/boundary_heat_flux/interface.h>
#include <aspect/boundary_convective_heating/interface.h>
#include <aspect/boundary_composition/interface.h>
#include <aspect/initial_temperature/interface.h>
#include <aspect/initial_composition/interface.h>
#include <aspect/prescribed_solution/interface.h>
#include <aspect/prescribed_stokes_solution/interface.h>
#include <aspect/boundary_velocity/interface.h>
#include <aspect/boundary_fluid_pressure/interface.h>
#include <aspect/boundary_traction/interface.h>
#include <aspect/mesh_refinement/interface.h>
#include <aspect/time_stepping/interface.h>
#include <aspect/postprocess/interface.h>
#include <aspect/adiabatic_conditions/interface.h>
#include <aspect/particle/manager.h>
#include <aspect/advection_field.h>

#include <boost/iostreams/tee.hpp>
#include <boost/iostreams/stream.hpp>

#include <memory>
#include <thread>

namespace WorldBuilder
{
  class World;
}


namespace aspect
{
  template <int dim>
  class MeltHandler;

  template <int dim>
  class NewtonHandler;

  template <int dim>
  class StokesMatrixFreeHandler;

  namespace StokesSolver
  {
    template <int dim>
    class Interface;

    template <int dim>
    class Direct;

    template <int dim>
    class MatrixBased;
  }

  template <int dim, int velocity_degree>
  class StokesMatrixFreeHandlerLocalSmoothingImplementation;

  namespace MeshDeformation
  {
    template <int dim>
    class MeshDeformationHandler;
  }

  template <int dim>
  class VolumeOfFluidHandler;

  namespace internal
  {
    namespace Assembly
    {
      namespace Scratch
      {
        template <int dim>      struct StokesPreconditioner;
        template <int dim>      struct StokesSystem;
        template <int dim>      struct AdvectionSystem;
      }

      namespace CopyData
      {
        template <int dim>      struct StokesPreconditioner;
        template <int dim>      struct StokesSystem;
        template <int dim>      struct AdvectionSystem;
      }
    }
  }

  namespace Assemblers
  {
    template <int dim>      class Interface;
    template <int dim>      class Manager;
  }

  struct DefectCorrectionResiduals
  {
    double initial_residual;
    double velocity_residual;
    double pressure_residual;
    double residual;
    double residual_old;
    double switch_initial_residual;
    double newton_residual_for_derivative_scaling_factor;
    std::pair<double,double> stokes_residuals;
  };

  template <int dim>
  struct RotationProperties
  {
    RotationProperties()
      :
      scalar_moment_of_inertia(numbers::signaling_nan<double>()),
      scalar_angular_momentum(numbers::signaling_nan<double>()),
      scalar_rotation(numbers::signaling_nan<double>()),
      tensor_moment_of_inertia(numbers::signaling_nan<SymmetricTensor<2,dim>>()),
      tensor_angular_momentum(numbers::signaling_nan<Tensor<1,dim>>()),
      tensor_rotation(numbers::signaling_nan<Tensor<1,dim>>())
    {};

    double scalar_moment_of_inertia;
    double scalar_angular_momentum;
    double scalar_rotation;

    SymmetricTensor<2,dim> tensor_moment_of_inertia;
    Tensor<1,dim> tensor_angular_momentum;
    Tensor<1,dim> tensor_rotation;
  };

  DeclExceptionMsg(ExcNonlinearSolverNoConvergence,
                   "Nonlinear solver failed to converge in the prescribed number of steps. "
                   "Consider changing `Max nonlinear iterations` or `Nonlinear solver failure "
                   "strategy`.");

  template <int dim>
  class Simulator : public EnableObserverPointer
  {
    public:
      Simulator (const MPI_Comm mpi_communicator,
                 ParameterHandler &prm);

      ~Simulator ();

      static
      void declare_parameters (ParameterHandler &prm, const unsigned int mpi_rank);

      void run ();

      void
      write_plugin_graph (std::ostream &output_stream) const;

      using NonlinearSolver = typename Parameters<dim>::NonlinearSolver;

      using NullspaceRemoval = typename Parameters<dim>::NullspaceRemoval;

      using AdvectionField = aspect::AdvectionField;

    private:

      bool simulator_is_past_initialization;

      struct IntermediaryConstructorAction
      {
        IntermediaryConstructorAction (const std::function<void ()> &action);
      };

      void setup_dofs ();

      void setup_introspection ();

      void set_initial_temperature_and_compositional_fields ();

      void compute_initial_pressure_field ();

      void compute_initial_velocity_boundary_constraints (AffineConstraints<double> &constraints);

      void compute_current_velocity_boundary_constraints (AffineConstraints<double> &constraints);

      void compute_current_constraints ();

      double compute_pressure_scaling_factor () const;

      void start_timestep ();

      void solve_timestep ();

      void solve_no_advection_no_stokes ();

      void solve_no_advection_single_stokes ();

      void solve_no_advection_single_stokes_first_timestep_only ();

      void solve_no_advection_iterated_stokes ();

      void solve_no_advection_iterated_defect_correction_stokes ();

      void solve_single_advection_no_stokes ();

      void solve_single_advection_single_stokes ();

      void solve_single_advection_iterated_stokes ();

      void solve_single_advection_iterated_defect_correction_stokes ();

      void solve_single_advection_iterated_newton_stokes (bool use_newton_iterations);

      void solve_iterated_advection_no_stokes ();

      void solve_iterated_advection_and_stokes ();

      void solve_iterated_advection_and_defect_correction_stokes ();

      void solve_iterated_advection_and_newton_stokes (bool use_newton_iterations);

      void build_stokes_preconditioner ();

      void build_advection_preconditioner (const AdvectionField &advection_field,
                                           aspect::LinearAlgebra::PreconditionILU &preconditioner,
                                           const double diagonal_strengthening);

      void assemble_stokes_system ();

      double assemble_and_solve_temperature (const double &initial_residual = 0,
                                             double *residual = nullptr);

      std::vector<double> assemble_and_solve_composition (const std::vector<double> &initial_residual = {},
                                                          const unsigned int nonlinear_iteration = 0,
                                                          std::vector<double> *residual = nullptr);

      double assemble_and_solve_stokes (const double &initial_nonlinear_residual = 0,
                                        double *nonlinear_residual = nullptr);

      void do_one_defect_correction_Stokes_step(DefectCorrectionResiduals &dcr,
                                                const bool use_picard);

      void assemble_advection_system (const AdvectionField &advection_field);

      double solve_advection (const AdvectionField &advection_field);

      void interpolate_particle_properties (const std::vector<AdvectionField> &advection_fields);

      std::pair<double,double>
      solve_stokes (LinearAlgebra::BlockVector &solution_vector);

      void postprocess ();

      void refine_mesh (const unsigned int max_grid_level);

      unsigned int determine_last_good_snapshot() const;

      unsigned int determine_resume_snapshot() const;

      void create_snapshot(const bool is_additional_checkpoint = false);

      void resume_from_snapshot();

      template <class Archive>
      void serialize (Archive &ar, const unsigned int version);
      Table<2,DoFTools::Coupling>
      setup_system_matrix_coupling () const;

      void setup_system_matrix (const std::vector<IndexSet> &system_partitioning);

      void setup_system_preconditioner (const std::vector<IndexSet> &system_partitioning);

      std::unique_ptr<Assemblers::Manager<dim>> assemblers;

      void set_assemblers ();

      void set_advection_assemblers ();

      void set_stokes_assemblers ();

      void assemble_stokes_preconditioner ();

      void
      local_assemble_stokes_preconditioner (const typename DoFHandler<dim>::active_cell_iterator &cell,
                                            internal::Assembly::Scratch::StokesPreconditioner<dim> &scratch,
                                            internal::Assembly::CopyData::StokesPreconditioner<dim> &data);

      void
      copy_local_to_global_stokes_preconditioner (const internal::Assembly::CopyData::StokesPreconditioner<dim> &data);

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

      void
      copy_local_to_global_stokes_system (const internal::Assembly::CopyData::StokesSystem<dim> &data);

      void
      local_assemble_advection_face_terms(const AdvectionField &advection_field,
                                          const typename DoFHandler<dim>::active_cell_iterator &cell,
                                          internal::Assembly::Scratch::AdvectionSystem<dim> &scratch,
                                          internal::Assembly::CopyData::AdvectionSystem<dim> &data);
      void
      local_assemble_advection_system (const AdvectionField &advection_field,
                                       const Vector<double>           &viscosity_per_cell,
                                       const typename DoFHandler<dim>::active_cell_iterator &cell,
                                       internal::Assembly::Scratch::AdvectionSystem<dim>  &scratch,
                                       internal::Assembly::CopyData::AdvectionSystem<dim> &data);

      void
      copy_local_to_global_advection_system (const AdvectionField &advection_field,
                                             const internal::Assembly::CopyData::AdvectionSystem<dim> &data);

      void make_pressure_rhs_compatible(LinearAlgebra::BlockVector &vector);

      template <typename T>
      void get_artificial_viscosity (Vector<T> &viscosity_per_cell,
                                     const AdvectionField &advection_field,
                                     const bool skip_interior_cells = false) const;

      double normalize_pressure(LinearAlgebra::BlockVector &vector) const;

      void denormalize_pressure(const double                      pressure_adjustment,
                                LinearAlgebra::BlockVector       &vector) const;

      void apply_limiter_to_dg_solutions (const AdvectionField &advection_field);

      void compute_unique_advection_support_points (const std::vector<AdvectionField> &advection_fields,
                                                    std::vector<Point<dim>> &unique_support_points,
                                                    std::vector<std::vector<unsigned int>> &support_point_index_by_field) const;

      void compute_reactions (const std::vector<AdvectionField> &advection_fields_with_reactions);

      void update_solution_vectors_with_reaction_results (const unsigned int block_index,
                                                          const LinearAlgebra::BlockVector &distributed_vector,
                                                          const LinearAlgebra::BlockVector &distributed_reaction_vector);

      void initialize_current_linearization_point ();

      void interpolate_material_output_into_advection_field (const std::vector<AdvectionField> &adv_field);


      void interpolate_onto_velocity_system(const TensorFunction<1,dim> &func,
                                            LinearAlgebra::Vector &vec) const;

      double perform_line_search(const DefectCorrectionResiduals &dcr,
                                 const bool use_picard,
                                 const LinearAlgebra::BlockVector &search_direction);

      void setup_nullspace_constraints(AffineConstraints<double> &constraints);


      void remove_nullspace(LinearAlgebra::BlockVector &solution,
                            LinearAlgebra::BlockVector &distributed_stokes_solution) const;

      RotationProperties<dim>
      compute_net_angular_momentum(const bool use_constant_density,
                                   const LinearAlgebra::BlockVector &solution,
                                   const bool limit_to_top_faces = false) const;

      void remove_net_angular_momentum(const bool use_constant_density,
                                       LinearAlgebra::BlockVector &relevant_dst,
                                       LinearAlgebra::BlockVector &tmp_distributed_stokes,
                                       const bool limit_to_top_faces = false) const;

      void replace_outflow_boundary_ids(const unsigned int boundary_id_offset,
                                        const bool is_composition,
                                        const unsigned int composition_index);

      void restore_outflow_boundary_ids(const unsigned int boundary_id_offset);

      void remove_net_linear_momentum(const bool use_constant_density,
                                      LinearAlgebra::BlockVector &relevant_dst,
                                      LinearAlgebra::BlockVector &tmp_distributed_stokes) const;

      double get_maximal_velocity (const LinearAlgebra::BlockVector &solution) const;

      double get_entropy_variation (const double average_field,
                                    const AdvectionField &advection_field) const;

      std::pair<double,double>
      get_extrapolated_advection_field_range (const AdvectionField &advection_field) const;

      void exchange_refinement_flags();


      bool synchronize_periodic_refinement_flags();


      void prepare_periodic_refinement();


      void maybe_write_timing_output () const;

      bool maybe_write_checkpoint (const std::time_t last_checkpoint_time,
                                   const bool        force_writing_checkpoint);

      bool maybe_do_initial_refinement (const unsigned int max_refinement_level);

      void maybe_refine_mesh (const double new_time_step,
                              unsigned int &max_refinement_level);

      void advance_time (const double step_size);

      double
      compute_viscosity(internal::Assembly::Scratch::AdvectionSystem<dim> &scratch,
                        const double                        global_u_infty,
                        const double                        global_field_variation,
                        const double                        average_field,
                        const double                        global_entropy_variation,
                        const double                        cell_diameter,
                        const AdvectionField     &advection_field) const;

      void
      compute_advection_system_residual(internal::Assembly::Scratch::AdvectionSystem<dim> &scratch,
                                        const double                        average_field,
                                        const AdvectionField               &advection_field,
                                        double                             &max_residual,
                                        double                             &max_velocity,
                                        double                             &max_density,
                                        double                             &max_specific_heat,
                                        double                             &conductivity) const;

      bool
      stokes_matrix_depends_on_solution () const;

      bool
      stokes_A_block_is_symmetric () const;

      bool
      is_stokes_matrix_free() const;

      void
      check_consistency_of_formulation ();

      void
      select_default_solver_and_averaging ();

      void
      check_consistency_of_boundary_conditions () const;

      double
      compute_initial_newton_residual ();

      double
      compute_Eisenstat_Walker_linear_tolerance(const bool EisenstatWalkerChoiceOne,
                                                const double maximum_linear_stokes_solver_tolerance,
                                                const double linear_stokes_solver_tolerance,
                                                const double stokes_residual,
                                                const double newton_residual,
                                                const double newton_residual_old);

      void output_statistics();

      double
      compute_initial_stokes_residual();

      Parameters<dim>                     parameters;

      std::unique_ptr<MeltHandler<dim>> melt_handler;

      std::unique_ptr<NewtonHandler<dim>> newton_handler;

      SimulatorSignals<dim>               signals;

      const IntermediaryConstructorAction post_signal_creation;

      std::unique_ptr<VolumeOfFluidHandler<dim>> volume_of_fluid_handler;

      Introspection<dim>                  introspection;


      MPI_Comm                            mpi_communicator;

      std::ofstream log_file_stream;

      using TeeDevice = boost::iostreams::tee_device<std::ostream, std::ofstream>;
      using TeeStream = boost::iostreams::stream<TeeDevice>;

      TeeDevice iostream_tee_device;
      TeeStream iostream_tee_stream;

      ConditionalOStream                  pcout;

      TableHandler                        statistics;

      std::size_t                         statistics_last_write_size;
      std::size_t                         statistics_last_hash;

      mutable TimerOutput                 computing_timer;

      Timer wall_timer;

      double total_walltime_until_last_snapshot;

      unsigned int last_checkpoint_id;

      unsigned int last_regular_checkpoint_id;

      unsigned int last_additional_checkpoint_id;

      std::thread                         output_statistics_thread;

      const std::shared_ptr<InitialTopographyModel::Interface<dim>>          initial_topography_model;
      const std::unique_ptr<GeometryModel::Interface<dim>>                   geometry_model;
      const IntermediaryConstructorAction                                    post_geometry_model_creation_action;
      const std::unique_ptr<MaterialModel::Interface<dim>>                   material_model;
      const std::unique_ptr<GravityModel::Interface<dim>>                    gravity_model;

      BoundaryTemperature::Manager<dim>                                      boundary_temperature_manager;
      BoundaryConvectiveHeating::Manager<dim>                                boundary_convective_heating_manager;
      BoundaryComposition::Manager<dim>                                      boundary_composition_manager;
      PrescribedSolution::Manager<dim>                                       prescribed_solution_manager;
      const std::unique_ptr<PrescribedStokesSolution::Interface<dim>>        prescribed_stokes_solution;

      std::shared_ptr<InitialTemperature::Manager<dim>>                      initial_temperature_manager;
      std::shared_ptr<InitialComposition::Manager<dim>>                      initial_composition_manager;

      const std::unique_ptr<AdiabaticConditions::Interface<dim>>             adiabatic_conditions;
#ifdef ASPECT_WITH_WORLD_BUILDER

      std::shared_ptr<WorldBuilder::World>                      world_builder;
#endif
      BoundaryVelocity::Manager<dim>                            boundary_velocity_manager;
      BoundaryTraction::Manager<dim>                            boundary_traction_manager;
      const std::unique_ptr<BoundaryHeatFlux::Interface<dim>>   boundary_heat_flux;

      double                                                    time;
      double                                                    time_step;
      double                                                    old_time_step;
      unsigned int                                              timestep_number;
      unsigned int                                              pre_refinement_step;
      unsigned int                                              nonlinear_iteration;
      unsigned int                                              nonlinear_solver_failures;
      unsigned int                                              linear_solver_failures;
      TimeStepping::Manager<dim>                                time_stepping_manager;
      LateralAveraging<dim>                                     lateral_averaging;
      parallel::distributed::Triangulation<dim>                 triangulation;
      double                                                    global_Omega_diameter;
      double                                                    global_volume;

      MeshRefinement::Manager<dim>                              mesh_refinement_manager;
      HeatingModel::Manager<dim>                                heating_model_manager;
      PrescribedDilation::Manager<dim>                          prescribed_dilation_manager;

      std::unique_ptr<Mapping<dim>>                             mapping;

      const FESystem<dim>                                       finite_element;

      DoFHandler<dim>                                           dof_handler;

      Postprocess::Manager<dim>                                 postprocess_manager;

      std::vector<Particle::Manager<dim>>                       particle_managers;

      AffineConstraints<double>                                 constraints;
      AffineConstraints<double>                                 current_constraints;

      double                                                    last_pressure_normalization_adjustment;

      double                                                    pressure_scaling;

      bool                           do_pressure_rhs_compatibility_modification;

      LinearAlgebra::BlockSparseMatrix                          system_matrix;

      LinearAlgebra::BlockVector                                inverse_lumped_mass_matrix;

      LinearAlgebra::BlockSparseMatrix                          system_preconditioner_matrix;

      LinearAlgebra::BlockVector                                solution;
      LinearAlgebra::BlockVector                                old_solution;
      LinearAlgebra::BlockVector                                old_old_solution;
      LinearAlgebra::BlockVector                                system_rhs;

      LinearAlgebra::BlockVector                                current_linearization_point;

      // only used if is_compressible()
      LinearAlgebra::BlockVector                                pressure_shape_function_integrals;

      // only used if operator split is enabled
      LinearAlgebra::BlockVector                                operator_split_reaction_vector;



      std::unique_ptr<LinearAlgebra::PreconditionAMG>           Amg_preconditioner;
      std::unique_ptr<LinearAlgebra::PreconditionBase>          Mp_preconditioner;

      bool                                                      rebuild_sparsity_and_matrices;

      bool                                                      rebuild_stokes_matrix;

      bool                                                      assemble_newton_stokes_matrix;

      bool                                                      assemble_newton_stokes_system;

      bool                                                      rebuild_stokes_preconditioner;

      std::unique_ptr<MeshDeformation::MeshDeformationHandler<dim>> mesh_deformation;

      std::unique_ptr<StokesSolver::Interface<dim>> stokes_solver;

      friend class boost::serialization::access;
      friend class SimulatorAccess<dim>;
      friend class MeshDeformation::MeshDeformationHandler<dim>;
      friend class VolumeOfFluidHandler<dim>;
      friend class StokesMatrixFreeHandler<dim>;
      friend class StokesSolver::MatrixBased<dim>;
      template <int dimension, int velocity_degree> friend class StokesMatrixFreeHandlerLocalSmoothingImplementation;
      friend struct Parameters<dim>;
  };
}


#endif