|
SpECTRE
v2024.04.12
|
A conjugate gradient solver for linear systems of equations \(Ax=b\) where the operator \(A\) is symmetric. More...
#include <ConjugateGradient.hpp>
Public Types | |
| using | fields_tag = FieldsTag |
| using | options_group = OptionsGroup |
| using | source_tag = SourceTag |
| using | operand_tag = db::add_tag_prefix< LinearSolver::Tags::Operand, fields_tag > |
| Apply the linear operator to this tag in each iteration. | |
| using | component_list = tmpl::list< detail::ResidualMonitor< Metavariables, FieldsTag, OptionsGroup > > |
| The parallel components used by the conjugate gradient linear solver. | |
| using | initialize_element = detail::InitializeElement< FieldsTag, OptionsGroup > |
| using | register_element = tmpl::list<> |
| template<typename ApplyOperatorActions , typename Label = OptionsGroup> | |
| using | solve = tmpl::list< detail::PrepareSolve< FieldsTag, OptionsGroup, Label, SourceTag >, detail::InitializeHasConverged< FieldsTag, OptionsGroup, Label >, ApplyOperatorActions, detail::PerformStep< FieldsTag, OptionsGroup, Label >, detail::UpdateFieldValues< FieldsTag, OptionsGroup, Label >, detail::UpdateOperand< FieldsTag, OptionsGroup, Label > > |
A conjugate gradient solver for linear systems of equations \(Ax=b\) where the operator \(A\) is symmetric.
The only operation we need to supply to the algorithm is the result of the operation \(A(p)\) (see Linear Solver) that in the case of the conjugate gradient algorithm must be symmetric. Each step of the algorithm expects that \(A(p)\) is computed and stored in the DataBox as db::add_tag_prefix<LinearSolver::Tags::OperatorAppliedTo, db::add_tag_prefix<LinearSolver::Tags::Operand, FieldsTag>>. To perform a solve, add the solve action list to an array parallel component. Pass the actions that compute \(A(q)\), as well as any further actions you wish to run in each step of the algorithm, as the first template parameter to solve. If you add the solve action list multiple times, use the second template parameter to label each solve with a different type.
Note that the operand \(p\) for which \(A(p)\) needs to be computed is not the field \(x\) we are solving for but db::add_tag_prefix<LinearSolver::Tags::Operand, FieldsTag>. This field is initially set to the residual \(r = b - A(x_0)\) where \(x_0\) is the initial value of the FieldsTag.
When the algorithm step is performed after the operator action \(A(p)\) has been computed and stored in the DataBox, the conjugate gradient algorithm implemented here will converge the field \(x\) towards the solution and update the operand \(p\) in the process. This requires two reductions over all elements that are received by a ResidualMonitor singleton parallel component, processed, and then broadcast back to all elements. The actions are implemented in the cg::detail namespace and constitute the full algorithm in the following order:
PerformStep (on elements): Compute the inner product \(\langle p, A(p)\rangle\) and reduce.ComputeAlpha (on ResidualMonitor): Compute \(\alpha=\frac{r^2}{\langle p, A(p)\rangle}\) and broadcast.UpdateFieldValues (on elements): Update \(x\) and \(r\), then compute the inner product \(\langle r, r\rangle\) and reduce to find the new \(r^2\).UpdateResidual (on ResidualMonitor): Store the new \(r^2\) and broadcast the ratio of the new and old \(r^2\), as well as a termination flag if the Convergence::Tags::Criteria are met.UpdateOperand (on elements): Update \(p\).