From eb6acf68adee091ae3d0148858c89737d7eccc2d Mon Sep 17 00:00:00 2001
From: Michael Innerberger <michael.innerberger@asc.tuwien.ac.at>
Date: Wed, 19 Jul 2023 07:21:54 -0400
Subject: [PATCH] Port 1->1->p cacade

---
 lib/solvers/LocalMgLowOrderVcycle.m           | 179 ------------------
 lib/solvers/chooseIterativeSolver.m           |   2 +-
 .../smoothers/JacobiLowOrderCascadeSmoother.m |  88 +++++++++
 3 files changed, 89 insertions(+), 180 deletions(-)
 delete mode 100644 lib/solvers/LocalMgLowOrderVcycle.m
 create mode 100644 lib/solvers/smoothers/JacobiLowOrderCascadeSmoother.m

diff --git a/lib/solvers/LocalMgLowOrderVcycle.m b/lib/solvers/LocalMgLowOrderVcycle.m
deleted file mode 100644
index f8920c0..0000000
--- 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/chooseIterativeSolver.m b/lib/solvers/chooseIterativeSolver.m
index fa30db2..9bc2cf9 100644
--- a/lib/solvers/chooseIterativeSolver.m
+++ b/lib/solvers/chooseIterativeSolver.m
@@ -55,7 +55,7 @@ switch class
                 if order == 1
                     solver = OptimalVcycleMultigridSolver(P1JacobiSmoother(fes, blf, P));
                 else
-                    solver = LocalMgLowOrderVcycle(fes, blf);
+                    solver = OptimalVcycleMultigridSolver(JacobiLowOrderCascadeSmoother(fes, blf));
                 end
             case "highOrderVcycle"
                 if order == 1
diff --git a/lib/solvers/smoothers/JacobiLowOrderCascadeSmoother.m b/lib/solvers/smoothers/JacobiLowOrderCascadeSmoother.m
new file mode 100644
index 0000000..ad5dd1e
--- /dev/null
+++ b/lib/solvers/smoothers/JacobiLowOrderCascadeSmoother.m
@@ -0,0 +1,88 @@
+% JacobiLowOrderCascadeSmoother (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 JacobiLowOrderCascadeSmoother < MultilevelSmoother
+    properties (Access=protected)
+        P
+        loFes
+        p1Smoother
+        hoSmoother
+        inclusionMatrix
+        intergridMatrix
+        freeVertices
+        freeVerticesOld
+    end
+
+    %% methods
+    methods (Access=public)
+        function obj = JacobiLowOrderCascadeSmoother(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
+                % intermediate level smoothers and transfer operators
+                obj.intergridMatrix{L} = obj.P.matrix(obj.freeVertices, obj.freeVerticesOld);
+                obj.changedPatches{L} = find(obj.changedPatches{L}(obj.freeVertices));
+                obj.p1Smoother{L} = obj.p1Smoother{L}(obj.changedPatches{L});
+                % finest level smoother and transfer operator
+                obj.hoSmoother = assemblePatchwise(obj.blf, obj.fes, ':');
+                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
+
+        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
-- 
GitLab