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 
36 #include <deal.II/fe/fe_system.h>
37 #include <deal.II/fe/mapping.h>
38 #include <deal.II/base/tensor_function.h>
39 
40 DEAL_II_ENABLE_EXTRA_DIAGNOSTICS
41 
42 #if !DEAL_II_VERSION_GTE(9,1,0)
43 # include <deal.II/lac/constraint_matrix.h>
44 #else
45 # include <deal.II/lac/affine_constraints.h>
46 #endif
47 
48 #include <aspect/global.h>
70 #include <aspect/particle/world.h>
71 
72 #include <boost/iostreams/tee.hpp>
73 #include <boost/iostreams/stream.hpp>
74 #include <memory>
75 
76 namespace aspect
77 {
78  using namespace dealii;
79 
80 #if DEAL_II_VERSION_GTE(9,1,0)
81 
89  using ConstraintMatrix = class ::AffineConstraints<double>;
90 #endif
91 
92  template <int dim>
93  class MeltHandler;
94 
95  template <int dim>
96  class NewtonHandler;
97 
98  template <int dim>
100 
101  template <int dim, int velocity_degree>
103 
104  namespace MeshDeformation
105  {
106  template <int dim>
108  }
109 
110  template <int dim>
112 
113  namespace internal
114  {
115  namespace Assembly
116  {
117  namespace Scratch
118  {
119  template <int dim> struct StokesPreconditioner;
120  template <int dim> struct StokesSystem;
121  template <int dim> struct AdvectionSystem;
122  }
123 
124  namespace CopyData
125  {
126  template <int dim> struct StokesPreconditioner;
127  template <int dim> struct StokesSystem;
128  template <int dim> struct AdvectionSystem;
129  }
130  }
131  }
132 
133  namespace Assemblers
134  {
135  template <int dim> class Interface;
136  template <int dim> class Manager;
137  }
138 
140  {
144  double residual;
145  double residual_old;
148  std::pair<double,double> stokes_residuals;
149  };
150 
158  template <int dim>
159  class Simulator
160  {
161  public:
175  Simulator (const MPI_Comm mpi_communicator,
176  ParameterHandler &prm);
177 
182  ~Simulator ();
183 
195  static
197 
206  void run ();
207 
219  void
220  write_plugin_graph (std::ostream &output_stream) const;
221 
226 
231 
232 
240  {
245  enum FieldType { temperature_field, compositional_field };
246 
252 
258  const unsigned int compositional_variable;
259 
271  AdvectionField (const FieldType field_type,
272  const unsigned int compositional_variable = numbers::invalid_unsigned_int);
273 
281  static
283 
291  static
292  AdvectionField composition (const unsigned int compositional_variable);
293 
297  bool
298  is_temperature () const;
299 
304  bool
305  is_discontinuous (const Introspection<dim> &introspection) const;
306 
312  advection_method (const Introspection<dim> &introspection) const;
313 
318  unsigned int component_index(const Introspection<dim> &introspection) const;
319 
324  unsigned int block_index(const Introspection<dim> &introspection) const;
325 
332  unsigned int field_index() const;
333 
339  unsigned int base_element(const Introspection<dim> &introspection) const;
340 
347  FEValuesExtractors::Scalar scalar_extractor(const Introspection<dim> &introspection) const;
348 
353  unsigned int polynomial_degree(const Introspection<dim> &introspection) const;
354  };
355 
356  private:
357 
377 
395  {
396  IntermediaryConstructorAction (const std::function<void ()> &action);
397  };
398 
413  void setup_dofs ();
414 
423  void setup_introspection ();
424 
436  void set_initial_temperature_and_compositional_fields ();
437 
453  void compute_initial_pressure_field ();
454 
460  void compute_initial_velocity_boundary_constraints (ConstraintMatrix &constraints);
461 
467  void compute_current_velocity_boundary_constraints (ConstraintMatrix &constraints);
468 
481  void compute_current_constraints ();
482 
495  void compute_pressure_scaling_factor ();
496 
508  void start_timestep ();
509 
517  void solve_timestep ();
518 
530  void solve_single_advection_single_stokes ();
531 
544  void solve_no_advection_iterated_stokes ();
545 
556  void solve_no_advection_single_stokes ();
557 
569  void solve_first_timestep_only_single_stokes ();
570 
583  void solve_iterated_advection_and_stokes ();
584 
597  void solve_single_advection_iterated_stokes ();
598 
613  void solve_iterated_advection_and_newton_stokes ();
614 
630  void solve_single_advection_iterated_newton_stokes ();
631 
643  void solve_single_advection_no_stokes ();
644 
656  void solve_no_advection_no_stokes ();
657 
666  void build_stokes_preconditioner ();
667 
675  void build_advection_preconditioner (const AdvectionField &advection_field,
677 
684  void assemble_stokes_system ();
685 
695  double assemble_and_solve_temperature (const bool compute_initial_residual = false,
696  double *initial_residual = nullptr);
697 
710  std::vector<double> assemble_and_solve_composition (const bool compute_initial_residual = false,
711  std::vector<double> *initial_residual = nullptr);
712 
730  double assemble_and_solve_stokes (const bool compute_initial_residual = false,
731  double *initial_nonlinear_residual = nullptr);
732 
744  void assemble_and_solve_defect_correction_Stokes(DefectCorrectionResiduals &dcr,
745  bool use_picard);
746 
754  void assemble_advection_system (const AdvectionField &advection_field);
755 
766  double solve_advection (const AdvectionField &advection_field);
767 
771  void interpolate_particle_properties (const AdvectionField &advection_field);
772 
843  std::pair<double,double>
844  solve_stokes ();
845 
849  std::pair<double,double>
850  solve_stokes_block_gmg ();
851 
864  void postprocess ();
865 
881  void refine_mesh (const unsigned int max_grid_level);
882 
903  void create_snapshot();
904 
918  void resume_from_snapshot();
919 
926  template <class Archive>
927  void serialize (Archive &ar, const unsigned int version);
943  void setup_system_matrix (const std::vector<IndexSet> &system_partitioning);
944 
956  void setup_system_preconditioner (const std::vector<IndexSet> &system_partitioning);
957 
979  std::unique_ptr<Assemblers::Manager<dim> > assemblers;
980 
989  void set_assemblers ();
990 
1001  void set_default_assemblers ();
1002 
1009  void assemble_stokes_preconditioner ();
1010 
1018  void
1019  local_assemble_stokes_preconditioner (const typename DoFHandler<dim>::active_cell_iterator &cell,
1022 
1030  void
1031  copy_local_to_global_stokes_preconditioner (const internal::Assembly::CopyData::StokesPreconditioner<dim> &data);
1032 
1040  void
1041  local_assemble_stokes_system (const typename DoFHandler<dim>::active_cell_iterator &cell,
1044 
1052  void
1053  copy_local_to_global_stokes_system (const internal::Assembly::CopyData::StokesSystem<dim> &data);
1054 
1062  void
1063  local_assemble_advection_face_terms(const AdvectionField &advection_field,
1064  const typename DoFHandler<dim>::active_cell_iterator &cell,
1074  void
1075  local_assemble_advection_system (const AdvectionField &advection_field,
1076  const Vector<double> &viscosity_per_cell,
1077  const typename DoFHandler<dim>::active_cell_iterator &cell,
1080 
1088  void
1089  copy_local_to_global_advection_system (const AdvectionField &advection_field,
1091 
1116  void make_pressure_rhs_compatible(LinearAlgebra::BlockVector &vector);
1117 
1127  template <typename T>
1128  void get_artificial_viscosity (Vector<T> &viscosity_per_cell,
1129  const AdvectionField &advection_field,
1130  const bool skip_interior_cells = false) const;
1131 
1143  void compute_Vs_anomaly(Vector<float> &values) const;
1144 
1159  void compute_Vp_anomaly(Vector<float> &values) const;
1160 
1201  double normalize_pressure(LinearAlgebra::BlockVector &vector) const;
1202 
1247  void denormalize_pressure(const double pressure_adjustment,
1249  const LinearAlgebra::BlockVector &relevant_vector) const;
1250 
1259  void apply_limiter_to_dg_solutions (const AdvectionField &advection_field);
1260 
1261 
1284  void compute_reactions ();
1285 
1296  void initialize_current_linearization_point ();
1297 
1298 
1318  void interpolate_material_output_into_advection_field (const AdvectionField &adv_field);
1319 
1320 
1328  void interpolate_onto_velocity_system(const TensorFunction<1,dim> &func,
1329  LinearAlgebra::Vector &vec);
1330 
1331 
1346  void setup_nullspace_constraints(ConstraintMatrix &constraints);
1347 
1348 
1360  void remove_nullspace(LinearAlgebra::BlockVector &relevant_dst,
1361  LinearAlgebra::BlockVector &tmp_distributed_stokes);
1362 
1376  void remove_net_angular_momentum( const bool use_constant_density,
1377  LinearAlgebra::BlockVector &relevant_dst,
1378  LinearAlgebra::BlockVector &tmp_distributed_stokes);
1379 
1387  void replace_outflow_boundary_ids(const unsigned int boundary_id_offset);
1388 
1396  void restore_outflow_boundary_ids(const unsigned int boundary_id_offset);
1397 
1411  void remove_net_linear_momentum( const bool use_constant_density,
1412  LinearAlgebra::BlockVector &relevant_dst,
1413  LinearAlgebra::BlockVector &tmp_distributed_stokes);
1414 
1422  double get_maximal_velocity (const LinearAlgebra::BlockVector &solution) const;
1423 
1436  double get_entropy_variation (const double average_field,
1437  const AdvectionField &advection_field) const;
1438 
1447  std::pair<double,double>
1448  get_extrapolated_advection_field_range (const AdvectionField &advection_field) const;
1449 
1458  void maybe_write_timing_output () const;
1459 
1467  bool maybe_write_checkpoint (const time_t last_checkpoint_time,
1468  const std::pair<bool,bool> termination_output);
1469 
1486  bool maybe_do_initial_refinement (const unsigned int max_refinement_level);
1487 
1496  void maybe_refine_mesh (const double new_time_step,
1497  unsigned int &max_refinement_level);
1498 
1509  double compute_time_step () const;
1510 
1518  void advance_time (const double step_size);
1519 
1527  double
1528  compute_viscosity(internal::Assembly::Scratch::AdvectionSystem<dim> &scratch,
1529  const double global_u_infty,
1530  const double global_field_variation,
1531  const double average_field,
1532  const double global_entropy_variation,
1533  const double cell_diameter,
1534  const AdvectionField &advection_field) const;
1535 
1544  void
1545  compute_advection_system_residual(internal::Assembly::Scratch::AdvectionSystem<dim> &scratch,
1546  const double average_field,
1547  const AdvectionField &advection_field,
1548  double &max_residual,
1549  double &max_velocity,
1550  double &max_density,
1551  double &max_specific_heat,
1552  double &conductivity) const;
1553 
1581  void
1582  compute_material_model_input_values (const LinearAlgebra::BlockVector &input_solution,
1583  const FEValuesBase<dim,dim> &input_finite_element_values,
1584  const typename DoFHandler<dim>::active_cell_iterator &cell,
1585  const bool compute_strainrate,
1586  MaterialModel::MaterialModelInputs<dim> &material_model_inputs) const;
1587 
1588 
1603  bool
1604  stokes_matrix_depends_on_solution () const;
1605 
1618  void
1619  check_consistency_of_formulation ();
1620 
1631  void
1632  check_consistency_of_boundary_conditions () const;
1633 
1637  double
1638  compute_initial_newton_residual (const LinearAlgebra::BlockVector &linearized_stokes_initial_guess);
1639 
1648  double
1649  compute_Eisenstat_Walker_linear_tolerance(const bool EisenstatWalkerChoiceOne,
1650  const double maximum_linear_stokes_solver_tolerance,
1651  const double linear_stokes_solver_tolerance,
1652  const double stokes_residual,
1653  const double newton_residual,
1654  const double newton_residual_old);
1655 
1665  void output_statistics();
1666 
1678  double
1679  compute_initial_stokes_residual();
1680 
1691 
1697  std::unique_ptr<MeltHandler<dim> > melt_handler;
1698 
1704  std::unique_ptr<NewtonHandler<dim> > newton_handler;
1705 
1707 
1709 
1717  std::unique_ptr<VolumeOfFluidHandler<dim> > volume_of_fluid_handler;
1718 
1720 
1721 
1723 
1728  std::ofstream log_file_stream;
1729 
1730  typedef boost::iostreams::tee_device<std::ostream, std::ofstream> TeeDevice;
1731  typedef boost::iostreams::stream< TeeDevice > TeeStream;
1732 
1735 
1741 
1752 
1754 
1771  const std::unique_ptr<InitialTopographyModel::Interface<dim> > initial_topography_model;
1772  const std::unique_ptr<GeometryModel::Interface<dim> > geometry_model;
1774  const std::unique_ptr<MaterialModel::Interface<dim> > material_model;
1775  const std::unique_ptr<GravityModel::Interface<dim> > gravity_model;
1778  const std::unique_ptr<PrescribedStokesSolution::Interface<dim> > prescribed_stokes_solution;
1781  const std::unique_ptr<AdiabaticConditions::Interface<dim> > adiabatic_conditions;
1782 #ifdef ASPECT_WITH_WORLD_BUILDER
1783  const std::unique_ptr<WorldBuilder::World> world_builder;
1784 #endif
1786  std::map<types::boundary_id,std::unique_ptr<BoundaryTraction::Interface<dim> > > boundary_traction;
1787  const std::unique_ptr<BoundaryHeatFlux::Interface<dim> > boundary_heat_flux;
1788 
1792  std::unique_ptr<Particle::World<dim> > particle_world;
1793 
1801  double time;
1802  double time_step;
1804  unsigned int timestep_number;
1805  unsigned int pre_refinement_step;
1806  unsigned int nonlinear_iteration;
1836 
1839 
1848  std::unique_ptr<Mapping<dim> > mapping;
1849 
1851 
1853 
1855 
1867  ConstraintMatrix constraints;
1868  ConstraintMatrix current_constraints;
1869 
1875 
1882 
1890 
1910 
1924 
1929 
1931 
1932  // only used if is_compressible()
1934 
1935  // only used if operator split is enabled
1937 
1938 
1939 
1940  std::unique_ptr<LinearAlgebra::PreconditionAMG> Amg_preconditioner;
1941  std::unique_ptr<LinearAlgebra::PreconditionBase> Mp_preconditioner;
1942 
1948 
1953  private:
1954 
1960  std::unique_ptr<MeshDeformation::MeshDeformationHandler<dim> > mesh_deformation;
1961 
1965  std::unique_ptr<StokesMatrixFreeHandler<dim> > stokes_matrix_free;
1966 
1967  friend class boost::serialization::access;
1968  friend class SimulatorAccess<dim>;
1969  friend class MeshDeformation::MeshDeformationHandler<dim>; // MeshDeformationHandler needs access to the internals of the Simulator
1970  friend class VolumeOfFluidHandler<dim>; // VolumeOfFluidHandler needs access to the internals of the Simulator
1971  friend class StokesMatrixFreeHandler<dim>;
1972  template <int dimension, int velocity_degree>
1974  friend struct Parameters<dim>;
1975  };
1976 }
1977 
1978 
1979 #endif
The NullspaceRemoval struct.
Definition: parameters.h:83
unsigned int nonlinear_iteration
Definition: simulator.h:1806
BoundaryVelocity::Manager< dim > boundary_velocity_manager
Definition: simulator.h:1785
const std::unique_ptr< AdiabaticConditions::Interface< dim > > adiabatic_conditions
Definition: simulator.h:1781
static const unsigned int invalid_unsigned_int
TerminationCriteria::Manager< dim > termination_manager
Definition: simulator.h:1815
void write_plugin_graph(std::ostream &output_stream)
parallel::distributed::Triangulation< dim > triangulation
Definition: simulator.h:1833
TimerOutput computing_timer
Definition: simulator.h:1753
BoundaryTemperature::Manager< dim > boundary_temperature_manager
Definition: simulator.h:1776
double pressure_scaling
Definition: simulator.h:1881
double old_time_step
Definition: simulator.h:1803
bool assemble_newton_stokes_system
Definition: simulator.h:1946
LinearAlgebra::BlockVector old_solution
Definition: simulator.h:1926
const IntermediaryConstructorAction post_geometry_model_creation_action
Definition: simulator.h:1773
std::unique_ptr< VolumeOfFluidHandler< dim > > volume_of_fluid_handler
Definition: simulator.h:1717
bool rebuild_stokes_matrix
Definition: simulator.h:1944
Parameters< dim > parameters
Definition: simulator.h:1690
std::map< types::boundary_id, std::unique_ptr< BoundaryTraction::Interface< dim > > > boundary_traction
Definition: simulator.h:1786
boost::iostreams::tee_device< std::ostream, std::ofstream > TeeDevice
Definition: simulator.h:1730
InitialComposition::Manager< dim > initial_composition_manager
Definition: simulator.h:1779
std::unique_ptr< MeshDeformation::MeshDeformationHandler< dim > > mesh_deformation
Definition: simulator.h:1960
LinearAlgebra::BlockVector system_rhs
Definition: simulator.h:1928
std::pair< double, double > stokes_residuals
Definition: simulator.h:148
boost::iostreams::stream< TeeDevice > TeeStream
Definition: simulator.h:1731
LateralAveraging< dim > lateral_averaging
Definition: simulator.h:1824
TeeStream iostream_tee_stream
Definition: simulator.h:1734
ConstraintMatrix constraints
Definition: simulator.h:1867
MeshRefinement::Manager< dim > mesh_refinement_manager
Definition: simulator.h:1837
const std::unique_ptr< MaterialModel::Interface< dim > > material_model
Definition: simulator.h:1774
double last_pressure_normalization_adjustment
Definition: simulator.h:1874
void declare_parameters(ParameterHandler &prm)
const std::unique_ptr< GeometryModel::Interface< dim > > geometry_model
Definition: simulator.h:1772
LinearAlgebra::BlockVector operator_split_reaction_vector
Definition: simulator.h:1936
Parameters< dim >::NonlinearSolver NonlinearSolver
Definition: simulator.h:225
DoFHandler< dim > dof_handler
Definition: simulator.h:1852
std::unique_ptr< StokesMatrixFreeHandler< dim > > stokes_matrix_free
Definition: simulator.h:1965
std::unique_ptr< Assemblers::Manager< dim > > assemblers
Definition: simulator.h:979
bool rebuild_sparsity_and_matrices
Definition: simulator.h:1943
const IntermediaryConstructorAction post_signal_creation
Definition: simulator.h:1708
double global_Omega_diameter
Definition: simulator.h:1834
unsigned int timestep_number
Definition: simulator.h:1804
bool assemble_newton_stokes_matrix
Definition: simulator.h:1945
ConditionalOStream pcout
Definition: simulator.h:1740
unsigned int pre_refinement_step
Definition: simulator.h:1805
TableHandler statistics
Definition: simulator.h:1751
Introspection< dim > introspection
Definition: simulator.h:1719
MPI_Comm mpi_communicator
Definition: simulator.h:1722
LinearAlgebra::BlockSparseMatrix system_preconditioner_matrix
Definition: simulator.h:1923
bool do_pressure_rhs_compatibility_modification
Definition: simulator.h:1889
LinearAlgebra::BlockVector current_linearization_point
Definition: simulator.h:1930
InitialTemperature::Manager< dim > initial_temperature_manager
Definition: simulator.h:1780
LinearAlgebra::BlockVector old_old_solution
Definition: simulator.h:1927
Definition: compat.h:53
LinearAlgebra::BlockVector pressure_shape_function_integrals
Definition: simulator.h:1933
std::ofstream log_file_stream
Definition: simulator.h:1728
std::unique_ptr< LinearAlgebra::PreconditionAMG > Amg_preconditioner
Definition: simulator.h:1940
ConstraintMatrix current_constraints
Definition: simulator.h:1868
TeeDevice iostream_tee_device
Definition: simulator.h:1733
LinearAlgebra::BlockSparseMatrix system_matrix
Definition: simulator.h:1909
std::unique_ptr< MeltHandler< dim > > melt_handler
Definition: simulator.h:1697
bool rebuild_stokes_preconditioner
Definition: simulator.h:1947
const unsigned int compositional_variable
Definition: simulator.h:258
std::unique_ptr< Mapping< dim > > mapping
Definition: simulator.h:1848
const std::unique_ptr< PrescribedStokesSolution::Interface< dim > > prescribed_stokes_solution
Definition: simulator.h:1778
HeatingModel::Manager< dim > heating_model_manager
Definition: simulator.h:1838
SimulatorSignals< dim > signals
Definition: simulator.h:1706
const std::unique_ptr< BoundaryHeatFlux::Interface< dim > > boundary_heat_flux
Definition: simulator.h:1787
const FESystem< dim > finite_element
Definition: simulator.h:1850
bool simulator_is_past_initialization
Definition: simulator.h:376
Threads::Thread output_statistics_thread
Definition: simulator.h:1762
std::unique_ptr< Particle::World< dim > > particle_world
Definition: simulator.h:1792
const std::unique_ptr< InitialTopographyModel::Interface< dim > > initial_topography_model
Definition: simulator.h:1771
LinearAlgebra::BlockVector solution
Definition: simulator.h:1925
std::unique_ptr< NewtonHandler< dim > > newton_handler
Definition: simulator.h:1704
Postprocess::Manager< dim > postprocess_manager
Definition: simulator.h:1854
const std::unique_ptr< GravityModel::Interface< dim > > gravity_model
Definition: simulator.h:1775
double newton_residual_for_derivative_scaling_factor
Definition: simulator.h:147
double global_volume
Definition: simulator.h:1835
BoundaryComposition::Manager< dim > boundary_composition_manager
Definition: simulator.h:1777
std::unique_ptr< LinearAlgebra::PreconditionBase > Mp_preconditioner
Definition: simulator.h:1941
Parameters< dim >::NullspaceRemoval NullspaceRemoval
Definition: simulator.h:230