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