ASPECT
simulator.h
Go to the documentation of this file.
1 /*
2  Copyright (C) 2011 - 2019 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>
26 #include <deal.II/base/parameter_handler.h>
27 #include <deal.II/base/conditional_ostream.h>
29 
30 DEAL_II_DISABLE_EXTRA_DIAGNOSTICS
31 
32 #include <deal.II/distributed/tria.h>
33 
34 #include <deal.II/dofs/dof_handler.h>
35 #include <deal.II/dofs/dof_tools.h>
36 
37 #include <deal.II/fe/fe_system.h>
38 #include <deal.II/fe/mapping.h>
39 #include <deal.II/base/tensor_function.h>
40 
41 DEAL_II_ENABLE_EXTRA_DIAGNOSTICS
42 
43 #if !DEAL_II_VERSION_GTE(9,1,0)
44 # include <deal.II/lac/constraint_matrix.h>
45 #else
46 # include <deal.II/lac/affine_constraints.h>
47 #endif
48 
49 #include <aspect/global.h>
71 #include <aspect/particle/world.h>
72 
73 #include <boost/iostreams/tee.hpp>
74 #include <boost/iostreams/stream.hpp>
75 #include <memory>
76 
77 namespace aspect
78 {
79  using namespace dealii;
80 
81 #if DEAL_II_VERSION_GTE(9,1,0)
82 
90  using ConstraintMatrix = class ::AffineConstraints<double>;
91 #endif
92 
93  template <int dim>
94  class MeltHandler;
95 
96  template <int dim>
97  class NewtonHandler;
98 
99  template <int dim>
101 
102  template <int dim, int velocity_degree>
104 
105  namespace MeshDeformation
106  {
107  template <int dim>
109  }
110 
111  template <int dim>
113 
114  namespace internal
115  {
116  namespace Assembly
117  {
118  namespace Scratch
119  {
120  template <int dim> struct StokesPreconditioner;
121  template <int dim> struct StokesSystem;
122  template <int dim> struct AdvectionSystem;
123  }
124 
125  namespace CopyData
126  {
127  template <int dim> struct StokesPreconditioner;
128  template <int dim> struct StokesSystem;
129  template <int dim> struct AdvectionSystem;
130  }
131  }
132  }
133 
134  namespace Assemblers
135  {
136  template <int dim> class Interface;
137  template <int dim> class Manager;
138  }
139 
141  {
145  double residual;
146  double residual_old;
149  std::pair<double,double> stokes_residuals;
150  };
151 
159  template <int dim>
160  class Simulator
161  {
162  public:
176  Simulator (const MPI_Comm mpi_communicator,
177  ParameterHandler &prm);
178 
183  ~Simulator ();
184 
196  static
198 
207  void run ();
208 
220  void
221  write_plugin_graph (std::ostream &output_stream) const;
222 
227 
232 
233 
241  {
246  enum FieldType { temperature_field, compositional_field };
247 
253 
259  const unsigned int compositional_variable;
260 
272  AdvectionField (const FieldType field_type,
273  const unsigned int compositional_variable = numbers::invalid_unsigned_int);
274 
282  static
284 
292  static
293  AdvectionField composition (const unsigned int compositional_variable);
294 
298  bool
299  is_temperature () const;
300 
305  bool
306  is_discontinuous (const Introspection<dim> &introspection) const;
307 
313  advection_method (const Introspection<dim> &introspection) const;
314 
319  unsigned int component_index(const Introspection<dim> &introspection) const;
320 
325  unsigned int block_index(const Introspection<dim> &introspection) const;
326 
333  unsigned int field_index() const;
334 
340  unsigned int base_element(const Introspection<dim> &introspection) const;
341 
348  FEValuesExtractors::Scalar scalar_extractor(const Introspection<dim> &introspection) const;
349 
354  unsigned int polynomial_degree(const Introspection<dim> &introspection) const;
355  };
356 
357  private:
358 
378 
396  {
397  IntermediaryConstructorAction (const std::function<void ()> &action);
398  };
399 
414  void setup_dofs ();
415 
424  void setup_introspection ();
425 
437  void set_initial_temperature_and_compositional_fields ();
438 
454  void compute_initial_pressure_field ();
455 
461  void compute_initial_velocity_boundary_constraints (ConstraintMatrix &constraints);
462 
468  void compute_current_velocity_boundary_constraints (ConstraintMatrix &constraints);
469 
482  void compute_current_constraints ();
483 
496  void compute_pressure_scaling_factor ();
497 
509  void start_timestep ();
510 
518  void solve_timestep ();
519 
531  void solve_single_advection_single_stokes ();
532 
545  void solve_no_advection_iterated_stokes ();
546 
557  void solve_no_advection_single_stokes ();
558 
570  void solve_first_timestep_only_single_stokes ();
571 
584  void solve_iterated_advection_and_stokes ();
585 
598  void solve_single_advection_iterated_stokes ();
599 
614  void solve_iterated_advection_and_newton_stokes ();
615 
631  void solve_single_advection_iterated_newton_stokes ();
632 
644  void solve_single_advection_no_stokes ();
645 
657  void solve_no_advection_no_stokes ();
658 
667  void build_stokes_preconditioner ();
668 
676  void build_advection_preconditioner (const AdvectionField &advection_field,
678 
685  void assemble_stokes_system ();
686 
696  double assemble_and_solve_temperature (const bool compute_initial_residual = false,
697  double *initial_residual = nullptr);
698 
711  std::vector<double> assemble_and_solve_composition (const bool compute_initial_residual = false,
712  std::vector<double> *initial_residual = nullptr);
713 
731  double assemble_and_solve_stokes (const bool compute_initial_residual = false,
732  double *initial_nonlinear_residual = nullptr);
733 
745  void assemble_and_solve_defect_correction_Stokes(DefectCorrectionResiduals &dcr,
746  bool use_picard);
747 
755  void assemble_advection_system (const AdvectionField &advection_field);
756 
767  double solve_advection (const AdvectionField &advection_field);
768 
772  void interpolate_particle_properties (const AdvectionField &advection_field);
773 
844  std::pair<double,double>
845  solve_stokes ();
846 
850  std::pair<double,double>
851  solve_stokes_block_gmg ();
852 
865  void postprocess ();
866 
882  void refine_mesh (const unsigned int max_grid_level);
883 
904  void create_snapshot();
905 
919  void resume_from_snapshot();
920 
927  template <class Archive>
928  void serialize (Archive &ar, const unsigned int version);
946  setup_system_matrix_coupling () const;
947 
955  void setup_system_matrix (const std::vector<IndexSet> &system_partitioning);
956 
968  void setup_system_preconditioner (const std::vector<IndexSet> &system_partitioning);
969 
991  std::unique_ptr<Assemblers::Manager<dim> > assemblers;
992 
1001  void set_assemblers ();
1002 
1013  void set_default_assemblers ();
1014 
1021  void assemble_stokes_preconditioner ();
1022 
1030  void
1031  local_assemble_stokes_preconditioner (const typename DoFHandler<dim>::active_cell_iterator &cell,
1034 
1042  void
1043  copy_local_to_global_stokes_preconditioner (const internal::Assembly::CopyData::StokesPreconditioner<dim> &data);
1044 
1052  void
1053  local_assemble_stokes_system (const typename DoFHandler<dim>::active_cell_iterator &cell,
1056 
1064  void
1065  copy_local_to_global_stokes_system (const internal::Assembly::CopyData::StokesSystem<dim> &data);
1066 
1074  void
1075  local_assemble_advection_face_terms(const AdvectionField &advection_field,
1076  const typename DoFHandler<dim>::active_cell_iterator &cell,
1086  void
1087  local_assemble_advection_system (const AdvectionField &advection_field,
1088  const Vector<double> &viscosity_per_cell,
1089  const typename DoFHandler<dim>::active_cell_iterator &cell,
1092 
1100  void
1101  copy_local_to_global_advection_system (const AdvectionField &advection_field,
1103 
1128  void make_pressure_rhs_compatible(LinearAlgebra::BlockVector &vector);
1129 
1139  template <typename T>
1140  void get_artificial_viscosity (Vector<T> &viscosity_per_cell,
1141  const AdvectionField &advection_field,
1142  const bool skip_interior_cells = false) const;
1143 
1155  void compute_Vs_anomaly(Vector<float> &values) const;
1156 
1171  void compute_Vp_anomaly(Vector<float> &values) const;
1172 
1213  double normalize_pressure(LinearAlgebra::BlockVector &vector) const;
1214 
1259  void denormalize_pressure(const double pressure_adjustment,
1261  const LinearAlgebra::BlockVector &relevant_vector) const;
1262 
1271  void apply_limiter_to_dg_solutions (const AdvectionField &advection_field);
1272 
1273 
1296  void compute_reactions ();
1297 
1307  void update_solution_vectors_with_reaction_results (const unsigned int block_index,
1308  const LinearAlgebra::BlockVector &distributed_vector,
1309  const LinearAlgebra::BlockVector &distributed_reaction_vector);
1310 
1321  void initialize_current_linearization_point ();
1322 
1342  void interpolate_material_output_into_advection_field (const AdvectionField &adv_field);
1343 
1344 
1352  void interpolate_onto_velocity_system(const TensorFunction<1,dim> &func,
1353  LinearAlgebra::Vector &vec);
1354 
1355 
1370  void setup_nullspace_constraints(ConstraintMatrix &constraints);
1371 
1372 
1384  void remove_nullspace(LinearAlgebra::BlockVector &relevant_dst,
1385  LinearAlgebra::BlockVector &tmp_distributed_stokes);
1386 
1400  void remove_net_angular_momentum( const bool use_constant_density,
1401  LinearAlgebra::BlockVector &relevant_dst,
1402  LinearAlgebra::BlockVector &tmp_distributed_stokes);
1403 
1411  void replace_outflow_boundary_ids(const unsigned int boundary_id_offset);
1412 
1420  void restore_outflow_boundary_ids(const unsigned int boundary_id_offset);
1421 
1435  void remove_net_linear_momentum( const bool use_constant_density,
1436  LinearAlgebra::BlockVector &relevant_dst,
1437  LinearAlgebra::BlockVector &tmp_distributed_stokes);
1438 
1446  double get_maximal_velocity (const LinearAlgebra::BlockVector &solution) const;
1447 
1460  double get_entropy_variation (const double average_field,
1461  const AdvectionField &advection_field) const;
1462 
1471  std::pair<double,double>
1472  get_extrapolated_advection_field_range (const AdvectionField &advection_field) const;
1473 
1482  void maybe_write_timing_output () const;
1483 
1491  bool maybe_write_checkpoint (const time_t last_checkpoint_time,
1492  const std::pair<bool,bool> termination_output);
1493 
1510  bool maybe_do_initial_refinement (const unsigned int max_refinement_level);
1511 
1520  void maybe_refine_mesh (const double new_time_step,
1521  unsigned int &max_refinement_level);
1522 
1533  double compute_time_step () const;
1534 
1542  void advance_time (const double step_size);
1543 
1551  double
1552  compute_viscosity(internal::Assembly::Scratch::AdvectionSystem<dim> &scratch,
1553  const double global_u_infty,
1554  const double global_field_variation,
1555  const double average_field,
1556  const double global_entropy_variation,
1557  const double cell_diameter,
1558  const AdvectionField &advection_field) const;
1559 
1568  void
1569  compute_advection_system_residual(internal::Assembly::Scratch::AdvectionSystem<dim> &scratch,
1570  const double average_field,
1571  const AdvectionField &advection_field,
1572  double &max_residual,
1573  double &max_velocity,
1574  double &max_density,
1575  double &max_specific_heat,
1576  double &conductivity) const;
1577 
1605  void
1606  compute_material_model_input_values (const LinearAlgebra::BlockVector &input_solution,
1607  const FEValuesBase<dim,dim> &input_finite_element_values,
1608  const typename DoFHandler<dim>::active_cell_iterator &cell,
1609  const bool compute_strainrate,
1610  MaterialModel::MaterialModelInputs<dim> &material_model_inputs) const;
1611 
1612 
1627  bool
1628  stokes_matrix_depends_on_solution () const;
1629 
1642  void
1643  check_consistency_of_formulation ();
1644 
1655  void
1656  check_consistency_of_boundary_conditions () const;
1657 
1661  double
1662  compute_initial_newton_residual (const LinearAlgebra::BlockVector &linearized_stokes_initial_guess);
1663 
1672  double
1673  compute_Eisenstat_Walker_linear_tolerance(const bool EisenstatWalkerChoiceOne,
1674  const double maximum_linear_stokes_solver_tolerance,
1675  const double linear_stokes_solver_tolerance,
1676  const double stokes_residual,
1677  const double newton_residual,
1678  const double newton_residual_old);
1679 
1689  void output_statistics();
1690 
1702  double
1703  compute_initial_stokes_residual();
1704 
1715 
1721  std::unique_ptr<MeltHandler<dim> > melt_handler;
1722 
1728  std::unique_ptr<NewtonHandler<dim> > newton_handler;
1729 
1731 
1733 
1741  std::unique_ptr<VolumeOfFluidHandler<dim> > volume_of_fluid_handler;
1742 
1744 
1745 
1747 
1752  std::ofstream log_file_stream;
1753 
1754  typedef boost::iostreams::tee_device<std::ostream, std::ofstream> TeeDevice;
1755  typedef boost::iostreams::stream< TeeDevice > TeeStream;
1756 
1759 
1765 
1776 
1793 
1795 
1812  const std::unique_ptr<InitialTopographyModel::Interface<dim> > initial_topography_model;
1813  const std::unique_ptr<GeometryModel::Interface<dim> > geometry_model;
1815  const std::unique_ptr<MaterialModel::Interface<dim> > material_model;
1816  const std::unique_ptr<GravityModel::Interface<dim> > gravity_model;
1819  const std::unique_ptr<PrescribedStokesSolution::Interface<dim> > prescribed_stokes_solution;
1822  const std::unique_ptr<AdiabaticConditions::Interface<dim> > adiabatic_conditions;
1823 #ifdef ASPECT_WITH_WORLD_BUILDER
1824  const std::unique_ptr<WorldBuilder::World> world_builder;
1825 #endif
1827  std::map<types::boundary_id,std::unique_ptr<BoundaryTraction::Interface<dim> > > boundary_traction;
1828  const std::unique_ptr<BoundaryHeatFlux::Interface<dim> > boundary_heat_flux;
1829 
1833  std::unique_ptr<Particle::World<dim> > particle_world;
1834 
1842  double time;
1843  double time_step;
1845  unsigned int timestep_number;
1846  unsigned int pre_refinement_step;
1847  unsigned int nonlinear_iteration;
1877 
1880 
1889  std::unique_ptr<Mapping<dim> > mapping;
1890 
1892 
1894 
1896 
1908  ConstraintMatrix constraints;
1909  ConstraintMatrix current_constraints;
1910 
1916 
1923 
1931 
1951 
1965 
1970 
1972 
1973  // only used if is_compressible()
1975 
1976  // only used if operator split is enabled
1978 
1979 
1980 
1981  std::unique_ptr<LinearAlgebra::PreconditionAMG> Amg_preconditioner;
1982  std::unique_ptr<LinearAlgebra::PreconditionBase> Mp_preconditioner;
1983 
1989 
1994  private:
1995 
2001  std::unique_ptr<MeshDeformation::MeshDeformationHandler<dim> > mesh_deformation;
2002 
2006  std::unique_ptr<StokesMatrixFreeHandler<dim> > stokes_matrix_free;
2007 
2008  friend class boost::serialization::access;
2009  friend class SimulatorAccess<dim>;
2010  friend class MeshDeformation::MeshDeformationHandler<dim>; // MeshDeformationHandler needs access to the internals of the Simulator
2011  friend class VolumeOfFluidHandler<dim>; // VolumeOfFluidHandler needs access to the internals of the Simulator
2012  friend class StokesMatrixFreeHandler<dim>;
2013  template <int dimension, int velocity_degree>
2015  friend struct Parameters<dim>;
2016  };
2017 }
2018 
2019 
2020 #endif
The NullspaceRemoval struct.
Definition: parameters.h:83
unsigned int nonlinear_iteration
Definition: simulator.h:1847
BoundaryVelocity::Manager< dim > boundary_velocity_manager
Definition: simulator.h:1826
const std::unique_ptr< AdiabaticConditions::Interface< dim > > adiabatic_conditions
Definition: simulator.h:1822
static const unsigned int invalid_unsigned_int
TerminationCriteria::Manager< dim > termination_manager
Definition: simulator.h:1856
void write_plugin_graph(std::ostream &output_stream)
parallel::distributed::Triangulation< dim > triangulation
Definition: simulator.h:1874
TimerOutput computing_timer
Definition: simulator.h:1794
BoundaryTemperature::Manager< dim > boundary_temperature_manager
Definition: simulator.h:1817
double pressure_scaling
Definition: simulator.h:1922
double old_time_step
Definition: simulator.h:1844
bool assemble_newton_stokes_system
Definition: simulator.h:1987
LinearAlgebra::BlockVector old_solution
Definition: simulator.h:1967
const IntermediaryConstructorAction post_geometry_model_creation_action
Definition: simulator.h:1814
std::unique_ptr< VolumeOfFluidHandler< dim > > volume_of_fluid_handler
Definition: simulator.h:1741
bool rebuild_stokes_matrix
Definition: simulator.h:1985
Parameters< dim > parameters
Definition: simulator.h:1714
std::map< types::boundary_id, std::unique_ptr< BoundaryTraction::Interface< dim > > > boundary_traction
Definition: simulator.h:1827
boost::iostreams::tee_device< std::ostream, std::ofstream > TeeDevice
Definition: simulator.h:1754
InitialComposition::Manager< dim > initial_composition_manager
Definition: simulator.h:1820
std::unique_ptr< MeshDeformation::MeshDeformationHandler< dim > > mesh_deformation
Definition: simulator.h:2001
LinearAlgebra::BlockVector system_rhs
Definition: simulator.h:1969
std::pair< double, double > stokes_residuals
Definition: simulator.h:149
boost::iostreams::stream< TeeDevice > TeeStream
Definition: simulator.h:1755
LateralAveraging< dim > lateral_averaging
Definition: simulator.h:1865
TeeStream iostream_tee_stream
Definition: simulator.h:1758
ConstraintMatrix constraints
Definition: simulator.h:1908
MeshRefinement::Manager< dim > mesh_refinement_manager
Definition: simulator.h:1878
std::size_t statistics_last_write_size
Definition: simulator.h:1791
const std::unique_ptr< MaterialModel::Interface< dim > > material_model
Definition: simulator.h:1815
double last_pressure_normalization_adjustment
Definition: simulator.h:1915
void declare_parameters(ParameterHandler &prm)
const std::unique_ptr< GeometryModel::Interface< dim > > geometry_model
Definition: simulator.h:1813
LinearAlgebra::BlockVector operator_split_reaction_vector
Definition: simulator.h:1977
Parameters< dim >::NonlinearSolver NonlinearSolver
Definition: simulator.h:226
DoFHandler< dim > dof_handler
Definition: simulator.h:1893
std::unique_ptr< StokesMatrixFreeHandler< dim > > stokes_matrix_free
Definition: simulator.h:2006
std::unique_ptr< Assemblers::Manager< dim > > assemblers
Definition: simulator.h:991
bool rebuild_sparsity_and_matrices
Definition: simulator.h:1984
const IntermediaryConstructorAction post_signal_creation
Definition: simulator.h:1732
double global_Omega_diameter
Definition: simulator.h:1875
unsigned int timestep_number
Definition: simulator.h:1845
bool assemble_newton_stokes_matrix
Definition: simulator.h:1986
ConditionalOStream pcout
Definition: simulator.h:1764
unsigned int pre_refinement_step
Definition: simulator.h:1846
std::size_t statistics_last_hash
Definition: simulator.h:1792
TableHandler statistics
Definition: simulator.h:1775
Introspection< dim > introspection
Definition: simulator.h:1743
MPI_Comm mpi_communicator
Definition: simulator.h:1746
LinearAlgebra::BlockSparseMatrix system_preconditioner_matrix
Definition: simulator.h:1964
bool do_pressure_rhs_compatibility_modification
Definition: simulator.h:1930
LinearAlgebra::BlockVector current_linearization_point
Definition: simulator.h:1971
InitialTemperature::Manager< dim > initial_temperature_manager
Definition: simulator.h:1821
LinearAlgebra::BlockVector old_old_solution
Definition: simulator.h:1968
Definition: compat.h:53
LinearAlgebra::BlockVector pressure_shape_function_integrals
Definition: simulator.h:1974
std::ofstream log_file_stream
Definition: simulator.h:1752
std::unique_ptr< LinearAlgebra::PreconditionAMG > Amg_preconditioner
Definition: simulator.h:1981
ConstraintMatrix current_constraints
Definition: simulator.h:1909
TeeDevice iostream_tee_device
Definition: simulator.h:1757
LinearAlgebra::BlockSparseMatrix system_matrix
Definition: simulator.h:1950
std::unique_ptr< MeltHandler< dim > > melt_handler
Definition: simulator.h:1721
bool rebuild_stokes_preconditioner
Definition: simulator.h:1988
const unsigned int compositional_variable
Definition: simulator.h:259
std::unique_ptr< Mapping< dim > > mapping
Definition: simulator.h:1889
const std::unique_ptr< PrescribedStokesSolution::Interface< dim > > prescribed_stokes_solution
Definition: simulator.h:1819
HeatingModel::Manager< dim > heating_model_manager
Definition: simulator.h:1879
SimulatorSignals< dim > signals
Definition: simulator.h:1730
const std::unique_ptr< BoundaryHeatFlux::Interface< dim > > boundary_heat_flux
Definition: simulator.h:1828
const FESystem< dim > finite_element
Definition: simulator.h:1891
bool simulator_is_past_initialization
Definition: simulator.h:377
Threads::Thread output_statistics_thread
Definition: simulator.h:1803
std::unique_ptr< Particle::World< dim > > particle_world
Definition: simulator.h:1833
const std::unique_ptr< InitialTopographyModel::Interface< dim > > initial_topography_model
Definition: simulator.h:1812
LinearAlgebra::BlockVector solution
Definition: simulator.h:1966
std::unique_ptr< NewtonHandler< dim > > newton_handler
Definition: simulator.h:1728
Postprocess::Manager< dim > postprocess_manager
Definition: simulator.h:1895
const std::unique_ptr< GravityModel::Interface< dim > > gravity_model
Definition: simulator.h:1816
double newton_residual_for_derivative_scaling_factor
Definition: simulator.h:148
double global_volume
Definition: simulator.h:1876
BoundaryComposition::Manager< dim > boundary_composition_manager
Definition: simulator.h:1818
std::unique_ptr< LinearAlgebra::PreconditionBase > Mp_preconditioner
Definition: simulator.h:1982
Parameters< dim >::NullspaceRemoval NullspaceRemoval
Definition: simulator.h:231