diff --git a/lib/solvers/LocalMgHighOrderVcycle.m b/lib/solvers/LocalMgHighOrderVcycle.m
deleted file mode 100644
index 1ac07b7c4b4b43fd3bf136c2fd9d11c58ed03dd3..0000000000000000000000000000000000000000
--- a/lib/solvers/LocalMgHighOrderVcycle.m
+++ /dev/null
@@ -1,172 +0,0 @@
-% OptimalLocalMGSolver (subclass of MGSolver) Solves linear equations
-% iteratively using high-order p-optimal multilevel local smoothing (additive
-% Schwarz/block Jacobi).
-% One iteration = Vcycle(0,1) (no pre-/one post-smoothing step).
-%
-% With the solver comes an
-
-classdef LocalMgHighOrderVcycle < MGSolver
- %% properties
- properties (GetAccess=public,SetAccess=protected)
- nLevels
- end
-
- properties (Access=protected)
- blf
- hoFes
- loFes
- P
- intergridMatrix
- inclusionMatrix
- freeDofs
- freeDofsOld
- changedPatches
- patchwiseA
- Acoarse
- Afine
- coarseLoFes
- coarseHoFes
- end
-
- %% event data
- properties (Access=private)
- listenerHandle
- end
-
- %% methods
- methods (Access=public)
- function obj = LocalMgHighOrderVcycle(fes, blf, P)
- arguments
- fes FeSpace
- blf BilinearForm
- P Prolongation
- end
- obj = obj@MGSolver();
-
- assert(fes.finiteElement.order > 1, ...
- 'LocalMgLowOrderVcycle only works for higher order finite elements.')
- assert(isempty(blf.b), ...
- 'Multigrid solvers only tested for symmetric problems.')
-
- obj.nLevels = 0;
-
- mesh = fes.mesh;
- obj.P = P;
- obj.loFes = FeSpace(mesh, LowestOrderH1Fe(), 'dirichlet', fes.bnd.dirichlet);
- obj.hoFes = fes;
- obj.blf = blf;
- obj.Afine = cell(1,1);
- obj.Acoarse = assemble(obj.blf, obj.loFes);
- obj.Acoarse = obj.Acoarse(getFreeDofs(obj.loFes), getFreeDofs(obj.loFes));
-
- % set up FE spaces on coarsest level for projection
- coarseMesh = clone(fes.mesh);
- obj.coarseLoFes = FeSpace(coarseMesh, LowestOrderH1Fe, 'dirichlet', fes.bnd.dirichlet);
- obj.coarseHoFes = FeSpace(coarseMesh, HigherOrderH1Fe(fes.finiteElement.order), 'dirichlet', fes.bnd.dirichlet);
-
- obj.listenerHandle = mesh.listener('IsAboutToRefine', @obj.getChangedPatches);
- end
-
- function setupSystemMatrix(obj, A)
- obj.nLevels = obj.nLevels + 1;
-
- L = obj.nLevels;
- freeVertices = getFreeDofs(obj.loFes);
-
- obj.freeDofsOld = obj.freeDofs;
- obj.freeDofs = getFreeDofs(obj.hoFes);
-
- if L >= 2
- obj.intergridMatrix{L} = obj.P.matrix(obj.freeDofs, obj.freeDofsOld);
- obj.changedPatches{L} = freeVertices(obj.changedPatches{L}(freeVertices));
- obj.inclusionMatrix = SpaceProlongation(obj.coarseLoFes, obj.coarseHoFes) ...
- .matrix(getFreeDofs(obj.coarseLoFes), getFreeDofs(obj.coarseHoFes));
- if L == 2
- obj.patchwiseA{L} = assemblePatchwise(obj.blf, obj.hoFes, ':');
- else
- obj.patchwiseA{L} = assemblePatchwise(obj.blf, obj.hoFes, obj.changedPatches{L});
- end
- end
-
- setupSystemMatrix@MGSolver(obj, A);
- obj.Afine{L} = obj.A;
- end
-
- function setupRhs(obj, b, varargin)
- setupRhs@MGSolver(obj, b, varargin{:});
- end
-
- % Geometric MultiGrid
- function [Cx, algError2] = Vcycle(obj, res)
- assert(isequal(size(res, 1), size(obj.A, 1)), ...
- 'Setup for multilevel iteration not complete!')
-
- % if there is only one coarse level: exact solve
- if obj.nLevels == 1
- Cx = obj.A \ res;
- algError2 = scalarProduct(Cx, obj.A*Cx);
- return
- end
-
- L = obj.nLevels;
- algError2 = zeros(1, size(res, 2)); % built-in estimator of the algebraic error
-
- % descending cascade in p: no smoothing
- residual{L} = res;
- for k = L:-1:2
- residual{k-1} = obj.intergridMatrix{k}'*residual{k};
- end
-
- % coarse level in P1
- residual{1} = obj.inclusionMatrix * residual{1};
- sigma = obj.Acoarse \ residual{1};
- algError2 = algError2 + scalarProduct(sigma, obj.Acoarse*sigma);
-
- % ascending cascade in Pp
- sigma = obj.inclusionMatrix' * sigma;
-
- for k = 2:L
- sigma = obj.intergridMatrix{k}*sigma;
- uptres = residual{k} - obj.Afine{k}*sigma; %updated residual
- rho = obj.patchwiseA{k} \ uptres;
- [aeUpd, sUpd] = computeOptimalUpdate(obj.Afine{k}, uptres, rho);
- sigma = sigma + sUpd;
- algError2 = algError2 + aeUpd;
- end
-
- Cx = sigma;
- 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
-
-%% auxiliary functions
-% error correction with optimal stepsize
-function [etaUpdate, sigmaUpdate] = computeOptimalUpdate(A, res, rho)
- rhoArho = scalarProduct(rho, A*rho);
- lambda = scalarProduct(res, rho) ./ rhoArho;
- sigmaUpdate = lambda.*rho;
- etaUpdate = 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
-
-function a = scalarProduct(x, y)
- a = sum(x.*y, 1);
-end
\ No newline at end of file
diff --git a/lib/solvers/MGSolver.m b/lib/solvers/MGSolver.m
deleted file mode 100644
index 036213e537cb8e4ec83c1aaa521c959c2eb66f4a..0000000000000000000000000000000000000000
--- a/lib/solvers/MGSolver.m
+++ /dev/null
@@ -1,74 +0,0 @@
-% MGSolver (abstract subclass of IterativeSolver) Solves linear equations
-% iteratively, using a Vcycle geometric multigrid solver.
-%
-% Abstract methods to be implemented by subclasses:
-%
-% solver.Vcycle(x) performs one Vcycle iteration.
-%
-% See also IterativeSolver
-
-classdef MGSolver < IterativeSolver
- %% properties
- properties (SetAccess=protected, GetAccess=public)
- residualCNorm
- algEstimator
- end
-
- properties (Access=protected)
- residual
- Cresidual
- normb
- end
-
- %% methods
- methods (Access=public)
- % initialize iteration
- function setupSystemMatrix(obj, A, varargin)
- setupSystemMatrix@IterativeSolver(obj, 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 = sum(obj.residual.*obj.Cresidual, 1);
- obj.normb = sqrt(sum(b.^2, 1));
- end
-
- % stopping criterion
- function tf = isConverged(obj)
- tf = ((obj.iterationCount >= obj.maxIter) ...
- | (sqrt(obj.residualCNorm) < obj.tol) ...
- | (sqrt(obj.residualCNorm)./obj.normb < 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) = sum(obj.residual(:,idx).*obj.Cresidual(:,idx), 1);
- end
- end
-
- methods (Abstract,Access=public)
- Vcycle(obj, x)
- end
-end
-
-
-
-
-
-
-
-
diff --git a/lib/solvers/chooseIterativeSolver.m b/lib/solvers/chooseIterativeSolver.m
index 9bc2cf95231e39f832f7b21faad1791314702795..5058d490f878596961034050428649b48499a063 100644
--- a/lib/solvers/chooseIterativeSolver.m
+++ b/lib/solvers/chooseIterativeSolver.m
@@ -53,19 +53,20 @@ switch class
switch variant
case {"", "lowOrderVcycle"}
if order == 1
- solver = OptimalVcycleMultigridSolver(P1JacobiSmoother(fes, blf, P));
+ smoother = P1JacobiSmoother(fes, blf, P);
else
- solver = OptimalVcycleMultigridSolver(JacobiLowOrderCascadeSmoother(fes, blf));
+ smoother = JacobiLowOrderCascadeSmoother(fes, blf);
end
case "highOrderVcycle"
if order == 1
- solver = OptimalVcycleMultigridSolver(P1JacobiSmoother(fes, blf, P));
+ smoother = P1JacobiSmoother(fes, blf, P);
else
- solver = LocalMgHighOrderVcycle(fes, blf, P);
+ smoother = JacobiHighOrderCascadeSmoother(fes, blf, P);
end
otherwise
error('No multigrid variant %s!', variant)
end
+ solver = OptimalVcycleMultigridSolver(smoother);
% direct solve (for testing purposes)
case "direct"
diff --git a/lib/solvers/smoothers/JacobiHighOrderCascadeSmoother.m b/lib/solvers/smoothers/JacobiHighOrderCascadeSmoother.m
new file mode 100644
index 0000000000000000000000000000000000000000..793261b86d04ada5761991f725058557e73274cf
--- /dev/null
+++ b/lib/solvers/smoothers/JacobiHighOrderCascadeSmoother.m
@@ -0,0 +1,85 @@
+% JacobiHighOrderCascadeSmoother (subclass of MultilevelSmoother) multilevel
+% Jacobi smoother for higher order finite elements: smooth locally with
+% patchwise stiffness matrix on changed patches on every level, except for
+% coarsest level, where global p1-solving is done. Because of the lifting from
+% p1 to higher order between the coarsest and the second level, all patches
+% are assumed to be changed on the second level
+%
+% Reference: [Innerberger, Miraçi, Praetorius, Streitberger; Optimal
+% computational costs of AFEM with optimal local hp-robust multigrid solver].
+%
+% See also: MultilevelSmoother, OptimalVcycleMultigridSolver
+
+
+classdef JacobiHighOrderCascadeSmoother < MultilevelSmoother
+ properties (Access=protected)
+ P
+ loFes
+ patchwiseA
+ inclusionMatrix
+ intergridMatrix
+ freeDofs
+ freeDofsOld
+ end
+
+ %% methods
+ methods (Access=public)
+ function obj = JacobiHighOrderCascadeSmoother(fes, blf, P)
+ obj = obj@MultilevelSmoother(fes, blf);
+
+ assert(fes.finiteElement.order > 1, ...
+ 'JacobiHighOrderCascadeSmoother only works for higher order finite elements.')
+
+ obj.P = P;
+ obj.loFes = FeSpace(fes.mesh, LowestOrderH1Fe(), 'dirichlet', fes.bnd.dirichlet);
+ end
+
+ function nonInvertedSmootherMatrix = setup(obj, A)
+ setup@MultilevelSmoother(obj, A);
+
+ L = obj.nLevels;
+ obj.freeDofsOld = obj.freeDofs;
+ obj.freeDofs = getFreeDofs(obj.fes);
+ freeVertices = getFreeDofs(obj.loFes);
+
+ hoFes = obj.fes;
+ if L == 1
+ % set up FE inclusion and lowest order matrix on coarsest level
+ obj.inclusionMatrix = SpaceProlongation(obj.loFes, hoFes) ...
+ .matrix(getFreeDofs(obj.loFes), getFreeDofs(hoFes));
+ nonInvertedSmootherMatrix = assemble(obj.blf, obj.loFes);
+ nonInvertedSmootherMatrix = nonInvertedSmootherMatrix(freeVertices, freeVertices);
+ else
+ % note: patchwiseA is not already indexed by free vertices -> use freeVertices here
+ obj.changedPatches{L} = freeVertices(obj.changedPatches{L}(freeVertices));
+ obj.intergridMatrix{L} = obj.P.matrix(obj.freeDofs, obj.freeDofsOld);
+ nonInvertedSmootherMatrix = A;
+ if L == 2
+ obj.patchwiseA{L} = assemblePatchwise(obj.blf, hoFes, ':');
+ else
+ obj.patchwiseA{L} = assemblePatchwise(obj.blf, hoFes, obj.changedPatches{L});
+ end
+ end
+ end
+
+ function Cx = smooth(obj, x, k)
+ % higher order global (k=2) or local (else) patchwise smooting
+ Cx = obj.patchwiseA{k} \ x;
+ end
+
+ % order of FESpace per level: 1 -> p -> ... -> p
+ function Px = prolongate(obj, x, k)
+ if k == 2
+ x = obj.inclusionMatrix' * x;
+ end
+ Px = obj.intergridMatrix{k} * x;
+ end
+
+ function Px = restrict(obj, x, k)
+ Px = obj.intergridMatrix{k}' * x;
+ if k == 2
+ Px = obj.inclusionMatrix * Px;
+ end
+ end
+ end
+end
\ No newline at end of file
diff --git a/lib/solvers/smoothers/JacobiLowOrderCascadeSmoother.m b/lib/solvers/smoothers/JacobiLowOrderCascadeSmoother.m
index ad5dd1e038e06072b08b35a03d78f9514ef3af4b..10541f2f7db960a80826b9f75230c74d7172f3e8 100644
--- a/lib/solvers/smoothers/JacobiLowOrderCascadeSmoother.m
+++ b/lib/solvers/smoothers/JacobiLowOrderCascadeSmoother.m
@@ -48,9 +48,10 @@ classdef JacobiLowOrderCascadeSmoother < MultilevelSmoother
nonInvertedSmootherMatrix = p1Matrix;
if L >= 2
+ % note: p1 smoother is already indexed by free vertices -> use find here
+ obj.changedPatches{L} = find(obj.changedPatches{L}(obj.freeVertices));
% intermediate level smoothers and transfer operators
obj.intergridMatrix{L} = obj.P.matrix(obj.freeVertices, obj.freeVerticesOld);
- obj.changedPatches{L} = find(obj.changedPatches{L}(obj.freeVertices));
obj.p1Smoother{L} = obj.p1Smoother{L}(obj.changedPatches{L});
% finest level smoother and transfer operator
obj.hoSmoother = assemblePatchwise(obj.blf, obj.fes, ':');
@@ -71,6 +72,7 @@ classdef JacobiLowOrderCascadeSmoother < MultilevelSmoother
end
end
+ % order of FESpace per level: 1 -> ... -> 1 -> p
function Px = prolongate(obj, x, k)
Px = obj.intergridMatrix{k} * x;
if k == obj.nLevels