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> 87 namespace StokesSolver
93 template <
int dim,
int velocity_degree>
96 namespace MeshDeformation
151 scalar_moment_of_inertia(numbers::signaling_nan<double>()),
152 scalar_angular_momentum(numbers::signaling_nan<double>()),
153 scalar_rotation(numbers::signaling_nan<double>()),
154 tensor_moment_of_inertia(numbers::signaling_nan<SymmetricTensor<2,dim>>()),
155 tensor_angular_momentum(numbers::signaling_nan<Tensor<1,dim>>()),
156 tensor_rotation(numbers::signaling_nan<Tensor<1,dim>>())
179 "Nonlinear solver failed to converge in the prescribed number of steps. " 180 "Consider changing `Max nonlinear iterations` or `Nonlinear solver failure " 207 Simulator (
const MPI_Comm mpi_communicator,
208 ParameterHandler &prm);
304 const unsigned int compositional_variable = numbers::invalid_unsigned_int);
324 AdvectionField composition (
const unsigned int compositional_variable);
330 is_temperature ()
const;
365 unsigned int sparsity_pattern_block_index(
const Introspection<dim> &introspection)
const;
373 unsigned int field_index()
const;
388 FEValuesExtractors::Scalar scalar_extractor(
const Introspection<dim> &introspection)
const;
471 void setup_introspection ();
484 void set_initial_temperature_and_compositional_fields ();
501 void compute_initial_pressure_field ();
508 void compute_initial_velocity_boundary_constraints (AffineConstraints<double> &constraints);
515 void compute_current_velocity_boundary_constraints (AffineConstraints<double> &constraints);
529 void compute_current_constraints ();
544 double compute_pressure_scaling_factor ()
const;
557 void start_timestep ();
566 void solve_timestep ();
579 void solve_single_advection_single_stokes ();
593 void solve_no_advection_iterated_stokes ();
605 void solve_no_advection_single_stokes ();
618 void solve_first_timestep_only_single_stokes ();
632 void solve_iterated_advection_and_stokes ();
646 void solve_single_advection_iterated_stokes ();
660 void solve_no_advection_iterated_defect_correction_stokes ();
674 void solve_single_advection_iterated_defect_correction_stokes ();
689 void solve_iterated_advection_and_defect_correction_stokes ();
709 void solve_iterated_advection_and_newton_stokes (
bool use_newton_iterations);
730 void solve_single_advection_and_iterated_newton_stokes (
bool use_newton_iterations);
743 void solve_single_advection_no_stokes ();
756 void solve_no_advection_no_stokes ();
766 void build_stokes_preconditioner ();
775 void build_advection_preconditioner (
const AdvectionField &advection_field,
777 const double diagonal_strengthening);
785 void assemble_stokes_system ();
801 double assemble_and_solve_temperature (
const double &initial_residual = 0,
802 double *residual =
nullptr);
819 std::vector<double> assemble_and_solve_composition (
const std::vector<double> &initial_residual = {},
820 std::vector<double> *residual =
nullptr);
836 double assemble_and_solve_stokes (
const double &initial_nonlinear_residual = 0,
837 double *nonlinear_residual =
nullptr);
852 const bool use_picard);
861 void assemble_advection_system (
const AdvectionField &advection_field);
879 void interpolate_particle_properties (
const std::vector<AdvectionField> &advection_fields);
957 std::pair<double,double>
989 void refine_mesh (
const unsigned int max_grid_level);
1011 void create_snapshot();
1026 void resume_from_snapshot();
1034 template <
class Archive>
1035 void serialize (Archive &ar,
const unsigned int version);
1052 Table<2,DoFTools::Coupling>
1053 setup_system_matrix_coupling ()
const;
1062 void setup_system_matrix (
const std::vector<IndexSet> &system_partitioning);
1075 void setup_system_preconditioner (
const std::vector<IndexSet> &system_partitioning);
1108 void set_assemblers ();
1120 void set_advection_assemblers ();
1131 void set_stokes_assemblers ();
1139 void assemble_stokes_preconditioner ();
1149 local_assemble_stokes_preconditioner (
const typename DoFHandler<dim>::active_cell_iterator &cell,
1171 local_assemble_stokes_system (
const typename DoFHandler<dim>::active_cell_iterator &cell,
1193 local_assemble_advection_face_terms(
const AdvectionField &advection_field,
1194 const typename DoFHandler<dim>::active_cell_iterator &cell,
1205 local_assemble_advection_system (
const AdvectionField &advection_field,
1206 const Vector<double> &viscosity_per_cell,
1207 const typename DoFHandler<dim>::active_cell_iterator &cell,
1219 copy_local_to_global_advection_system (
const AdvectionField &advection_field,
1257 template <
typename T>
1258 void get_artificial_viscosity (Vector<T> &viscosity_per_cell,
1260 const bool skip_interior_cells =
false)
const;
1273 void compute_Vs_anomaly(Vector<float> &values)
const;
1289 void compute_Vp_anomaly(Vector<float> &values)
const;
1364 void denormalize_pressure(
const double pressure_adjustment,
1375 void apply_limiter_to_dg_solutions (
const AdvectionField &advection_field);
1393 void compute_unique_advection_support_points (
const std::vector<AdvectionField> &advection_fields,
1394 std::vector<Point<dim>> &unique_support_points,
1395 std::vector<std::vector<unsigned int>> &support_point_index_by_field)
const;
1419 void compute_reactions ();
1430 void update_solution_vectors_with_reaction_results (
const unsigned int block_index,
1444 void initialize_current_linearization_point ();
1465 void interpolate_material_output_into_advection_field (
const std::vector<AdvectionField> &adv_field);
1475 void interpolate_onto_velocity_system(
const TensorFunction<1,dim> &func,
1493 void setup_nullspace_constraints(AffineConstraints<double> &constraints);
1526 compute_net_angular_momentum(
const bool use_constant_density,
1528 const bool limit_to_top_faces =
false)
const;
1546 void remove_net_angular_momentum(
const bool use_constant_density,
1549 const bool limit_to_top_faces =
false)
const;
1558 void replace_outflow_boundary_ids(
const unsigned int boundary_id_offset);
1567 void restore_outflow_boundary_ids(
const unsigned int boundary_id_offset);
1582 void remove_net_linear_momentum(
const bool use_constant_density,
1607 double get_entropy_variation (
const double average_field,
1618 std::pair<double,double>
1619 get_extrapolated_advection_field_range (
const AdvectionField &advection_field)
const;
1629 void exchange_refinement_flags();
1640 void maybe_write_timing_output ()
const;
1649 bool maybe_write_checkpoint (
const time_t last_checkpoint_time,
1650 const bool force_writing_checkpoint);
1668 bool maybe_do_initial_refinement (
const unsigned int max_refinement_level);
1678 void maybe_refine_mesh (
const double new_time_step,
1679 unsigned int &max_refinement_level);
1688 void advance_time (
const double step_size);
1699 const double global_u_infty,
1700 const double global_field_variation,
1701 const double average_field,
1702 const double global_entropy_variation,
1703 const double cell_diameter,
1716 const double average_field,
1718 double &max_residual,
1719 double &max_velocity,
1720 double &max_density,
1721 double &max_specific_heat,
1722 double &conductivity)
const;
1739 stokes_matrix_depends_on_solution ()
const;
1751 stokes_A_block_is_symmetric ()
const;
1766 check_consistency_of_formulation ();
1774 select_default_solver_and_averaging ();
1787 check_consistency_of_boundary_conditions ()
const;
1793 compute_initial_newton_residual ();
1804 compute_Eisenstat_Walker_linear_tolerance(
const bool EisenstatWalkerChoiceOne,
1805 const double maximum_linear_stokes_solver_tolerance,
1806 const double linear_stokes_solver_tolerance,
1807 const double stokes_residual,
1808 const double newton_residual,
1809 const double newton_residual_old);
1820 void output_statistics();
1834 compute_initial_stokes_residual();
1885 using TeeDevice = boost::iostreams::tee_device<std::ostream, std::ofstream>;
1979 #ifdef ASPECT_WITH_WORLD_BUILDER 1991 std::shared_ptr<WorldBuilder::World> world_builder;
2188 friend class boost::serialization::access;
2193 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
std::unique_ptr< StokesSolver::Direct< dim > > stokes_direct
::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
std::unique_ptr< MeshDeformation::MeshDeformationHandler< dim > > mesh_deformation
LinearAlgebra::BlockVector system_rhs
std::pair< double, double > stokes_residuals
std::vector< Particle::Manager< dim > > particle_managers
std::thread output_statistics_thread
LateralAveraging< dim > lateral_averaging
DeclExceptionMsg(ExcNonlinearSolverNoConvergence, "Nonlinear solver failed to converge in the prescribed number of steps. " "Consider changing `Max nonlinear iterations` or `Nonlinear solver failure " "strategy`.")
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
unsigned int nonlinear_solver_failures
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
BoundaryTraction::Manager< dim > boundary_traction_manager
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
LinearAlgebra::BlockVector inverse_lumped_mass_matrix
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
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