Draft CR solver with quasi-interpolation as restriction

lib/solvers/@IterativeSolver/choose.m

@@ -21,7 +21,8 @@ arguments
         "iChol", "jacobi", ...
         "additiveSchwarzLowOrder", "additiveSchwarzHighOrder" ...
         "lowOrderVcycle", "highOrderVcycle", ...
-        "multiplicativeCR", "vcycleCR", "additiveCR"])} = ""
+        "multiplicativeCR", "vcycleCR", "additiveCR", ...
+        "multiplicativeCRQI", "vcycleCRQI"])} = ""
 end
 
 order = fes.finiteElement.order;
@@ -85,10 +86,18 @@ switch method
                 Av = LoMeshAveraging(fes);  % TODO: use proper choose function
                 smoother = CRAveragingJacobiCascade(fes, blf, Av, P);
                 solver = LocalMultiplicativeMultigridCRSolver(smoother);
+            case "multiplicativeCRQI"
+                P = LoMeshQuasiInterpolation(fes);
+                smoother = CRQuasiInterpolationJacobiCascade(fes, blf, P);
+                solver = LocalMultiplicativeMultigridCRQISolver(smoother);
             case "vcycleCR"
                 Av = LoMeshAveraging(fes);  % TODO: use proper choose function
                 smoother = CRAveragingJacobiCascade(fes, blf, Av, P);
                 solver = LocalMultigridVcycleCRSolver(smoother);
+            case "vcycleCRQI"
+                P = LoMeshQuasiInterpolation(fes);
+                smoother = CRQuasiInterpolationJacobiCascade(fes, blf, P);
+                solver = OptimalVcycleMultigridSolver(smoother);
             otherwise
                 error('No multigrid variant %s!', variant)
         end

lib/solvers/LocalMultiplicativeMultigridCRQISolver.m

+% LocalMultiplicativeMultigridCRSolver (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, MultigridSolver, MultiLevelSmoother, DiagonalJacobi,
+%   LowOrderBlockJacobi, HighOrderBlockJacobi
+
+classdef LocalMultiplicativeMultigridCRQISolver < MultigridSolver
+    %% methods
+    methods
+        function obj = LocalMultiplicativeMultigridCRQISolver(smoother)
+            obj = obj@MultigridSolver(smoother);
+        end
+    end
+        
+    methods (Access=protected)
+        function [Cx, algEstimator] = cycle(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
+
+            % sequential multiplicative V-cycle update
+            algEstimator = 0;
+            res = x;
+            v = obj.x;
+            for ell = 1:L
+                % Apply V-cycle down to level ell
+                xUpdate = obj.Vcycle(res, ell);
+                % A priori step size
+                if ell == L
+                    rho = 1;
+                else
+                    rho = 0.8;
+                end
+                v = v + rho * xUpdate;
+                % Compute optimal step size
+                % updatedResidual = x - obj.A * Cx;
+                % rho = obj.smoother.smooth(updatedResidual, ell+1);
+                % [etaUpdate2, xUpdate] = MultigridSolver.computeOptimalUpdate(obj.A, updatedResidual, rho);
+                % Cx = Cx + xUpdate;
+                % Update estimator
+                % algEstimator = algEstimator + etaUpdate2;
+                % Update residual
+                res = res - obj.A * rho * xUpdate;
+            end
+            Cx = v;
+        end
+
+        function Cx = Vcycle(obj, x, ell) 
+            % one step of the Vcycle
+            
+            L = obj.nLevels;
+
+            if ell == L
+                Cx = obj.smoother.smooth(x, ell);
+                return
+            end
+
+            % descending cascade: no smoothing
+            res = x;
+            for k = L:-1:ell+1
+                res = obj.smoother.restrict(res, k);
+            end
+
+            % smoothing step on level ell
+            if ell == 1
+                Cx = obj.projectedMatrix{1} \ res;
+            else
+                Cx = obj.smoother.smooth(res, ell);
+            end
+
+            % ascending cascade: no smoothing
+            for k = ell+1:L
+                Cx = obj.smoother.prolongate(Cx, k);
+            end
+        end
+    end
+end
\ No newline at end of file

lib/solvers/smoothers/CRQuasiInterpolationJacobiCascade.m

+% CRQuasiInterpolationJacobiCascade (subclass of MultilevelSmoother) multilevel Jacobi smoother
+%   for lowest order nonconforming Crouzeix-Raviart finite elements: smooth 
+%   with diagonal of stiffness matrix on changed patches on every level.
+%
+% See also: MultilevelSmoother, OptimalVcycleMultigridSolver
+
+
+classdef CRQuasiInterpolationJacobiCascade < MultilevelSmoother
+    properties (Access=protected)
+        P
+        patchwiseA
+        intergridMatrix
+        freeDofs
+        freeDofsOld
+        confFes
+    end
+    
+    %% event data
+    properties (Access=private)
+        listenerHandle
+    end
+
+    %% methods
+    methods (Access=public)
+        function obj = CRQuasiInterpolationJacobiCascade(fes, blf, P)
+            arguments
+                fes
+                blf
+                P Prolongation
+            end
+
+            assert(isa(fes.finiteElement, 'LowestOrderCRFe'), ...
+                ['CRQuasiInterpolationJacobiCascade only works for lowest order nonconforming', ... 
+                'Crouzeix-Raviart finite elements.'])
+            assert(fes.finiteElement.order == 1, ...
+                'CRQuasiInterpolationJacobiCascade only works for lowest order finite elements.')
+            
+            obj = obj@MultilevelSmoother(fes, blf);
+            obj.P = P;
+            obj.confFes = FeSpace(fes.mesh, LowestOrderH1Fe);
+        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.confFes);
+            nonInvertedSmootherMatrix = A;
+            
+            if L > 1
+                obj.intergridMatrix{L} = obj.P.matrix(obj.freeDofs, obj.freeDofsOld);
+                obj.changedPatches{L} = find(obj.changedPatches{L}(freeVertices));
+                obj.patchwiseA{L} = assemblePatchwise(obj.blf, obj.fes, obj.changedPatches{L});
+            end
+        end
+
+        function Cx = smooth(obj, x, k, nSteps)
+            if nargin < 4
+                nSteps = 1;
+            end
+            Cx = obj.patchwiseA{k} \ x;
+            for j = 2:nSteps
+                Cx = obj.patchwiseA{k} \ Cx;
+            end
+        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
+
+        % Debugging
+        function plot(obj)
+            plot(obj.fes.mesh)
+        end
+    end
+end

tests/TestCRIterativeSolver.m

@@ -4,8 +4,10 @@ properties (TestParameter)
     variant = struct('CG', ["cg", ""],...
         'direct', ["direct", ""], ...
         'additiveCR', ["pcg", "additiveCR"], ...
-        ... 'MG', ["multigrid", "multiplicativeCR"]), ...
-        'MGvcycle', ["multigrid", "vcycleCR"]);
+        'MG', ["multigrid", "multiplicativeCR"], ...
+        'MGQI', ["multigrid", "multiplicativeCRQI"], ...
+        'MGvcycle', ["multigrid", "vcycleCR"], ...
+        'MGvcycleQI', ["multigrid", "vcycleCRQI"]);
 end
 
 methods (Test)