Split additive Schwarz PCG into p=1 and higher order

There is still some error in the higher order version! The solver passes
the tests only if there is no smoothing on intermediate levels!

lib/solvers/AdditiveSchwartzPcg.m

@@ -12,7 +12,7 @@ classdef AdditiveSchwartzPcg < PcgSolver
         blf
         hoFes
         loFes
-        p1Matrix
+        p1Acoarse
         p1Smoother
         P
         intergridMatrix
@@ -29,21 +29,27 @@ classdef AdditiveSchwartzPcg < PcgSolver
     
     %% methods
     methods (Access=public)
-        function obj = AdditiveSchwartzPcg(fes, blf)
+        function obj = AdditiveSchwartzPcg(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.highestOrderIsOne = (fes.finiteElement.order == 1);
             obj.nLevels = 0;
             
             mesh = fes.mesh;
-            obj.loFes = FeSpace(mesh, LowestOrderH1Fe(), 'dirichlet', ':');
+            obj.loFes = FeSpace(mesh, LowestOrderH1Fe, 'dirichlet', ':');
             obj.hoFes = fes;
-            obj.P = LoFeProlongation(obj.loFes);
-            obj.blf = blf; % TODO: check coefficients of blf for compatibility with theoretical results!
+            obj.P = P;
+            obj.blf = blf;
             
             obj.listenerHandle = mesh.listener('IsAboutToRefine', @obj.getChangedPatches);
         end
@@ -55,21 +61,16 @@ classdef AdditiveSchwartzPcg < PcgSolver
             obj.freeDofsOld = obj.freeDofs;
             obj.freeDofs = getFreeDofs(obj.loFes);
             
-            if obj.highestOrderIsOne
-                obj.p1Matrix{L} = A;
-            else
-                obj.p1Matrix{L} = assemble(obj.blf, obj.loFes);
-                obj.p1Matrix{L} = obj.p1Matrix{L}(obj.freeDofs, obj.freeDofs);
-            end
-            obj.p1Smoother{L} = full(diag(obj.p1Matrix{L})).^(-1);
+            p1Matrix = assemble(obj.blf, obj.loFes);
+            p1Matrix = p1Matrix(obj.freeDofs, obj.freeDofs);
             
-            if L >= 2
+            if L == 1
+                obj.p1Acoarse = p1Matrix;
+            else
+                obj.p1Smoother{L} = full(diag(p1Matrix)).^(-1);
                 obj.intergridMatrix{L} = obj.P.matrix(obj.freeDofs, obj.freeDofsOld);
                 obj.changedPatches{L} = find(obj.changedPatches{L}(obj.freeDofs));
-
-                if ~obj.highestOrderIsOne
-                    obj.patchwiseA = assemblePatchwise(obj.blf, obj.hoFes, ':');
-                end
+                obj.patchwiseA = assemblePatchwise(obj.blf, obj.hoFes, ':');
             end
             
             setupSystemMatrix@PcgSolver(obj, A);
@@ -92,14 +93,8 @@ classdef AdditiveSchwartzPcg < PcgSolver
             end
 
             % descending cascade
-            if obj.highestOrderIsOne
-                rho{L} = p1LocalSmoothing(obj, L, res);
-            else
-                %finest level high-order patch-problems ONLY for p > 1
-                rho{L} = hoGlobalSmoothing(obj, res);
-            end
-
-            residual = projectFromPto1(obj, res);
+            rho{L} = hoGlobalSmoothing(obj, res);
+            residual = interpolateFreeData(res, obj.hoFes, obj.loFes);
             for k = L-1:-1:2
                 residual = obj.intergridMatrix{k+1}'*residual;
                 rho{k} = p1LocalSmoothing(obj, k, residual);
@@ -107,7 +102,7 @@ classdef AdditiveSchwartzPcg < PcgSolver
             
             % exact solve on coarsest level
             residual = obj.intergridMatrix{2}'*residual;
-            sigma = obj.p1Matrix{1} \ residual;
+            sigma = obj.p1Acoarse \ residual;
             
             % ascending cascade
             for k = 2:L-1
@@ -115,7 +110,7 @@ classdef AdditiveSchwartzPcg < PcgSolver
                 sigma = sigma + rho{k};
             end
             sigma = obj.intergridMatrix{L}*sigma;
-            sigma = projectFrom1toP(obj, sigma);
+            sigma = interpolateFreeData(sigma, obj.loFes, obj.hoFes);
             sigma = sigma + rho{L};
            
             Cx = sigma;
@@ -145,22 +140,6 @@ classdef AdditiveSchwartzPcg < PcgSolver
         function rho = hoGlobalSmoothing(obj, res)
             rho = obj.patchwiseA \ res;
         end
-        
-        function y = projectFrom1toP(obj, x)
-            if obj.highestOrderIsOne
-                y = x;
-            else
-                y = interpolateFreeData(x, obj.loFes, obj.hoFes);
-            end
-        end
-        
-        function y = projectFromPto1(obj, x)
-            if obj.highestOrderIsOne
-                y = x;
-            else
-                y = interpolateFreeData(x, obj.hoFes, obj.loFes);
-            end
-        end
     end
 end
 

lib/solvers/LowestOrderAdditiveSchwartzPcg.m

+% LowestOrderAdditiveSchwartzPcg (subclass of PcgSolver) Solves linear
+%   equations iteratively using the CG method with optimal multilevel
+%   additive Schwartz preconditioner for lowest order finite elements.
+
+classdef LowestOrderAdditiveSchwartzPcg < 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 = LowestOrderAdditiveSchwartzPcg(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

lib/solvers/LowestOrderLocalMg.m

@@ -34,10 +34,11 @@ classdef LowestOrderLocalMg < MGSolver
 
     %% methods
     methods (Access=public)
-        function obj = LowestOrderLocalMg(fes, blf)
+        function obj = LowestOrderLocalMg(fes, blf, P)
             arguments
                 fes FeSpace
                 blf BilinearForm
+                P Prolongation
             end
             obj = obj@MGSolver();
             
@@ -50,7 +51,7 @@ classdef LowestOrderLocalMg < MGSolver
             
             mesh = fes.mesh;
             obj.fes = fes;
-            obj.P = Prolongation.chooseFor(obj.fes);
+            obj.P = P;
             obj.blf = blf;
             
             obj.listenerHandle = mesh.listener('IsAboutToRefine', @obj.getChangedPatches);

lib/solvers/chooseIterativeSolver.m

@@ -26,7 +26,8 @@ switch class
                 if order == 1, solver = JacobiPcgSolver(fes);
                 else, solver = BlockJacobiPcgSolver(fes, blf); end
             case {"", "additiveSchwarz"}
-                solver = AdditiveSchwartzPcg(fes, blf);
+                if order == 1, solver = LowestOrderAdditiveSchwartzPcg(fes, blf, P);
+                else, solver = AdditiveSchwartzPcg(fes, blf, P); end
             otherwise
                 error('No PCG variant %s!', variant)
         end
@@ -36,10 +37,10 @@ switch class
     case "multigrid"
         switch variant
             case {"", "lowOrderVcycle"}
-                if order == 1, solver = LowestOrderLocalMg(fes, blf);
+                if order == 1, solver = LowestOrderLocalMg(fes, blf, P);
                 else, solver = LocalMgLowOrderVcycle(fes, blf); end
             case "highOrderVcycle"
-                if order == 1, solver = LowestOrderLocalMg(fes, blf);
+                if order == 1, solver = LowestOrderLocalMg(fes, blf, P);
                 else, solver = LocalMgHighOrderVcycle(fes, blf, P); end
             otherwise
                 error('No multigrid variant %s!', variant)