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>
98  class FreeSurfaceHandler;
99 
100  template <int dim>
102 
103  namespace internal
104  {
105  namespace Assembly
106  {
107  namespace Scratch
108  {
109  template <int dim> struct StokesPreconditioner;
110  template <int dim> struct StokesSystem;
111  template <int dim> struct AdvectionSystem;
112  }
113 
114  namespace CopyData
115  {
116  template <int dim> struct StokesPreconditioner;
117  template <int dim> struct StokesSystem;
118  template <int dim> struct AdvectionSystem;
119  }
120  }
121  }
122 
123  namespace Assemblers
124  {
125  template <int dim> class Interface;
126  template <int dim> class Manager;
127  }
128 
136  template <int dim>
137  class Simulator
138  {
139  public:
153  Simulator (const MPI_Comm mpi_communicator,
154  ParameterHandler &prm);
155 
160  ~Simulator ();
161 
173  static
175 
184  void run ();
185 
197  void
198  write_plugin_graph (std::ostream &output_stream) const;
199 
204 
209 
210 
218  {
223  enum FieldType { temperature_field, compositional_field };
224 
230 
236  const unsigned int compositional_variable;
237 
249  AdvectionField (const FieldType field_type,
250  const unsigned int compositional_variable = numbers::invalid_unsigned_int);
251 
259  static
261 
269  static
270  AdvectionField composition (const unsigned int compositional_variable);
271 
275  bool
276  is_temperature () const;
277 
282  bool
283  is_discontinuous (const Introspection<dim> &introspection) const;
284 
290  advection_method (const Introspection<dim> &introspection) const;
291 
296  unsigned int component_index(const Introspection<dim> &introspection) const;
297 
302  unsigned int block_index(const Introspection<dim> &introspection) const;
303 
310  unsigned int field_index() const;
311 
317  unsigned int base_element(const Introspection<dim> &introspection) const;
318 
325  FEValuesExtractors::Scalar scalar_extractor(const Introspection<dim> &introspection) const;
326 
331  unsigned int polynomial_degree(const Introspection<dim> &introspection) const;
332  };
333 
334  private:
335 
336 
354  {
355  IntermediaryConstructorAction (const std::function<void ()> &action);
356  };
357 
372  void setup_dofs ();
373 
382  void setup_introspection ();
383 
395  void set_initial_temperature_and_compositional_fields ();
396 
412  void compute_initial_pressure_field ();
413 
419  void compute_initial_velocity_boundary_constraints (ConstraintMatrix &constraints);
420 
426  void compute_current_velocity_boundary_constraints (ConstraintMatrix &constraints);
427 
440  void compute_current_constraints ();
441 
453  void start_timestep ();
454 
462  void solve_timestep ();
463 
475  void solve_single_advection_single_stokes ();
476 
489  void solve_no_advection_iterated_stokes ();
490 
502  void solve_first_timestep_only_single_stokes ();
503 
516  void solve_iterated_advection_and_stokes ();
517 
530  void solve_single_advection_iterated_stokes ();
531 
546  void solve_iterated_advection_and_newton_stokes ();
547 
559  void solve_single_advection_no_stokes ();
560 
572  void solve_no_advection_no_stokes ();
573 
582  void build_stokes_preconditioner ();
583 
591  void build_advection_preconditioner (const AdvectionField &advection_field,
593 
600  void assemble_stokes_system ();
601 
611  double assemble_and_solve_temperature (const bool compute_initial_residual = false,
612  double *initial_residual = nullptr);
613 
626  std::vector<double> assemble_and_solve_composition (const bool compute_initial_residual = false,
627  std::vector<double> *initial_residual = nullptr);
628 
646  double assemble_and_solve_stokes (const bool compute_initial_residual = false,
647  double *initial_nonlinear_residual = nullptr);
648 
656  void assemble_advection_system (const AdvectionField &advection_field);
657 
668  double solve_advection (const AdvectionField &advection_field);
669 
673  void interpolate_particle_properties (const AdvectionField &advection_field);
674 
745  std::pair<double,double>
746  solve_stokes ();
747 
760  void postprocess ();
761 
777  void refine_mesh (const unsigned int max_grid_level);
778 
799  void create_snapshot();
800 
814  void resume_from_snapshot();
815 
822  template <class Archive>
823  void serialize (Archive &ar, const unsigned int version);
839  void setup_system_matrix (const std::vector<IndexSet> &system_partitioning);
840 
852  void setup_system_preconditioner (const std::vector<IndexSet> &system_partitioning);
853 
875  std::unique_ptr<Assemblers::Manager<dim> > assemblers;
876 
885  void set_assemblers ();
886 
897  void set_default_assemblers ();
898 
905  void assemble_stokes_preconditioner ();
906 
914  void
915  local_assemble_stokes_preconditioner (const typename DoFHandler<dim>::active_cell_iterator &cell,
918 
926  void
927  copy_local_to_global_stokes_preconditioner (const internal::Assembly::CopyData::StokesPreconditioner<dim> &data);
928 
936  void
937  local_assemble_stokes_system (const typename DoFHandler<dim>::active_cell_iterator &cell,
940 
948  void
949  copy_local_to_global_stokes_system (const internal::Assembly::CopyData::StokesSystem<dim> &data);
950 
958  void
959  local_assemble_advection_face_terms(const AdvectionField &advection_field,
960  const typename DoFHandler<dim>::active_cell_iterator &cell,
970  void
971  local_assemble_advection_system (const AdvectionField &advection_field,
972  const Vector<double> &viscosity_per_cell,
973  const typename DoFHandler<dim>::active_cell_iterator &cell,
976 
984  void
985  copy_local_to_global_advection_system (const AdvectionField &advection_field,
987 
1012  void make_pressure_rhs_compatible(LinearAlgebra::BlockVector &vector);
1013 
1023  template <typename T>
1024  void get_artificial_viscosity (Vector<T> &viscosity_per_cell,
1025  const AdvectionField &advection_field,
1026  const bool skip_interior_cells = false) const;
1027 
1039  void compute_Vs_anomaly(Vector<float> &values) const;
1040 
1055  void compute_Vp_anomaly(Vector<float> &values) const;
1056 
1097  double normalize_pressure(LinearAlgebra::BlockVector &vector) const;
1098 
1143  void denormalize_pressure(const double pressure_adjustment,
1145  const LinearAlgebra::BlockVector &relevant_vector) const;
1146 
1155  void apply_limiter_to_dg_solutions (const AdvectionField &advection_field);
1156 
1157 
1180  void compute_reactions ();
1181 
1182 
1205  void interpolate_material_output_into_compositional_field (const unsigned int compositional_index);
1206 
1207 
1215  void interpolate_onto_velocity_system(const TensorFunction<1,dim> &func,
1216  LinearAlgebra::Vector &vec);
1217 
1218 
1233  void setup_nullspace_constraints(ConstraintMatrix &constraints);
1234 
1235 
1247  void remove_nullspace(LinearAlgebra::BlockVector &relevant_dst,
1248  LinearAlgebra::BlockVector &tmp_distributed_stokes);
1249 
1263  void remove_net_angular_momentum( const bool use_constant_density,
1264  LinearAlgebra::BlockVector &relevant_dst,
1265  LinearAlgebra::BlockVector &tmp_distributed_stokes);
1266 
1274  void replace_outflow_boundary_ids(const unsigned int boundary_id_offset);
1275 
1283  void restore_outflow_boundary_ids(const unsigned int boundary_id_offset);
1284 
1298  void remove_net_linear_momentum( const bool use_constant_density,
1299  LinearAlgebra::BlockVector &relevant_dst,
1300  LinearAlgebra::BlockVector &tmp_distributed_stokes);
1301 
1309  double get_maximal_velocity (const LinearAlgebra::BlockVector &solution) const;
1310 
1323  double get_entropy_variation (const double average_field,
1324  const AdvectionField &advection_field) const;
1325 
1334  std::pair<double,double>
1335  get_extrapolated_advection_field_range (const AdvectionField &advection_field) const;
1336 
1345  void maybe_write_timing_output () const;
1346 
1354  bool maybe_write_checkpoint (const time_t last_checkpoint_time,
1355  const std::pair<bool,bool> termination_output);
1356 
1373  bool maybe_do_initial_refinement (const unsigned int max_refinement_level);
1374 
1383  void maybe_refine_mesh (const double new_time_step,
1384  unsigned int &max_refinement_level);
1385 
1396  double compute_time_step () const;
1397 
1405  double
1406  compute_viscosity(internal::Assembly::Scratch::AdvectionSystem<dim> &scratch,
1407  const double global_u_infty,
1408  const double global_field_variation,
1409  const double average_field,
1410  const double global_entropy_variation,
1411  const double cell_diameter,
1412  const AdvectionField &advection_field) const;
1413 
1422  void
1423  compute_advection_system_residual(internal::Assembly::Scratch::AdvectionSystem<dim> &scratch,
1424  const double average_field,
1425  const AdvectionField &advection_field,
1426  double &max_residual,
1427  double &max_velocity,
1428  double &max_density,
1429  double &max_specific_heat,
1430  double &conductivity) const;
1431 
1459  void
1460  compute_material_model_input_values (const LinearAlgebra::BlockVector &input_solution,
1461  const FEValuesBase<dim,dim> &input_finite_element_values,
1462  const typename DoFHandler<dim>::active_cell_iterator &cell,
1463  const bool compute_strainrate,
1464  MaterialModel::MaterialModelInputs<dim> &material_model_inputs) const;
1465 
1466 
1481  bool
1482  stokes_matrix_depends_on_solution () const;
1483 
1496  void
1497  check_consistency_of_formulation ();
1498 
1509  void
1510  check_consistency_of_boundary_conditions () const;
1511 
1515  double
1516  compute_initial_newton_residual (const LinearAlgebra::BlockVector &linearized_stokes_initial_guess);
1517 
1526  double
1527  compute_Eisenstat_Walker_linear_tolerance(const bool EisenstatWalkerChoiceOne,
1528  const double maximum_linear_stokes_solver_tolerance,
1529  const double linear_stokes_solver_tolerance,
1530  const double stokes_residual,
1531  const double newton_residual,
1532  const double newton_residual_old);
1533 
1543  void output_statistics();
1544 
1556  double
1557  compute_initial_stokes_residual();
1558 
1569 
1575  std::unique_ptr<MeltHandler<dim> > melt_handler;
1576 
1582  std::unique_ptr<NewtonHandler<dim> > newton_handler;
1583 
1585 
1587 
1595  std::unique_ptr<VolumeOfFluidHandler<dim> > volume_of_fluid_handler;
1596 
1598 
1599 
1601 
1606  std::ofstream log_file_stream;
1607 
1608  typedef boost::iostreams::tee_device<std::ostream, std::ofstream> TeeDevice;
1609  typedef boost::iostreams::stream< TeeDevice > TeeStream;
1610 
1613 
1619 
1630 
1632 
1649  const std::unique_ptr<InitialTopographyModel::Interface<dim> > initial_topography_model;
1650  const std::unique_ptr<GeometryModel::Interface<dim> > geometry_model;
1652  const std::unique_ptr<MaterialModel::Interface<dim> > material_model;
1653  const std::unique_ptr<GravityModel::Interface<dim> > gravity_model;
1656  const std::unique_ptr<PrescribedStokesSolution::Interface<dim> > prescribed_stokes_solution;
1659  const std::unique_ptr<AdiabaticConditions::Interface<dim> > adiabatic_conditions;
1660  const std::unique_ptr<WorldBuilder::World> world_builder;
1662  std::map<types::boundary_id,std::unique_ptr<BoundaryTraction::Interface<dim> > > boundary_traction;
1663  const std::unique_ptr<BoundaryHeatFlux::Interface<dim> > boundary_heat_flux;
1664 
1672  double time;
1673  double time_step;
1675  unsigned int timestep_number;
1676  unsigned int pre_refinement_step;
1677  unsigned int nonlinear_iteration;
1707 
1710 
1719  std::unique_ptr<Mapping<dim> > mapping;
1720 
1722 
1724 
1726 
1738  ConstraintMatrix constraints;
1739  ConstraintMatrix current_constraints;
1740 
1746 
1753 
1761 
1781 
1795 
1800 
1802 
1803  // only used if is_compressible()
1805 
1806  // only used if operator split is enabled
1808 
1809 
1810 
1811  std::unique_ptr<LinearAlgebra::PreconditionAMG> Amg_preconditioner;
1812  std::unique_ptr<LinearAlgebra::PreconditionBase> Mp_preconditioner;
1813 
1819 
1824  private:
1825 
1831  std::unique_ptr<FreeSurfaceHandler<dim> > free_surface;
1832 
1833  friend class boost::serialization::access;
1834  friend class SimulatorAccess<dim>;
1835  friend class FreeSurfaceHandler<dim>; // FreeSurfaceHandler needs access to the internals of the Simulator
1836  friend class VolumeOfFluidHandler<dim>; // VolumeOfFluidHandler needs access to the internals of the Simulator
1837  friend struct Parameters<dim>;
1838  };
1839 }
1840 
1841 
1842 #endif
The NullspaceRemoval struct.
Definition: parameters.h:81
unsigned int nonlinear_iteration
Definition: simulator.h:1677
BoundaryVelocity::Manager< dim > boundary_velocity_manager
Definition: simulator.h:1661
const std::unique_ptr< AdiabaticConditions::Interface< dim > > adiabatic_conditions
Definition: simulator.h:1659
static const unsigned int invalid_unsigned_int
TerminationCriteria::Manager< dim > termination_manager
Definition: simulator.h:1686
void write_plugin_graph(std::ostream &output_stream)
parallel::distributed::Triangulation< dim > triangulation
Definition: simulator.h:1704
TimerOutput computing_timer
Definition: simulator.h:1631
BoundaryTemperature::Manager< dim > boundary_temperature_manager
Definition: simulator.h:1654
double pressure_scaling
Definition: simulator.h:1752
double old_time_step
Definition: simulator.h:1674
bool assemble_newton_stokes_system
Definition: simulator.h:1817
LinearAlgebra::BlockVector old_solution
Definition: simulator.h:1797
const IntermediaryConstructorAction post_geometry_model_creation_action
Definition: simulator.h:1651
std::unique_ptr< VolumeOfFluidHandler< dim > > volume_of_fluid_handler
Definition: simulator.h:1595
bool rebuild_stokes_matrix
Definition: simulator.h:1815
Parameters< dim > parameters
Definition: simulator.h:1568
std::unique_ptr< FreeSurfaceHandler< dim > > free_surface
Definition: simulator.h:1831
std::map< types::boundary_id, std::unique_ptr< BoundaryTraction::Interface< dim > > > boundary_traction
Definition: simulator.h:1662
boost::iostreams::tee_device< std::ostream, std::ofstream > TeeDevice
Definition: simulator.h:1608
InitialComposition::Manager< dim > initial_composition_manager
Definition: simulator.h:1657
LinearAlgebra::BlockVector system_rhs
Definition: simulator.h:1799
boost::iostreams::stream< TeeDevice > TeeStream
Definition: simulator.h:1609
LateralAveraging< dim > lateral_averaging
Definition: simulator.h:1695
TeeStream iostream_tee_stream
Definition: simulator.h:1612
ConstraintMatrix constraints
Definition: simulator.h:1738
MeshRefinement::Manager< dim > mesh_refinement_manager
Definition: simulator.h:1708
const std::unique_ptr< MaterialModel::Interface< dim > > material_model
Definition: simulator.h:1652
double last_pressure_normalization_adjustment
Definition: simulator.h:1745
void declare_parameters(ParameterHandler &prm)
const std::unique_ptr< GeometryModel::Interface< dim > > geometry_model
Definition: simulator.h:1650
LinearAlgebra::BlockVector operator_split_reaction_vector
Definition: simulator.h:1807
Parameters< dim >::NonlinearSolver NonlinearSolver
Definition: simulator.h:203
DoFHandler< dim > dof_handler
Definition: simulator.h:1723
std::unique_ptr< Assemblers::Manager< dim > > assemblers
Definition: simulator.h:875
bool rebuild_sparsity_and_matrices
Definition: simulator.h:1814
const IntermediaryConstructorAction post_signal_creation
Definition: simulator.h:1586
double global_Omega_diameter
Definition: simulator.h:1705
unsigned int timestep_number
Definition: simulator.h:1675
bool assemble_newton_stokes_matrix
Definition: simulator.h:1816
ConditionalOStream pcout
Definition: simulator.h:1618
unsigned int pre_refinement_step
Definition: simulator.h:1676
TableHandler statistics
Definition: simulator.h:1629
Introspection< dim > introspection
Definition: simulator.h:1597
MPI_Comm mpi_communicator
Definition: simulator.h:1600
LinearAlgebra::BlockSparseMatrix system_preconditioner_matrix
Definition: simulator.h:1794
bool do_pressure_rhs_compatibility_modification
Definition: simulator.h:1760
LinearAlgebra::BlockVector current_linearization_point
Definition: simulator.h:1801
InitialTemperature::Manager< dim > initial_temperature_manager
Definition: simulator.h:1658
LinearAlgebra::BlockVector old_old_solution
Definition: simulator.h:1798
Definition: compat.h:38
LinearAlgebra::BlockVector pressure_shape_function_integrals
Definition: simulator.h:1804
std::ofstream log_file_stream
Definition: simulator.h:1606
std::unique_ptr< LinearAlgebra::PreconditionAMG > Amg_preconditioner
Definition: simulator.h:1811
ConstraintMatrix current_constraints
Definition: simulator.h:1739
TeeDevice iostream_tee_device
Definition: simulator.h:1611
LinearAlgebra::BlockSparseMatrix system_matrix
Definition: simulator.h:1780
std::unique_ptr< MeltHandler< dim > > melt_handler
Definition: simulator.h:1575
bool rebuild_stokes_preconditioner
Definition: simulator.h:1818
const unsigned int compositional_variable
Definition: simulator.h:236
std::unique_ptr< Mapping< dim > > mapping
Definition: simulator.h:1719
const std::unique_ptr< PrescribedStokesSolution::Interface< dim > > prescribed_stokes_solution
Definition: simulator.h:1656
HeatingModel::Manager< dim > heating_model_manager
Definition: simulator.h:1709
SimulatorSignals< dim > signals
Definition: simulator.h:1584
const std::unique_ptr< BoundaryHeatFlux::Interface< dim > > boundary_heat_flux
Definition: simulator.h:1663
const FESystem< dim > finite_element
Definition: simulator.h:1721
Threads::Thread output_statistics_thread
Definition: simulator.h:1640
const std::unique_ptr< InitialTopographyModel::Interface< dim > > initial_topography_model
Definition: simulator.h:1649
LinearAlgebra::BlockVector solution
Definition: simulator.h:1796
std::unique_ptr< NewtonHandler< dim > > newton_handler
Definition: simulator.h:1582
Postprocess::Manager< dim > postprocess_manager
Definition: simulator.h:1725
const std::unique_ptr< GravityModel::Interface< dim > > gravity_model
Definition: simulator.h:1653
double global_volume
Definition: simulator.h:1706
BoundaryComposition::Manager< dim > boundary_composition_manager
Definition: simulator.h:1655
const std::unique_ptr< WorldBuilder::World > world_builder
Definition: simulator.h:1660
std::unique_ptr< LinearAlgebra::PreconditionBase > Mp_preconditioner
Definition: simulator.h:1812
Parameters< dim >::NullspaceRemoval NullspaceRemoval
Definition: simulator.h:208