239 *
using namespace dealii::LinearAlgebraPETSc;
242 *
using namespace dealii::LinearAlgebraTrilinos;
292 * <a name=
"step_55-Linearsolversandpreconditioners"></a>
319 *
template <
class Matrix,
class Preconditioner>
325 *
template <
typename VectorType>
326 *
void vmult(VectorType &dst,
const VectorType &src)
const;
334 *
template <
class Matrix,
class Preconditioner>
339 *
, preconditioner(preconditioner)
344 *
template <
class Matrix,
class Preconditioner>
345 *
template <
typename VectorType>
348 *
const VectorType &src)
const
350 *
SolverControl solver_control(src.size(), 1e-8 * src.l2_norm());
356 *
cg.solve(*matrix, dst, src, preconditioner);
358 *
catch (std::exception &e)
360 *
Assert(
false, ExcMessage(
e.what()));
371 *
template <
class PreconditionerA,
class PreconditionerS>
378 *
void vmult(LA::MPI::BlockVector &dst,
379 *
const LA::MPI::BlockVector &src)
const;
386 *
template <
class PreconditionerA,
class PreconditionerS>
395 *
template <
class PreconditionerA,
class PreconditionerS>
397 *
LA::MPI::BlockVector &dst,
398 *
const LA::MPI::BlockVector &src)
const
409 * <a name=
"step_55-Problemsetup"></a>
429 *
virtual void vector_value(
const Point<dim> &p,
438 *
const double R_x = p[0];
439 *
const double R_y = p[1];
454 *
Utilities::fixed_power<2>(-
std::sqrt(25.0 + 4 *
pi2) + 5.0) *
461 *
Utilities::fixed_power<3>(-
std::sqrt(25.0 + 4 *
pi2) + 5.0) *
482 *
virtual void vector_value(
const Point<dim> &p,
490 *
const double R_x = p[0];
491 *
const double R_y = p[1];
516 *
(-6538034.74494422 +
535 * <a name=
"step_55-Themainprogram"></a>
560 *
void refine_grid();
564 *
const double viscosity;
576 *
LA::MPI::BlockSparseMatrix system_matrix;
623 * <a name=
"step_55-SystemSetup"></a>
638 *
dof_handler.distribute_dofs(fe);
657 *
pcout <<
" Number of degrees of freedom: " << dof_handler.n_dofs() <<
" ("
658 *
<<
n_u <<
'+' <<
n_p <<
')' << std::endl;
667 *
const IndexSet &locally_owned_dofs = dof_handler.locally_owned_dofs();
669 *
locally_owned_dofs.get_view(
n_u,
n_u +
n_p)};
680 * because we only perform global refinement, it is still a good idea
681 * to put this function call in, in case adaptive refinement gets
686 * constraints.reinit(locally_owned_dofs, locally_relevant_dofs);
688 * const FEValuesExtractors::Vector velocities(0);
689 * DoFTools::make_hanging_node_constraints(dof_handler, constraints);
690 * VectorTools::interpolate_boundary_values(dof_handler,
692 * ExactSolution<dim>(),
694 * fe.component_mask(velocities));
695 * constraints.close();
700 * Now we create the system matrix based on a BlockDynamicSparsityPattern.
709 *
system_matrix.
clear();
712 *
for (
unsigned int c = 0; c < dim + 1; ++c)
713 *
for (
unsigned int d = 0;
d < dim + 1; ++
d)
714 *
if (c == dim && d == dim)
716 *
else if (c == dim || d == dim || c == d)
728 *
dof_handler.locally_owned_dofs(),
746 *
for (
unsigned int c = 0; c < dim + 1; ++c)
747 *
for (
unsigned int d = 0;
d < dim + 1; ++
d)
748 *
if (c == dim && d == dim)
760 *
dof_handler.locally_owned_dofs()),
769 * function call
with two std::vector<IndexSet>
will create a ghosted
784 * <a name=
"step_55-Assembly"></a>
809 *
const unsigned int dofs_per_cell = fe.n_dofs_per_cell();
819 *
std::vector<Tensor<2, dim>>
grad_phi_u(dofs_per_cell);
820 *
std::vector<double>
div_phi_u(dofs_per_cell);
821 *
std::vector<double>
phi_p(dofs_per_cell);
823 *
std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
827 *
for (
const auto &cell : dof_handler.active_cell_iterators())
828 *
if (cell->is_locally_owned())
834 *
fe_values.reinit(cell);
837 *
for (
unsigned int q = 0;
q < n_q_points; ++
q)
839 *
for (
unsigned int k = 0;
k < dofs_per_cell; ++
k)
846 *
for (
unsigned int i = 0; i < dofs_per_cell; ++i)
848 *
for (
unsigned int j = 0;
j < dofs_per_cell; ++
j)
862 *
fe.system_to_component_index(i).first;
863 *
cell_rhs(i) += fe_values.shape_value(i,
q) *
869 *
cell->get_dof_indices(local_dof_indices);
891 * <a name=
"step_55-Solving"></a>
907 *
LA::MPI::PreconditionAMG
prec_A;
909 *
LA::MPI::PreconditionAMG::AdditionalData
data;
912 *
data.symmetric_operator =
true;
914 *
prec_A.initialize(system_matrix.block(0, 0),
data);
917 *
LA::MPI::PreconditionAMG
prec_S;
919 *
LA::MPI::PreconditionAMG::AdditionalData
data;
922 *
data.symmetric_operator =
true;
933 *
LA::MPI::PreconditionAMG>;
961 *
solver.solve(system_matrix,
966 *
pcout <<
" Solved in " << solver_control.last_step() <<
" iterations."
979 *
const double mean_pressure =
993 * <a name=
"step_55-Therest"></a>
1002 *
template <
int dim>
1012 *
template <
int dim>
1058 *
std::vector<DataComponentInterpretation::DataComponentInterpretation>
1059 *
data_component_interpretation(
1061 *
data_component_interpretation.push_back(
1065 *
data_out.attach_dof_handler(dof_handler);
1069 *
data_component_interpretation);
1085 *
data_component_interpretation);
1090 *
for (
unsigned int i = 0; i <
subdomain.size(); ++i)
1092 *
data_out.add_data_vector(
subdomain,
"subdomain");
1094 *
data_out.build_patches();
1096 *
data_out.write_vtu_with_pvtu_record(
1097 *
"./",
"solution", cycle, mpi_communicator, 2);
1102 *
template <
int dim>
1106 *
pcout <<
"Running using PETSc." << std::endl;
1108 *
pcout <<
"Running using Trilinos." << std::endl;
1110 *
const unsigned int n_cycles = 5;
1111 *
for (
unsigned int cycle = 0; cycle < n_cycles; ++cycle)
1113 *
pcout <<
"Cycle " << cycle <<
':' << std::endl;
1131 *
pcout << std::endl;
1142 *
using namespace dealii;
1143 *
using namespace Step55;
1150 *
catch (std::exception &exc)
1152 *
std::cerr << std::endl
1154 *
<<
"----------------------------------------------------"
1156 *
std::cerr <<
"Exception on processing: " << std::endl
1157 *
<< exc.what() << std::endl
1158 *
<<
"Aborting!" << std::endl
1159 *
<<
"----------------------------------------------------"
1166 *
std::cerr << std::endl
1168 *
<<
"----------------------------------------------------"
1170 *
std::cerr <<
"Unknown exception!" << std::endl
1171 *
<<
"Aborting!" << std::endl
1172 *
<<
"----------------------------------------------------"
1400<a name=
"step-55-extensions"></a>
1401<a name=
"step_55-Possibilitiesforextensions"></a><
h3>Possibilities
for extensions</
h3>
1431 QGauss<dim - 1>(fe.degree + 1),
std::conditional_t< rank_==1, std::array< Number, dim >, std::array< Tensor< rank_ - 1, dim, Number >, dim > > values
#define Assert(cond, exc)
void loop(IteratorType begin, std_cxx20::type_identity_t< IteratorType > end, DOFINFO &dinfo, INFOBOX &info, const std::function< void(std_cxx20::type_identity_t< DOFINFO > &, typename INFOBOX::CellInfo &)> &cell_worker, const std::function< void(std_cxx20::type_identity_t< DOFINFO > &, typename INFOBOX::CellInfo &)> &boundary_worker, const std::function< void(std_cxx20::type_identity_t< DOFINFO > &, std_cxx20::type_identity_t< DOFINFO > &, typename INFOBOX::CellInfo &, typename INFOBOX::CellInfo &)> &face_worker, AssemblerType &assembler, const LoopControl &lctrl=LoopControl())
void make_sparsity_pattern(const DoFHandler< dim, spacedim > &dof_handler, SparsityPatternBase &sparsity_pattern, const AffineConstraints< number > &constraints={}, const bool keep_constrained_dofs=true, const types::subdomain_id subdomain_id=numbers::invalid_subdomain_id)
@ update_values
Shape function values.
@ update_JxW_values
Transformed quadrature weights.
@ update_gradients
Shape function gradients.
@ update_quadrature_points
Transformed quadrature points.
std::vector< index_type > data
std::vector< value_type > split(const typename ::Triangulation< dim, spacedim >::cell_iterator &parent, const value_type parent_value)
@ component_is_part_of_vector
void component_wise(DoFHandler< dim, spacedim > &dof_handler, const std::vector< unsigned int > &target_component=std::vector< unsigned int >())
void hyper_cube(Triangulation< dim, spacedim > &tria, const double left=0., const double right=1., const bool colorize=false)
@ matrix
Contents is actually a matrix.
@ symmetric
Matrix is symmetric.
@ diagonal
Matrix is diagonal.
constexpr types::blas_int zero
constexpr types::blas_int one
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
std::vector< unsigned int > serial(const std::vector< unsigned int > &targets, const std::function< RequestType(const unsigned int)> &create_request, const std::function< AnswerType(const unsigned int, const RequestType &)> &answer_request, const std::function< void(const unsigned int, const AnswerType &)> &process_answer, const MPI_Comm comm)
unsigned int n_mpi_processes(const MPI_Comm mpi_communicator)
std::vector< T > all_gather(const MPI_Comm comm, const T &object_to_send)
unsigned int this_mpi_process(const MPI_Comm mpi_communicator)
void run(const Iterator &begin, const std_cxx20::type_identity_t< Iterator > &end, Worker worker, Copier copier, const ScratchData &sample_scratch_data, const CopyData &sample_copy_data, const unsigned int queue_length, const unsigned int chunk_size)
int(&) functions(const void *v1, const void *v2)
::VectorizedArray< Number, width > exp(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > cos(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > sin(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > sqrt(const ::VectorizedArray< Number, width > &)
const ::parallel::distributed::Triangulation< dim, spacedim > * triangulation