From e7bec66b54ac1df2960a9507e2eade85ef8b990f Mon Sep 17 00:00:00 2001
From: Michael Innerberger <michael.innerberger@asc.tuwien.ac.at>
Date: Mon, 4 Jul 2022 15:56:12 +0200
Subject: [PATCH] Add modified Zarantonello iteration (allows for larger delta)
---
experiments/iterativeModifiedZarantonello.m | 159 ++++++++++++++++++++
1 file changed, 159 insertions(+)
create mode 100644 experiments/iterativeModifiedZarantonello.m
diff --git a/experiments/iterativeModifiedZarantonello.m b/experiments/iterativeModifiedZarantonello.m
new file mode 100644
index 0000000..6ab3536
--- /dev/null
+++ b/experiments/iterativeModifiedZarantonello.m
@@ -0,0 +1,159 @@
+% ******************************************************************************
+% Adaptive FEM algorithm for a non-symmetric/indefinite problem using a
+% modified Zarantonello iteration (where the energy scalar product is used
+% instead of the Laplace operator).
+% ******************************************************************************
+
+%% paramters
+nEmax = 1e4;
+theta = 0.5;
+coeff = 150;
+lambda = 0.1;
+delta = 0.03;
+
+%% initialization
+mesh = Mesh.loadFromGeometry('Lshape');
+mesh.refineUniform(1);
+fes = FeSpace(mesh, LowestOrderH1Fe());
+u = FeFunction(fes);
+u.setData(0);
+uref = FeFunction(fes);
+uold = FeFunction(fes);
+update = FeFunction(fes);
+
+%% problem data
+blf = BilinearForm(fes);
+blf.a = Constant(mesh, 5);
+
+blfZ = BilinearForm(fes);
+blfZ.a = Constant(mesh, 1);
+
+[lf, perturbation] = setUpConvectionProblem(fes, uold, coeff);
+% [lf, perturbation] = setUpHelmholtzProblem(fes, uold, coeff);
+
+lfRef = LinearForm(fes);
+lfRef.f = Constant(mesh, 1);
+
+P = LoFeProlongation(fes);
+
+%% adaptive loop
+ell = 0;
+[nElements, estimator, exactDiff, nIter] = deal(zeros(1000,1));
+meshSufficientlyFine = false;
+while ~meshSufficientlyFine
+ ell = ell + 1;
+ freeDofs = getFreeDofs(fes);
+ A = assemble(blf);
+ A = A(freeDofs,freeDofs);
+ Z = assemble(blfZ);
+ Z = Z(freeDofs,freeDofs);
+
+ K = assemble(perturbation);
+ Aref = A + K(freeDofs,freeDofs);
+ Fref = assemble(lfRef);
+ uref.setFreeData(Aref \ Fref(freeDofs));
+
+ % inner iteration: sufficiently many kacanov iterations
+ n = 0;
+ solutionAccurateEnough = false;
+ while ~solutionAccurateEnough
+ n = n + 1;
+ uold.setData(u.data);
+ F = assemble(lf);
+ res = A*u.data(freeDofs)' - F(freeDofs);
+ update.setFreeData(A \ res);
+ u.setData(u.data - delta*update.data);
+ eta2 = estimate(blf, lf, u);
+
+ v = u.data(freeDofs) - uold.data(freeDofs);
+ iterationDiff = sqrt(v * A * v');
+
+ solutionAccurateEnough = (iterationDiff < delta*lambda*sqrt(sum(eta2)));
+ end
+
+ nElements(ell) = mesh.nElements;
+ estimator(ell) = sqrt(sum(eta2));
+ v = u.data(freeDofs) - uref.data(freeDofs);
+ exactDiff(ell) = sqrt(v * A * v');
+ nIter(ell) = n;
+ printLogMessage('number of elements: %d, estimator: %.2e after %d iterations', ...
+ mesh.nElements, estimator(ell), n);
+
+ meshSufficientlyFine = (mesh.nElements > nEmax);
+
+ if ~meshSufficientlyFine
+ marked = markDoerflerBinning(eta2, theta);
+ mesh.refineLocally(addLargestElement(marked, mesh), 'NVB');
+ u.setData(prolongate(P, u));
+ end
+end
+
+%% plot solution and estimator
+if mesh.nElements < 2e5
+ plot(u);
+end
+
+figure()
+idx = 1:ell;
+yyaxis left
+ylabel('error / estimator')
+loglog(nElements(idx), estimator(idx), '-o', 'LineWidth', 2)
+hold on
+loglog(nElements(idx), exactDiff(idx), '-x', 'LineWidth', 2)
+shift = estimator(1)*nElements(1);
+loglog(nElements(idx), shift*nElements(idx).^(-1/2), '--', 'LineWidth', 2, 'Color', 'k')
+yyaxis right
+semilogx(nElements(idx), nIter(idx), '-x', 'LineWidth', 2)
+ylim([0, max(nIter)+1])
+legend('\eta_H', '|u_H - u^{k}_H|_{H^1}', '#T_H^{-1/2}', 'iteration count')
+ylabel('number of iterations')
+xlabel('#T_H')
+title('Kacanov iteration for non-symmetric / indefinite problem')
+
+
+%% local function for residual a posteriori error estimation
+% \eta(T)^2 = h_T^2 * || f ||_{L^2(T)}^2 + h_T * || [A*Du] ||_{L^2(E)}^2
+function indicators = estimate(blf, lf, u)
+ mesh = blf.fes.mesh;
+
+ f = CompositeFunction(@(f) f.^2, lf.f);
+ qrTri = QuadratureRule.ofOrder(2);
+ volumeRes = integrateElement(f, qrTri);
+
+ qrEdge = QuadratureRule.ofOrder(1, '1D');
+ f = CompositeFunction(@(A,Du) A.* Du, blf.a, Gradient(u));
+ edgeRes = integrateNormalJump(f, qrEdge, @(j) j.^2, {}, ':');
+ dirichlet = vertcat(mesh.boundaries{:});
+ edgeRes(dirichlet) = 0;
+
+ %% combine the resdiuals suitably
+ hT = sqrt(getAffineTransformation(mesh).area);
+ indicators = hT.^2 .* volumeRes + ...
+ hT .* sum(edgeRes(mesh.element2edges), Dim.Vector);
+end
+
+% -\Delta u + [c;c]*\nabla u = 1
+function [lf, perturbation] = setUpConvectionProblem(fes, uold, coeff)
+ lf = LinearForm(fes);
+ perturbation = BilinearForm(fes);
+
+ lf.f = CompositeFunction(@(Du) 1 - vectorProduct(coeff*[1;1], Du), Gradient(uold));
+ perturbation.b = Constant(fes.mesh, coeff*[1;1]);
+end
+
+% -\Delta u - c^2*u = 1
+function [lf, perturbation] = setUpHelmholtzProblem(fes, uold, coeff)
+ lf = LinearForm(fes);
+ perturbation = BilinearForm(fes);
+
+ lf.f = CompositeFunction(@(u) 1 + coeff^2*u, uold);
+ lf.qrf = QuadratureRule.ofOrder(2);
+ perturbation.c = Constant(fes.mesh, -coeff^2);
+ perturbation.qrc = QuadratureRule.ofOrder(2);
+end
+
+function newMarked = addLargestElement(marked, mesh)
+ trafo = getAffineTransformation(mesh);
+ [~,maxElement] = max(trafo.area);
+ newMarked = union(marked, maxElement);
+end
--
GitLab