diff --git a/lib/solvers/OptimalVcycleMultigridSolver.m b/lib/solvers/OptimalVcycleMultigridSolver.m
new file mode 100644
index 0000000000000000000000000000000000000000..9c3cc0baa13945cc83fef787961aa38e967c5358
--- /dev/null
+++ b/lib/solvers/OptimalVcycleMultigridSolver.m
@@ -0,0 +1,137 @@
+% OptimalVcycleMultigridSolver (subclass of IterativeSolver) Solves linear
+% equations iteratively, using a geometric multigrid solver with V-cycle and
+% specified MultiLevelSmoother.
+%
+% To obtain the solver described in [Innerberger, Miraçi, Praetorius,
+% Streitberger; Optimal computational costs of AFEM with optimal local
+% hp-robust multigrid solver], choose the smoother as DiagonalJacobi (p=1),
+% LowOrderBlockJacobi (p > 1, 1->1->p iteration), or HighOrderBlockJacobi
+% (p > 1, 1->p->p iteration).
+% This solver also comes with built-in estimator for the algebraic error,
+% which is also returned by this multigrid solver.
+%
+% See also IterativeSolver, MultiLevelSmoother, DiagonalJacobi,
+% LowOrderBlockJacobi, HighOrderBlockJacobi
+
+classdef OptimalVcycleMultigridSolver < IterativeSolver
+ %% properties
+ properties (SetAccess=protected, GetAccess=public)
+ residualCNorm
+ algEstimator
+ residual
+ Cresidual
+ end
+
+ properties (Dependent, Access=public)
+ nLevels
+ end
+
+ properties (Access=protected)
+ smoother
+ matrix
+ lhsNorm
+ end
+
+ methods
+ function obj = OptimalVcycleMultigridSolver(smoother)
+ obj = obj@IterativeSolver();
+ obj.smoother = smoother;
+ end
+
+ function n = get.nLevels(obj)
+ n = obj.smoother.nLevels;
+ end
+
+ function setupSystemMatrix(obj, A, varargin)
+ % override superclass method, because ALL matrices are needed
+ setupSystemMatrix@IterativeSolver(obj, A, varargin{:});
+
+ obj.smoother.setup(A, varargin{:});
+ obj.matrix{obj.nLevels} = A;
+ end
+
+ function setupRhs(obj, b, varargin)
+ setupRhs@IterativeSolver(obj, b, varargin{:});
+ % initialize residual & hierarchy
+ obj.residual = b - obj.A * obj.x;
+ [obj.Cresidual, obj.algEstimator] = obj.Vcycle(obj.residual);
+ obj.residualCNorm = sqrt(dot(obj.residual, obj.Cresidual, Dim.Vector));
+ obj.lhsNorm = vecnorm(b, 2, Dim.Vector);
+ end
+
+ % stopping criterion
+ function tf = isConverged(obj)
+ tf = ((obj.iterationCount >= obj.maxIter) ...
+ | (obj.residualCNorm < obj.tol) ...
+ | (obj.residualCNorm ./ obj.lhsNorm < obj.tol));
+ end
+ end
+
+ methods (Access=protected)
+ % one step
+ function computeUpdate(obj)
+ % only iterate on active components
+ idx = find(obj.activeComponents);
+
+ % update solution & residual
+ obj.x(:,idx) = obj.x(:,idx) + obj.Cresidual(:,idx);
+ obj.residual(:,idx) = obj.residual(:,idx) - obj.A * obj.Cresidual(:,idx);
+
+ % residual & update error correction
+ [obj.Cresidual(:,idx), obj.algEstimator(idx)] = obj.Vcycle(obj.residual(:,idx));
+ obj.residualCNorm(:,idx) = sqrt(dot(obj.residual(:,idx), obj.Cresidual(:,idx), Dim.Vector));
+ end
+
+ function [Cx, algEstimator] = Vcycle(obj, x)
+ assert(isequal(size(x, 1), size(obj.A, 1)), ...
+ 'Setup for multilevel iteration not complete!')
+
+ % if there is only one coarse level: exact solve
+ L = obj.nLevels;
+ if L == 1
+ Cx = obj.A \ x;
+ algEstimator = sqrt(dot(Cx, obj.A * Cx, Dim.Vector));
+ return
+ end
+
+ % descending cascade: no smoothing
+ res{L} = x;
+ for k = L:-1:2
+ res{k-1} = obj.smoother.restrict(res{k}, k);
+ end
+
+ % exact solve on coarsest level to compute accumulative lifting of x
+ Cx = obj.matrix{1} \ res{1};
+ eta2 = dot(Cx, obj.matrix{1} * Cx, Dim.Vector);
+
+ % ascending cascade: (local) smoothing and compute error estimator
+ for k = 2:L
+ Cx = obj.smoother.prolongate(Cx, k);
+ updatedResidual = res{k} - obj.matrix{k} * Cx;
+ rho = obj.smoother.smooth(updatedResidual, k);
+ [etaUpdate2, xUpdate] = computeOptimalUpdate(obj.matrix{k}, updatedResidual, rho);
+ Cx = Cx + xUpdate;
+ eta2 = eta2 + etaUpdate2;
+ end
+
+ algEstimator = sqrt(eta2);
+ end
+ end
+end
+
+%% auxiliary functions
+% error correction with optimal stepsize
+function [etaUpdate2, resUpdate] = computeOptimalUpdate(A, res, rho)
+ rhoArho = dot(rho, A*rho, Dim.Vector);
+ if max(abs(rho)) < eps
+ lambda = 1;
+ else
+ lambda = dot(res, rho, Dim.Vector) ./ rhoArho;
+ end
+ resUpdate = lambda .* rho;
+ etaUpdate2 = rhoArho .* lambda.^2;
+
+ if any(lambda > 3)
+ warning('MG step-sizes no longer bound by d+1. Optimality of step size cannot be guaranteed!')
+ end
+end
diff --git a/lib/solvers/chooseIterativeSolver.m b/lib/solvers/chooseIterativeSolver.m
index 02d735c3b8fb17f82f0f0c2bae9c2a5658f8383a..fa30db22d2f6701c3fc1c43c78c4b31d43c63937 100644
--- a/lib/solvers/chooseIterativeSolver.m
+++ b/lib/solvers/chooseIterativeSolver.m
@@ -53,13 +53,13 @@ switch class
switch variant
case {"", "lowOrderVcycle"}
if order == 1
- solver = LowestOrderLocalMg(fes, blf, P);
+ solver = OptimalVcycleMultigridSolver(P1JacobiSmoother(fes, blf, P));
else
solver = LocalMgLowOrderVcycle(fes, blf);
end
case "highOrderVcycle"
if order == 1
- solver = LowestOrderLocalMg(fes, blf, P);
+ solver = OptimalVcycleMultigridSolver(P1JacobiSmoother(fes, blf, P));
else
solver = LocalMgHighOrderVcycle(fes, blf, P);
end
diff --git a/lib/solvers/smoothers/MultilevelSmoother.m b/lib/solvers/smoothers/MultilevelSmoother.m
new file mode 100644
index 0000000000000000000000000000000000000000..45ff623060c116210a886c73c0587de8027eb30a
--- /dev/null
+++ b/lib/solvers/smoothers/MultilevelSmoother.m
@@ -0,0 +1,40 @@
+% MultilevelSmoother (abstract handle class) Interface for smoothers to use in a
+% multigrid solver.
+%
+% Abstract methods to implement:
+%
+% smoother.setup(A) hook to set up the smoother for the matrix A.
+%
+% Cx = smoother.smooth(x, k) applies the smoother to the vector x on level k.
+%
+% Px = smoother.prolongate(x, k) prolongate vector x from level k to level k+1.
+%
+% Px = smoother.restrict(x, k) restrict vector x from level k+1 to level k.
+%
+% See also: OptimalVcycleMultigridSolver
+
+
+classdef MultilevelSmoother < handle
+ %% properties
+ properties (GetAccess=public, SetAccess=protected)
+ nLevels
+ end
+
+ %% methods
+ methods (Access=public)
+ function obj = MultilevelSmoother()
+ obj.nLevels = 0;
+ end
+
+ function setup(obj, ~)
+ obj.nLevels = obj.nLevels + 1;
+ end
+ end
+
+ %% abstract methods
+ methods (Abstract, Access=public)
+ smooth(obj, x, k)
+ prolongate(obj, x, k)
+ restrict(obj, x, k)
+ end
+end
\ No newline at end of file
diff --git a/lib/solvers/smoothers/P1JacobiSmoother.m b/lib/solvers/smoothers/P1JacobiSmoother.m
new file mode 100644
index 0000000000000000000000000000000000000000..59d4311240877cbe1245f1d580ae5e4aa8e43a86
--- /dev/null
+++ b/lib/solvers/smoothers/P1JacobiSmoother.m
@@ -0,0 +1,91 @@
+% P1JacobiSmoother (subclass of MultilevelSmoother) multilevel Jacobi smoother
+% for lowest order finite elements: smooth with diagonal of stiffness matrix
+% on changed patches on every level.
+%
+% See also: MultilevelSmoother, OptimalVcycleMultigridSolver
+
+
+classdef P1JacobiSmoother < MultilevelSmoother
+ properties (Access=protected)
+ fes
+ blf
+ inverseDiagonal
+ P
+ intergridMatrix
+ freeVertices
+ freeVerticesOld
+ changedPatches
+ end
+
+ %% event data
+ properties (Access=private)
+ listenerHandle
+ end
+
+ %% methods
+ methods (Access=public)
+ function obj = P1JacobiSmoother(fes, blf, P)
+ arguments
+ fes FeSpace
+ blf BilinearForm
+ P Prolongation
+ end
+ obj = obj@MultilevelSmoother();
+
+ assert(fes.finiteElement.order == 1, ...
+ 'P1JacobiSmoother only works for lowest order finite elements.')
+
+ obj.fes = fes;
+ obj.P = P;
+ obj.blf = blf;
+
+ mesh = fes.mesh;
+ obj.listenerHandle = mesh.listener('IsAboutToRefine', @obj.getChangedPatches);
+ end
+
+ function setup(obj, A)
+ setup@MultilevelSmoother(obj, A);
+
+ L = obj.nLevels;
+ obj.freeVerticesOld = obj.freeVertices;
+ obj.freeVertices = getFreeDofs(obj.fes);
+
+ obj.inverseDiagonal{L} = full(diag(A)).^(-1);
+
+ if L >= 2
+ obj.intergridMatrix{L} = obj.P.matrix(obj.freeVertices, obj.freeVerticesOld);
+ obj.changedPatches{L} = find(obj.changedPatches{L}(obj.freeVertices));
+ obj.inverseDiagonal{L} = obj.inverseDiagonal{L}(obj.changedPatches{L},:);
+ end
+ end
+
+ function Cx = smooth(obj, x, k)
+ idx = obj.changedPatches{k};
+ Cx = zeros(size(x));
+ Cx(idx,:) = obj.inverseDiagonal{k} .* x(idx,:);
+ end
+
+ function Px = prolongate(obj, x, k)
+ Px = obj.intergridMatrix{k} * x;
+ end
+
+ function Px = restrict(obj, x, k)
+ Px = obj.intergridMatrix{k}' * x;
+ end
+ end
+
+ methods (Access=protected)
+ % callback function to be executed before mesh refinement:
+ % -> patches change for all vertices adjacent to refined edges and
+ % all new vertices
+ function getChangedPatches(obj, mesh, bisecData)
+ nCOld = mesh.nCoordinates;
+ nCNew = mesh.nCoordinates + nnz(bisecData.bisectedEdges) + ...
+ bisecData.nRefinedElements'*bisecData.nInnerNodes;
+ bisectedEdgeNodes = unique(mesh.edges(:,bisecData.bisectedEdges));
+ obj.changedPatches{obj.nLevels+1} = false(nCNew, 1);
+ idx = [unique(bisectedEdgeNodes); ((nCOld+1):nCNew)'];
+ obj.changedPatches{obj.nLevels+1}(idx) = true;
+ end
+ end
+end
\ No newline at end of file