*** Wartungsfenster jeden ersten Mittwoch vormittag im Monat ***

Skip to content
Snippets Groups Projects
Commit eb6acf68 authored by Innerberger, Michael's avatar Innerberger, Michael
Browse files

Port 1->1->p cacade

parent 2c143659
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
lib/solvers/LocalMgLowOrderVcycle.m deleted 100644 → 0
+ 0
− 179
View file @ 2c143659
% 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
lib/solvers/chooseIterativeSolver.m
+ 1
− 1
View file @ eb6acf68
......@@ -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
......
lib/solvers/smoothers/JacobiLowOrderCascadeSmoother.m 0 → 100644
+ 88
− 0
View file @ eb6acf68
% 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
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment

[ 18.Apr 0h : 3,749 users | 915 groups | 5,096 projects ]