From e2e40aee142ad6edb523e5a30143431efb3eafb0 Mon Sep 17 00:00:00 2001
From: Michael Innerberger <michael.innerberger@asc.tuwien.ac.at>
Date: Wed, 19 Jul 2023 04:40:48 -0400
Subject: [PATCH] Refactor MultilevelSmoother to facilitate higher-order
variants
---
lib/solvers/OptimalVcycleMultigridSolver.m | 29 ++++++++++++++--------
lib/solvers/smoothers/MultilevelSmoother.m | 25 +++++++++++++++++--
lib/solvers/smoothers/P1JacobiSmoother.m | 16 +++++-------
3 files changed, 47 insertions(+), 23 deletions(-)
diff --git a/lib/solvers/OptimalVcycleMultigridSolver.m b/lib/solvers/OptimalVcycleMultigridSolver.m
index 9c3cc0b..496c2a7 100644
--- a/lib/solvers/OptimalVcycleMultigridSolver.m
+++ b/lib/solvers/OptimalVcycleMultigridSolver.m
@@ -1,6 +1,7 @@
% OptimalVcycleMultigridSolver (subclass of IterativeSolver) Solves linear
-% equations iteratively, using a geometric multigrid solver with V-cycle and
-% specified MultiLevelSmoother.
+% equations iteratively, using a geometric multigrid solver with specified
+% MultiLevelSmoother.
+% One iteration = Vcycle(0,1) (no pre-/one post-smoothing step).
%
% To obtain the solver described in [Innerberger, Miraçi, Praetorius,
% Streitberger; Optimal computational costs of AFEM with optimal local
@@ -28,7 +29,7 @@ classdef OptimalVcycleMultigridSolver < IterativeSolver
properties (Access=protected)
smoother
- matrix
+ projectedMatrix
lhsNorm
end
@@ -46,8 +47,8 @@ classdef OptimalVcycleMultigridSolver < IterativeSolver
% override superclass method, because ALL matrices are needed
setupSystemMatrix@IterativeSolver(obj, A, varargin{:});
- obj.smoother.setup(A, varargin{:});
- obj.matrix{obj.nLevels} = A;
+ projectedA = obj.smoother.setup(A, varargin{:});
+ obj.projectedMatrix{obj.nLevels} = projectedA;
end
function setupRhs(obj, b, varargin)
@@ -101,20 +102,26 @@ classdef OptimalVcycleMultigridSolver < IterativeSolver
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);
+ Cx = obj.projectedMatrix{1} \ res{1};
+ eta2 = dot(Cx, obj.projectedMatrix{1} * Cx, Dim.Vector);
% ascending cascade: (local) smoothing and compute error estimator
- for k = 2:L
+ for k = 2:(L-1)
Cx = obj.smoother.prolongate(Cx, k);
- updatedResidual = res{k} - obj.matrix{k} * Cx;
+ updatedResidual = res{k} - obj.projectedMatrix{k} * Cx;
rho = obj.smoother.smooth(updatedResidual, k);
- [etaUpdate2, xUpdate] = computeOptimalUpdate(obj.matrix{k}, updatedResidual, rho);
+ [etaUpdate2, xUpdate] = computeOptimalUpdate(obj.projectedMatrix{k}, updatedResidual, rho);
Cx = Cx + xUpdate;
eta2 = eta2 + etaUpdate2;
end
- algEstimator = sqrt(eta2);
+ % final smoothing and residual update (with non-projected matrix)
+ Cx = obj.smoother.prolongate(Cx, L);
+ updatedResidual = res{L} - obj.A * Cx;
+ rho = obj.smoother.smooth(updatedResidual, L);
+ [etaUpdate2, xUpdate] = computeOptimalUpdate(obj.A, updatedResidual, rho);
+ Cx = Cx + xUpdate;
+ algEstimator = sqrt(eta2 + etaUpdate2);
end
end
end
diff --git a/lib/solvers/smoothers/MultilevelSmoother.m b/lib/solvers/smoothers/MultilevelSmoother.m
index 45ff623..0616a03 100644
--- a/lib/solvers/smoothers/MultilevelSmoother.m
+++ b/lib/solvers/smoothers/MultilevelSmoother.m
@@ -5,6 +5,9 @@
%
% smoother.setup(A) hook to set up the smoother for the matrix A.
%
+% nonInvertedSmootherMatrix = smoother.setup(A) set up the smoother and also
+% return the non-inverted matrix, based on which the smoother is computed.
+%
% 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.
@@ -20,14 +23,32 @@ classdef MultilevelSmoother < handle
nLevels
end
+ properties (Access=protected)
+ fes
+ blf
+ end
+
%% methods
methods (Access=public)
- function obj = MultilevelSmoother()
+ function obj = MultilevelSmoother(fes, blf)
+ arguments
+ fes FeSpace
+ blf BilinearForm
+ end
+
+ assert(isempty(blf.b), ...
+ 'Multilevel smoothers only tested for symmetric problems.')
+
obj.nLevels = 0;
+ obj.fes = fes;
+ obj.blf = blf;
end
- function setup(obj, ~)
+ function nonInvertedSmootherMatrix = setup(obj, ~)
obj.nLevels = obj.nLevels + 1;
+
+ % should be implemented by subclasses
+ nonInvertedSmootherMatrix = [];
end
end
diff --git a/lib/solvers/smoothers/P1JacobiSmoother.m b/lib/solvers/smoothers/P1JacobiSmoother.m
index 59d4311..91a82f9 100644
--- a/lib/solvers/smoothers/P1JacobiSmoother.m
+++ b/lib/solvers/smoothers/P1JacobiSmoother.m
@@ -7,10 +7,8 @@
classdef P1JacobiSmoother < MultilevelSmoother
properties (Access=protected)
- fes
- blf
- inverseDiagonal
P
+ inverseDiagonal
intergridMatrix
freeVertices
freeVerticesOld
@@ -26,30 +24,28 @@ classdef P1JacobiSmoother < MultilevelSmoother
methods (Access=public)
function obj = P1JacobiSmoother(fes, blf, P)
arguments
- fes FeSpace
- blf BilinearForm
+ fes
+ blf
P Prolongation
end
- obj = obj@MultilevelSmoother();
+ obj = obj@MultilevelSmoother(fes, blf);
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)
+ function nonInvertedSmootherMatrix = setup(obj, A)
setup@MultilevelSmoother(obj, A);
L = obj.nLevels;
obj.freeVerticesOld = obj.freeVertices;
obj.freeVertices = getFreeDofs(obj.fes);
+ nonInvertedSmootherMatrix = A;
obj.inverseDiagonal{L} = full(diag(A)).^(-1);
if L >= 2
--
GitLab