ASPECT
simulator.h
Go to the documentation of this file.
1 /*
2  Copyright (C) 2011 - 2018 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 
30 DEAL_II_DISABLE_EXTRA_DIAGNOSTICS
31 
33 
35 
36 #include <deal.II/fe/fe_system.h>
37 #include <deal.II/fe/mapping.h>
39 
40 DEAL_II_ENABLE_EXTRA_DIAGNOSTICS
41 
42 #if !DEAL_II_VERSION_GTE(9,1,0)
44 #else
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  namespace internal
101  {
102  namespace Assembly
103  {
104  namespace Scratch
105  {
106  template <int dim> struct StokesPreconditioner;
107  template <int dim> struct StokesSystem;
108  template <int dim> struct AdvectionSystem;
109  }
110 
111  namespace CopyData
112  {
113  template <int dim> struct StokesPreconditioner;
114  template <int dim> struct StokesSystem;
115  template <int dim> struct AdvectionSystem;
116  }
117  }
118  }
119 
120  namespace Assemblers
121  {
122  template <int dim> class Interface;
123  template <int dim> class Manager;
124  }
125 
133  template <int dim>
134  class Simulator
135  {
136  public:
150  Simulator (const MPI_Comm mpi_communicator,
151  ParameterHandler &prm);
152 
157  ~Simulator ();
158 
170  static
172 
181  void run ();
182 
194  void
195  write_plugin_graph (std::ostream &output_stream) const;
196 
201 
206 
207 
215  {
220  enum FieldType { temperature_field, compositional_field };
221 
227 
233  const unsigned int compositional_variable;
234 
246  AdvectionField (const FieldType field_type,
247  const unsigned int compositional_variable = numbers::invalid_unsigned_int);
248 
256  static
258 
266  static
267  AdvectionField composition (const unsigned int compositional_variable);
268 
272  bool
273  is_temperature () const;
274 
279  bool
280  is_discontinuous (const Introspection<dim> &introspection) const;
281 
287  advection_method (const Introspection<dim> &introspection) const;
288 
293  unsigned int component_index(const Introspection<dim> &introspection) const;
294 
299  unsigned int block_index(const Introspection<dim> &introspection) const;
300 
307  unsigned int field_index() const;
308 
314  unsigned int base_element(const Introspection<dim> &introspection) const;
315 
322  FEValuesExtractors::Scalar scalar_extractor(const Introspection<dim> &introspection) const;
323 
328  unsigned int polynomial_degree(const Introspection<dim> &introspection) const;
329  };
330 
331 
332  private:
333 
334 
352  {
353  IntermediaryConstructorAction (const std::function<void ()> &action);
354  };
355 
370  void setup_dofs ();
371 
380  void setup_introspection ();
381 
393  void set_initial_temperature_and_compositional_fields ();
394 
410  void compute_initial_pressure_field ();
411 
417  void compute_initial_velocity_boundary_constraints (ConstraintMatrix &constraints);
418 
424  void compute_current_velocity_boundary_constraints (ConstraintMatrix &constraints);
425 
438  void compute_current_constraints ();
439 
451  void start_timestep ();
452 
460  void solve_timestep ();
461 
473  void solve_single_advection_single_stokes ();
474 
487  void solve_no_advection_iterated_stokes ();
488 
500  void solve_first_timestep_only_single_stokes ();
501 
514  void solve_iterated_advection_and_stokes ();
515 
528  void solve_single_advection_iterated_stokes ();
529 
544  void solve_iterated_advection_and_newton_stokes ();
545 
557  void solve_single_advection_no_stokes ();
558 
570  void solve_no_advection_no_stokes ();
571 
580  void build_stokes_preconditioner ();
581 
589  void build_advection_preconditioner (const AdvectionField &advection_field,
591 
598  void assemble_stokes_system ();
599 
609  double assemble_and_solve_temperature (const bool compute_initial_residual = false,
610  double *initial_residual = nullptr);
611 
624  std::vector<double> assemble_and_solve_composition (const bool compute_initial_residual = false,
625  std::vector<double> *initial_residual = nullptr);
626 
644  double assemble_and_solve_stokes (const bool compute_initial_residual = false,
645  double *initial_residual = nullptr);
646 
654  void assemble_advection_system (const AdvectionField &advection_field);
655 
666  double solve_advection (const AdvectionField &advection_field);
667 
671  void interpolate_particle_properties (const AdvectionField &advection_field);
672 
743  std::pair<double,double>
744  solve_stokes ();
745 
758  void postprocess ();
759 
775  void refine_mesh (const unsigned int max_grid_level);
776 
797  void create_snapshot();
798 
812  void resume_from_snapshot();
813 
820  template <class Archive>
821  void serialize (Archive &ar, const unsigned int version);
837  void setup_system_matrix (const std::vector<IndexSet> &system_partitioning);
838 
850  void setup_system_preconditioner (const std::vector<IndexSet> &system_partitioning);
851 
873  std::unique_ptr<Assemblers::Manager<dim> > assemblers;
874 
883  void set_assemblers ();
884 
895  void set_default_assemblers ();
896 
903  void assemble_stokes_preconditioner ();
904 
912  void
913  local_assemble_stokes_preconditioner (const typename DoFHandler<dim>::active_cell_iterator &cell,
916 
924  void
925  copy_local_to_global_stokes_preconditioner (const internal::Assembly::CopyData::StokesPreconditioner<dim> &data);
926 
934  void
935  local_assemble_stokes_system (const typename DoFHandler<dim>::active_cell_iterator &cell,
938 
946  void
947  copy_local_to_global_stokes_system (const internal::Assembly::CopyData::StokesSystem<dim> &data);
948 
956  void
957  local_assemble_advection_face_terms(const AdvectionField &advection_field,
958  const typename DoFHandler<dim>::active_cell_iterator &cell,
968  void
969  local_assemble_advection_system (const AdvectionField &advection_field,
970  const Vector<double> &viscosity_per_cell,
971  const typename DoFHandler<dim>::active_cell_iterator &cell,
974 
982  void
983  copy_local_to_global_advection_system (const AdvectionField &advection_field,
985 
1010  void make_pressure_rhs_compatible(LinearAlgebra::BlockVector &vector);
1011 
1021  template <typename T>
1022  void get_artificial_viscosity (Vector<T> &viscosity_per_cell,
1023  const AdvectionField &advection_field,
1024  const bool skip_interior_cells = false) const;
1025 
1037  void compute_Vs_anomaly(Vector<float> &values) const;
1038 
1053  void compute_Vp_anomaly(Vector<float> &values) const;
1054 
1095  double normalize_pressure(LinearAlgebra::BlockVector &vector) const;
1096 
1141  void denormalize_pressure(const double pressure_adjustment,
1143  const LinearAlgebra::BlockVector &relevant_vector) const;
1144 
1153  void apply_limiter_to_dg_solutions (const AdvectionField &advection_field);
1154 
1155 
1178  void compute_reactions ();
1179 
1180 
1203  void interpolate_material_output_into_compositional_field (const unsigned int compositional_index);
1204 
1205 
1213  void interpolate_onto_velocity_system(const TensorFunction<1,dim> &func,
1214  LinearAlgebra::Vector &vec);
1215 
1216 
1231  void setup_nullspace_constraints(ConstraintMatrix &constraints);
1232 
1233 
1245  void remove_nullspace(LinearAlgebra::BlockVector &relevant_dst,
1246  LinearAlgebra::BlockVector &tmp_distributed_stokes);
1247 
1261  void remove_net_angular_momentum( const bool use_constant_density,
1262  LinearAlgebra::BlockVector &relevant_dst,
1263  LinearAlgebra::BlockVector &tmp_distributed_stokes);
1264 
1272  void replace_outflow_boundary_ids(const unsigned int boundary_id_offset);
1273 
1281  void restore_outflow_boundary_ids(const unsigned int boundary_id_offset);
1282 
1296  void remove_net_linear_momentum( const bool use_constant_density,
1297  LinearAlgebra::BlockVector &relevant_dst,
1298  LinearAlgebra::BlockVector &tmp_distributed_stokes);
1299 
1307  double get_maximal_velocity (const LinearAlgebra::BlockVector &solution) const;
1308 
1321  double get_entropy_variation (const double average_value,
1322  const AdvectionField &advection_field) const;
1323 
1332  std::pair<double,double>
1333  get_extrapolated_advection_field_range (const AdvectionField &advection_field) const;
1334 
1343  void maybe_write_timing_output () const;
1344 
1352  bool maybe_write_checkpoint (const time_t last_checkpoint_time,
1353  const std::pair<bool,bool> termination_output);
1354 
1371  bool maybe_do_initial_refinement (const unsigned int max_refinement_level);
1372 
1381  void maybe_refine_mesh (const double new_time_step,
1382  unsigned int &max_refinement_level);
1383 
1394  double compute_time_step () const;
1395 
1403  double
1404  compute_viscosity(internal::Assembly::Scratch::AdvectionSystem<dim> &scratch,
1405  const double global_u_infty,
1406  const double global_T_variation,
1407  const double average_temperature,
1408  const double global_entropy_variation,
1409  const double cell_diameter,
1410  const AdvectionField &advection_field) const;
1411 
1420  void
1421  compute_advection_system_residual(internal::Assembly::Scratch::AdvectionSystem<dim> &scratch,
1422  const double average_field,
1423  const AdvectionField &advection_field,
1424  double &max_residual,
1425  double &max_velocity,
1426  double &max_density,
1427  double &max_specific_heat,
1428  double &conductivity) const;
1429 
1457  void
1458  compute_material_model_input_values (const LinearAlgebra::BlockVector &input_solution,
1459  const FEValuesBase<dim,dim> &input_finite_element_values,
1460  const typename DoFHandler<dim>::active_cell_iterator &cell,
1461  const bool compute_strainrate,
1462  MaterialModel::MaterialModelInputs<dim> &material_model_inputs) const;
1463 
1464 
1479  bool
1480  stokes_matrix_depends_on_solution () const;
1481 
1494  void
1495  check_consistency_of_formulation ();
1496 
1507  void
1508  check_consistency_of_boundary_conditions () const;
1509 
1513  double
1514  compute_initial_newton_residual (const LinearAlgebra::BlockVector &linearized_stokes_initial_guess);
1515 
1524  double
1525  compute_Eisenstat_Walker_linear_tolerance(const bool EisenstatWalkerChoiceOne,
1526  const double maximum_linear_stokes_solver_tolerance,
1527  const double linear_stokes_solver_tolerance,
1528  const double stokes_residual,
1529  const double newton_residual,
1530  const double newton_residual_old);
1531 
1541  void output_statistics();
1542 
1554  double
1555  compute_initial_stokes_residual();
1556 
1567 
1573  std::unique_ptr<MeltHandler<dim> > melt_handler;
1574 
1580  std::unique_ptr<NewtonHandler<dim> > newton_handler;
1581 
1585 
1586 
1588 
1593  std::ofstream log_file_stream;
1594 
1595  typedef boost::iostreams::tee_device<std::ostream, std::ofstream> TeeDevice;
1596  typedef boost::iostreams::stream< TeeDevice > TeeStream;
1597 
1600 
1606 
1617 
1619 
1636  const std::unique_ptr<InitialTopographyModel::Interface<dim> > initial_topography_model;
1637  const std::unique_ptr<GeometryModel::Interface<dim> > geometry_model;
1639  const std::unique_ptr<MaterialModel::Interface<dim> > material_model;
1640  const std::unique_ptr<GravityModel::Interface<dim> > gravity_model;
1643  const std::unique_ptr<PrescribedStokesSolution::Interface<dim> > prescribed_stokes_solution;
1646  const std::unique_ptr<AdiabaticConditions::Interface<dim> > adiabatic_conditions;
1647  const std::unique_ptr<WorldBuilder::World> world_builder;
1649  std::map<types::boundary_id,std::shared_ptr<BoundaryTraction::Interface<dim> > > boundary_traction;
1650  const std::unique_ptr<BoundaryHeatFlux::Interface<dim> > boundary_heat_flux;
1651 
1659  double time;
1660  double time_step;
1662  unsigned int timestep_number;
1663  unsigned int pre_refinement_step;
1664  unsigned int nonlinear_iteration;
1694 
1697 
1706  std::unique_ptr<Mapping<dim> > mapping;
1707 
1709 
1711 
1713 
1725  ConstraintMatrix constraints;
1726  ConstraintMatrix current_constraints;
1727 
1733 
1740 
1748 
1768 
1782 
1787 
1789 
1790  // only used if is_compressible()
1792 
1793  // only used if operator split is enabled
1795 
1796 
1797 
1798  std::unique_ptr<LinearAlgebra::PreconditionAMG> Amg_preconditioner;
1799  std::unique_ptr<LinearAlgebra::PreconditionBase> Mp_preconditioner;
1800 
1806 
1811  private:
1812 
1818  std::unique_ptr<FreeSurfaceHandler<dim> > free_surface;
1819 
1820  friend class boost::serialization::access;
1821  friend class SimulatorAccess<dim>;
1822  friend class FreeSurfaceHandler<dim>; // FreeSurfaceHandler needs access to the internals of the Simulator
1823  friend struct Parameters<dim>;
1824  };
1825 }
1826 
1827 
1828 #endif
The NullspaceRemoval struct.
Definition: parameters.h:81
unsigned int nonlinear_iteration
Definition: simulator.h:1664
BoundaryVelocity::Manager< dim > boundary_velocity_manager
Definition: simulator.h:1648
const std::unique_ptr< AdiabaticConditions::Interface< dim > > adiabatic_conditions
Definition: simulator.h:1646
static const unsigned int invalid_unsigned_int
TerminationCriteria::Manager< dim > termination_manager
Definition: simulator.h:1673
void write_plugin_graph(std::ostream &output_stream)
parallel::distributed::Triangulation< dim > triangulation
Definition: simulator.h:1691
TimerOutput computing_timer
Definition: simulator.h:1618
BoundaryTemperature::Manager< dim > boundary_temperature_manager
Definition: simulator.h:1641
double pressure_scaling
Definition: simulator.h:1739
double old_time_step
Definition: simulator.h:1661
bool assemble_newton_stokes_system
Definition: simulator.h:1804
LinearAlgebra::BlockVector old_solution
Definition: simulator.h:1784
const IntermediaryConstructorAction post_geometry_model_creation_action
Definition: simulator.h:1638
bool rebuild_stokes_matrix
Definition: simulator.h:1802
Parameters< dim > parameters
Definition: simulator.h:1566
std::unique_ptr< FreeSurfaceHandler< dim > > free_surface
Definition: simulator.h:1818
boost::iostreams::tee_device< std::ostream, std::ofstream > TeeDevice
Definition: simulator.h:1595
InitialComposition::Manager< dim > initial_composition_manager
Definition: simulator.h:1644
values
LinearAlgebra::BlockVector system_rhs
Definition: simulator.h:1786
boost::iostreams::stream< TeeDevice > TeeStream
Definition: simulator.h:1596
LateralAveraging< dim > lateral_averaging
Definition: simulator.h:1682
TeeStream iostream_tee_stream
Definition: simulator.h:1599
ConstraintMatrix constraints
Definition: simulator.h:1725
MeshRefinement::Manager< dim > mesh_refinement_manager
Definition: simulator.h:1695
const std::unique_ptr< MaterialModel::Interface< dim > > material_model
Definition: simulator.h:1639
double last_pressure_normalization_adjustment
Definition: simulator.h:1732
void declare_parameters(ParameterHandler &prm)
const std::unique_ptr< GeometryModel::Interface< dim > > geometry_model
Definition: simulator.h:1637
LinearAlgebra::BlockVector operator_split_reaction_vector
Definition: simulator.h:1794
Parameters< dim >::NonlinearSolver NonlinearSolver
Definition: simulator.h:200
DoFHandler< dim > dof_handler
Definition: simulator.h:1710
double average_value(const std::vector< double > &volume_fractions, const std::vector< double > &parameter_values, const CompositionalAveragingOperation &average_type)
std::unique_ptr< Assemblers::Manager< dim > > assemblers
Definition: simulator.h:873
bool rebuild_sparsity_and_matrices
Definition: simulator.h:1801
const IntermediaryConstructorAction post_signal_creation
Definition: simulator.h:1583
double global_Omega_diameter
Definition: simulator.h:1692
unsigned int timestep_number
Definition: simulator.h:1662
bool assemble_newton_stokes_matrix
Definition: simulator.h:1803
ConditionalOStream pcout
Definition: simulator.h:1605
unsigned int pre_refinement_step
Definition: simulator.h:1663
TableHandler statistics
Definition: simulator.h:1616
Introspection< dim > introspection
Definition: simulator.h:1584
MPI_Comm mpi_communicator
Definition: simulator.h:1587
LinearAlgebra::BlockSparseMatrix system_preconditioner_matrix
Definition: simulator.h:1781
bool do_pressure_rhs_compatibility_modification
Definition: simulator.h:1747
LinearAlgebra::BlockVector current_linearization_point
Definition: simulator.h:1788
InitialTemperature::Manager< dim > initial_temperature_manager
Definition: simulator.h:1645
LinearAlgebra::BlockVector old_old_solution
Definition: simulator.h:1785
Definition: compat.h:38
LinearAlgebra::BlockVector pressure_shape_function_integrals
Definition: simulator.h:1791
std::ofstream log_file_stream
Definition: simulator.h:1593
std::unique_ptr< LinearAlgebra::PreconditionAMG > Amg_preconditioner
Definition: simulator.h:1798
ConstraintMatrix current_constraints
Definition: simulator.h:1726
TeeDevice iostream_tee_device
Definition: simulator.h:1598
LinearAlgebra::BlockSparseMatrix system_matrix
Definition: simulator.h:1767
std::unique_ptr< MeltHandler< dim > > melt_handler
Definition: simulator.h:1573
bool rebuild_stokes_preconditioner
Definition: simulator.h:1805
const unsigned int compositional_variable
Definition: simulator.h:233
std::unique_ptr< Mapping< dim > > mapping
Definition: simulator.h:1706
const std::unique_ptr< PrescribedStokesSolution::Interface< dim > > prescribed_stokes_solution
Definition: simulator.h:1643
HeatingModel::Manager< dim > heating_model_manager
Definition: simulator.h:1696
void run(const std::vector< std::vector< Iterator >> &colored_iterators, Worker worker, Copier copier, const ScratchData &sample_scratch_data, const CopyData &sample_copy_data, const unsigned int queue_length=2 *MultithreadInfo::n_threads(), const unsigned int chunk_size=8)
std::map< types::boundary_id, std::shared_ptr< BoundaryTraction::Interface< dim > > > boundary_traction
Definition: simulator.h:1649
SimulatorSignals< dim > signals
Definition: simulator.h:1582
const std::unique_ptr< BoundaryHeatFlux::Interface< dim > > boundary_heat_flux
Definition: simulator.h:1650
const FESystem< dim > finite_element
Definition: simulator.h:1708
Threads::Thread output_statistics_thread
Definition: simulator.h:1627
const std::unique_ptr< InitialTopographyModel::Interface< dim > > initial_topography_model
Definition: simulator.h:1636
LinearAlgebra::BlockVector solution
Definition: simulator.h:1783
std::unique_ptr< NewtonHandler< dim > > newton_handler
Definition: simulator.h:1580
Postprocess::Manager< dim > postprocess_manager
Definition: simulator.h:1712
const std::unique_ptr< GravityModel::Interface< dim > > gravity_model
Definition: simulator.h:1640
double global_volume
Definition: simulator.h:1693
BoundaryComposition::Manager< dim > boundary_composition_manager
Definition: simulator.h:1642
const std::unique_ptr< WorldBuilder::World > world_builder
Definition: simulator.h:1647
std::unique_ptr< LinearAlgebra::PreconditionBase > Mp_preconditioner
Definition: simulator.h:1799
Parameters< dim >::NullspaceRemoval NullspaceRemoval
Definition: simulator.h:205