22 #ifndef _aspect_simulator_h 23 #define _aspect_simulator_h 25 #include <deal.II/base/timer.h> 26 #include <deal.II/base/parameter_handler.h> 27 #include <deal.II/base/conditional_ostream.h> 28 #include <deal.II/base/symmetric_tensor.h> 30 DEAL_II_DISABLE_EXTRA_DIAGNOSTICS
32 #include <deal.II/lac/affine_constraints.h> 34 #include <deal.II/distributed/tria.h> 36 #include <deal.II/dofs/dof_handler.h> 37 #include <deal.II/dofs/dof_tools.h> 39 #include <deal.II/fe/fe_system.h> 40 #include <deal.II/fe/mapping.h> 41 #include <deal.II/base/tensor_function.h> 43 DEAL_II_ENABLE_EXTRA_DIAGNOSTICS
69 #include <boost/iostreams/tee.hpp> 70 #include <boost/iostreams/stream.hpp> 89 template <
int dim,
int velocity_degree>
92 namespace MeshDeformation
147 scalar_moment_of_inertia(numbers::signaling_nan<double>()),
148 scalar_angular_momentum(numbers::signaling_nan<double>()),
149 scalar_rotation(numbers::signaling_nan<double>()),
150 tensor_moment_of_inertia(numbers::signaling_nan<SymmetricTensor<2,dim>>()),
151 tensor_angular_momentum(numbers::signaling_nan<Tensor<1,dim>>()),
152 tensor_rotation(numbers::signaling_nan<Tensor<1,dim>>())
195 Simulator (
const MPI_Comm mpi_communicator,
196 ParameterHandler &prm);
292 const unsigned int compositional_variable = numbers::invalid_unsigned_int);
312 AdvectionField composition (
const unsigned int compositional_variable);
318 is_temperature ()
const;
352 unsigned int field_index()
const;
367 FEValuesExtractors::Scalar scalar_extractor(
const Introspection<dim> &introspection)
const;
443 void setup_introspection ();
456 void set_initial_temperature_and_compositional_fields ();
473 void compute_initial_pressure_field ();
480 void compute_initial_velocity_boundary_constraints (AffineConstraints<double> &constraints);
487 void compute_current_velocity_boundary_constraints (AffineConstraints<double> &constraints);
501 void compute_current_constraints ();
516 double compute_pressure_scaling_factor ()
const;
529 void start_timestep ();
538 void solve_timestep ();
551 void solve_single_advection_single_stokes ();
565 void solve_no_advection_iterated_stokes ();
577 void solve_no_advection_single_stokes ();
590 void solve_first_timestep_only_single_stokes ();
604 void solve_iterated_advection_and_stokes ();
618 void solve_single_advection_iterated_stokes ();
632 void solve_no_advection_iterated_defect_correction_stokes ();
646 void solve_single_advection_iterated_defect_correction_stokes ();
661 void solve_iterated_advection_and_defect_correction_stokes ();
677 void solve_iterated_advection_and_newton_stokes ();
694 void solve_single_advection_iterated_newton_stokes ();
707 void solve_single_advection_no_stokes ();
720 void solve_no_advection_no_stokes ();
730 void build_stokes_preconditioner ();
739 void build_advection_preconditioner (
const AdvectionField &advection_field,
741 const double diagonal_strengthening);
749 void assemble_stokes_system ();
760 double assemble_and_solve_temperature (
const bool compute_initial_residual =
false,
761 double *initial_residual =
nullptr);
775 std::vector<double> assemble_and_solve_composition (
const bool compute_initial_residual =
false,
776 std::vector<double> *initial_residual =
nullptr);
795 double assemble_and_solve_stokes (
const bool compute_initial_residual =
false,
796 double *initial_nonlinear_residual =
nullptr);
810 const bool use_picard);
819 void assemble_advection_system (
const AdvectionField &advection_field);
839 void interpolate_particle_properties (
const AdvectionField &advection_field);
845 void interpolate_particle_properties (
const std::vector<AdvectionField> &advection_fields);
917 std::pair<double,double>
923 std::pair<double,double>
924 solve_stokes_block_gmg ();
955 void refine_mesh (
const unsigned int max_grid_level);
977 void create_snapshot();
992 void resume_from_snapshot();
1000 template <
class Archive>
1001 void serialize (Archive &ar,
const unsigned int version);
1018 Table<2,DoFTools::Coupling>
1019 setup_system_matrix_coupling ()
const;
1028 void setup_system_matrix (
const std::vector<IndexSet> &system_partitioning);
1041 void setup_system_preconditioner (
const std::vector<IndexSet> &system_partitioning);
1074 void set_assemblers ();
1086 void set_advection_assemblers ();
1097 void set_stokes_assemblers ();
1105 void assemble_stokes_preconditioner ();
1115 local_assemble_stokes_preconditioner (
const typename DoFHandler<dim>::active_cell_iterator &cell,
1137 local_assemble_stokes_system (
const typename DoFHandler<dim>::active_cell_iterator &cell,
1159 local_assemble_advection_face_terms(
const AdvectionField &advection_field,
1160 const typename DoFHandler<dim>::active_cell_iterator &cell,
1171 local_assemble_advection_system (
const AdvectionField &advection_field,
1172 const Vector<double> &viscosity_per_cell,
1173 const typename DoFHandler<dim>::active_cell_iterator &cell,
1185 copy_local_to_global_advection_system (
const AdvectionField &advection_field,
1223 template <
typename T>
1224 void get_artificial_viscosity (Vector<T> &viscosity_per_cell,
1226 const bool skip_interior_cells =
false)
const;
1239 void compute_Vs_anomaly(Vector<float> &values)
const;
1255 void compute_Vp_anomaly(Vector<float> &values)
const;
1343 void denormalize_pressure(
const double pressure_adjustment,
1355 void apply_limiter_to_dg_solutions (
const AdvectionField &advection_field);
1380 void compute_reactions ();
1391 void update_solution_vectors_with_reaction_results (
const unsigned int block_index,
1405 void initialize_current_linearization_point ();
1426 void interpolate_material_output_into_advection_field (
const AdvectionField &adv_field);
1436 void interpolate_onto_velocity_system(
const TensorFunction<1,dim> &func,
1454 void setup_nullspace_constraints(AffineConstraints<double> &constraints);
1486 compute_net_angular_momentum(
const bool use_constant_density,
1488 const bool limit_to_top_faces =
false)
const;
1506 void remove_net_angular_momentum(
const bool use_constant_density,
1509 const bool limit_to_top_faces =
false);
1518 void replace_outflow_boundary_ids(
const unsigned int boundary_id_offset);
1527 void restore_outflow_boundary_ids(
const unsigned int boundary_id_offset);
1542 void remove_net_linear_momentum(
const bool use_constant_density,
1567 double get_entropy_variation (
const double average_field,
1578 std::pair<double,double>
1579 get_extrapolated_advection_field_range (
const AdvectionField &advection_field)
const;
1589 void exchange_refinement_flags();
1600 void maybe_write_timing_output ()
const;
1609 bool maybe_write_checkpoint (
const time_t last_checkpoint_time,
1610 const bool force_writing_checkpoint);
1628 bool maybe_do_initial_refinement (
const unsigned int max_refinement_level);
1638 void maybe_refine_mesh (
const double new_time_step,
1639 unsigned int &max_refinement_level);
1648 void advance_time (
const double step_size);
1659 const double global_u_infty,
1660 const double global_field_variation,
1661 const double average_field,
1662 const double global_entropy_variation,
1663 const double cell_diameter,
1676 const double average_field,
1678 double &max_residual,
1679 double &max_velocity,
1680 double &max_density,
1681 double &max_specific_heat,
1682 double &conductivity)
const;
1699 stokes_matrix_depends_on_solution ()
const;
1714 check_consistency_of_formulation ();
1727 check_consistency_of_boundary_conditions ()
const;
1744 compute_Eisenstat_Walker_linear_tolerance(
const bool EisenstatWalkerChoiceOne,
1745 const double maximum_linear_stokes_solver_tolerance,
1746 const double linear_stokes_solver_tolerance,
1747 const double stokes_residual,
1748 const double newton_residual,
1749 const double newton_residual_old);
1760 void output_statistics();
1774 compute_initial_stokes_residual();
1825 using TeeDevice = boost::iostreams::tee_device<std::ostream, std::ofstream>;
1919 #ifdef ASPECT_WITH_WORLD_BUILDER 1931 std::shared_ptr<WorldBuilder::World> world_builder;
2121 friend class boost::serialization::access;
2126 template <
int dimension,
int velocity_degree>
The NullspaceRemoval struct.
unsigned int nonlinear_iteration
BoundaryVelocity::Manager< dim > boundary_velocity_manager
const std::unique_ptr< AdiabaticConditions::Interface< dim > > adiabatic_conditions
std::shared_ptr< InitialComposition::Manager< dim > > initial_composition_manager
void write_plugin_graph(std::ostream &output_stream)
parallel::distributed::Triangulation< dim > triangulation
TimerOutput computing_timer
BoundaryTemperature::Manager< dim > boundary_temperature_manager
boost::iostreams::tee_device< std::ostream, std::ofstream > TeeDevice
::TrilinosWrappers::MPI::BlockVector BlockVector
double scalar_moment_of_inertia
Tensor< 1, dim > tensor_rotation
Tensor< 1, dim > tensor_angular_momentum
AffineConstraints< double > constraints
bool assemble_newton_stokes_system
double scalar_angular_momentum
LinearAlgebra::BlockVector old_solution
const IntermediaryConstructorAction post_geometry_model_creation_action
std::unique_ptr< VolumeOfFluidHandler< dim > > volume_of_fluid_handler
bool rebuild_stokes_matrix
Parameters< dim > parameters
SymmetricTensor< 2, dim > tensor_moment_of_inertia
AffineConstraints< double > current_constraints
StructuredDataLookup< dim > DEAL_II_DEPRECATED
std::map< types::boundary_id, std::unique_ptr< BoundaryTraction::Interface< dim > > > boundary_traction
std::unique_ptr< MeshDeformation::MeshDeformationHandler< dim > > mesh_deformation
LinearAlgebra::BlockVector system_rhs
std::pair< double, double > stokes_residuals
std::thread output_statistics_thread
LateralAveraging< dim > lateral_averaging
TeeStream iostream_tee_stream
::TrilinosWrappers::MPI::Vector Vector
MeshRefinement::Manager< dim > mesh_refinement_manager
std::size_t statistics_last_write_size
const std::unique_ptr< MaterialModel::Interface< dim > > material_model
::Particles::ParticleHandler< dim > particle_handler_copy
double last_pressure_normalization_adjustment
::TrilinosWrappers::BlockSparseMatrix BlockSparseMatrix
double total_walltime_until_last_snapshot
void declare_parameters(ParameterHandler &prm)
const std::unique_ptr< GeometryModel::Interface< dim > > geometry_model
LinearAlgebra::BlockVector operator_split_reaction_vector
DoFHandler< dim > dof_handler
std::unique_ptr< StokesMatrixFreeHandler< dim > > stokes_matrix_free
std::unique_ptr< Assemblers::Manager< dim > > assemblers
bool rebuild_sparsity_and_matrices
const IntermediaryConstructorAction post_signal_creation
double global_Omega_diameter
unsigned int timestep_number
bool assemble_newton_stokes_matrix
unsigned int pre_refinement_step
TimeStepping::Manager< dim > time_stepping_manager
std::size_t statistics_last_hash
Introspection< dim > introspection
MPI_Comm mpi_communicator
typename Parameters< dim >::NullspaceRemoval NullspaceRemoval
LinearAlgebra::BlockSparseMatrix system_preconditioner_matrix
bool do_pressure_rhs_compatibility_modification
LinearAlgebra::BlockVector current_linearization_point
LinearAlgebra::BlockVector old_old_solution
LinearAlgebra::BlockVector pressure_shape_function_integrals
std::ofstream log_file_stream
std::unique_ptr< LinearAlgebra::PreconditionAMG > Amg_preconditioner
const FieldType field_type
std::shared_ptr< InitialTemperature::Manager< dim > > initial_temperature_manager
TeeDevice iostream_tee_device
boost::iostreams::stream< TeeDevice > TeeStream
LinearAlgebra::BlockSparseMatrix system_matrix
std::unique_ptr< MeltHandler< dim > > melt_handler
bool rebuild_stokes_preconditioner
const unsigned int compositional_variable
std::unique_ptr< Mapping< dim > > mapping
const std::unique_ptr< PrescribedStokesSolution::Interface< dim > > prescribed_stokes_solution
double switch_initial_residual
HeatingModel::Manager< dim > heating_model_manager
SimulatorSignals< dim > signals
const std::unique_ptr< BoundaryHeatFlux::Interface< dim > > boundary_heat_flux
const FESystem< dim > finite_element
bool simulator_is_past_initialization
std::unique_ptr< Particle::World< dim > > particle_world
const std::unique_ptr< InitialTopographyModel::Interface< dim > > initial_topography_model
LinearAlgebra::BlockVector solution
std::unique_ptr< NewtonHandler< dim > > newton_handler
typename Parameters< dim >::NonlinearSolver NonlinearSolver
Postprocess::Manager< dim > postprocess_manager
const std::unique_ptr< GravityModel::Interface< dim > > gravity_model
double newton_residual_for_derivative_scaling_factor
BoundaryComposition::Manager< dim > boundary_composition_manager
std::unique_ptr< LinearAlgebra::PreconditionBase > Mp_preconditioner
::TrilinosWrappers::PreconditionILU PreconditionILU