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 
143  template <int dim>
144  class Simulator
145  {
146  public:
160  Simulator (const MPI_Comm mpi_communicator,
161  ParameterHandler &prm);
162 
167  ~Simulator ();
168 
180  static
182 
191  void run ();
192 
204  void
205  write_plugin_graph (std::ostream &output_stream) const;
206 
211 
216 
217 
225  {
230  enum FieldType { temperature_field, compositional_field };
231 
237 
243  const unsigned int compositional_variable;
244 
256  AdvectionField (const FieldType field_type,
257  const unsigned int compositional_variable = numbers::invalid_unsigned_int);
258 
266  static
268 
276  static
277  AdvectionField composition (const unsigned int compositional_variable);
278 
282  bool
283  is_temperature () const;
284 
289  bool
290  is_discontinuous (const Introspection<dim> &introspection) const;
291 
297  advection_method (const Introspection<dim> &introspection) const;
298 
303  unsigned int component_index(const Introspection<dim> &introspection) const;
304 
309  unsigned int block_index(const Introspection<dim> &introspection) const;
310 
317  unsigned int field_index() const;
318 
324  unsigned int base_element(const Introspection<dim> &introspection) const;
325 
332  FEValuesExtractors::Scalar scalar_extractor(const Introspection<dim> &introspection) const;
333 
338  unsigned int polynomial_degree(const Introspection<dim> &introspection) const;
339  };
340 
341  private:
342 
362 
380  {
381  IntermediaryConstructorAction (const std::function<void ()> &action);
382  };
383 
398  void setup_dofs ();
399 
408  void setup_introspection ();
409 
421  void set_initial_temperature_and_compositional_fields ();
422 
438  void compute_initial_pressure_field ();
439 
445  void compute_initial_velocity_boundary_constraints (AffineConstraints<double> &constraints);
446 
452  void compute_current_velocity_boundary_constraints (AffineConstraints<double> &constraints);
453 
466  void compute_current_constraints ();
467 
480  void compute_pressure_scaling_factor ();
481 
493  void start_timestep ();
494 
502  void solve_timestep ();
503 
515  void solve_single_advection_single_stokes ();
516 
529  void solve_no_advection_iterated_stokes ();
530 
541  void solve_no_advection_single_stokes ();
542 
554  void solve_first_timestep_only_single_stokes ();
555 
568  void solve_iterated_advection_and_stokes ();
569 
582  void solve_single_advection_iterated_stokes ();
583 
598  void solve_iterated_advection_and_newton_stokes ();
599 
615  void solve_single_advection_iterated_newton_stokes ();
616 
628  void solve_single_advection_no_stokes ();
629 
641  void solve_no_advection_no_stokes ();
642 
651  void build_stokes_preconditioner ();
652 
660  void build_advection_preconditioner (const AdvectionField &advection_field,
662  const double diagonal_strengthening);
663 
670  void assemble_stokes_system ();
671 
681  double assemble_and_solve_temperature (const bool compute_initial_residual = false,
682  double *initial_residual = nullptr);
683 
696  std::vector<double> assemble_and_solve_composition (const bool compute_initial_residual = false,
697  std::vector<double> *initial_residual = nullptr);
698 
716  double assemble_and_solve_stokes (const bool compute_initial_residual = false,
717  double *initial_nonlinear_residual = nullptr);
718 
730  void assemble_and_solve_defect_correction_Stokes(DefectCorrectionResiduals &dcr,
731  const bool use_picard);
732 
740  void assemble_advection_system (const AdvectionField &advection_field);
741 
752  double solve_advection (const AdvectionField &advection_field);
753 
757  void interpolate_particle_properties (const AdvectionField &advection_field);
758 
829  std::pair<double,double>
830  solve_stokes ();
831 
835  std::pair<double,double>
836  solve_stokes_block_gmg ();
837 
850  void postprocess ();
851 
867  void refine_mesh (const unsigned int max_grid_level);
868 
889  void create_snapshot();
890 
904  void resume_from_snapshot();
905 
912  template <class Archive>
913  void serialize (Archive &ar, const unsigned int version);
931  setup_system_matrix_coupling () const;
932 
940  void setup_system_matrix (const std::vector<IndexSet> &system_partitioning);
941 
953  void setup_system_preconditioner (const std::vector<IndexSet> &system_partitioning);
954 
976  std::unique_ptr<Assemblers::Manager<dim> > assemblers;
977 
986  void set_assemblers ();
987 
998  void set_default_assemblers ();
999 
1006  void assemble_stokes_preconditioner ();
1007 
1015  void
1016  local_assemble_stokes_preconditioner (const typename DoFHandler<dim>::active_cell_iterator &cell,
1019 
1027  void
1028  copy_local_to_global_stokes_preconditioner (const internal::Assembly::CopyData::StokesPreconditioner<dim> &data);
1029 
1037  void
1038  local_assemble_stokes_system (const typename DoFHandler<dim>::active_cell_iterator &cell,
1041 
1049  void
1050  copy_local_to_global_stokes_system (const internal::Assembly::CopyData::StokesSystem<dim> &data);
1051 
1059  void
1060  local_assemble_advection_face_terms(const AdvectionField &advection_field,
1061  const typename DoFHandler<dim>::active_cell_iterator &cell,
1071  void
1072  local_assemble_advection_system (const AdvectionField &advection_field,
1073  const Vector<double> &viscosity_per_cell,
1074  const typename DoFHandler<dim>::active_cell_iterator &cell,
1077 
1085  void
1086  copy_local_to_global_advection_system (const AdvectionField &advection_field,
1088 
1113  void make_pressure_rhs_compatible(LinearAlgebra::BlockVector &vector);
1114 
1124  template <typename T>
1125  void get_artificial_viscosity (Vector<T> &viscosity_per_cell,
1126  const AdvectionField &advection_field,
1127  const bool skip_interior_cells = false) const;
1128 
1140  void compute_Vs_anomaly(Vector<float> &values) const;
1141 
1156  void compute_Vp_anomaly(Vector<float> &values) const;
1157 
1198  double normalize_pressure(LinearAlgebra::BlockVector &vector) const;
1199 
1244  void denormalize_pressure(const double pressure_adjustment,
1246  const LinearAlgebra::BlockVector &relevant_vector) const;
1247 
1256  void apply_limiter_to_dg_solutions (const AdvectionField &advection_field);
1257 
1258 
1281  void compute_reactions ();
1282 
1292  void update_solution_vectors_with_reaction_results (const unsigned int block_index,
1293  const LinearAlgebra::BlockVector &distributed_vector,
1294  const LinearAlgebra::BlockVector &distributed_reaction_vector);
1295 
1306  void initialize_current_linearization_point ();
1307 
1327  void interpolate_material_output_into_advection_field (const AdvectionField &adv_field);
1328 
1329 
1337  void interpolate_onto_velocity_system(const TensorFunction<1,dim> &func,
1338  LinearAlgebra::Vector &vec);
1339 
1340 
1355  void setup_nullspace_constraints(AffineConstraints<double> &constraints);
1356 
1357 
1369  void remove_nullspace(LinearAlgebra::BlockVector &relevant_dst,
1370  LinearAlgebra::BlockVector &tmp_distributed_stokes);
1371 
1385  void remove_net_angular_momentum( const bool use_constant_density,
1386  LinearAlgebra::BlockVector &relevant_dst,
1387  LinearAlgebra::BlockVector &tmp_distributed_stokes);
1388 
1396  void replace_outflow_boundary_ids(const unsigned int boundary_id_offset);
1397 
1405  void restore_outflow_boundary_ids(const unsigned int boundary_id_offset);
1406 
1420  void remove_net_linear_momentum( const bool use_constant_density,
1421  LinearAlgebra::BlockVector &relevant_dst,
1422  LinearAlgebra::BlockVector &tmp_distributed_stokes);
1423 
1431  double get_maximal_velocity (const LinearAlgebra::BlockVector &solution) const;
1432 
1445  double get_entropy_variation (const double average_field,
1446  const AdvectionField &advection_field) const;
1447 
1456  std::pair<double,double>
1457  get_extrapolated_advection_field_range (const AdvectionField &advection_field) const;
1458 
1467  void maybe_write_timing_output () const;
1468 
1476  bool maybe_write_checkpoint (const time_t last_checkpoint_time,
1477  const bool force_writing_checkpoint);
1478 
1495  bool maybe_do_initial_refinement (const unsigned int max_refinement_level);
1496 
1505  void maybe_refine_mesh (const double new_time_step,
1506  unsigned int &max_refinement_level);
1507 
1515  void advance_time (const double step_size);
1516 
1524  double
1525  compute_viscosity(internal::Assembly::Scratch::AdvectionSystem<dim> &scratch,
1526  const double global_u_infty,
1527  const double global_field_variation,
1528  const double average_field,
1529  const double global_entropy_variation,
1530  const double cell_diameter,
1531  const AdvectionField &advection_field) const;
1532 
1541  void
1542  compute_advection_system_residual(internal::Assembly::Scratch::AdvectionSystem<dim> &scratch,
1543  const double average_field,
1544  const AdvectionField &advection_field,
1545  double &max_residual,
1546  double &max_velocity,
1547  double &max_density,
1548  double &max_specific_heat,
1549  double &conductivity) const;
1550 
1578  void
1579  compute_material_model_input_values (const LinearAlgebra::BlockVector &input_solution,
1580  const FEValuesBase<dim,dim> &input_finite_element_values,
1581  const typename DoFHandler<dim>::active_cell_iterator &cell,
1582  const bool compute_strainrate,
1583  MaterialModel::MaterialModelInputs<dim> &material_model_inputs) const;
1584 
1585 
1600  bool
1601  stokes_matrix_depends_on_solution () const;
1602 
1615  void
1616  check_consistency_of_formulation ();
1617 
1628  void
1629  check_consistency_of_boundary_conditions () const;
1630 
1634  double
1635  compute_initial_newton_residual (const LinearAlgebra::BlockVector &linearized_stokes_initial_guess);
1636 
1645  double
1646  compute_Eisenstat_Walker_linear_tolerance(const bool EisenstatWalkerChoiceOne,
1647  const double maximum_linear_stokes_solver_tolerance,
1648  const double linear_stokes_solver_tolerance,
1649  const double stokes_residual,
1650  const double newton_residual,
1651  const double newton_residual_old);
1652 
1662  void output_statistics();
1663 
1675  double
1676  compute_initial_stokes_residual();
1677 
1688 
1694  std::unique_ptr<MeltHandler<dim> > melt_handler;
1695 
1701  std::unique_ptr<NewtonHandler<dim> > newton_handler;
1702 
1704 
1706 
1714  std::unique_ptr<VolumeOfFluidHandler<dim> > volume_of_fluid_handler;
1715 
1717 
1718 
1720 
1725  std::ofstream log_file_stream;
1726 
1727  using TeeDevice = boost::iostreams::tee_device<std::ostream, std::ofstream>;
1728  using TeeStream = boost::iostreams::stream< TeeDevice >;
1729 
1732 
1738 
1749 
1766 
1768 
1785  const std::unique_ptr<InitialTopographyModel::Interface<dim> > initial_topography_model;
1786  const std::unique_ptr<GeometryModel::Interface<dim> > geometry_model;
1788  const std::unique_ptr<MaterialModel::Interface<dim> > material_model;
1789  const std::unique_ptr<GravityModel::Interface<dim> > gravity_model;
1792  const std::unique_ptr<PrescribedStokesSolution::Interface<dim> > prescribed_stokes_solution;
1795  const std::unique_ptr<AdiabaticConditions::Interface<dim> > adiabatic_conditions;
1796 #ifdef ASPECT_WITH_WORLD_BUILDER
1797  const std::unique_ptr<WorldBuilder::World> world_builder;
1798 #endif
1800  std::map<types::boundary_id,std::unique_ptr<BoundaryTraction::Interface<dim> > > boundary_traction;
1801  const std::unique_ptr<BoundaryHeatFlux::Interface<dim> > boundary_heat_flux;
1802 
1806  std::unique_ptr<Particle::World<dim> > particle_world;
1807 
1815  double time;
1816  double time_step;
1818  unsigned int timestep_number;
1819  unsigned int pre_refinement_step;
1820  unsigned int nonlinear_iteration;
1850 
1853 
1862  std::unique_ptr<Mapping<dim> > mapping;
1863 
1865 
1867 
1869 
1883 
1889 
1896 
1904 
1924 
1938 
1943 
1945 
1946  // only used if is_compressible()
1948 
1949  // only used if operator split is enabled
1951 
1952 
1953 
1954  std::unique_ptr<LinearAlgebra::PreconditionAMG> Amg_preconditioner;
1955  std::unique_ptr<LinearAlgebra::PreconditionBase> Mp_preconditioner;
1956 
1962 
1967  private:
1968 
1974  std::unique_ptr<MeshDeformation::MeshDeformationHandler<dim> > mesh_deformation;
1975 
1979  std::unique_ptr<StokesMatrixFreeHandler<dim> > stokes_matrix_free;
1980 
1981  friend class boost::serialization::access;
1982  friend class SimulatorAccess<dim>;
1983  friend class MeshDeformation::MeshDeformationHandler<dim>; // MeshDeformationHandler needs access to the internals of the Simulator
1984  friend class VolumeOfFluidHandler<dim>; // VolumeOfFluidHandler needs access to the internals of the Simulator
1985  friend class StokesMatrixFreeHandler<dim>;
1986  template <int dimension, int velocity_degree>
1988  friend struct Parameters<dim>;
1989  };
1990 }
1991 
1992 
1993 #endif
The NullspaceRemoval struct.
Definition: parameters.h:83
unsigned int nonlinear_iteration
Definition: simulator.h:1820
BoundaryVelocity::Manager< dim > boundary_velocity_manager
Definition: simulator.h:1799
const std::unique_ptr< AdiabaticConditions::Interface< dim > > adiabatic_conditions
Definition: simulator.h:1795
static const unsigned int invalid_unsigned_int
void write_plugin_graph(std::ostream &output_stream)
parallel::distributed::Triangulation< dim > triangulation
Definition: simulator.h:1847
TimerOutput computing_timer
Definition: simulator.h:1767
BoundaryTemperature::Manager< dim > boundary_temperature_manager
Definition: simulator.h:1790
boost::iostreams::tee_device< std::ostream, std::ofstream > TeeDevice
Definition: simulator.h:1727
double pressure_scaling
Definition: simulator.h:1895
double old_time_step
Definition: simulator.h:1817
AffineConstraints< double > constraints
Definition: simulator.h:1881
bool assemble_newton_stokes_system
Definition: simulator.h:1960
LinearAlgebra::BlockVector old_solution
Definition: simulator.h:1940
const IntermediaryConstructorAction post_geometry_model_creation_action
Definition: simulator.h:1787
std::unique_ptr< VolumeOfFluidHandler< dim > > volume_of_fluid_handler
Definition: simulator.h:1714
bool rebuild_stokes_matrix
Definition: simulator.h:1958
Parameters< dim > parameters
Definition: simulator.h:1687
AffineConstraints< double > current_constraints
Definition: simulator.h:1882
std::map< types::boundary_id, std::unique_ptr< BoundaryTraction::Interface< dim > > > boundary_traction
Definition: simulator.h:1800
InitialComposition::Manager< dim > initial_composition_manager
Definition: simulator.h:1793
std::unique_ptr< MeshDeformation::MeshDeformationHandler< dim > > mesh_deformation
Definition: simulator.h:1974
LinearAlgebra::BlockVector system_rhs
Definition: simulator.h:1942
std::pair< double, double > stokes_residuals
Definition: simulator.h:133
LateralAveraging< dim > lateral_averaging
Definition: simulator.h:1838
TeeStream iostream_tee_stream
Definition: simulator.h:1731
MeshRefinement::Manager< dim > mesh_refinement_manager
Definition: simulator.h:1851
std::size_t statistics_last_write_size
Definition: simulator.h:1764
boost::iostreams::stream< TeeDevice > TeeStream
Definition: simulator.h:1728
const std::unique_ptr< MaterialModel::Interface< dim > > material_model
Definition: simulator.h:1788
#define DEAL_II_DISABLE_EXTRA_DIAGNOSTICS
double last_pressure_normalization_adjustment
Definition: simulator.h:1888
void declare_parameters(ParameterHandler &prm)
const std::unique_ptr< GeometryModel::Interface< dim > > geometry_model
Definition: simulator.h:1786
double initial_residual
LinearAlgebra::BlockVector operator_split_reaction_vector
Definition: simulator.h:1950
DoFHandler< dim > dof_handler
Definition: simulator.h:1866
std::unique_ptr< StokesMatrixFreeHandler< dim > > stokes_matrix_free
Definition: simulator.h:1979
std::unique_ptr< Assemblers::Manager< dim > > assemblers
Definition: simulator.h:976
bool rebuild_sparsity_and_matrices
Definition: simulator.h:1957
const IntermediaryConstructorAction post_signal_creation
Definition: simulator.h:1705
double global_Omega_diameter
Definition: simulator.h:1848
unsigned int timestep_number
Definition: simulator.h:1818
bool assemble_newton_stokes_matrix
Definition: simulator.h:1959
ConditionalOStream pcout
Definition: simulator.h:1737
unsigned int pre_refinement_step
Definition: simulator.h:1819
TimeStepping::Manager< dim > time_stepping_manager
Definition: simulator.h:1829
std::size_t statistics_last_hash
Definition: simulator.h:1765
TableHandler statistics
Definition: simulator.h:1748
Introspection< dim > introspection
Definition: simulator.h:1716
MPI_Comm mpi_communicator
Definition: simulator.h:1719
typename Parameters< dim >::NullspaceRemoval NullspaceRemoval
Definition: simulator.h:215
LinearAlgebra::BlockSparseMatrix system_preconditioner_matrix
Definition: simulator.h:1937
bool do_pressure_rhs_compatibility_modification
Definition: simulator.h:1903
LinearAlgebra::BlockVector current_linearization_point
Definition: simulator.h:1944
InitialTemperature::Manager< dim > initial_temperature_manager
Definition: simulator.h:1794
LinearAlgebra::BlockVector old_old_solution
Definition: simulator.h:1941
LinearAlgebra::BlockVector pressure_shape_function_integrals
Definition: simulator.h:1947
std::ofstream log_file_stream
Definition: simulator.h:1725
std::unique_ptr< LinearAlgebra::PreconditionAMG > Amg_preconditioner
Definition: simulator.h:1954
TeeDevice iostream_tee_device
Definition: simulator.h:1730
LinearAlgebra::BlockSparseMatrix system_matrix
Definition: simulator.h:1923
#define DEAL_II_ENABLE_EXTRA_DIAGNOSTICS
std::unique_ptr< MeltHandler< dim > > melt_handler
Definition: simulator.h:1694
bool rebuild_stokes_preconditioner
Definition: simulator.h:1961
const unsigned int compositional_variable
Definition: simulator.h:243
std::unique_ptr< Mapping< dim > > mapping
Definition: simulator.h:1862
const std::unique_ptr< PrescribedStokesSolution::Interface< dim > > prescribed_stokes_solution
Definition: simulator.h:1792
HeatingModel::Manager< dim > heating_model_manager
Definition: simulator.h:1852
SimulatorSignals< dim > signals
Definition: simulator.h:1703
const std::unique_ptr< BoundaryHeatFlux::Interface< dim > > boundary_heat_flux
Definition: simulator.h:1801
const FESystem< dim > finite_element
Definition: simulator.h:1864
bool simulator_is_past_initialization
Definition: simulator.h:361
Threads::Thread output_statistics_thread
Definition: simulator.h:1776
std::unique_ptr< Particle::World< dim > > particle_world
Definition: simulator.h:1806
const std::unique_ptr< InitialTopographyModel::Interface< dim > > initial_topography_model
Definition: simulator.h:1785
LinearAlgebra::BlockVector solution
Definition: simulator.h:1939
std::unique_ptr< NewtonHandler< dim > > newton_handler
Definition: simulator.h:1701
typename Parameters< dim >::NonlinearSolver NonlinearSolver
Definition: simulator.h:210
Postprocess::Manager< dim > postprocess_manager
Definition: simulator.h:1868
const std::unique_ptr< GravityModel::Interface< dim > > gravity_model
Definition: simulator.h:1789
double newton_residual_for_derivative_scaling_factor
Definition: simulator.h:132
double global_volume
Definition: simulator.h:1849
BoundaryComposition::Manager< dim > boundary_composition_manager
Definition: simulator.h:1791
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:1955