diff --git a/examples/algebraicSolver.m b/examples/algebraicSolver.m
index 15084f31b955b2fc67e5442a9a104c4c26881c7e..3ef04b1eb3a6f50650b74ff9ad1cf3380613f7b9 100644
--- a/examples/algebraicSolver.m
+++ b/examples/algebraicSolver.m
@@ -33,7 +33,7 @@ for p = 1:pmax
% choose the iterative solver of choice: multigrid (with the variants
% lowOrderVcycle and highOrderVcycle), pcg (with jacobi and additive
% Schwarz preconditioner), and the cg solver
- [solver, P] = chooseIterativeSolver(fes, blf, 'multigrid', 'lowOrderVcycle');
+ [solver, P] = IterativeSolver.choose(fes, blf, 'multigrid', 'lowOrderVcycle');
solver.tol = 1e-6;
solver.maxIter = 100;
diff --git a/examples/goafemIterativeSolver.m b/examples/goafemIterativeSolver.m
index 801bbfb0ba1caaa3660e4ee777db7969d4fcc842..50a7df3835d03592ed0f894e33081a62b8935177 100644
--- a/examples/goafemIterativeSolver.m
+++ b/examples/goafemIterativeSolver.m
@@ -37,7 +37,7 @@ for p = 1:pmax
lfG.qrf = QuadratureRule.ofOrder(2*p);
%% set up solver and lifting operator for nested iteration
- [solver, P] = chooseIterativeSolver(fes, blf, 'multigrid', 'lowOrderVcycle');
+ [solver, P] = IterativeSolver.choose(fes, blf, 'multigrid', 'lowOrderVcycle');
solver.tol = 1e-6;
solver.maxIter = 1000;
diff --git a/lib/solvers/IterativeSolver.m b/lib/solvers/@IterativeSolver/IterativeSolver.m
similarity index 92%
rename from lib/solvers/IterativeSolver.m
rename to lib/solvers/@IterativeSolver/IterativeSolver.m
index 191676e17144f503abcdefd31db89cafbba62aec..be5b90b45bbca63d72b0cb500295acc8cb234787 100644
--- a/lib/solvers/IterativeSolver.m
+++ b/lib/solvers/@IterativeSolver/IterativeSolver.m
@@ -25,6 +25,10 @@
%
% tf = solver.isConverged(solver) returns state of convergence of the solver
% for each column of the right-hand side.
+%
+% [solver, P] = IterativeSolver.choose(fes, blf, class, variant) convenience
+% factory method. Returns instances of IterativeSolver and Prolongation
+% accoding to the given parameters.
classdef IterativeSolver < handle
@@ -102,6 +106,11 @@ classdef IterativeSolver < handle
obj.activeComponents = obj.activeComponents & ~tf;
end
end
+
+ %% convenience factory function
+ methods (Static, Access=public)
+ [solver, P] = choose(fes, blf, class, variant)
+ end
%% validation functions to use within this class
methods (Static, Access=protected)
diff --git a/lib/solvers/@IterativeSolver/choose.m b/lib/solvers/@IterativeSolver/choose.m
new file mode 100644
index 0000000000000000000000000000000000000000..ef5d4fa6dbc32592f5c496f354ff839080dbd3a7
--- /dev/null
+++ b/lib/solvers/@IterativeSolver/choose.m
@@ -0,0 +1,91 @@
+% choose Return suitable solver (instance of subclass of IterativeSolver) and
+% suitable Prolongation P.
+%
+% solver = choose(fes, blf, class) returns a solver of class 'class' for given
+% finite element space fes and bilinear form blf. For each class, sensible
+% default variants are chosen.
+%
+% solver = choose(fes, blf, class, variant) returns a solver of class 'class'
+% of specific variant 'variant'.
+%
+% [solver, P] = choose(__) additionally to the solver, returns a suitable
+% Prolongation P, that can be used to prolongate the solution of the
+% solver to the next finer mesh.
+
+function [solver, P] = choose(fes, blf, class, variant)
+arguments
+ fes FeSpace
+ blf BilinearForm
+ class (1,1) string {mustBeMember(class, ["cg", "pcg", "multigrid", "direct"])}
+ variant (1,1) string {mustBeMember(variant, [ "", ...
+ "iChol", "jacobi", ...
+ "additiveSchwarzLowOrder", "additiveSchwarzHighOrder" ...
+ "lowOrderVcycle", "highOrderVcycle"])} = ""
+end
+
+order = fes.finiteElement.order;
+P = Prolongation.chooseFor(fes);
+
+switch class
+ % non-preconditioned CG
+ case "cg"
+ solver = PcgSolver(NoPreconditioner());
+
+ % preconditioned CG family
+ case "pcg"
+ switch variant
+ case "iChol"
+ preconditioner = IncompleteCholesky();
+ case "jacobi"
+ if order == 1
+ preconditioner = DiagonalJacobi();
+ else
+ preconditioner = BlockJacobi(fes, blf);
+ end
+ case {"", "additiveSchwarzLowOrder"}
+ if order == 1
+ preconditioner = OptimalMLAdditiveSchwarz(P1JacobiCascade(fes, blf, P));
+ else
+ preconditioner = OptimalMLAdditiveSchwarz(JacobiLowOrderCascade(fes, blf));
+ end
+ case "additiveSchwarzHighOrder"
+ if order == 1
+ preconditioner = OptimalMLAdditiveSchwarz(P1JacobiCascade(fes, blf, P));
+ else
+ preconditioner = OptimalMLAdditiveSchwarz(JacobiHighOrderCascade(fes, blf, P));
+ end
+ otherwise
+ error('No PCG variant %s!', variant)
+ end
+ solver = PcgSolver(preconditioner);
+
+ % geometric , Pmultigrid family
+ case "multigrid"
+ switch variant
+ case {"", "lowOrderVcycle"}
+ if order == 1
+ smoother = P1JacobiCascade(fes, blf, P);
+ else
+ smoother = JacobiLowOrderCascade(fes, blf);
+ end
+ case "highOrderVcycle"
+ if order == 1
+ smoother = P1JacobiCascade(fes, blf, P);
+ else
+ smoother = JacobiHighOrderCascade(fes, blf, P);
+ end
+ otherwise
+ error('No multigrid variant %s!', variant)
+ end
+ solver = OptimalVcycleMultigridSolver(smoother);
+
+ % direct solve (for testing purposes)
+ case "direct"
+ solver = DirectSolver();
+
+ % default
+ otherwise
+ error('No iterative solver of class %s!', class)
+end
+
+end
diff --git a/lib/solvers/AdditiveSchwarzHighOrderPcg.m b/lib/solvers/AdditiveSchwarzHighOrderPcg.m
deleted file mode 100644
index 076f2fdd50e58d4262e85366e13f98763f060221..0000000000000000000000000000000000000000
--- a/lib/solvers/AdditiveSchwarzHighOrderPcg.m
+++ /dev/null
@@ -1,161 +0,0 @@
-% AdditiveSchwarzHighOrderPcg (subclass of PcgSolver) Solves linear equations
-% iteratively using the CG method with optimal multilevel additive
-% Schwartz preconditioner.
-
-classdef AdditiveSchwarzHighOrderPcg < PcgSolver
- %% properties
- properties (GetAccess=public,SetAccess=protected)
- nLevels
- end
-
- properties (Access=protected)
- blf
- hoFes
- loFes
- coarseLoFes
- coarseHoFes
- p1Acoarse
- smoother
- P
- intergridMatrix
- freeVertices
- freeVerticesOld
- freeDofs
- freeDofsOld
- changedPatches
- patchwiseA
- end
-
- properties (Access=private)
- listenerHandle
- end
-
- %% methods
- methods (Access=public)
- function obj = AdditiveSchwarzHighOrderPcg(fes, blf, P)
- arguments
- fes FeSpace
- blf BilinearForm
- P Prolongation
- end
- obj = obj@PcgSolver();
-
- assert(fes.finiteElement.order > 1, ...
- 'AdditiveSchwartzPcg only works for higher order finite elements.')
- assert(isempty(blf.b), ...
- 'Additive Schwarz PCG solvers only tested for symmetric problems.')
-
- obj.nLevels = 0;
-
- mesh = fes.mesh;
- obj.loFes = FeSpace(mesh, LowestOrderH1Fe, 'dirichlet', ':');
- obj.hoFes = fes;
- obj.P = P;
- obj.blf = blf;
-
- obj.listenerHandle = mesh.listener('IsAboutToRefine', @obj.getChangedPatches);
- end
-
- function setupSystemMatrix(obj, A)
- obj.nLevels = obj.nLevels + 1;
-
- L = obj.nLevels;
- obj.freeVerticesOld = obj.freeVertices;
- obj.freeVertices = getFreeDofs(obj.loFes);
-
- obj.freeDofsOld = obj.freeDofs;
- obj.freeDofs = getFreeDofs(obj.hoFes);
-
- ppMatrix = assemble(obj.blf, obj.hoFes);
- ppMatrix = ppMatrix(getFreeDofs(obj.hoFes), getFreeDofs(obj.hoFes));
-
- p1Matrix = assemble(obj.blf, obj.loFes);
- p1Matrix = p1Matrix(obj.freeVertices, obj.freeVertices);
-
- if L == 1
- obj.p1Acoarse = p1Matrix;
- fes = obj.hoFes;
- 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);
- else
- obj.smoother{L} = full(diag(ppMatrix)).^(-1);
- obj.intergridMatrix{L} = obj.P.matrix(obj.freeDofs, obj.freeDofsOld);
- obj.changedPatches{L} = obj.changedPatches{L}(obj.freeVertices);
- obj.patchwiseA{L} = assemblePatchwise(obj.blf, obj.hoFes, obj.changedPatches{L});
- end
-
- setupSystemMatrix@PcgSolver(obj, A);
- end
-
- function setupRhs(obj, b, varargin)
- setupRhs@PcgSolver(obj, b, varargin{:});
- end
-
- % preconditioner: inverse of diagonal on each level
- function Cx = preconditionAction(obj, residual)
- assert(~isempty(obj.nLevels), 'Data for multilevel iteration not given!')
-
- L = obj.nLevels;
- rho = cell(L, 1);
-
- if L == 1
- Cx = obj.A \ residual;
- return
- end
-
- % descending cascade
- for k = L:-1:2
- rho{k} = localSmoothing(obj, k, residual);
- residual = obj.intergridMatrix{k}'*residual;
- end
-
- % exact solve on coarsest level
- residual = interpolateFreeData(residual, obj.coarseHoFes, obj.coarseLoFes);
- sigma = obj.p1Acoarse \ residual;
- sigma = interpolateFreeData(sigma, obj.coarseLoFes, obj.coarseHoFes);
-
- % ascending cascade
- for k = 2:L
- sigma = obj.intergridMatrix{k}*sigma;
- sigma = sigma + rho{k};
- 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)
- % 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 = localSmoothing(obj, k, res)
- rho = obj.patchwiseA{k} \ res;
- end
- end
-end
-
- %% auxiliary functions
-function interpolatedData = interpolateFreeData(data, fromFes, toFes)
- freeDofs = getFreeDofs(toFes);
- nComponents = size(data, 2);
- interpolatedData = zeros(numel(freeDofs), nComponents);
- feFunctionWrapper = FeFunction(fromFes);
- for k = 1:nComponents
- feFunctionWrapper.setFreeData(data(:,k));
- wholeData = nodalInterpolation(feFunctionWrapper, toFes);
- interpolatedData(:,k) = wholeData(freeDofs);
- end
-end
\ No newline at end of file
diff --git a/lib/solvers/AdditiveSchwarzLowOrderPcg.m b/lib/solvers/AdditiveSchwarzLowOrderPcg.m
deleted file mode 100644
index 448e069da4be9e6797e538c9896ae7c2107ba1de..0000000000000000000000000000000000000000
--- a/lib/solvers/AdditiveSchwarzLowOrderPcg.m
+++ /dev/null
@@ -1,152 +0,0 @@
-% AdditiveSchwartzPcg (subclass of PcgSolver) Solves linear equations
-% iteratively using the CG method with optimal multilevel additive
-% Schwartz preconditioner.
-
-classdef AdditiveSchwarzLowOrderPcg < PcgSolver
- %% properties
- properties (GetAccess=public,SetAccess=protected)
- nLevels
- end
-
- properties (Access=protected)
- blf
- hoFes
- loFes
- p1Acoarse
- p1Smoother
- P
- inclusionMatrix
- intergridMatrix
- freeVertices
- freeVerticesOld
- changedPatches
- patchwiseA
- patchwiseP1Matrix
- end
-
- properties (Access=private)
- listenerHandle
- end
-
- %% methods
- methods (Access=public)
- function obj = AdditiveSchwarzLowOrderPcg(fes, blf)
- arguments
- fes FeSpace
- blf BilinearForm
- end
- obj = obj@PcgSolver();
-
- assert(fes.finiteElement.order > 1, ...
- 'AdditiveSchwartzPcg only works for higher order finite elements.')
- assert(isempty(blf.b), ...
- 'Additive Schwarz PCG solvers only tested for symmetric problems.')
-
- obj.nLevels = 0;
-
- mesh = fes.mesh;
- obj.loFes = FeSpace(mesh, LowestOrderH1Fe, 'dirichlet', ':');
- obj.hoFes = fes;
- obj.P = Prolongation.chooseFor(obj.loFes);
- obj.blf = blf;
-
- obj.listenerHandle = mesh.listener('IsAboutToRefine', @obj.getChangedPatches);
- end
-
- function setupSystemMatrix(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;
- 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
-
- setupSystemMatrix@PcgSolver(obj, A);
- end
-
- function setupRhs(obj, b, varargin)
- setupRhs@PcgSolver(obj, b, varargin{:});
- end
-
- % preconditioner: inverse of diagonal on each level
- function Cx = preconditionAction(obj, residual)
- assert(~isempty(obj.nLevels), 'Data for multilevel iteration not given!')
-
- L = obj.nLevels;
- rho = cell(L, 1);
-
- if L == 1
- Cx = obj.A \ residual;
- return
- end
-
- % descending cascade
- rho{L} = hoGlobalSmoothing(obj, residual);
-
- residual = obj.inclusionMatrix * residual;
-
- for k = L:-1:3
- residual = obj.intergridMatrix{k}'*residual;
- rho{k-1} = p1LocalSmoothing(obj, k-1, residual);
- end
-
- % exact solve on coarsest level
- residual = obj.intergridMatrix{2}'*residual;
- sigma = obj.p1Acoarse \ residual;
- sigma = obj.intergridMatrix{2}*sigma;
-
- % ascending cascade
- for k = 3:L
- sigma = sigma + rho{k-1};
- sigma = obj.intergridMatrix{k}*sigma;
- end
-
- sigma = obj.inclusionMatrix' * sigma;
- sigma = sigma + rho{L};
-
- 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)
- % 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,:);
- % DEBUG: Local patchwise smoothing in P1
- % rho = obj.patchwiseP1Matrix{k} \ res;
- end
-
- function rho = hoGlobalSmoothing(obj, res)
- rho = obj.patchwiseA \ res;
- end
- end
-end
\ No newline at end of file
diff --git a/lib/solvers/BlockJacobiPcgSolver.m b/lib/solvers/BlockJacobiPcgSolver.m
deleted file mode 100644
index 71cbd08b19cc02f0f64c52d8be600c365408b470..0000000000000000000000000000000000000000
--- a/lib/solvers/BlockJacobiPcgSolver.m
+++ /dev/null
@@ -1,42 +0,0 @@
-% BlockJacobiPcgSolver (subclass of PcgSolver) Solves linear equations
-% iteratively using the CG method with block-diagonal preconditioner.
-
-classdef BlockJacobiPcgSolver < PcgSolver
- %% properties
- properties (GetAccess=public,SetAccess=protected)
- fes
- end
-
- properties (Access=private)
- C
- blf
- end
-
- %% methods
- methods (Access=public)
- % preconditioner action: patchwise inverse of diagonal
- function obj = BlockJacobiPcgSolver(fes, blf)
- arguments
- fes FeSpace
- blf BilinearForm
- end
-
- obj = obj@PcgSolver();
- obj.fes = fes;
- obj.blf = blf;
- end
-
- function setupSystemMatrix(obj, A)
- obj.C = assemblePatchwise(obj.blf, obj.fes, ':');
- setupSystemMatrix@PcgSolver(obj, A);
- end
-
- function setupRhs(obj, b, varargin)
- setupRhs@PcgSolver(obj, b, varargin{:});
- end
-
- function Cx = preconditionAction(obj, x)
- Cx = obj.C \ x;
- end
- end
-end
\ No newline at end of file
diff --git a/lib/solvers/ICholPcgSolver.m b/lib/solvers/ICholPcgSolver.m
deleted file mode 100644
index 58a300faaf4ad73d33336da2f8b4ccbfa093025c..0000000000000000000000000000000000000000
--- a/lib/solvers/ICholPcgSolver.m
+++ /dev/null
@@ -1,39 +0,0 @@
-% ICholPcgSolver (subclass of PcgSolver) Solves linear equations iteratively
-% using the CG method with incomplete Cholesky decomposition preconditioner.
-
-classdef ICholPcgSolver < PcgSolver
- %% properties
- properties (GetAccess=public,SetAccess=protected)
- fes
- end
-
- properties (Access=private)
- C
- end
-
- %% methods
- methods (Access=public)
- % preconditioner action: inverse of diagonal
- function obj = ICholPcgSolver(fes)
- arguments
- fes FeSpace
- end
-
- obj = obj@PcgSolver();
- obj.fes = fes;
- end
-
- function setupSystemMatrix(obj, A)
- obj.C = ichol(A);
- setupSystemMatrix@PcgSolver(obj, A);
- end
-
- function setupRhs(obj, b, varargin)
- setupRhs@PcgSolver(obj, b, varargin{:});
- end
-
- function Cx = preconditionAction(obj, x)
- Cx = obj.C' \ (obj.C \ x);
- end
- end
-end
\ No newline at end of file
diff --git a/lib/solvers/JacobiPcgSolver.m b/lib/solvers/JacobiPcgSolver.m
deleted file mode 100644
index 0b54610008a77e2cc73a6f0630da6addb03b760f..0000000000000000000000000000000000000000
--- a/lib/solvers/JacobiPcgSolver.m
+++ /dev/null
@@ -1,39 +0,0 @@
-% JacobiPcgSolver (subclass of PcgSolver) Solves linear equations
-% iteratively using the CG method with diagonal preconditioner.
-
-classdef JacobiPcgSolver < PcgSolver
- %% properties
- properties (GetAccess=public,SetAccess=protected)
- fes
- end
-
- properties (Access=private)
- C
- end
-
- %% methods
- methods (Access=public)
- % preconditioner action: inverse of diagonal
- function obj = JacobiPcgSolver(fes)
- arguments
- fes FeSpace
- end
-
- obj = obj@PcgSolver();
- obj.fes = fes;
- end
-
- function setupSystemMatrix(obj, A)
- obj.C = full(diag(A)).^(-1);
- setupSystemMatrix@PcgSolver(obj, A);
- end
-
- function setupRhs(obj, b, varargin)
- setupRhs@PcgSolver(obj, b, varargin{:});
- end
-
- function Cx = preconditionAction(obj, x)
- Cx = obj.C .* x;
- end
- end
-end
\ No newline at end of file
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/LocalMgLowOrderVcycle.m b/lib/solvers/LocalMgLowOrderVcycle.m
deleted file mode 100644
index f8920c03d1057f10ed21b3f411112f4e4c83abd5..0000000000000000000000000000000000000000
--- a/lib/solvers/LocalMgLowOrderVcycle.m
+++ /dev/null
@@ -1,179 +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 LocalMgLowOrderVcycle < MGSolver
- %% properties
- properties (GetAccess=public,SetAccess=protected)
- nLevels
- end
-
- properties (Access=protected)
- blf
- hoFes
- loFes
- p1Matrix
- p1Smoother
- P
- intergridMatrix
- inclusionMatrix
- freeVertices
- freeVerticesOld
- changedPatches
- patchwiseA
- end
-
- %% event data
- properties (Access=private)
- listenerHandle
- end
-
- %% methods
- methods (Access=public)
- function obj = LocalMgLowOrderVcycle(fes, blf)
- arguments
- fes FeSpace
- blf BilinearForm
- 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.loFes = FeSpace(mesh, LowestOrderH1Fe(), 'dirichlet', ':');
- obj.hoFes = fes;
- obj.P = Prolongation.chooseFor(obj.loFes);
- obj.blf = blf;
-
- obj.listenerHandle = mesh.listener('IsAboutToRefine', @obj.getChangedPatches);
- end
-
- function setupSystemMatrix(obj, A)
- obj.nLevels = obj.nLevels + 1;
-
- L = obj.nLevels;
- obj.freeVerticesOld = obj.freeVertices;
- obj.freeVertices = getFreeDofs(obj.loFes);
- freeVerticesHigher = getFreeDofs(obj.hoFes);
-
- obj.p1Matrix{L} = assemble(obj.blf, obj.loFes);
- obj.p1Matrix{L} = obj.p1Matrix{L}(obj.freeVertices, obj.freeVertices);
- obj.p1Smoother{L} = full(diag(obj.p1Matrix{L})).^(-1);
-
- if L >= 2
- obj.intergridMatrix{L} = obj.P.matrix(obj.freeVertices, obj.freeVerticesOld);
- obj.changedPatches{L} = find(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
- setupSystemMatrix@MGSolver(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 P1: no smoothing
- residual{L} = obj.inclusionMatrix * res; %interpolateFreeData(res, obj.hoFes, obj.loFes);
- for k = L:-1:2
- residual{k-1} = obj.intergridMatrix{k}'*residual{k};
- end
-
- % exact solve on coarsest level to compute accumulative lifting
- % of the residual (sigma)
- sigma = obj.p1Matrix{1} \ residual{1};
- algError2 = algError2 + scalarProduct(sigma, obj.p1Matrix{1}*sigma);
-
- % ascending cascade in P1 AND local smoothing
- for k = 2:(L-1)
- sigma = obj.intergridMatrix{k}*sigma;
- uptres = residual{k} - obj.p1Matrix{k}*sigma; %updated residual
- rho = p1LocalSmoothing(obj, k, uptres);
- [aeUpd, sUpd] = computeOptimalUpdate(obj.p1Matrix{k}, uptres, rho);
- sigma = sigma + sUpd;
- algError2 = algError2 + aeUpd;
- end
-
- % smoothing on finest level dependent on p
- sigma = obj.intergridMatrix{L}*sigma;
- sigma = obj.inclusionMatrix' * sigma; %interpolateFreeData(sigma, obj.loFes, obj.hoFes);
- uptres = res - obj.A*sigma;
- rho = hoGlobalSmoothing(obj, uptres);
- [aeUpd, sUpd] = computeOptimalUpdate(obj.A, uptres, rho);
- sigma = sigma + sUpd;
- algError2 = algError2 + aeUpd;
-
- 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
-
- 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
-
-%% auxiliary functions
-% error correction with optimal stepsize
-function [etaUpdate, sigmaUpdate] = computeOptimalUpdate(A, res, rho)
- rhoArho = scalarProduct(rho, A*rho);
- if max(abs(rho)) < eps
- lambda = 1;
- else
- lambda = scalarProduct(res, rho) ./ rhoArho;
- end
- 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
diff --git a/lib/solvers/LowestOrderAdditiveSchwarzPcg.m b/lib/solvers/LowestOrderAdditiveSchwarzPcg.m
deleted file mode 100644
index 9e19823e6c22b4141f296dd8ae4d903e33d143c3..0000000000000000000000000000000000000000
--- a/lib/solvers/LowestOrderAdditiveSchwarzPcg.m
+++ /dev/null
@@ -1,130 +0,0 @@
-% LowestOrderAdditiveSchwartzPcg (subclass of PcgSolver) Solves linear
-% equations iteratively using the CG method with optimal multilevel
-% additive Schwartz preconditioner for lowest order finite elements.
-
-classdef LowestOrderAdditiveSchwarzPcg < PcgSolver
- %% properties
- properties (GetAccess=public,SetAccess=protected)
- nLevels
- end
-
- properties (Access=protected)
- blf
- fes
- matrix
- smoother
- P
- intergridMatrix
- freeDofs
- freeDofsOld
- changedPatches
- highestOrderIsOne
- end
-
- properties (Access=private)
- listenerHandle
- end
-
- %% methods
- methods (Access=public)
- function obj = LowestOrderAdditiveSchwarzPcg(fes, blf, P)
- arguments
- fes FeSpace
- blf BilinearForm
- P Prolongation
- end
- obj = obj@PcgSolver();
-
- assert(fes.finiteElement.order == 1, ...
- 'LowestOrderAdditiveSchwartzPcg only works for lowest order finite elements.')
- assert(isempty(blf.b), ...
- 'Additive Schwarz PCG solvers only tested for symmetric problems.')
-
- obj.nLevels = 0;
-
- mesh = fes.mesh;
- obj.fes = fes;
- obj.P = P;
- obj.blf = blf;
-
- obj.listenerHandle = mesh.listener('IsAboutToRefine', @obj.getChangedPatches);
- end
-
- function setupSystemMatrix(obj, A)
- obj.nLevels = obj.nLevels + 1;
-
- L = obj.nLevels;
- obj.freeDofsOld = obj.freeDofs;
- obj.freeDofs = getFreeDofs(obj.fes);
-
- obj.matrix{L} = A;
- obj.smoother{L} = full(diag(obj.matrix{L})).^(-1);
-
- if L >= 2
- obj.intergridMatrix{L} = obj.P.matrix(obj.freeDofs, obj.freeDofsOld);
- obj.changedPatches{L} = find(obj.changedPatches{L}(obj.freeDofs));
- end
-
- setupSystemMatrix@PcgSolver(obj, A);
- end
-
- function setupRhs(obj, b, varargin)
- setupRhs@PcgSolver(obj, b, varargin{:});
- end
-
- % preconditioner: inverse of diagonal on each level
- function Cx = preconditionAction(obj, res)
- assert(~isempty(obj.nLevels), 'Data for multilevel iteration not given!')
-
- L = obj.nLevels;
- rho = cell(L, 1);
-
- if L == 1
- Cx = obj.A \ res;
- return
- end
-
- % descending cascade
- rho{L} = localSmoothing(obj, L, res);
- residual = res;
- 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;
- sigma = obj.matrix{1} \ residual;
-
- % ascending cascade
- for k = 2:L-1
- sigma = obj.intergridMatrix{k}*sigma;
- sigma = sigma + rho{k};
- end
- sigma = obj.intergridMatrix{L}*sigma + rho{L};
-
- 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
-
- 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/lib/solvers/LowestOrderLocalMg.m b/lib/solvers/LowestOrderLocalMg.m
deleted file mode 100644
index cde8fb1d573b76a5d2de1c0927f0222c1637b3a9..0000000000000000000000000000000000000000
--- a/lib/solvers/LowestOrderLocalMg.m
+++ /dev/null
@@ -1,162 +0,0 @@
-% LowestOrderLocalMg (subclass of MGSolver) Solves linear equations
-% stemming from lowest order AFEM computations iteratively using a
-% multigrid solver where one iteration is characterized by Vcycle(0,1)
-% (no pre-/one post-smoothing step).
-% On every level, local Jacobi smoothing is employed on changed patches
-% and the smoothed update is weighted by an optimal step size.
-%
-% See also [Innerberger, Miraçi, Praetorius, Streitberger; Optimal
-% computational costs of AFEM with optimal local hp-robust multigrid
-% solver].
-
-classdef LowestOrderLocalMg < MGSolver
- %% properties
- properties (GetAccess=public,SetAccess=protected)
- nLevels
- end
-
- properties (Access=protected)
- fes
- blf
- matrix
- smoother
- P
- intergridMatrix
- freeVertices
- freeVerticesOld
- changedPatches
- end
-
- %% event data
- properties (Access=private)
- listenerHandle
- end
-
- %% methods
- methods (Access=public)
- function obj = LowestOrderLocalMg(fes, blf, P)
- arguments
- fes FeSpace
- blf BilinearForm
- P Prolongation
- end
- obj = obj@MGSolver();
-
- assert(fes.finiteElement.order == 1, ...
- 'LowestOrderLocalMg only works for lowest order finite elements.')
- assert(isempty(blf.b), ...
- 'Multigrid solvers only tested for symmetric problems.')
-
- obj.nLevels = 0;
-
- mesh = fes.mesh;
- obj.fes = fes;
- obj.P = P;
- obj.blf = blf;
-
- obj.listenerHandle = mesh.listener('IsAboutToRefine', @obj.getChangedPatches);
- end
-
- function setupSystemMatrix(obj, A)
- obj.nLevels = obj.nLevels + 1;
-
- L = obj.nLevels;
- obj.freeVerticesOld = obj.freeVertices;
- obj.freeVertices = getFreeDofs(obj.fes);
-
- obj.matrix{L} = A;
- obj.smoother{L} = full(diag(obj.matrix{L})).^(-1);
-
- if L >= 2
- obj.intergridMatrix{L} = obj.P.matrix(obj.freeVertices, obj.freeVerticesOld);
- obj.changedPatches{L} = find(obj.changedPatches{L}(obj.freeVertices));
- end
- setupSystemMatrix@MGSolver(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 P1: no smoothing
- residual{L} = res;
- for k = L:-1:2
- residual{k-1} = obj.intergridMatrix{k}'*residual{k};
- end
-
- % exact solve on coarsest level to compute accumulative lifting
- % of the residual (sigma)
- sigma = obj.matrix{1} \ residual{1};
- algError2 = algError2 + scalarProduct(sigma, obj.matrix{1}*sigma);
-
- % ascending cascade in P1 AND local smoothing
- for k = 2:L
- sigma = obj.intergridMatrix{k}*sigma;
- uptres = residual{k} - obj.matrix{k}*sigma; %updated residual
- rho = localSmoothing(obj, k, uptres);
- [aeUpd, sUpd] = computeOptimalUpdate(obj.matrix{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
-
- 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
-
-%% auxiliary functions
-% error correction with optimal stepsize
-function [etaUpdate, sigmaUpdate] = computeOptimalUpdate(A, res, rho)
- rhoArho = scalarProduct(rho, A*rho);
- if max(abs(rho)) < eps
- lambda = 1;
- else
- lambda = scalarProduct(res, rho) ./ rhoArho;
- end
- 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
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/OptimalVcycleMultigridSolver.m b/lib/solvers/OptimalVcycleMultigridSolver.m
new file mode 100644
index 0000000000000000000000000000000000000000..496c2a73082af19af51919bf7ebb1c0ae04f8572
--- /dev/null
+++ b/lib/solvers/OptimalVcycleMultigridSolver.m
@@ -0,0 +1,144 @@
+% OptimalVcycleMultigridSolver (subclass of IterativeSolver) Solves linear
+% 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
+% 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
+ projectedMatrix
+ 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{:});
+
+ projectedA = obj.smoother.setup(A, varargin{:});
+ obj.projectedMatrix{obj.nLevels} = projectedA;
+ 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.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-1)
+ Cx = obj.smoother.prolongate(Cx, k);
+ updatedResidual = res{k} - obj.projectedMatrix{k} * Cx;
+ rho = obj.smoother.smooth(updatedResidual, k);
+ [etaUpdate2, xUpdate] = computeOptimalUpdate(obj.projectedMatrix{k}, updatedResidual, rho);
+ Cx = Cx + xUpdate;
+ eta2 = eta2 + etaUpdate2;
+ end
+
+ % 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
+
+%% 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/PcgSolver.m b/lib/solvers/PcgSolver.m
index 3a21ecc58737624493b4679cdb3dbfe9f211c153..1a278c82dda3a969218f6e1d6f730baa11c4bdcc 100644
--- a/lib/solvers/PcgSolver.m
+++ b/lib/solvers/PcgSolver.m
@@ -10,27 +10,41 @@ classdef PcgSolver < IterativeSolver
end
properties (Access=protected)
+ C
residual
Cresidual
searchDirection
normb
end
+
+ %% constructor
+ methods (Access=public)
+ function obj = PcgSolver(preconditioner)
+ arguments
+ preconditioner (1,1) Preconditioner
+ end
+
+ obj = obj@IterativeSolver();
+ obj.C = preconditioner;
+ end
+ end
%% extend superclass methods
methods (Access=public)
function setupSystemMatrix(obj, A)
setupSystemMatrix@IterativeSolver(obj, A);
+ obj.C.setup(A);
end
function setupRhs(obj, b, varargin)
setupRhs@IterativeSolver(obj, b, varargin{:});
% initialize residual & search direction
- obj.residual = b - obj.A*obj.x;
- obj.Cresidual = obj.preconditionAction(obj.residual);
+ obj.residual = b - obj.A * obj.x;
+ obj.Cresidual = obj.C.apply(obj.residual);
obj.searchDirection = obj.Cresidual;
- obj.residualCNorm = sum(obj.residual.*obj.Cresidual, 1);
- obj.normb = sqrt(sum(b.^2, 1));
+ obj.residualCNorm = sqrt(dot(obj.residual, obj.Cresidual, 1));
+ obj.normb = vecnorm(b, 2, 1);
end
end
@@ -38,8 +52,8 @@ classdef PcgSolver < IterativeSolver
methods (Access=public)
function tf = isConverged(obj)
tf = ((obj.iterationCount >= obj.maxIter) ...
- | (sqrt(obj.residualCNorm) < obj.tol) ...
- | (sqrt(obj.residualCNorm)./obj.normb < obj.tol));
+ | (obj.residualCNorm < obj.tol) ...
+ | (obj.residualCNorm ./ obj.normb < obj.tol));
end
end
@@ -50,30 +64,23 @@ classdef PcgSolver < IterativeSolver
% update solution
AsearchDirection = obj.A * obj.searchDirection(:,idx);
- if sum(obj.searchDirection(:,idx).*AsearchDirection, 1) < eps
+ dAd = dot(obj.searchDirection(:,idx), AsearchDirection, 1);
+ if dAd < eps
alpha = 1;
else
- alpha = obj.residualCNorm(:,idx) ./ sum(obj.searchDirection(:,idx).*AsearchDirection, 1);
+ alpha = obj.residualCNorm(:,idx).^2 ./ dAd;
end
obj.x(:,idx) = obj.x(:,idx) + alpha .* obj.searchDirection(:,idx);
- % DEBUG:
- % disp(['alpha = ', num2str(alpha)])
-
% update residual
obj.residual(:,idx) = obj.residual(:,idx) - alpha .* AsearchDirection;
- obj.Cresidual(:,idx) = obj.preconditionAction(obj.residual(:,idx));
+ obj.Cresidual(:,idx) = obj.C.apply(obj.residual(:,idx));
residualCNormOld = obj.residualCNorm(:,idx);
- obj.residualCNorm(:,idx) = sum(obj.residual(:,idx).*obj.Cresidual(:,idx), 1);
+ obj.residualCNorm(:,idx) = sqrt(dot(obj.residual(:,idx), obj.Cresidual(:,idx), 1));
% update search direction
- beta = obj.residualCNorm(:,idx) ./ residualCNormOld;
+ beta = (obj.residualCNorm(:,idx) ./ residualCNormOld).^2;
obj.searchDirection(:,idx) = obj.Cresidual(:,idx) + beta .* obj.searchDirection(:,idx);
end
end
-
- %% abstract preconditioning
- methods (Abstract, Access=public)
- preconditionAction(obj, x)
- end
end
\ No newline at end of file
diff --git a/lib/solvers/chooseIterativeSolver.m b/lib/solvers/chooseIterativeSolver.m
deleted file mode 100644
index c0d3bc50f1223bd9d51b41151084f8a2bb9f3aea..0000000000000000000000000000000000000000
--- a/lib/solvers/chooseIterativeSolver.m
+++ /dev/null
@@ -1,66 +0,0 @@
-% chooseIterativeSolver Return suitable solver (instance of subclass of
-% IterativeSolver) and suitable Prolongation P.
-
-function [solver, P] = chooseIterativeSolver(fes, blf, class, variant)
-arguments
- fes FeSpace
- blf BilinearForm
- class (1,1) string {mustBeMember(class, ["cg", "pcg", "multigrid", "direct"])}
- variant (1,1) string {mustBeMember(variant, [ "", ...
- "iChol", "jacobi", ...
- "additiveSchwarzLowOrder", "additiveSchwarzHighOrder" ...
- "lowOrderVcycle", "highOrderVcycle"])} = ""
-end
-
-order = fes.finiteElement.order;
-P = Prolongation.chooseFor(fes);
-
-switch class
- % non-preconditioned CG
- case "cg"
- solver = CgSolver();
-
- % preconditioned CG family
- case "pcg"
- switch variant
- case "iChol"
- solver = ICholPcgSolver(fes);
- case "jacobi"
- if order == 1, solver = JacobiPcgSolver(fes);
- else, solver = BlockJacobiPcgSolver(fes, blf); end
- case {"", "additiveSchwarzLowOrder"}
- if order == 1, solver = LowestOrderAdditiveSchwarzPcg(fes, blf, P);
- else
- solver = AdditiveSchwarzLowOrderPcg(fes, blf);
- end
- case "additiveSchwarzHighOrder"
- solver = AdditiveSchwarzHighOrderPcg(fes, blf, P);
- otherwise
- error('No PCG variant %s!', variant)
- end
-
- % geometric multigrid family
- case "multigrid"
- switch variant
- case {"", "lowOrderVcycle"}
- if order == 1, solver = LowestOrderLocalMg(fes, blf, P);
- else, solver = LocalMgLowOrderVcycle(fes, blf); end
- case "highOrderVcycle"
- if order == 1, solver = LowestOrderLocalMg(fes, blf, P);
- else, solver = LocalMgHighOrderVcycle(fes, blf, P); end
- otherwise
- error('No multigrid variant %s!', variant)
- end
-
- % direct solve (for testing purposes)
- case "direct"
- solver = DirectSolver();
-
- % default
- otherwise
- error('No iterative solver of class %s!', class)
-end
-
-end
-
-
diff --git a/lib/solvers/preconditioners/BlockJacobi.m b/lib/solvers/preconditioners/BlockJacobi.m
new file mode 100644
index 0000000000000000000000000000000000000000..d4bf3ec1cffaab76bcea9d48ed8d4f036e92a34c
--- /dev/null
+++ b/lib/solvers/preconditioners/BlockJacobi.m
@@ -0,0 +1,34 @@
+% BlockJacobi (subclass of Preconditioner) block-diagonal preconditioner.
+%
+% See also: Preconditioner, PcgSolver
+
+classdef BlockJacobi < Preconditioner
+ %% properties
+ properties (Access=private)
+ C
+ fes
+ blf
+ end
+
+ %% methods
+ methods (Access=public)
+ % preconditioner action: patchwise inverse of diagonal
+ function obj = BlockJacobi(fes, blf)
+ arguments
+ fes FeSpace
+ blf BilinearForm
+ end
+
+ obj.fes = fes;
+ obj.blf = blf;
+ end
+
+ function setup(obj, A)
+ obj.C = assemblePatchwise(obj.blf, obj.fes, ':');
+ end
+
+ function Cx = apply(obj, x)
+ Cx = obj.C \ x;
+ end
+ end
+end
\ No newline at end of file
diff --git a/lib/solvers/preconditioners/DiagonalJacobi.m b/lib/solvers/preconditioners/DiagonalJacobi.m
new file mode 100644
index 0000000000000000000000000000000000000000..9ade2832478bd8c204fcb6c3c99b89e5c1132500
--- /dev/null
+++ b/lib/solvers/preconditioners/DiagonalJacobi.m
@@ -0,0 +1,22 @@
+% DiagonalJacobi (subclass of Preconditioner) diagonal preconditioner.
+%
+% See also: Preconditioner, PcgSolver
+
+classdef DiagonalJacobi < Preconditioner
+ %% properties
+ properties (Access=private)
+ C
+ end
+
+ %% methods
+ methods (Access=public)
+ % preconditioner action: inverse of diagonal
+ function setup(obj, A)
+ obj.C = full(diag(A)).^(-1);
+ end
+
+ function Cx = apply(obj, x)
+ Cx = obj.C .* x;
+ end
+ end
+end
\ No newline at end of file
diff --git a/lib/solvers/preconditioners/IncompleteCholesky.m b/lib/solvers/preconditioners/IncompleteCholesky.m
new file mode 100644
index 0000000000000000000000000000000000000000..6a04a2577a04519a146d399cbe4ecbeb9b242583
--- /dev/null
+++ b/lib/solvers/preconditioners/IncompleteCholesky.m
@@ -0,0 +1,22 @@
+% IncompleteCholesky (subclass of Preconditioner) incomplete Cholesky
+% decomposition preconditioner.
+%
+% See also: Preconditioner, PcgSolver
+
+classdef IncompleteCholesky < Preconditioner
+ %% properties
+ properties (Access=private)
+ C
+ end
+
+ %% methods
+ methods (Access=public)
+ function setup(obj, A)
+ obj.C = ichol(A);
+ end
+
+ function Cx = apply(obj, x)
+ Cx = obj.C' \ (obj.C \ x);
+ end
+ end
+end
\ No newline at end of file
diff --git a/lib/solvers/preconditioners/NoPreconditioner.m b/lib/solvers/preconditioners/NoPreconditioner.m
new file mode 100644
index 0000000000000000000000000000000000000000..fa413b5e0f21d730a0c570740656d09f603c7d5d
--- /dev/null
+++ b/lib/solvers/preconditioners/NoPreconditioner.m
@@ -0,0 +1,12 @@
+% NoPreconditioner (subclass of Preconditioner) identity preconditioner (=no-op).
+%
+% See also: Preconditioner, PcgSolver
+
+
+classdef NoPreconditioner < Preconditioner
+ methods (Access=public)
+ function Cx = apply(~, x)
+ Cx = x;
+ 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 0000000000000000000000000000000000000000..c7499bdc6d56acac17c2531fd484ba6ac3c603b9
--- /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/Preconditioner.m b/lib/solvers/preconditioners/Preconditioner.m
new file mode 100644
index 0000000000000000000000000000000000000000..88833d8618a7f99c5c25dc8e83edf8e46f840a16
--- /dev/null
+++ b/lib/solvers/preconditioners/Preconditioner.m
@@ -0,0 +1,23 @@
+% Preconditioner (abstract handle class) Interface for preconditioners.
+%
+% pc.setup(A, ...) hook to set up the preconditioner for the matrix A.
+% Default implementation does nothing.
+%
+% Abstract methods to implement:
+%
+% Cx = pc.apply(x) applies the preconditioner to the vector x.
+%
+% See also: PcgSolver
+
+
+classdef Preconditioner < handle
+ methods (Access=public)
+ function setup(~, ~)
+ % Default implementation does nothing.
+ end
+ end
+
+ methods (Abstract, Access=public)
+ apply(obj, x)
+ end
+end
\ No newline at end of file
diff --git a/lib/solvers/smoothers/JacobiHighOrderCascade.m b/lib/solvers/smoothers/JacobiHighOrderCascade.m
new file mode 100644
index 0000000000000000000000000000000000000000..d58ab6995e87d15155a6a8aac5e58fe7a68779a2
--- /dev/null
+++ b/lib/solvers/smoothers/JacobiHighOrderCascade.m
@@ -0,0 +1,85 @@
+% JacobiHighOrderCascade (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 JacobiHighOrderCascade < MultilevelSmoother
+ properties (Access=protected)
+ P
+ loFes
+ patchwiseA
+ inclusionMatrix
+ intergridMatrix
+ freeDofs
+ freeDofsOld
+ end
+
+ %% methods
+ methods (Access=public)
+ function obj = JacobiHighOrderCascade(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/JacobiLowOrderCascade.m b/lib/solvers/smoothers/JacobiLowOrderCascade.m
new file mode 100644
index 0000000000000000000000000000000000000000..b0b9e5545ee9e79ae3e8420dedf94b822bca81f3
--- /dev/null
+++ b/lib/solvers/smoothers/JacobiLowOrderCascade.m
@@ -0,0 +1,90 @@
+% JacobiLowOrderCascade (subclass of MultilevelSmoother) multilevel
+% Jacobi smoother for higher order finite elements: smooth locally with
+% diagonal of P1 stiffness matrix on changed patches on every level, except
+% for finest level, where global patchwise higher order smoothing is done.
+%
+% Reference: [Innerberger, Miraçi, Praetorius, Streitberger; Optimal
+% computational costs of AFEM with optimal local hp-robust multigrid solver].
+%
+% See also: MultilevelSmoother, OptimalVcycleMultigridSolver
+
+
+classdef JacobiLowOrderCascade < MultilevelSmoother
+ properties (Access=protected)
+ P
+ loFes
+ p1Smoother
+ hoSmoother
+ inclusionMatrix
+ intergridMatrix
+ freeVertices
+ freeVerticesOld
+ end
+
+ %% methods
+ methods (Access=public)
+ function obj = JacobiLowOrderCascade(fes, blf)
+ obj = obj@MultilevelSmoother(fes, blf);
+
+ assert(fes.finiteElement.order > 1, ...
+ 'JacobiLowOrderCascadeSmoother only works for higher order finite elements.')
+
+ mesh = fes.mesh;
+ obj.loFes = FeSpace(mesh, LowestOrderH1Fe(), 'dirichlet', fes.bnd.dirichlet);
+ obj.P = Prolongation.chooseFor(obj.loFes);
+
+ end
+
+ function nonInvertedSmootherMatrix = setup(obj, A)
+ setup@MultilevelSmoother(obj, A);
+
+ L = obj.nLevels;
+ obj.freeVerticesOld = obj.freeVertices;
+ obj.freeVertices = getFreeDofs(obj.loFes);
+
+ p1Matrix = assemble(obj.blf, obj.loFes);
+ p1Matrix = p1Matrix(obj.freeVertices, obj.freeVertices);
+ obj.p1Smoother{L} = full(diag(p1Matrix)).^(-1);
+ 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.p1Smoother{L} = obj.p1Smoother{L}(obj.changedPatches{L});
+ % finest level smoother and transfer operator
+ obj.hoSmoother = assemblePatchwise(obj.blf, obj.fes, ':');
+ obj.inclusionMatrix = SpaceProlongation(obj.loFes, obj.fes) ...
+ .matrix(getFreeDofs(obj.loFes), getFreeDofs(obj.fes));
+ end
+ end
+
+ function Cx = smooth(obj, x, k)
+ if k == obj.nLevels
+ % higher order global patchwise smooting
+ Cx = obj.hoSmoother \ x;
+ else
+ % local p1 patchwise smoothing
+ idx = obj.changedPatches{k};
+ Cx = zeros(size(x));
+ Cx(idx,:) = obj.p1Smoother{k} .* x(idx,:);
+ 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
+ Px = obj.inclusionMatrix' * Px;
+ end
+ end
+
+ function Px = restrict(obj, x, k)
+ if k == obj.nLevels
+ x = obj.inclusionMatrix * x;
+ end
+ Px = obj.intergridMatrix{k}' * x;
+ end
+ end
+end
\ No newline at end of file
diff --git a/lib/solvers/smoothers/MultilevelSmoother.m b/lib/solvers/smoothers/MultilevelSmoother.m
new file mode 100644
index 0000000000000000000000000000000000000000..ca7271db059c51bb0e23f923b20bd0e54110dc7c
--- /dev/null
+++ b/lib/solvers/smoothers/MultilevelSmoother.m
@@ -0,0 +1,84 @@
+% 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.
+%
+% 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.
+%
+% 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
+
+ properties (Access=protected)
+ fes
+ blf
+ changedPatches
+ end
+
+ properties (Access=private)
+ listenerHandle
+ end
+
+ %% methods
+ methods (Access=public)
+ 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;
+
+ mesh = fes.mesh;
+ obj.listenerHandle = mesh.listener('IsAboutToRefine', @obj.getChangedPatches);
+ end
+
+ function nonInvertedSmootherMatrix = setup(obj, ~)
+ obj.nLevels = obj.nLevels + 1;
+
+ % should be implemented by subclasses
+ nonInvertedSmootherMatrix = [];
+ 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
+
+ %% 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/P1JacobiCascade.m b/lib/solvers/smoothers/P1JacobiCascade.m
new file mode 100644
index 0000000000000000000000000000000000000000..93d15999bc92e42a7ebcb05955c717d93cdd4844
--- /dev/null
+++ b/lib/solvers/smoothers/P1JacobiCascade.m
@@ -0,0 +1,69 @@
+% P1JacobiCascade (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 P1JacobiCascade < MultilevelSmoother
+ properties (Access=protected)
+ P
+ inverseDiagonal
+ intergridMatrix
+ freeVertices
+ freeVerticesOld
+ end
+
+ %% event data
+ properties (Access=private)
+ listenerHandle
+ end
+
+ %% methods
+ methods (Access=public)
+ function obj = P1JacobiCascade(fes, blf, P)
+ arguments
+ fes
+ blf
+ P Prolongation
+ end
+
+ assert(fes.finiteElement.order == 1, ...
+ 'P1JacobiSmoother only works for lowest order finite elements.')
+
+ obj = obj@MultilevelSmoother(fes, blf);
+ obj.P = P;
+ end
+
+ 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
+ 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
+end
\ No newline at end of file
diff --git a/runAdditiveSchwarzPreconditioner.m b/runAdditiveSchwarzPreconditioner.m
index 3e6d606bd1631a6eaad60923c6a2d12516e26ae9..971ffd8d1f8979c73c28c1216fafa76b00187343 100644
--- a/runAdditiveSchwarzPreconditioner.m
+++ b/runAdditiveSchwarzPreconditioner.m
@@ -32,10 +32,10 @@ function leveldata = runAdditiveSchwarzPreconditioner(pMax)
lf.f = Constant(mesh, 1);
% Choose iterative solver
- solver = chooseIterativeSolver(fes, blf, "pcg", "additiveSchwarzLowOrder");
- % solver = chooseIterativeSolver(fes, blf, "pcg", "additiveSchwarzHighOrder");
- % solver = chooseIterativeSolver(fes, blf, "pcg", "iChol");
- % solver = chooseIterativeSolver(fes, blf, "multigrid", "lowOrderVcycle");
+ solver = IterativeSolver.choose(fes, blf, "pcg", "additiveSchwarzLowOrder");
+ % solver = IterativeSolver.choose(fes, blf, "pcg", "additiveSchwarzHighOrder");
+ % solver = IterativeSolver.choose(fes, blf, "pcg", "iChol");
+ % solver = IterativeSolver.choose(fes, blf, "multigrid", "lowOrderVcycle");
for k = 1:nLevels-1
% Assemble matrices on intermediate meshes
diff --git a/runIterativeSolver.m b/runIterativeSolver.m
index 350c6b1f2801128319ed64ee03192f39facac566..ce2bcf8b6d8cd180066937d7c5843e07e05907f6 100644
--- a/runIterativeSolver.m
+++ b/runIterativeSolver.m
@@ -28,10 +28,10 @@ function leveldata = runIterativeSolver(maxNiter)
lf.f = Constant(mesh, 1);
% Choose iterative solver
- solver = chooseIterativeSolver(fes, blf, "pcg", "additiveSchwarzLowOrder");
- % solver = chooseIterativeSolver(fes, blf, "pcg", "additiveSchwarzHighOrder");
- % solver = chooseIterativeSolver(fes, blf, "pcg", "iChol");
- % solver = chooseIterativeSolver(fes, blf, "multigrid", "lowOrderVcycle");
+ solver = IterativeSolver.choose(fes, blf, "pcg", "additiveSchwarzLowOrder");
+ % solver = IterativeSolver.choose(fes, blf, "pcg", "additiveSchwarzHighOrder");
+ % solver = IterativeSolver.choose(fes, blf, "pcg", "iChol");
+ % solver = IterativeSolver.choose(fes, blf, "multigrid", "lowOrderVcycle");
% Initialize LevelData
leveldata = LevelData('results');
diff --git a/tests/TestIterativeSolver.m b/tests/TestIterativeSolver.m
index 043d4e7fb5ea8ed03f1e85ecd84491f71a9575fc..6a32e6af2f72e24d59a80be2ceb5d94690187682 100644
--- a/tests/TestIterativeSolver.m
+++ b/tests/TestIterativeSolver.m
@@ -6,16 +6,33 @@ properties (TestParameter)
'direct', ["direct", ""], ...
'jacobiPCG', ["pcg", "jacobi"], ...
'iChol', ["pcg", "iChol"], ...
- 'additiveSchwarzPCG', ["pcg", "additiveSchwarzLowOrder"], ...
+ 'loAdditiveSchwarzPCG', ["pcg", "additiveSchwarzLowOrder"], ...
+ 'hoAdditiveSchwarzPCG', ["pcg", "additiveSchwarzHighOrder"], ...
'lowMG', ["multigrid", "lowOrderVcycle"], ...
'highMG', ["multigrid", "highOrderVcycle"]);
end
methods (Test)
- function oneStepDecreasesNorm(testCase, p, variant)
+ function firstStepDecreasesNorm(testCase, p, variant)
+ [~, fes, blf, lf] = setupProblem(testCase, p);
+ s = variant(1); v = variant(2);
+ solver = IterativeSolver.choose(fes, blf, s, v);
+
+ [xstar, A, F] = assembleData(testCase, blf, lf, fes);
+ solver.setupSystemMatrix(A);
+ solver.setupRhs(F);
+
+ normBefore = sqrt((xstar - solver.x)' * A * (xstar - solver.x));
+ solver.step();
+ normAfterwards = sqrt((xstar - solver.x)' * A * (xstar - solver.x));
+
+ testCase.verifyGreaterThan(normBefore, normAfterwards);
+ end
+
+ function laterStepDecreasesNorm(testCase, p, variant)
[mesh, fes, blf, lf] = setupProblem(testCase, p);
s = variant(1); v = variant(2);
- solver = chooseIterativeSolver(fes, blf, s, v);
+ solver = IterativeSolver.choose(fes, blf, s, v);
for k = 1:3
[xstar, A, F] = assembleData(testCase, blf, lf, fes);
@@ -30,11 +47,11 @@ methods (Test)
testCase.verifyGreaterThan(normBefore, normAfterwards);
end
-
+
function solverIsLinear(testCase, p, variant)
[mesh, fes, blf, lf] = setupProblem(testCase, p);
s = variant(1); v = variant(2);
- solver = chooseIterativeSolver(fes, blf, s, v);
+ solver = IterativeSolver.choose(fes, blf, s, v);
for k = 1:3
[xstar, A, F] = assembleData(testCase, blf, lf, fes);
@@ -44,15 +61,9 @@ methods (Test)
solver.setupRhs([F, -pi*F], 0*[F, F]);
solver.solve();
-% for i = 1:60
-% solver.step();
-% end
- %xmat = cgs(A, F, 1e-16, 100);
-
- testCase.verifyEqual(-pi*solver.x(:,1), solver.x(:,2), 'RelTol', 2e-5);
- testCase.verifyEqual(solver.x(:,1), xstar, 'RelTol', 2e-5);
- %testCase.verifyEqual(solver.x(:, 1), xmat(:, 1), 'RelTol', 2e-5)
- %testCase.verifyEqual(xmat(:,1), xstar, 'RelTol', 2e-5);
+
+ testCase.verifyEqual(-pi*solver.x(:,1), solver.x(:,2), 'RelTol', 2e-5);
+ testCase.verifyEqual(solver.x(:,1), xstar, 'RelTol', 2e-5);
end
end