497 * <a name=
"step_72-ThecodeMinimalSurfaceProblemcodeclassimplementation"></a>
503 * <a name=
"step_72-MinimalSurfaceProblemMinimalSurfaceProblem"></a>
522 * <a name=
"step_72-MinimalSurfaceProblemsetup_system"></a>
535 *
dof_handler.distribute_dofs(fe);
558 *
sparsity_pattern.copy_from(
dsp);
559 *
system_matrix.reinit(sparsity_pattern);
565 * <a name=
"step_72-Assemblingthelinearsystem"></a>
571 * <a name=
"step_72-Manualassembly"></a>
604 *
const unsigned int dofs_per_cell = fe.n_dofs_per_cell();
633 *
using CellIteratorType =
decltype(dof_handler.begin_active());
643 *
const ScratchData sample_scratch_data(fe,
648 *
const CopyData sample_copy_data(dofs_per_cell);
667 * `scratch_data`
object.
679 *
ScratchData &scratch_data,
680 *
CopyData ©_data) {
681 *
const auto &fe_values = scratch_data.reinit(cell);
685 *
std::vector<types::global_dof_index> &local_dof_indices =
686 *
copy_data.local_dof_indices[0];
687 *
cell->get_dof_indices(local_dof_indices);
691 *
For Newton
's method, we require the gradient of the solution at the
692 * point about which the problem is being linearized.
696 * Once we have that, we can perform assembly for this cell in
697 * the usual way. One minor difference to @ref step_15 "step-15" is that we've
703 *
fe_values.n_quadrature_points);
707 *
for (
const unsigned int q : fe_values.quadrature_point_indices())
713 *
for (
const unsigned int i : fe_values.dof_indices())
717 *
(((fe_values.shape_grad(i,
q)
719 *
* fe_values.shape_grad(
j,
q))
721 *
(fe_values.shape_grad(i,
q)
723 *
* (fe_values.shape_grad(
j,
q)
726 *
* fe_values.JxW(
q));
728 *
cell_rhs(i) -= (fe_values.shape_grad(i,
q)
731 *
* fe_values.JxW(
q));
751 *
const std::vector<types::global_dof_index> &local_dof_indices =
752 *
copy_data.local_dof_indices[0];
770 *
MeshWorker::mesh_loop(dof_handler.active_cell_iterators(),
773 *
sample_scratch_data,
814 *
const unsigned int dofs_per_cell = fe.n_dofs_per_cell();
818 *
using CellIteratorType =
decltype(dof_handler.begin_active());
820 *
const ScratchData sample_scratch_data(fe,
825 *
const CopyData sample_copy_data(dofs_per_cell);
829 *
We'll define up front the AD data structures that we'll be using,
835 *
know that we'll only be linearizing the residual, which means that we
836 * only need to compute first-order derivatives. The return values from the
837 * calculations are to be of type `double`.
841 * We also need an extractor to retrieve some data related to the field
842 * solution to the problem.
845 * using ADHelper = Differentiation::AD::ResidualLinearization<
846 * Differentiation::AD::NumberTypes::sacado_dfad,
848 * using ADNumberType = typename ADHelper::ad_type;
850 * const FEValuesExtractors::Scalar u_fe(0);
854 * With this, let us define the lambda function that will be used
855 * to compute the cell contributions to the Jacobian matrix and
856 * the right hand side:
859 * const auto cell_worker = [&u_fe, this](const CellIteratorType &cell,
860 * ScratchData &scratch_data,
861 * CopyData ©_data) {
862 * const auto &fe_values = scratch_data.reinit(cell);
863 * const unsigned int dofs_per_cell = fe_values.get_fe().n_dofs_per_cell();
865 * FullMatrix<double> &cell_matrix = copy_data.matrices[0];
866 * Vector<double> &cell_rhs = copy_data.vectors[0];
867 * std::vector<types::global_dof_index> &local_dof_indices =
868 * copy_data.local_dof_indices[0];
869 * cell->get_dof_indices(local_dof_indices);
887 * local solution coefficients
matches the number
of local residual
891 *
const unsigned
int n_independent_variables = local_dof_indices.
size();
892 *
const unsigned int n_dependent_variables = dofs_per_cell;
934 *
fe_values.n_quadrature_points);
935 *
fe_values[
u_fe].get_function_gradients_from_local_dof_values(
950 * value explicitly. After that, apart from a sign change the residual
951 * assembly looks much the same as we saw for the cell RHS vector before:
952 * We loop over all quadrature points, ensure that the coefficient now
953 * encodes its dependence on the (sensitive) finite element DoF values by
954 * using the correct `ADNumberType`, and finally we assemble the
955 * components of the residual vector. For complete clarity, the finite
956 * element shape functions (and their gradients, etc.) as well as the
957 * "JxW" values remain scalar
958 * valued, but the @p coeff and the @p old_solution_gradients at each
959 * quadrature point are computed in terms of the independent
963 * std::vector<ADNumberType> residual_ad(n_dependent_variables,
964 * ADNumberType(0.0));
965 * for (const unsigned int q : fe_values.quadrature_point_indices())
967 * const ADNumberType coeff =
968 * 1.0 / std::sqrt(1.0 + old_solution_gradients[q] *
969 * old_solution_gradients[q]);
971 * for (const unsigned int i : fe_values.dof_indices())
973 * residual_ad[i] += (fe_values.shape_grad(i, q) // \nabla \phi_i
975 * * old_solution_gradients[q]) // * \nabla u_n
976 * * fe_values.JxW(q); // * dx
982 * Once we have the full cell residual vector computed, we can register
983 * it with the helper class.
987 * Thereafter, we compute the residual values (basically,
988 * extracting the real values from what we already computed) and
989 * their Jacobian (the linearization of each residual component
990 * with respect to all cell DoFs) at the evaluation point. For
991 * the purposes of assembly into the global linear system, we
992 * have to respect the sign difference between the residual and
993 * the RHS contribution: For Newton's method,
the right
hand
1003 *
ad_helper.compute_linearization(cell_matrix);
1011 *
const auto copier = [
this](
const CopyData ©_data) {
1014 *
const std::vector<types::global_dof_index> &local_dof_indices =
1015 *
copy_data.local_dof_indices[0];
1024 *
sample_scratch_data,
1032 * <a name=
"step_72-Assemblyviadifferentiationoftheenergyfunctional"></a>
1056 *
template <
int dim>
1059 *
system_matrix = 0;
1062 *
const unsigned int dofs_per_cell = fe.n_dofs_per_cell();
1066 *
using CellIteratorType =
decltype(dof_handler.begin_active());
1068 *
const ScratchData sample_scratch_data(fe,
1073 *
const CopyData sample_copy_data(dofs_per_cell);
1084 *
its gradient).
We'll then need to linearize the residual, implying that
1085 * second derivatives of the potential energy must be computed. You might
1086 * want to compare this with the definition of `ADHelper` used int
1087 * previous function, where we used
1088 * `Differentiation::AD::ResidualLinearization<Differentiation::AD::NumberTypes::sacado_dfad,double>`.
1091 * using ADHelper = Differentiation::AD::EnergyFunctional<
1092 * Differentiation::AD::NumberTypes::sacado_dfad_dfad,
1094 * using ADNumberType = typename ADHelper::ad_type;
1096 * const FEValuesExtractors::Scalar u_fe(0);
1100 * Let us then again define the lambda function that does the integration on
1105 * To initialize an instance of the helper class, we now only require
1106 * that the number of independent variables (that is, the number
1107 * of degrees of freedom associated with the element solution vector)
1108 * are known up front. This is because the second-derivative matrix that
1109 * results from an energy functional is necessarily square (and also,
1110 * incidentally, symmetric).
1113 * const auto cell_worker = [&u_fe, this](const CellIteratorType &cell,
1114 * ScratchData &scratch_data,
1115 * CopyData ©_data) {
1116 * const auto &fe_values = scratch_data.reinit(cell);
1118 * FullMatrix<double> &cell_matrix = copy_data.matrices[0];
1119 * Vector<double> &cell_rhs = copy_data.vectors[0];
1120 * std::vector<types::global_dof_index> &local_dof_indices =
1121 * copy_data.local_dof_indices[0];
1122 * cell->get_dof_indices(local_dof_indices);
1124 * const unsigned int n_independent_variables = local_dof_indices.size();
1125 * ADHelper ad_helper(n_independent_variables);
1129 * Once more, we register all cell DoFs values with the helper, followed
1130 * by extracting the "sensitive" variant of these values that are to be
1131 * used in subsequent operations that must be differentiated -- one of
1132 * those being the calculation of the solution gradients.
1135 * ad_helper.register_dof_values(current_solution, local_dof_indices);
1137 * const std::vector<ADNumberType> &dof_values_ad =
1138 * ad_helper.get_sensitive_dof_values();
1140 * std::vector<Tensor<1, dim, ADNumberType>> old_solution_gradients(
1141 * fe_values.n_quadrature_points);
1142 * fe_values[u_fe].get_function_gradients_from_local_dof_values(
1143 * dof_values_ad, old_solution_gradients);
1147 * We next create a variable that stores the cell total energy.
1148 * Once more we emphasize that we explicitly zero-initialize this value,
1149 * thereby ensuring the integrity of the data for this starting value.
1153 * The aim for our approach is then to compute the cell total
1154 * energy, which is the sum of the internal (due to right hand
1155 * side functions, typically linear in @f$U@f$) and external
1156 * energies. In this particular case, we have no external
1157 * energies (e.g., from source terms or Neumann boundary
1167 *
for (
const unsigned int q : fe_values.quadrature_point_indices())
1177 *
After we've computed the total energy on this cell, we'll
1190 *
ad_helper.compute_linearization(cell_matrix);
1199 *
const auto copier = [
this](
const CopyData ©_data) {
1202 *
const std::vector<types::global_dof_index> &local_dof_indices =
1203 *
copy_data.local_dof_indices[0];
1212 *
sample_scratch_data,
1221 * <a name=
"step_72-MinimalSurfaceProblemsolve"></a>
1222 * <
h4>MinimalSurfaceProblem::solve</
h4>
1229 *
template <
int dim>
1237 *
preconditioner.initialize(system_matrix, 1.2);
1251 * <a name=
"step_72-MinimalSurfaceProblemrefine_mesh"></a>
1260 *
template <
int dim>
1296 * <a name=
"step_72-MinimalSurfaceProblemcompute_residual"></a>
1297 * <
h4>MinimalSurfaceProblem::compute_residual</
h4>
1305 * @
ref step_15
"step-15".
1308 *
template <
int dim>
1322 *
const unsigned int dofs_per_cell = fe.n_dofs_per_cell();
1326 *
std::vector<Tensor<1, dim>>
gradients(n_q_points);
1328 *
std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
1330 *
for (
const auto &cell : dof_handler.active_cell_iterators())
1333 *
fe_values.reinit(cell);
1337 *
for (
unsigned int q = 0;
q < n_q_points; ++
q)
1339 *
const double coeff =
1340 *
1.0 /
std::sqrt(1.0 + gradients[
q] * gradients[
q]);
1342 *
for (
unsigned int i = 0; i < dofs_per_cell; ++i)
1346 *
* fe_values.JxW(
q));
1349 *
cell->get_dof_indices(local_dof_indices);
1351 *
local_dof_indices,
1355 *
return residual.l2_norm();
1363 * <a name=
"step_72-MinimalSurfaceProblemdetermine_step_length"></a>
1370 * @
ref step_15
"step-15".
1373 *
template <
int dim>
1384 * <a name=
"step_72-MinimalSurfaceProblemoutput_results"></a>
1400 *
data_out.attach_dof_handler(dof_handler);
1403 *
data_out.build_patches();
1408 *
data_out.write_vtu(output);
1415 * <a name=
"step_72-MinimalSurfaceProblemrun"></a>
1416 * <
h4>MinimalSurfaceProblem::run</
h4>
1426 * the time taken to assemble for each of the three methods, we've also
1433 *
template <
int dim>
1435 *
const double tolerance)
1437 *
std::cout <<
"******** Assembly approach ********" << std::endl;
1439 *
{
"Unassisted implementation (full hand linearization).",
1440 *
"Automated linearization of the finite element residual.",
1441 *
"Automated computation of finite element residual and linearization using a variational formulation."}};
1463 *
std::cout <<
" Initial residual: " << compute_residual(0) << std::endl;
1489 *
std::cout <<
" Residual: " << compute_residual(0) << std::endl;
1495 *
std::cout << std::endl;
1504 * <a name=
"step_72-Themainfunction"></a>
1523 *
using namespace Step72;
1538 *
parameters.tolerance);
1540 *
catch (std::exception &exc)
1542 *
std::cerr << std::endl
1544 *
<<
"----------------------------------------------------"
1546 *
std::cerr <<
"Exception on processing: " << std::endl
1547 *
<< exc.what() << std::endl
1548 *
<<
"Aborting!" << std::endl
1549 *
<<
"----------------------------------------------------"
1556 *
std::cerr << std::endl
1558 *
<<
"----------------------------------------------------"
1560 *
std::cerr <<
"Unknown exception!" << std::endl
1561 *
<<
"Aborting!" << std::endl
1562 *
<<
"----------------------------------------------------"
1604relative time between the different implementations):
1606******** Assembly approach ********
1607Unassisted implementation (full hand linearization).
1611+---------------------------------------------+------------+------------+
1612| Total wallclock time elapsed since start | 35.1s | |
1614| Section | no. calls | wall time | % of total |
1615+---------------------------------+-----------+------------+------------+
1616| Assemble | 50 | 1.56s | 4.5% |
1617| Solve | 50 | 30.8s | 88% |
1618+---------------------------------+-----------+------------+------------+
1620And for the implementation that linearizes the residual in an automated
1621manner using the Sacado dynamic forward AD number type:
1623******** Assembly approach ********
1624Automated linearization of the finite element residual.
1628+---------------------------------------------+------------+------------+
1629| Total wallclock time elapsed since start | 40.1s | |
1631| Section | no. calls | wall time | % of total |
1632+---------------------------------+-----------+------------+------------+
1633| Assemble | 50 | 8.8s | 22% |
1634| Solve | 50 | 28.6s | 71% |
1635+---------------------------------+-----------+------------+------------+
1637And, lastly, for the implementation that computes both the residual and
1638its linearization directly from an energy functional (using nested Sacado
1639dynamic forward AD numbers):
1641******** Assembly approach ********
1642Automated computation of finite element residual and linearization using a variational formulation.
1646+---------------------------------------------+------------+------------+
1647| Total wallclock time elapsed since start | 48.8s | |
1649| Section | no. calls | wall time | % of total |
1650+---------------------------------+-----------+------------+------------+
1651| Assemble | 50 | 16.7s | 34% |
1652| Solve | 50 | 29.3s | 60% |
1653+---------------------------------+-----------+------------+------------+
1681<a name=
"step_72-Possibilitiesforextensions"></a><
h3> Possibilities
for extensions </h3>
1711 each cell quadrature point).
1723<a name=
"step_72-PlainProg"></a>
const unsigned int dofs_per_cell
static void estimate(const Mapping< dim, spacedim > &mapping, const DoFHandler< dim, spacedim > &dof, const Quadrature< dim - 1 > &quadrature, const std::map< types::boundary_id, const Function< spacedim, Number > * > &neumann_bc, const ReadVector< Number > &solution, Vector< float > &error, const ComponentMask &component_mask={}, const Function< spacedim > *coefficients=nullptr, const unsigned int n_threads=numbers::invalid_unsigned_int, const types::subdomain_id subdomain_id=numbers::invalid_subdomain_id, const types::material_id material_id=numbers::invalid_material_id, const Strategy strategy=cell_diameter_over_24)
static void initialize(const std::string &filename="", const std::string &output_filename="", const ParameterHandler::OutputStyle output_style_for_output_filename=ParameterHandler::Short, ParameterHandler &prm=ParameterAcceptor::prm, const ParameterHandler::OutputStyle output_style_for_filename=ParameterHandler::DefaultStyle)
Tensor< rank, dim, Number > sum(const Tensor< rank, dim, Number > &local, const MPI_Comm mpi_communicator)
std::conditional_t< rank_==1, std::array< Number, dim >, std::array< Tensor< rank_ - 1, dim, Number >, dim > > values
static ::ExceptionBase & ExcNotImplemented()
#define AssertIndexRange(index, range)
#define AssertThrow(cond, exc)
void mesh_loop(const CellIteratorType &begin, const CellIteratorType &end, const CellWorkerFunctionType &cell_worker, const CopierType &copier, const ScratchData &sample_scratch_data, const CopyData &sample_copy_data, const AssembleFlags flags=assemble_own_cells, const BoundaryWorkerFunctionType &boundary_worker=BoundaryWorkerFunctionType(), const FaceWorkerFunctionType &face_worker=FaceWorkerFunctionType(), const unsigned int queue_length=2 *MultithreadInfo::n_threads(), const unsigned int chunk_size=8)
void make_hanging_node_constraints(const DoFHandler< dim, spacedim > &dof_handler, AffineConstraints< number > &constraints)
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_JxW_values
Transformed quadrature weights.
@ update_gradients
Shape function gradients.
@ update_quadrature_points
Transformed quadrature points.
std::vector< index_type > data
CGAL::Exact_predicates_exact_constructions_kernel_with_sqrt K
void hyper_ball(Triangulation< dim, spacedim > &tria, const Point< spacedim > ¢er={}, const double radius=1., const bool attach_spherical_manifold_on_boundary_cells=false)
void refine_and_coarsen_fixed_number(Triangulation< dim, spacedim > &triangulation, const Vector< Number > &criteria, const double top_fraction_of_cells, const double bottom_fraction_of_cells, const unsigned int max_n_cells=std::numeric_limits< unsigned int >::max())
@ matrix
Contents is actually a matrix.
constexpr types::blas_int zero
constexpr types::blas_int one
void cell_matrix(FullMatrix< double > &M, const FEValuesBase< dim > &fe, const FEValuesBase< dim > &fetest, const ArrayView< const std::vector< double > > &velocity, const double factor=1.)
void cell_residual(Vector< double > &result, const FEValuesBase< dim > &fe, const std::vector< Tensor< 1, dim > > &input, const ArrayView< const std::vector< double > > &velocity, double factor=1.)
Point< spacedim > point(const gp_Pnt &p, const double tolerance=1e-10)
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
SymmetricTensor< 2, dim, Number > b(const Tensor< 2, dim, Number > &F)
Tensor< 2, dim, Number > F(const Tensor< 2, dim, Number > &Grad_u)
constexpr const ReferenceCell Line
constexpr ReturnType< rank, T >::value_type & extract(T &t, const ArrayType &indices)
std::string int_to_string(const unsigned int value, const unsigned int digits=numbers::invalid_unsigned_int)
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)
void copy(const T *begin, const T *end, U *dest)
int(&) functions(const void *v1, const void *v2)
::VectorizedArray< Number, width > sin(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > sqrt(const ::VectorizedArray< Number, width > &)
const ::parallel::distributed::Triangulation< dim, spacedim > * triangulation