From 1bca69aeb1ca22c745b617bc2e5304ed527ef6bd Mon Sep 17 00:00:00 2001
From: Julian Streitberger <julian.streitberger@tuwien.ac.at>
Date: Mon, 20 Jan 2025 14:43:53 +0100
Subject: [PATCH] Clean up
---
examples/afemWithGmres.m | 46 +++++---------------------------------
lib/solvers/PgmresSolver.m | 4 ++--
2 files changed, 7 insertions(+), 43 deletions(-)
diff --git a/examples/afemWithGmres.m b/examples/afemWithGmres.m
index 056452f..0a19331 100644
--- a/examples/afemWithGmres.m
+++ b/examples/afemWithGmres.m
@@ -4,10 +4,10 @@
% ******************************************************************************
%% parameters
-nDofsMax = 2e5;
+nDofsMax = 5e5;
theta = 0.5;
pmax = 4;
-lamalg = 0.5;
+lamalg = 0.1;
%% initialization for every polynomial degree
[nElem, nDofs, errEst, iterations, time, directSolveTime] = deal(zeros(pmax, 1000));
@@ -43,14 +43,9 @@ for p = 4:pmax
lf.f = Constant(mesh, 1);
%% set up solver and operator for nested iteration
- % 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] = IterativeSolver.choose(fes, blf, "pcg", "iChol");
- %[solver, P] = IterativeSolver.choose(fes, blf, "pcg", "additiveSchwarzHighOrder");
[solver, P] = IterativeSolver.choose(fes, symmetricBlf, "gmres", "additiveSchwarzLowOrder");
solver.tol = 1e-6;
- solver.maxIter = 100;
+ solver.maxIter = 10;
%% adaptive loop
ell = 0;
@@ -72,30 +67,7 @@ for p = 4:pmax
solver.setupSystemMatrix(A, principalPart);
solver.setupRhs(rhs(freeDofs), u.data(freeDofs)');
- % exact solution as reference
- % x_ex = A \ rhs(freeDofs);
- % uEllExact.setFreeData(x_ex);
- % directSolveTime(p, ell) = toc;
-
- % %% iterative solve
- % % solver.solve() implements the algebraic solver loop and
- % % terminates automatically once the wanted accuracy is reached
- % % Another option is to manually implement a loop and implement one
- % % solver step with solver.step() and stop by setting a stopping
- % % criterion
- % solverIsConverged = false;
- % while ~solverIsConverged
- % solver.solve();
- % u.setFreeData(solver.x)
- % eta2 = estimate(blf, lf, u);
- % estimator = sqrt(sum(eta2));
- % solver.tol = lamalg*estimator;
- % solverIsConverged = solver.isConverged;
- % end
- %
-
solverIsConverged = false;
- solver.maxIter = 10;
k = 0;
while ~solverIsConverged
k = k + 1;
@@ -103,7 +75,8 @@ for p = 4:pmax
u.setFreeData(solver.x)
eta2 = estimate(blf, lf, u);
estimator = sqrt(sum(eta2));
- solverIsConverged = (solver.residualCNorm < lamalg*estimator);
+ solver.tol = lamalg*estimator;
+ solverIsConverged = solver.applyStoppingCriterion();
% store level-oriented data
leveldata.append('ell', int32(ell), 'k', int32(k), ...
'nDofs', int32(length(freeDofs)), ...
@@ -121,15 +94,6 @@ for p = 4:pmax
end
end
-
- %% estimate error and store data
- % nDofs(p,ell) = getDofs(fes).nDofs;
- % errEst(p,ell) = estimator;
- % iterations(p, ell) = solver.iterationCount();
- % nElem(p,ell) = mesh.nElements;
- % printLogMessage('number of dofs: %d, estimator: %.2e after %d iterations', ...
- % nDofs(p,ell), errEst(p,ell), solver.iterationCount);
-
%% stoping criterion
meshSufficientlyFine = (length(freeDofs) > nDofsMax);
diff --git a/lib/solvers/PgmresSolver.m b/lib/solvers/PgmresSolver.m
index 72694b8..d519fb2 100644
--- a/lib/solvers/PgmresSolver.m
+++ b/lib/solvers/PgmresSolver.m
@@ -111,13 +111,13 @@ classdef PgmresSolver < IterativeSolver
obj.H(i, k) = temp;
end
- % Neue Givens-Rotation berechnen
+ % Apply Givens rotation
beta = sqrt( obj.H(k, k)^2 + obj.H(k + 1, k)^2);
obj.cosine(k) = obj.H(k, k) / beta;
obj.sine(k) = obj.H(k + 1, k) / beta;
obj.H(k, k) = beta;
- % Residualvektor aktualisieren
+ % Update residual vector
obj.g(k + 1) = -obj.sine(k) * obj.g(k);
obj.g(k) = obj.cosine(k) * obj.g(k);
--
GitLab