ASPECT
simulator.h
Go to the documentation of this file.
1 /*
2  Copyright (C) 2011 - 2020 by the authors of the ASPECT code.
3 
4  This file is part of ASPECT.
5 
6  ASPECT is free software; you can redistribute it and/or modify
7  it under the terms of the GNU General Public License as published by
8  the Free Software Foundation; either version 2, or (at your option)
9  any later version.
10 
11  ASPECT is distributed in the hope that it will be useful,
12  but WITHOUT ANY WARRANTY; without even the implied warranty of
13  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14  GNU General Public License for more details.
15 
16  You should have received a copy of the GNU General Public License
17  along with ASPECT; see the file LICENSE. If not see
18  <http://www.gnu.org/licenses/>.
19 */
20 
21 
22 #ifndef _aspect_simulator_h
23 #define _aspect_simulator_h
24 
25 #include <deal.II/base/timer.h>
29 
31 
33 
35 
37 #include <deal.II/dofs/dof_tools.h>
38 
39 #include <deal.II/fe/fe_system.h>
40 #include <deal.II/fe/mapping.h>
42 
44 
45 #include <aspect/global.h>
67 #include <aspect/particle/world.h>
68 
69 #include <boost/iostreams/tee.hpp>
70 #include <boost/iostreams/stream.hpp>
71 #include <memory>
72 
73 namespace aspect
74 {
75  using namespace dealii;
76 
77  template <int dim>
78  class MeltHandler;
79 
80  template <int dim>
81  class NewtonHandler;
82 
83  template <int dim>
85 
86  template <int dim, int velocity_degree>
88 
89  namespace MeshDeformation
90  {
91  template <int dim>
93  }
94 
95  template <int dim>
97 
98  namespace internal
99  {
100  namespace Assembly
101  {
102  namespace Scratch
103  {
104  template <int dim> struct StokesPreconditioner;
105  template <int dim> struct StokesSystem;
106  template <int dim> struct AdvectionSystem;
107  }
108 
109  namespace CopyData
110  {
111  template <int dim> struct StokesPreconditioner;
112  template <int dim> struct StokesSystem;
113  template <int dim> struct AdvectionSystem;
114  }
115  }
116  }
117 
118  namespace Assemblers
119  {
120  template <int dim> class Interface;
121  template <int dim> class Manager;
122  }
123 
125  {
129  double residual;
130  double residual_old;
133  std::pair<double,double> stokes_residuals;
134  };
135 
139  template <int dim>
141  {
143  :
144  scalar_moment_of_inertia(numbers::signaling_nan<double>()),
145  scalar_angular_momentum(numbers::signaling_nan<double>()),
146  scalar_rotation(numbers::signaling_nan<double>()),
147  tensor_moment_of_inertia(numbers::signaling_nan<SymmetricTensor<2,dim>>()),
148  tensor_angular_momentum(numbers::signaling_nan<Tensor<1,dim>>()),
149  tensor_rotation(numbers::signaling_nan<Tensor<1,dim>>())
150  {};
151 
159 
166  };
167 
175  template <int dim>
176  class Simulator
177  {
178  public:
192  Simulator (const MPI_Comm mpi_communicator,
193  ParameterHandler &prm);
194 
199  ~Simulator ();
200 
212  static
214 
223  void run ();
224 
236  void
237  write_plugin_graph (std::ostream &output_stream) const;
238 
243 
248 
249 
257  {
262  enum FieldType { temperature_field, compositional_field };
263 
269 
275  const unsigned int compositional_variable;
276 
288  AdvectionField (const FieldType field_type,
289  const unsigned int compositional_variable = numbers::invalid_unsigned_int);
290 
298  static
300 
308  static
309  AdvectionField composition (const unsigned int compositional_variable);
310 
314  bool
315  is_temperature () const;
316 
321  bool
322  is_discontinuous (const Introspection<dim> &introspection) const;
323 
329  advection_method (const Introspection<dim> &introspection) const;
330 
335  unsigned int component_index(const Introspection<dim> &introspection) const;
336 
341  unsigned int block_index(const Introspection<dim> &introspection) const;
342 
349  unsigned int field_index() const;
350 
356  unsigned int base_element(const Introspection<dim> &introspection) const;
357 
364  FEValuesExtractors::Scalar scalar_extractor(const Introspection<dim> &introspection) const;
365 
370  unsigned int polynomial_degree(const Introspection<dim> &introspection) const;
371  };
372 
373  private:
374 
394 
412  {
413  IntermediaryConstructorAction (const std::function<void ()> &action);
414  };
415 
430  void setup_dofs ();
431 
440  void setup_introspection ();
441 
453  void set_initial_temperature_and_compositional_fields ();
454 
470  void compute_initial_pressure_field ();
471 
477  void compute_initial_velocity_boundary_constraints (AffineConstraints<double> &constraints);
478 
484  void compute_current_velocity_boundary_constraints (AffineConstraints<double> &constraints);
485 
498  void compute_current_constraints ();
499 
512  void compute_pressure_scaling_factor ();
513 
525  void start_timestep ();
526 
534  void solve_timestep ();
535 
547  void solve_single_advection_single_stokes ();
548 
561  void solve_no_advection_iterated_stokes ();
562 
573  void solve_no_advection_single_stokes ();
574 
586  void solve_first_timestep_only_single_stokes ();
587 
600  void solve_iterated_advection_and_stokes ();
601 
614  void solve_single_advection_iterated_stokes ();
615 
628  void solve_no_advection_iterated_defect_correction_stokes ();
629 
642  void solve_single_advection_iterated_defect_correction_stokes ();
643 
657  void solve_iterated_advection_and_defect_correction_stokes ();
658 
673  void solve_iterated_advection_and_newton_stokes ();
674 
690  void solve_single_advection_iterated_newton_stokes ();
691 
703  void solve_single_advection_no_stokes ();
704 
716  void solve_no_advection_no_stokes ();
717 
726  void build_stokes_preconditioner ();
727 
735  void build_advection_preconditioner (const AdvectionField &advection_field,
737  const double diagonal_strengthening);
738 
745  void assemble_stokes_system ();
746 
756  double assemble_and_solve_temperature (const bool compute_initial_residual = false,
757  double *initial_residual = nullptr);
758 
771  std::vector<double> assemble_and_solve_composition (const bool compute_initial_residual = false,
772  std::vector<double> *initial_residual = nullptr);
773 
791  double assemble_and_solve_stokes (const bool compute_initial_residual = false,
792  double *initial_nonlinear_residual = nullptr);
793 
805  void assemble_and_solve_defect_correction_Stokes(DefectCorrectionResiduals &dcr,
806  const bool use_picard);
807 
815  void assemble_advection_system (const AdvectionField &advection_field);
816 
827  double solve_advection (const AdvectionField &advection_field);
828 
832  void interpolate_particle_properties (const AdvectionField &advection_field);
833 
904  std::pair<double,double>
905  solve_stokes ();
906 
910  std::pair<double,double>
911  solve_stokes_block_gmg ();
912 
925  void postprocess ();
926 
942  void refine_mesh (const unsigned int max_grid_level);
943 
964  void create_snapshot();
965 
979  void resume_from_snapshot();
980 
987  template <class Archive>
988  void serialize (Archive &ar, const unsigned int version);
1006  setup_system_matrix_coupling () const;
1007 
1015  void setup_system_matrix (const std::vector<IndexSet> &system_partitioning);
1016 
1028  void setup_system_preconditioner (const std::vector<IndexSet> &system_partitioning);
1029 
1051  std::unique_ptr<Assemblers::Manager<dim> > assemblers;
1052 
1061  void set_assemblers ();
1062 
1073  void set_advection_assemblers ();
1074 
1084  void set_stokes_assemblers ();
1085 
1092  void assemble_stokes_preconditioner ();
1093 
1101  void
1102  local_assemble_stokes_preconditioner (const typename DoFHandler<dim>::active_cell_iterator &cell,
1105 
1113  void
1114  copy_local_to_global_stokes_preconditioner (const internal::Assembly::CopyData::StokesPreconditioner<dim> &data);
1115 
1123  void
1124  local_assemble_stokes_system (const typename DoFHandler<dim>::active_cell_iterator &cell,
1127 
1135  void
1136  copy_local_to_global_stokes_system (const internal::Assembly::CopyData::StokesSystem<dim> &data);
1137 
1145  void
1146  local_assemble_advection_face_terms(const AdvectionField &advection_field,
1147  const typename DoFHandler<dim>::active_cell_iterator &cell,
1157  void
1158  local_assemble_advection_system (const AdvectionField &advection_field,
1159  const Vector<double> &viscosity_per_cell,
1160  const typename DoFHandler<dim>::active_cell_iterator &cell,
1163 
1171  void
1172  copy_local_to_global_advection_system (const AdvectionField &advection_field,
1174 
1199  void make_pressure_rhs_compatible(LinearAlgebra::BlockVector &vector);
1200 
1210  template <typename T>
1211  void get_artificial_viscosity (Vector<T> &viscosity_per_cell,
1212  const AdvectionField &advection_field,
1213  const bool skip_interior_cells = false) const;
1214 
1226  void compute_Vs_anomaly(Vector<float> &values) const;
1227 
1242  void compute_Vp_anomaly(Vector<float> &values) const;
1243 
1284  double normalize_pressure(LinearAlgebra::BlockVector &vector) const;
1285 
1330  void denormalize_pressure(const double pressure_adjustment,
1332  const LinearAlgebra::BlockVector &relevant_vector) const;
1333 
1342  void apply_limiter_to_dg_solutions (const AdvectionField &advection_field);
1343 
1344 
1367  void compute_reactions ();
1368 
1378  void update_solution_vectors_with_reaction_results (const unsigned int block_index,
1379  const LinearAlgebra::BlockVector &distributed_vector,
1380  const LinearAlgebra::BlockVector &distributed_reaction_vector);
1381 
1392  void initialize_current_linearization_point ();
1393 
1413  void interpolate_material_output_into_advection_field (const AdvectionField &adv_field);
1414 
1415 
1423  void interpolate_onto_velocity_system(const TensorFunction<1,dim> &func,
1424  LinearAlgebra::Vector &vec);
1425 
1426 
1441  void setup_nullspace_constraints(AffineConstraints<double> &constraints);
1442 
1443 
1455  void remove_nullspace(LinearAlgebra::BlockVector &relevant_dst,
1456  LinearAlgebra::BlockVector &tmp_distributed_stokes);
1457 
1473  compute_net_angular_momentum(const bool use_constant_density,
1474  const LinearAlgebra::BlockVector &solution,
1475  const bool limit_to_top_faces = false) const;
1476 
1493  void remove_net_angular_momentum( const bool use_constant_density,
1494  LinearAlgebra::BlockVector &relevant_dst,
1495  LinearAlgebra::BlockVector &tmp_distributed_stokes,
1496  const bool limit_to_top_faces = false);
1497 
1505  void replace_outflow_boundary_ids(const unsigned int boundary_id_offset);
1506 
1514  void restore_outflow_boundary_ids(const unsigned int boundary_id_offset);
1515 
1529  void remove_net_linear_momentum( const bool use_constant_density,
1530  LinearAlgebra::BlockVector &relevant_dst,
1531  LinearAlgebra::BlockVector &tmp_distributed_stokes);
1532 
1540  double get_maximal_velocity (const LinearAlgebra::BlockVector &solution) const;
1541 
1554  double get_entropy_variation (const double average_field,
1555  const AdvectionField &advection_field) const;
1556 
1565  std::pair<double,double>
1566  get_extrapolated_advection_field_range (const AdvectionField &advection_field) const;
1567 
1576  void maybe_write_timing_output () const;
1577 
1585  bool maybe_write_checkpoint (const time_t last_checkpoint_time,
1586  const bool force_writing_checkpoint);
1587 
1604  bool maybe_do_initial_refinement (const unsigned int max_refinement_level);
1605 
1614  void maybe_refine_mesh (const double new_time_step,
1615  unsigned int &max_refinement_level);
1616 
1624  void advance_time (const double step_size);
1625 
1633  double
1634  compute_viscosity(internal::Assembly::Scratch::AdvectionSystem<dim> &scratch,
1635  const double global_u_infty,
1636  const double global_field_variation,
1637  const double average_field,
1638  const double global_entropy_variation,
1639  const double cell_diameter,
1640  const AdvectionField &advection_field) const;
1641 
1650  void
1651  compute_advection_system_residual(internal::Assembly::Scratch::AdvectionSystem<dim> &scratch,
1652  const double average_field,
1653  const AdvectionField &advection_field,
1654  double &max_residual,
1655  double &max_velocity,
1656  double &max_density,
1657  double &max_specific_heat,
1658  double &conductivity) const;
1659 
1687  void
1688  compute_material_model_input_values (const LinearAlgebra::BlockVector &input_solution,
1689  const FEValuesBase<dim,dim> &input_finite_element_values,
1690  const typename DoFHandler<dim>::active_cell_iterator &cell,
1691  const bool compute_strainrate,
1692  MaterialModel::MaterialModelInputs<dim> &material_model_inputs) const;
1693 
1694 
1709  bool
1710  stokes_matrix_depends_on_solution () const;
1711 
1724  void
1725  check_consistency_of_formulation ();
1726 
1737  void
1738  check_consistency_of_boundary_conditions () const;
1739 
1743  double
1744  compute_initial_newton_residual (const LinearAlgebra::BlockVector &linearized_stokes_initial_guess);
1745 
1754  double
1755  compute_Eisenstat_Walker_linear_tolerance(const bool EisenstatWalkerChoiceOne,
1756  const double maximum_linear_stokes_solver_tolerance,
1757  const double linear_stokes_solver_tolerance,
1758  const double stokes_residual,
1759  const double newton_residual,
1760  const double newton_residual_old);
1761 
1771  void output_statistics();
1772 
1784  double
1785  compute_initial_stokes_residual();
1786 
1797 
1803  std::unique_ptr<MeltHandler<dim> > melt_handler;
1804 
1810  std::unique_ptr<NewtonHandler<dim> > newton_handler;
1811 
1813 
1815 
1823  std::unique_ptr<VolumeOfFluidHandler<dim> > volume_of_fluid_handler;
1824 
1826 
1827 
1829 
1834  std::ofstream log_file_stream;
1835 
1836  using TeeDevice = boost::iostreams::tee_device<std::ostream, std::ofstream>;
1837  using TeeStream = boost::iostreams::stream< TeeDevice >;
1838 
1841 
1847 
1858 
1875 
1877 
1894  const std::unique_ptr<InitialTopographyModel::Interface<dim> > initial_topography_model;
1895  const std::unique_ptr<GeometryModel::Interface<dim> > geometry_model;
1897  const std::unique_ptr<MaterialModel::Interface<dim> > material_model;
1898  const std::unique_ptr<GravityModel::Interface<dim> > gravity_model;
1901  const std::unique_ptr<PrescribedStokesSolution::Interface<dim> > prescribed_stokes_solution;
1904  const std::unique_ptr<AdiabaticConditions::Interface<dim> > adiabatic_conditions;
1905 #ifdef ASPECT_WITH_WORLD_BUILDER
1906  const std::unique_ptr<WorldBuilder::World> world_builder;
1907 #endif
1909  std::map<types::boundary_id,std::unique_ptr<BoundaryTraction::Interface<dim> > > boundary_traction;
1910  const std::unique_ptr<BoundaryHeatFlux::Interface<dim> > boundary_heat_flux;
1911 
1915  std::unique_ptr<Particle::World<dim> > particle_world;
1916 
1922 
1930  double time;
1931  double time_step;
1933  unsigned int timestep_number;
1934  unsigned int pre_refinement_step;
1935  unsigned int nonlinear_iteration;
1965 
1968 
1977  std::unique_ptr<Mapping<dim> > mapping;
1978 
1980 
1982 
1984 
1998 
2004 
2011 
2019 
2039 
2053 
2058 
2060 
2061  // only used if is_compressible()
2063 
2064  // only used if operator split is enabled
2066 
2067 
2068 
2069  std::unique_ptr<LinearAlgebra::PreconditionAMG> Amg_preconditioner;
2070  std::unique_ptr<LinearAlgebra::PreconditionBase> Mp_preconditioner;
2071 
2077 
2082  private:
2083 
2089  std::unique_ptr<MeshDeformation::MeshDeformationHandler<dim> > mesh_deformation;
2090 
2094  std::unique_ptr<StokesMatrixFreeHandler<dim> > stokes_matrix_free;
2095 
2096  friend class boost::serialization::access;
2097  friend class SimulatorAccess<dim>;
2098  friend class MeshDeformation::MeshDeformationHandler<dim>; // MeshDeformationHandler needs access to the internals of the Simulator
2099  friend class VolumeOfFluidHandler<dim>; // VolumeOfFluidHandler needs access to the internals of the Simulator
2100  friend class StokesMatrixFreeHandler<dim>;
2101  template <int dimension, int velocity_degree>
2103  friend struct Parameters<dim>;
2104  };
2105 }
2106 
2107 
2108 #endif
The NullspaceRemoval struct.
Definition: parameters.h:101
unsigned int nonlinear_iteration
Definition: simulator.h:1935
BoundaryVelocity::Manager< dim > boundary_velocity_manager
Definition: simulator.h:1908
const std::unique_ptr< AdiabaticConditions::Interface< dim > > adiabatic_conditions
Definition: simulator.h:1904
static const unsigned int invalid_unsigned_int
void write_plugin_graph(std::ostream &output_stream)
parallel::distributed::Triangulation< dim > triangulation
Definition: simulator.h:1962
TimerOutput computing_timer
Definition: simulator.h:1876
BoundaryTemperature::Manager< dim > boundary_temperature_manager
Definition: simulator.h:1899
boost::iostreams::tee_device< std::ostream, std::ofstream > TeeDevice
Definition: simulator.h:1836
double pressure_scaling
Definition: simulator.h:2010
Tensor< 1, dim > tensor_rotation
Definition: simulator.h:165
double old_time_step
Definition: simulator.h:1932
Tensor< 1, dim > tensor_angular_momentum
Definition: simulator.h:164
AffineConstraints< double > constraints
Definition: simulator.h:1996
bool assemble_newton_stokes_system
Definition: simulator.h:2075
LinearAlgebra::BlockVector old_solution
Definition: simulator.h:2055
const IntermediaryConstructorAction post_geometry_model_creation_action
Definition: simulator.h:1896
std::unique_ptr< VolumeOfFluidHandler< dim > > volume_of_fluid_handler
Definition: simulator.h:1823
bool rebuild_stokes_matrix
Definition: simulator.h:2073
Parameters< dim > parameters
Definition: simulator.h:1796
SymmetricTensor< 2, dim > tensor_moment_of_inertia
Definition: simulator.h:163
AffineConstraints< double > current_constraints
Definition: simulator.h:1997
std::map< types::boundary_id, std::unique_ptr< BoundaryTraction::Interface< dim > > > boundary_traction
Definition: simulator.h:1909
InitialComposition::Manager< dim > initial_composition_manager
Definition: simulator.h:1902
std::unique_ptr< MeshDeformation::MeshDeformationHandler< dim > > mesh_deformation
Definition: simulator.h:2089
LinearAlgebra::BlockVector system_rhs
Definition: simulator.h:2057
std::pair< double, double > stokes_residuals
Definition: simulator.h:133
LateralAveraging< dim > lateral_averaging
Definition: simulator.h:1953
TeeStream iostream_tee_stream
Definition: simulator.h:1840
MeshRefinement::Manager< dim > mesh_refinement_manager
Definition: simulator.h:1966
std::size_t statistics_last_write_size
Definition: simulator.h:1873
boost::iostreams::stream< TeeDevice > TeeStream
Definition: simulator.h:1837
const std::unique_ptr< MaterialModel::Interface< dim > > material_model
Definition: simulator.h:1897
#define DEAL_II_DISABLE_EXTRA_DIAGNOSTICS
::Particles::ParticleHandler< dim > particle_handler_copy
Definition: simulator.h:1921
double last_pressure_normalization_adjustment
Definition: simulator.h:2003
void declare_parameters(ParameterHandler &prm)
const std::unique_ptr< GeometryModel::Interface< dim > > geometry_model
Definition: simulator.h:1895
double initial_residual
LinearAlgebra::BlockVector operator_split_reaction_vector
Definition: simulator.h:2065
DoFHandler< dim > dof_handler
Definition: simulator.h:1981
std::unique_ptr< StokesMatrixFreeHandler< dim > > stokes_matrix_free
Definition: simulator.h:2094
std::unique_ptr< Assemblers::Manager< dim > > assemblers
Definition: simulator.h:1051
bool rebuild_sparsity_and_matrices
Definition: simulator.h:2072
const IntermediaryConstructorAction post_signal_creation
Definition: simulator.h:1814
double global_Omega_diameter
Definition: simulator.h:1963
unsigned int timestep_number
Definition: simulator.h:1933
bool assemble_newton_stokes_matrix
Definition: simulator.h:2074
ConditionalOStream pcout
Definition: simulator.h:1846
unsigned int pre_refinement_step
Definition: simulator.h:1934
TimeStepping::Manager< dim > time_stepping_manager
Definition: simulator.h:1944
std::size_t statistics_last_hash
Definition: simulator.h:1874
TableHandler statistics
Definition: simulator.h:1857
Introspection< dim > introspection
Definition: simulator.h:1825
MPI_Comm mpi_communicator
Definition: simulator.h:1828
typename Parameters< dim >::NullspaceRemoval NullspaceRemoval
Definition: simulator.h:247
LinearAlgebra::BlockSparseMatrix system_preconditioner_matrix
Definition: simulator.h:2052
bool do_pressure_rhs_compatibility_modification
Definition: simulator.h:2018
LinearAlgebra::BlockVector current_linearization_point
Definition: simulator.h:2059
InitialTemperature::Manager< dim > initial_temperature_manager
Definition: simulator.h:1903
LinearAlgebra::BlockVector old_old_solution
Definition: simulator.h:2056
LinearAlgebra::BlockVector pressure_shape_function_integrals
Definition: simulator.h:2062
std::ofstream log_file_stream
Definition: simulator.h:1834
std::unique_ptr< LinearAlgebra::PreconditionAMG > Amg_preconditioner
Definition: simulator.h:2069
TeeDevice iostream_tee_device
Definition: simulator.h:1839
LinearAlgebra::BlockSparseMatrix system_matrix
Definition: simulator.h:2038
#define DEAL_II_ENABLE_EXTRA_DIAGNOSTICS
T signaling_nan()
std::unique_ptr< MeltHandler< dim > > melt_handler
Definition: simulator.h:1803
bool rebuild_stokes_preconditioner
Definition: simulator.h:2076
const unsigned int compositional_variable
Definition: simulator.h:275
std::unique_ptr< Mapping< dim > > mapping
Definition: simulator.h:1977
const std::unique_ptr< PrescribedStokesSolution::Interface< dim > > prescribed_stokes_solution
Definition: simulator.h:1901
HeatingModel::Manager< dim > heating_model_manager
Definition: simulator.h:1967
SimulatorSignals< dim > signals
Definition: simulator.h:1812
const std::unique_ptr< BoundaryHeatFlux::Interface< dim > > boundary_heat_flux
Definition: simulator.h:1910
const FESystem< dim > finite_element
Definition: simulator.h:1979
bool simulator_is_past_initialization
Definition: simulator.h:393
Threads::Thread output_statistics_thread
Definition: simulator.h:1885
std::unique_ptr< Particle::World< dim > > particle_world
Definition: simulator.h:1915
const std::unique_ptr< InitialTopographyModel::Interface< dim > > initial_topography_model
Definition: simulator.h:1894
LinearAlgebra::BlockVector solution
Definition: simulator.h:2054
std::unique_ptr< NewtonHandler< dim > > newton_handler
Definition: simulator.h:1810
void serialize(Archive &archive, const unsigned int version)
typename Parameters< dim >::NonlinearSolver NonlinearSolver
Definition: simulator.h:242
Postprocess::Manager< dim > postprocess_manager
Definition: simulator.h:1983
const std::unique_ptr< GravityModel::Interface< dim > > gravity_model
Definition: simulator.h:1898
double newton_residual_for_derivative_scaling_factor
Definition: simulator.h:132
double global_volume
Definition: simulator.h:1964
BoundaryComposition::Manager< dim > boundary_composition_manager
Definition: simulator.h:1900
void run(const Iterator &begin, const typename identity< Iterator >::type &end, Worker worker, Copier copier, const ScratchData &sample_scratch_data, const CopyData &sample_copy_data, const unsigned int queue_length, const unsigned int chunk_size)
std::unique_ptr< LinearAlgebra::PreconditionBase > Mp_preconditioner
Definition: simulator.h:2070