From a4f622140a92865f6a4e847d2b039eca57ae96e5 Mon Sep 17 00:00:00 2001
From: Michael Innerberger <michael.innerberger@asc.tuwien.ac.at>
Date: Wed, 19 Jul 2023 09:33:02 -0400
Subject: [PATCH] Switch addtive Schwarz implementation and delete old code
This somehow also fixes the problems with the higher order version of
the multilevel additive Schwarz preconditioner (this is tested now!)
---
lib/solvers/chooseIterativeSolver.m | 8 +-
.../AdditiveSchwarzHighOrderCascade.m | 141 -----------------
.../AdditiveSchwarzLowOrderCascade.m | 146 ------------------
.../OptimalMLAdditiveSchwarz.m | 76 +++++++++
.../preconditioners/P1AdditiveSchwarz.m | 120 --------------
tests/TestIterativeSolver.m | 1 +
6 files changed, 81 insertions(+), 411 deletions(-)
delete mode 100644 lib/solvers/preconditioners/AdditiveSchwarzHighOrderCascade.m
delete mode 100644 lib/solvers/preconditioners/AdditiveSchwarzLowOrderCascade.m
create mode 100644 lib/solvers/preconditioners/OptimalMLAdditiveSchwarz.m
delete mode 100644 lib/solvers/preconditioners/P1AdditiveSchwarz.m
diff --git a/lib/solvers/chooseIterativeSolver.m b/lib/solvers/chooseIterativeSolver.m
index 5058d49..55dbc42 100644
--- a/lib/solvers/chooseIterativeSolver.m
+++ b/lib/solvers/chooseIterativeSolver.m
@@ -33,15 +33,15 @@ switch class
end
case {"", "additiveSchwarzLowOrder"}
if order == 1
- preconditioner = P1AdditiveSchwarz(fes, blf, P);
+ preconditioner = OptimalMLAdditiveSchwarz(P1JacobiSmoother(fes, blf, P));
else
- preconditioner = AdditiveSchwarzLowOrderCascade(fes, blf);
+ preconditioner = OptimalMLAdditiveSchwarz(JacobiLowOrderCascadeSmoother(fes, blf));
end
case "additiveSchwarzHighOrder"
if order == 1
- preconditioner = P1AdditiveSchwarz(fes, blf, P);
+ preconditioner = OptimalMLAdditiveSchwarz(P1JacobiSmoother(fes, blf, P));
else
- preconditioner = AdditiveSchwarzHighOrderCascade(fes, blf, P);
+ preconditioner = OptimalMLAdditiveSchwarz(JacobiHighOrderCascadeSmoother(fes, blf, P));
end
otherwise
error('No PCG variant %s!', variant)
diff --git a/lib/solvers/preconditioners/AdditiveSchwarzHighOrderCascade.m b/lib/solvers/preconditioners/AdditiveSchwarzHighOrderCascade.m
deleted file mode 100644
index c0c4e1c..0000000
--- a/lib/solvers/preconditioners/AdditiveSchwarzHighOrderCascade.m
+++ /dev/null
@@ -1,141 +0,0 @@
-% AdditiveSchwarzHighOrderCascade (subclass of Preconditioner) optimal
-% multilevel additive Schwarz preconditioner for arbitrary order finite
-% elements: smooth with patchwise higher order matrix on changed patches on
-% all levels (local higher order smoothing) and globally with patchwise higher
-% order matrix on finest level (global higher order smoothing).
-%
-% See also: Preconditioner, PcgSolver
-
-% TODO: this does not work as intended! Test and fix it!
-classdef AdditiveSchwarzHighOrderCascade < Preconditioner
- %% properties
- properties (GetAccess=public,SetAccess=protected)
- nLevels
- end
-
- properties (Access=protected)
- blf
- hoFes
- loFes
- p1Acoarse
- hoAcoarse
- P
- inclusionMatrix
- intergridMatrix
- freeVertices
- freeDofs
- freeDofsOld
- changedPatches
- patchwiseA
- end
-
- properties (Access=private)
- listenerHandle
- end
-
- %% methods
- methods (Access=public)
- function obj = AdditiveSchwarzHighOrderCascade(fes, blf, P)
- arguments
- fes FeSpace
- blf BilinearForm
- P Prolongation
- end
-
- assert(fes.finiteElement.order > 1, ...
- 'Higher order additive Schwarz preconditioner only works for higher order finite elements.')
- assert(isempty(blf.b), ...
- 'Additive Schwarz preconditioner only tested for symmetric problems.')
-
- obj.nLevels = 0;
-
- mesh = fes.mesh;
- obj.loFes = FeSpace(mesh, LowestOrderH1Fe, 'dirichlet', fes.bnd.dirichlet);
- obj.hoFes = fes;
- obj.P = P;
- obj.blf = blf;
-
- obj.listenerHandle = mesh.listener('IsAboutToRefine', @obj.getChangedPatches);
- end
-
- function setup(obj, A)
- obj.nLevels = obj.nLevels + 1;
-
- L = obj.nLevels;
- obj.freeVertices = getFreeDofs(obj.loFes);
-
- obj.freeDofsOld = obj.freeDofs;
- obj.freeDofs = getFreeDofs(obj.hoFes);
-
- p1Matrix = assemble(obj.blf, obj.loFes);
- p1Matrix = p1Matrix(obj.freeVertices, obj.freeVertices);
-
- if L == 1
- obj.p1Acoarse = p1Matrix;
- obj.hoAcoarse = A;
- obj.inclusionMatrix = SpaceProlongation(obj.loFes, obj.hoFes) ...
- .matrix(getFreeDofs(obj.loFes), getFreeDofs(obj.hoFes));
- else
- obj.intergridMatrix{L} = obj.P.matrix(obj.freeDofs, obj.freeDofsOld);
- obj.changedPatches = obj.changedPatches(obj.freeVertices);
- if L == 2
- % since order changes to 2nd mesh, we have to do global smoothing
- obj.patchwiseA{L} = assemblePatchwise(obj.blf, obj.hoFes, ':');
- else
- obj.patchwiseA{L} = assemblePatchwise(obj.blf, obj.hoFes, obj.changedPatches);
- end
- end
- end
-
- % preconditioner: patchwise higher order solve on changed patches
- function Cx = apply(obj, x)
- assert(~isempty(obj.nLevels), 'Data for multilevel iteration not given!')
-
- L = obj.nLevels;
- rho = cell(L, 1);
-
- if L == 1
- Cx = obj.hoAcoarse \ x;
- return
- end
-
- % descending cascade
- for k = L:-1:2
- rho{k} = localSmoothing(obj, k, x);
- x = obj.intergridMatrix{k}' * x;
- end
-
- % exact solve on coarsest level
- x = obj.inclusionMatrix * x;
- Cx = obj.p1Acoarse \ x;
- Cx = obj.inclusionMatrix' * Cx;
-
- % ascending cascade
- for k = 2:L
- Cx = obj.intergridMatrix{k} * Cx;
- Cx = Cx + rho{k};
- end
- 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)
- % Selects vertices of bisected edges and all newly created
- % vertices
- nCOld = mesh.nCoordinates;
- nCNew = mesh.nCoordinates + nnz(bisecData.bisectedEdges) + ...
- bisecData.nRefinedElements' * bisecData.nInnerNodes;
- bisectedEdgeNodes = unique(mesh.edges(:,bisecData.bisectedEdges));
- obj.changedPatches = false(nCNew, 1);
- idx = [bisectedEdgeNodes; ((nCOld+1):nCNew)'];
- obj.changedPatches(idx) = true;
- end
-
- function rho = localSmoothing(obj, k, res)
- rho = obj.patchwiseA{k} \ res;
- end
- end
-end
\ No newline at end of file
diff --git a/lib/solvers/preconditioners/AdditiveSchwarzLowOrderCascade.m b/lib/solvers/preconditioners/AdditiveSchwarzLowOrderCascade.m
deleted file mode 100644
index f65917c..0000000
--- a/lib/solvers/preconditioners/AdditiveSchwarzLowOrderCascade.m
+++ /dev/null
@@ -1,146 +0,0 @@
-% AdditiveSchwarzLowOrderCascade (subclass of Preconditioner) optimal multilevel
-% additive Schwarz preconditioner for arbitrary order finite elements: smooth
-% with diagonal of P1 matrix on changed patches on every level (local P1
-% smoothing) and with patchwise higher order matrix on finest level (global
-% higher order smoothing).
-%
-% See also: Preconditioner, PcgSolver
-
-classdef AdditiveSchwarzLowOrderCascade < Preconditioner
- %% properties
- properties (GetAccess=public,SetAccess=protected)
- end
-
- properties (Access=protected)
- blf
- hoFes
- loFes
- P
- p1Acoarse
- hoAcoarse
- p1Smoother
- nLevels
- inclusionMatrix
- intergridMatrix
- freeVertices
- freeVerticesOld
- changedPatches
- patchwiseA
- patchwiseP1Matrix
- end
-
- properties (Access=private)
- listenerHandle
- end
-
- %% methods
- methods (Access=public)
- function obj = AdditiveSchwarzLowOrderCascade(fes, blf)
- arguments
- fes FeSpace
- blf BilinearForm
- end
-
- assert(fes.finiteElement.order > 1, ...
- 'Higher order additive Schwarz preconditioner only works for higher order finite elements.')
- assert(isempty(blf.b), ...
- 'Additive Schwarz preconditioner only tested for symmetric problems.')
-
- obj.nLevels = 0;
-
- mesh = fes.mesh;
- obj.loFes = FeSpace(mesh, LowestOrderH1Fe, 'dirichlet', fes.bnd.dirichlet);
- obj.hoFes = fes;
- obj.P = Prolongation.chooseFor(obj.loFes);
- obj.blf = blf;
-
- obj.listenerHandle = mesh.listener('IsAboutToRefine', @obj.getChangedPatches);
- end
-
- function setup(obj, A)
- obj.nLevels = obj.nLevels + 1;
-
- L = obj.nLevels;
- obj.freeVerticesOld = obj.freeVertices;
- obj.freeVertices = getFreeDofs(obj.loFes);
-
- p1Matrix = assemble(obj.blf, obj.loFes);
- p1Matrix = p1Matrix(obj.freeVertices, obj.freeVertices);
-
- if L == 1
- obj.p1Acoarse = p1Matrix;
- obj.hoAcoarse = A;
- else
- obj.p1Smoother{L} = full(diag(p1Matrix)).^(-1);
- obj.intergridMatrix{L} = obj.P.matrix(obj.freeVertices, obj.freeVerticesOld);
- obj.changedPatches{L} = obj.changedPatches{L}(obj.freeVertices);
- obj.patchwiseA = assemblePatchwise(obj.blf, obj.hoFes);
- obj.inclusionMatrix = SpaceProlongation(obj.loFes, obj.hoFes) ...
- .matrix(getFreeDofs(obj.loFes), getFreeDofs(obj.hoFes));
- end
- end
-
- % preconditioner: inverse of diagonal on each level
- function Cx = apply(obj, x)
- assert(~isempty(obj.nLevels), 'Data for multilevel iteration not given!')
-
- L = obj.nLevels;
- rho = cell(L, 1);
-
- if L == 1
- Cx = obj.hoAcoarse \ x;
- return
- end
-
- % descending cascade
- rho{L} = hoGlobalSmoothing(obj, x);
- x = obj.inclusionMatrix * x;
-
- for k = L:-1:3
- x = obj.intergridMatrix{k}' * x;
- rho{k-1} = p1LocalSmoothing(obj, k-1, x);
- end
-
- % exact solve on coarsest level
- x = obj.intergridMatrix{2}' * x;
- Cx = obj.p1Acoarse \ x;
- Cx = obj.intergridMatrix{2}*Cx;
-
- % ascending cascade
- for k = 3:L
- Cx = Cx + rho{k-1};
- Cx = obj.intergridMatrix{k} * Cx;
- end
-
- Cx = obj.inclusionMatrix' * Cx;
- Cx = Cx + rho{L};
- 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)
- % Selects vertices of bisected edges and all newly created
- % vertices
- 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 = [bisectedEdgeNodes; ((nCOld+1):nCNew)'];
- obj.changedPatches{obj.nLevels+1}(idx) = true;
- end
-
- function rho = p1LocalSmoothing(obj, k, res)
- idx = obj.changedPatches{k};
- rho = zeros(size(res));
- rho(idx,:) = obj.p1Smoother{k}(idx).*res(idx,:);
- end
-
- function rho = hoGlobalSmoothing(obj, res)
- rho = obj.patchwiseA \ res;
- end
- end
-end
\ No newline at end of file
diff --git a/lib/solvers/preconditioners/OptimalMLAdditiveSchwarz.m b/lib/solvers/preconditioners/OptimalMLAdditiveSchwarz.m
new file mode 100644
index 0000000..c7499bd
--- /dev/null
+++ b/lib/solvers/preconditioners/OptimalMLAdditiveSchwarz.m
@@ -0,0 +1,76 @@
+% OptimalMLAdditiveSchwarz (subclass of Preconditioner) optimal multilevel
+% additive Schwarz preconditioner with given smoother.
+%
+% See also: Preconditioner, PcgSolver, MultilevelSmoother
+
+classdef OptimalMLAdditiveSchwarz < Preconditioner
+ %% properties
+ properties (Access=protected)
+ smoother
+ Acoarse
+ projectedAcoarse
+ end
+
+ properties (Dependent)
+ nLevels
+ end
+
+ properties (Access=private)
+ listenerHandle
+ end
+
+ %% methods
+ methods (Access=public)
+ function obj = OptimalMLAdditiveSchwarz(smoother)
+ arguments
+ smoother MultilevelSmoother
+ end
+
+ obj = obj@Preconditioner();
+ obj.smoother = smoother;
+ end
+
+ function setup(obj, A)
+ if obj.nLevels == 0
+ obj.Acoarse = A;
+ obj.projectedAcoarse = obj.smoother.setup(A);
+ else
+ obj.smoother.setup(A);
+ end
+ end
+
+ function Cx = apply(obj, x)
+ assert(~isempty(obj.nLevels), 'Data for multilevel iteration not given!')
+
+ L = obj.nLevels;
+ rho = cell(L, 1);
+
+ if L == 1
+ Cx = obj.Acoarse \ x;
+ return
+ end
+
+ % descending cascade
+ residual = x;
+ for k = L:-1:2
+ rho{k} = obj.smoother.smooth(residual, k);
+ residual = obj.smoother.restrict(residual, k);
+ end
+
+ % exact solve on coarsest level
+ Cx = obj.projectedAcoarse \ residual;
+
+ % ascending cascade
+ for k = 2:L
+ Cx = obj.smoother.prolongate(Cx, k);
+ Cx = Cx + rho{k};
+ end
+ end
+ end
+
+ methods
+ function n = get.nLevels(obj)
+ n = obj.smoother.nLevels;
+ end
+ end
+end
\ No newline at end of file
diff --git a/lib/solvers/preconditioners/P1AdditiveSchwarz.m b/lib/solvers/preconditioners/P1AdditiveSchwarz.m
deleted file mode 100644
index a001f79..0000000
--- a/lib/solvers/preconditioners/P1AdditiveSchwarz.m
+++ /dev/null
@@ -1,120 +0,0 @@
-% P1AdditiveSchwarz (subclass of Preconditioner) optimal multilevel
-% additive Schwarz preconditioner for lowest order finite elements:
-% smooth with inverse diagonal on changed patches on every level.
-%
-% See also: Preconditioner, PcgSolver
-
-classdef P1AdditiveSchwarz < Preconditioner
- %% properties
- properties (Access=protected)
- fes
- blf
- P
- nLevels
- Acoarse
- smoother
- intergridMatrix
- freeDofs
- freeDofsOld
- changedPatches
- highestOrderIsOne
- end
-
- properties (Access=private)
- listenerHandle
- end
-
- %% methods
- methods (Access=public)
- function obj = P1AdditiveSchwarz(fes, blf, P)
- arguments
- fes FeSpace
- blf BilinearForm
- P Prolongation
- end
-
- assert(fes.finiteElement.order == 1, ...
- 'P1 Additive Schwarz preconditioner only works for lowest order finite elements.')
- assert(isempty(blf.b), ...
- 'Additive Schwarz preconditioner only tested for symmetric problems.')
-
- obj.nLevels = 0;
-
- mesh = fes.mesh;
- obj.fes = fes;
- obj.blf = blf;
- obj.P = P;
-
- obj.listenerHandle = mesh.listener('IsAboutToRefine', @obj.getChangedPatches);
- end
-
- function setup(obj, A)
- obj.nLevels = obj.nLevels + 1;
-
- L = obj.nLevels;
- obj.freeDofsOld = obj.freeDofs;
- obj.freeDofs = getFreeDofs(obj.fes);
-
- if L == 1
- obj.Acoarse = A;
- else
- obj.smoother{L} = full(diag(A)).^(-1);
- obj.intergridMatrix{L} = obj.P.matrix(obj.freeDofs, obj.freeDofsOld);
- obj.changedPatches{L} = find(obj.changedPatches{L}(obj.freeDofs));
- end
- end
-
- % action: inverse of diagonal on each level
- function Cx = apply(obj, x)
- assert(~isempty(obj.nLevels), 'Data for multilevel iteration not given!')
-
- L = obj.nLevels;
- rho = cell(L, 1);
-
- if L == 1
- Cx = obj.Acoarse \ x;
- return
- end
-
- % descending cascade
- rho{L} = localSmoothing(obj, L, x);
- residual = x;
- for k = L-1:-1:2
- residual = obj.intergridMatrix{k+1}' * residual;
- rho{k} = localSmoothing(obj, k, residual);
- end
-
- % exact solve on coarsest level
- residual = obj.intergridMatrix{2}' * residual;
- Cx = obj.Acoarse \ residual;
-
- % ascending cascade
- for k = 2:L-1
- Cx = obj.intergridMatrix{k} * Cx;
- Cx = Cx + rho{k};
- end
- Cx = obj.intergridMatrix{L} * Cx + rho{L};
- 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
-
- function rho = localSmoothing(obj, k, res)
- idx = obj.changedPatches{k};
- rho = zeros(size(res));
- rho(idx,:) = obj.smoother{k}(idx).*res(idx,:);
- end
- end
-end
\ No newline at end of file
diff --git a/tests/TestIterativeSolver.m b/tests/TestIterativeSolver.m
index a4f516c..cbde884 100644
--- a/tests/TestIterativeSolver.m
+++ b/tests/TestIterativeSolver.m
@@ -7,6 +7,7 @@ properties (TestParameter)
'jacobiPCG', ["pcg", "jacobi"], ...
'iChol', ["pcg", "iChol"], ...
'loAdditiveSchwarzPCG', ["pcg", "additiveSchwarzLowOrder"], ...
+ 'hoAdditiveSchwarzPCG', ["pcg", "additiveSchwarzHighOrder"], ...
'lowMG', ["multigrid", "lowOrderVcycle"], ...
'highMG', ["multigrid", "highOrderVcycle"]);
end
--
GitLab