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 
71 #include <boost/iostreams/tee.hpp>
72 #include <boost/iostreams/stream.hpp>
73 #include <memory>
74 
75 namespace aspect
76 {
77  using namespace dealii;
78 
79 #if DEAL_II_VERSION_GTE(9,1,0)
80 
88  using ConstraintMatrix = class ::AffineConstraints<double>;
89 #endif
90 
91  template <int dim>
92  class MeltHandler;
93 
94  template <int dim>
95  class NewtonHandler;
96 
97  template <int dim>
99 
100  namespace MeshDeformation
101  {
102  template <int dim>
104  }
105 
106  template <int dim>
108 
109  namespace internal
110  {
111  namespace Assembly
112  {
113  namespace Scratch
114  {
115  template <int dim> struct StokesPreconditioner;
116  template <int dim> struct StokesSystem;
117  template <int dim> struct AdvectionSystem;
118  }
119 
120  namespace CopyData
121  {
122  template <int dim> struct StokesPreconditioner;
123  template <int dim> struct StokesSystem;
124  template <int dim> struct AdvectionSystem;
125  }
126  }
127  }
128 
129  namespace Assemblers
130  {
131  template <int dim> class Interface;
132  template <int dim> class Manager;
133  }
134 
142  template <int dim>
143  class Simulator
144  {
145  public:
159  Simulator (const MPI_Comm mpi_communicator,
160  ParameterHandler &prm);
161 
166  ~Simulator ();
167 
179  static
181 
190  void run ();
191 
203  void
204  write_plugin_graph (std::ostream &output_stream) const;
205 
210 
215 
216 
224  {
229  enum FieldType { temperature_field, compositional_field };
230 
236 
242  const unsigned int compositional_variable;
243 
255  AdvectionField (const FieldType field_type,
256  const unsigned int compositional_variable = numbers::invalid_unsigned_int);
257 
265  static
267 
275  static
276  AdvectionField composition (const unsigned int compositional_variable);
277 
281  bool
282  is_temperature () const;
283 
288  bool
289  is_discontinuous (const Introspection<dim> &introspection) const;
290 
296  advection_method (const Introspection<dim> &introspection) const;
297 
302  unsigned int component_index(const Introspection<dim> &introspection) const;
303 
308  unsigned int block_index(const Introspection<dim> &introspection) const;
309 
316  unsigned int field_index() const;
317 
323  unsigned int base_element(const Introspection<dim> &introspection) const;
324 
331  FEValuesExtractors::Scalar scalar_extractor(const Introspection<dim> &introspection) const;
332 
337  unsigned int polynomial_degree(const Introspection<dim> &introspection) const;
338  };
339 
340  private:
341 
361 
379  {
380  IntermediaryConstructorAction (const std::function<void ()> &action);
381  };
382 
397  void setup_dofs ();
398 
407  void setup_introspection ();
408 
420  void set_initial_temperature_and_compositional_fields ();
421 
437  void compute_initial_pressure_field ();
438 
444  void compute_initial_velocity_boundary_constraints (ConstraintMatrix &constraints);
445 
451  void compute_current_velocity_boundary_constraints (ConstraintMatrix &constraints);
452 
465  void compute_current_constraints ();
466 
479  void compute_pressure_scaling_factor ();
480 
492  void start_timestep ();
493 
501  void solve_timestep ();
502 
514  void solve_single_advection_single_stokes ();
515 
528  void solve_no_advection_iterated_stokes ();
529 
541  void solve_first_timestep_only_single_stokes ();
542 
555  void solve_iterated_advection_and_stokes ();
556 
569  void solve_single_advection_iterated_stokes ();
570 
585  void solve_iterated_advection_and_newton_stokes ();
586 
602  void solve_single_advection_iterated_newton_stokes ();
603 
615  void solve_single_advection_no_stokes ();
616 
628  void solve_no_advection_no_stokes ();
629 
638  void build_stokes_preconditioner ();
639 
647  void build_advection_preconditioner (const AdvectionField &advection_field,
649 
656  void assemble_stokes_system ();
657 
667  double assemble_and_solve_temperature (const bool compute_initial_residual = false,
668  double *initial_residual = nullptr);
669 
682  std::vector<double> assemble_and_solve_composition (const bool compute_initial_residual = false,
683  std::vector<double> *initial_residual = nullptr);
684 
702  double assemble_and_solve_stokes (const bool compute_initial_residual = false,
703  double *initial_nonlinear_residual = nullptr);
704 
712  void assemble_advection_system (const AdvectionField &advection_field);
713 
724  double solve_advection (const AdvectionField &advection_field);
725 
729  void interpolate_particle_properties (const AdvectionField &advection_field);
730 
801  std::pair<double,double>
802  solve_stokes ();
803 
807  std::pair<double,double>
808  solve_stokes_block_gmg ();
809 
822  void postprocess ();
823 
839  void refine_mesh (const unsigned int max_grid_level);
840 
861  void create_snapshot();
862 
876  void resume_from_snapshot();
877 
884  template <class Archive>
885  void serialize (Archive &ar, const unsigned int version);
901  void setup_system_matrix (const std::vector<IndexSet> &system_partitioning);
902 
914  void setup_system_preconditioner (const std::vector<IndexSet> &system_partitioning);
915 
937  std::unique_ptr<Assemblers::Manager<dim> > assemblers;
938 
947  void set_assemblers ();
948 
959  void set_default_assemblers ();
960 
967  void assemble_stokes_preconditioner ();
968 
976  void
977  local_assemble_stokes_preconditioner (const typename DoFHandler<dim>::active_cell_iterator &cell,
980 
988  void
989  copy_local_to_global_stokes_preconditioner (const internal::Assembly::CopyData::StokesPreconditioner<dim> &data);
990 
998  void
999  local_assemble_stokes_system (const typename DoFHandler<dim>::active_cell_iterator &cell,
1002 
1010  void
1011  copy_local_to_global_stokes_system (const internal::Assembly::CopyData::StokesSystem<dim> &data);
1012 
1020  void
1021  local_assemble_advection_face_terms(const AdvectionField &advection_field,
1022  const typename DoFHandler<dim>::active_cell_iterator &cell,
1032  void
1033  local_assemble_advection_system (const AdvectionField &advection_field,
1034  const Vector<double> &viscosity_per_cell,
1035  const typename DoFHandler<dim>::active_cell_iterator &cell,
1038 
1046  void
1047  copy_local_to_global_advection_system (const AdvectionField &advection_field,
1049 
1074  void make_pressure_rhs_compatible(LinearAlgebra::BlockVector &vector);
1075 
1085  template <typename T>
1086  void get_artificial_viscosity (Vector<T> &viscosity_per_cell,
1087  const AdvectionField &advection_field,
1088  const bool skip_interior_cells = false) const;
1089 
1101  void compute_Vs_anomaly(Vector<float> &values) const;
1102 
1117  void compute_Vp_anomaly(Vector<float> &values) const;
1118 
1159  double normalize_pressure(LinearAlgebra::BlockVector &vector) const;
1160 
1205  void denormalize_pressure(const double pressure_adjustment,
1207  const LinearAlgebra::BlockVector &relevant_vector) const;
1208 
1217  void apply_limiter_to_dg_solutions (const AdvectionField &advection_field);
1218 
1219 
1242  void compute_reactions ();
1243 
1254  void initialize_current_linearization_point ();
1255 
1256 
1274  void interpolate_material_output_into_compositional_field (const unsigned int compositional_index);
1275 
1276 
1284  void interpolate_onto_velocity_system(const TensorFunction<1,dim> &func,
1285  LinearAlgebra::Vector &vec);
1286 
1287 
1302  void setup_nullspace_constraints(ConstraintMatrix &constraints);
1303 
1304 
1316  void remove_nullspace(LinearAlgebra::BlockVector &relevant_dst,
1317  LinearAlgebra::BlockVector &tmp_distributed_stokes);
1318 
1332  void remove_net_angular_momentum( const bool use_constant_density,
1333  LinearAlgebra::BlockVector &relevant_dst,
1334  LinearAlgebra::BlockVector &tmp_distributed_stokes);
1335 
1343  void replace_outflow_boundary_ids(const unsigned int boundary_id_offset);
1344 
1352  void restore_outflow_boundary_ids(const unsigned int boundary_id_offset);
1353 
1367  void remove_net_linear_momentum( const bool use_constant_density,
1368  LinearAlgebra::BlockVector &relevant_dst,
1369  LinearAlgebra::BlockVector &tmp_distributed_stokes);
1370 
1378  double get_maximal_velocity (const LinearAlgebra::BlockVector &solution) const;
1379 
1392  double get_entropy_variation (const double average_field,
1393  const AdvectionField &advection_field) const;
1394 
1403  std::pair<double,double>
1404  get_extrapolated_advection_field_range (const AdvectionField &advection_field) const;
1405 
1414  void maybe_write_timing_output () const;
1415 
1423  bool maybe_write_checkpoint (const time_t last_checkpoint_time,
1424  const std::pair<bool,bool> termination_output);
1425 
1442  bool maybe_do_initial_refinement (const unsigned int max_refinement_level);
1443 
1452  void maybe_refine_mesh (const double new_time_step,
1453  unsigned int &max_refinement_level);
1454 
1465  double compute_time_step () const;
1466 
1474  double
1475  compute_viscosity(internal::Assembly::Scratch::AdvectionSystem<dim> &scratch,
1476  const double global_u_infty,
1477  const double global_field_variation,
1478  const double average_field,
1479  const double global_entropy_variation,
1480  const double cell_diameter,
1481  const AdvectionField &advection_field) const;
1482 
1491  void
1492  compute_advection_system_residual(internal::Assembly::Scratch::AdvectionSystem<dim> &scratch,
1493  const double average_field,
1494  const AdvectionField &advection_field,
1495  double &max_residual,
1496  double &max_velocity,
1497  double &max_density,
1498  double &max_specific_heat,
1499  double &conductivity) const;
1500 
1528  void
1529  compute_material_model_input_values (const LinearAlgebra::BlockVector &input_solution,
1530  const FEValuesBase<dim,dim> &input_finite_element_values,
1531  const typename DoFHandler<dim>::active_cell_iterator &cell,
1532  const bool compute_strainrate,
1533  MaterialModel::MaterialModelInputs<dim> &material_model_inputs) const;
1534 
1535 
1550  bool
1551  stokes_matrix_depends_on_solution () const;
1552 
1565  void
1566  check_consistency_of_formulation ();
1567 
1578  void
1579  check_consistency_of_boundary_conditions () const;
1580 
1584  double
1585  compute_initial_newton_residual (const LinearAlgebra::BlockVector &linearized_stokes_initial_guess);
1586 
1595  double
1596  compute_Eisenstat_Walker_linear_tolerance(const bool EisenstatWalkerChoiceOne,
1597  const double maximum_linear_stokes_solver_tolerance,
1598  const double linear_stokes_solver_tolerance,
1599  const double stokes_residual,
1600  const double newton_residual,
1601  const double newton_residual_old);
1602 
1612  void output_statistics();
1613 
1625  double
1626  compute_initial_stokes_residual();
1627 
1638 
1644  std::unique_ptr<MeltHandler<dim> > melt_handler;
1645 
1651  std::unique_ptr<NewtonHandler<dim> > newton_handler;
1652 
1654 
1656 
1664  std::unique_ptr<VolumeOfFluidHandler<dim> > volume_of_fluid_handler;
1665 
1667 
1668 
1670 
1675  std::ofstream log_file_stream;
1676 
1677  typedef boost::iostreams::tee_device<std::ostream, std::ofstream> TeeDevice;
1678  typedef boost::iostreams::stream< TeeDevice > TeeStream;
1679 
1682 
1688 
1699 
1701 
1718  const std::unique_ptr<InitialTopographyModel::Interface<dim> > initial_topography_model;
1719  const std::unique_ptr<GeometryModel::Interface<dim> > geometry_model;
1721  const std::unique_ptr<MaterialModel::Interface<dim> > material_model;
1722  const std::unique_ptr<GravityModel::Interface<dim> > gravity_model;
1725  const std::unique_ptr<PrescribedStokesSolution::Interface<dim> > prescribed_stokes_solution;
1728  const std::unique_ptr<AdiabaticConditions::Interface<dim> > adiabatic_conditions;
1729  const std::unique_ptr<WorldBuilder::World> world_builder;
1731  std::map<types::boundary_id,std::unique_ptr<BoundaryTraction::Interface<dim> > > boundary_traction;
1732  const std::unique_ptr<BoundaryHeatFlux::Interface<dim> > boundary_heat_flux;
1733 
1741  double time;
1742  double time_step;
1744  unsigned int timestep_number;
1745  unsigned int pre_refinement_step;
1746  unsigned int nonlinear_iteration;
1776 
1779 
1788  std::unique_ptr<Mapping<dim> > mapping;
1789 
1791 
1793 
1795 
1807  ConstraintMatrix constraints;
1808  ConstraintMatrix current_constraints;
1809 
1815 
1822 
1830 
1850 
1864 
1869 
1871 
1872  // only used if is_compressible()
1874 
1875  // only used if operator split is enabled
1877 
1878 
1879 
1880  std::unique_ptr<LinearAlgebra::PreconditionAMG> Amg_preconditioner;
1881  std::unique_ptr<LinearAlgebra::PreconditionBase> Mp_preconditioner;
1882 
1888 
1893  private:
1894 
1900  std::unique_ptr<MeshDeformation::MeshDeformationHandler<dim> > mesh_deformation;
1901 
1905  std::unique_ptr<StokesMatrixFreeHandler<dim> > stokes_matrix_free;
1906 
1907  friend class boost::serialization::access;
1908  friend class SimulatorAccess<dim>;
1909  friend class MeshDeformation::MeshDeformationHandler<dim>; // MeshDeformationHandler needs access to the internals of the Simulator
1910  friend class VolumeOfFluidHandler<dim>; // VolumeOfFluidHandler needs access to the internals of the Simulator
1911  friend class StokesMatrixFreeHandler<dim>;
1912  friend struct Parameters<dim>;
1913  };
1914 }
1915 
1916 
1917 #endif
The NullspaceRemoval struct.
Definition: parameters.h:82
unsigned int nonlinear_iteration
Definition: simulator.h:1746
BoundaryVelocity::Manager< dim > boundary_velocity_manager
Definition: simulator.h:1730
const std::unique_ptr< AdiabaticConditions::Interface< dim > > adiabatic_conditions
Definition: simulator.h:1728
static const unsigned int invalid_unsigned_int
TerminationCriteria::Manager< dim > termination_manager
Definition: simulator.h:1755
void write_plugin_graph(std::ostream &output_stream)
parallel::distributed::Triangulation< dim > triangulation
Definition: simulator.h:1773
TimerOutput computing_timer
Definition: simulator.h:1700
BoundaryTemperature::Manager< dim > boundary_temperature_manager
Definition: simulator.h:1723
double pressure_scaling
Definition: simulator.h:1821
double old_time_step
Definition: simulator.h:1743
bool assemble_newton_stokes_system
Definition: simulator.h:1886
LinearAlgebra::BlockVector old_solution
Definition: simulator.h:1866
const IntermediaryConstructorAction post_geometry_model_creation_action
Definition: simulator.h:1720
std::unique_ptr< VolumeOfFluidHandler< dim > > volume_of_fluid_handler
Definition: simulator.h:1664
bool rebuild_stokes_matrix
Definition: simulator.h:1884
Parameters< dim > parameters
Definition: simulator.h:1637
std::map< types::boundary_id, std::unique_ptr< BoundaryTraction::Interface< dim > > > boundary_traction
Definition: simulator.h:1731
boost::iostreams::tee_device< std::ostream, std::ofstream > TeeDevice
Definition: simulator.h:1677
InitialComposition::Manager< dim > initial_composition_manager
Definition: simulator.h:1726
std::unique_ptr< MeshDeformation::MeshDeformationHandler< dim > > mesh_deformation
Definition: simulator.h:1900
LinearAlgebra::BlockVector system_rhs
Definition: simulator.h:1868
boost::iostreams::stream< TeeDevice > TeeStream
Definition: simulator.h:1678
LateralAveraging< dim > lateral_averaging
Definition: simulator.h:1764
TeeStream iostream_tee_stream
Definition: simulator.h:1681
ConstraintMatrix constraints
Definition: simulator.h:1807
MeshRefinement::Manager< dim > mesh_refinement_manager
Definition: simulator.h:1777
const std::unique_ptr< MaterialModel::Interface< dim > > material_model
Definition: simulator.h:1721
double last_pressure_normalization_adjustment
Definition: simulator.h:1814
void declare_parameters(ParameterHandler &prm)
const std::unique_ptr< GeometryModel::Interface< dim > > geometry_model
Definition: simulator.h:1719
LinearAlgebra::BlockVector operator_split_reaction_vector
Definition: simulator.h:1876
Parameters< dim >::NonlinearSolver NonlinearSolver
Definition: simulator.h:209
DoFHandler< dim > dof_handler
Definition: simulator.h:1792
std::unique_ptr< StokesMatrixFreeHandler< dim > > stokes_matrix_free
Definition: simulator.h:1905
std::unique_ptr< Assemblers::Manager< dim > > assemblers
Definition: simulator.h:937
bool rebuild_sparsity_and_matrices
Definition: simulator.h:1883
const IntermediaryConstructorAction post_signal_creation
Definition: simulator.h:1655
double global_Omega_diameter
Definition: simulator.h:1774
unsigned int timestep_number
Definition: simulator.h:1744
bool assemble_newton_stokes_matrix
Definition: simulator.h:1885
ConditionalOStream pcout
Definition: simulator.h:1687
unsigned int pre_refinement_step
Definition: simulator.h:1745
TableHandler statistics
Definition: simulator.h:1698
Introspection< dim > introspection
Definition: simulator.h:1666
MPI_Comm mpi_communicator
Definition: simulator.h:1669
LinearAlgebra::BlockSparseMatrix system_preconditioner_matrix
Definition: simulator.h:1863
bool do_pressure_rhs_compatibility_modification
Definition: simulator.h:1829
LinearAlgebra::BlockVector current_linearization_point
Definition: simulator.h:1870
InitialTemperature::Manager< dim > initial_temperature_manager
Definition: simulator.h:1727
LinearAlgebra::BlockVector old_old_solution
Definition: simulator.h:1867
Definition: compat.h:37
LinearAlgebra::BlockVector pressure_shape_function_integrals
Definition: simulator.h:1873
std::ofstream log_file_stream
Definition: simulator.h:1675
std::unique_ptr< LinearAlgebra::PreconditionAMG > Amg_preconditioner
Definition: simulator.h:1880
ConstraintMatrix current_constraints
Definition: simulator.h:1808
TeeDevice iostream_tee_device
Definition: simulator.h:1680
LinearAlgebra::BlockSparseMatrix system_matrix
Definition: simulator.h:1849
std::unique_ptr< MeltHandler< dim > > melt_handler
Definition: simulator.h:1644
bool rebuild_stokes_preconditioner
Definition: simulator.h:1887
const unsigned int compositional_variable
Definition: simulator.h:242
std::unique_ptr< Mapping< dim > > mapping
Definition: simulator.h:1788
const std::unique_ptr< PrescribedStokesSolution::Interface< dim > > prescribed_stokes_solution
Definition: simulator.h:1725
HeatingModel::Manager< dim > heating_model_manager
Definition: simulator.h:1778
SimulatorSignals< dim > signals
Definition: simulator.h:1653
const std::unique_ptr< BoundaryHeatFlux::Interface< dim > > boundary_heat_flux
Definition: simulator.h:1732
const FESystem< dim > finite_element
Definition: simulator.h:1790
bool simulator_is_past_initialization
Definition: simulator.h:360
Threads::Thread output_statistics_thread
Definition: simulator.h:1709
const std::unique_ptr< InitialTopographyModel::Interface< dim > > initial_topography_model
Definition: simulator.h:1718
LinearAlgebra::BlockVector solution
Definition: simulator.h:1865
std::unique_ptr< NewtonHandler< dim > > newton_handler
Definition: simulator.h:1651
Postprocess::Manager< dim > postprocess_manager
Definition: simulator.h:1794
const std::unique_ptr< GravityModel::Interface< dim > > gravity_model
Definition: simulator.h:1722
double global_volume
Definition: simulator.h:1775
BoundaryComposition::Manager< dim > boundary_composition_manager
Definition: simulator.h:1724
const std::unique_ptr< WorldBuilder::World > world_builder
Definition: simulator.h:1729
std::unique_ptr< LinearAlgebra::PreconditionBase > Mp_preconditioner
Definition: simulator.h:1881
Parameters< dim >::NullspaceRemoval NullspaceRemoval
Definition: simulator.h:214