From 6ccee37169781a4e5272553ab9c92f7aa95d8fce Mon Sep 17 00:00:00 2001 From: Andrey A Popov Date: Wed, 10 May 2023 13:46:24 -0500 Subject: [PATCH 1/3] initial commit. THIS DOES NOT WORK --- src/+otp/+rangangermann/+presets/Canonical.m | 40 +++++++++ .../+rangangermann/RangAngermannParameters.m | 19 +++++ .../+rangangermann/RangAngermannProblem.m | 57 +++++++++++++ src/+otp/+rangangermann/f.m | 37 ++++++++ src/+otp/+rangangermann/jacobian.m | 19 +++++ src/+otp/+rangangermann/mass.m | 5 ++ src/+otp/+rangangermann/test.m | 84 +++++++++++++++++++ 7 files changed, 261 insertions(+) create mode 100644 src/+otp/+rangangermann/+presets/Canonical.m create mode 100644 src/+otp/+rangangermann/RangAngermannParameters.m create mode 100644 src/+otp/+rangangermann/RangAngermannProblem.m create mode 100644 src/+otp/+rangangermann/f.m create mode 100644 src/+otp/+rangangermann/jacobian.m create mode 100644 src/+otp/+rangangermann/mass.m create mode 100644 src/+otp/+rangangermann/test.m diff --git a/src/+otp/+rangangermann/+presets/Canonical.m b/src/+otp/+rangangermann/+presets/Canonical.m new file mode 100644 index 00000000..11868e49 --- /dev/null +++ b/src/+otp/+rangangermann/+presets/Canonical.m @@ -0,0 +1,40 @@ +classdef Canonical < otp.rangangermann.RangAngermannProblem + %CANONICAL + + methods + function obj = Canonical(varargin) + + p = inputParser; + p.addParameter('Size', 32); + + exactu = @(t, x, y) (2*x + y).*sin(t); + exactv = @(t, x, y) (x + 3*y).*cos(t); + + p.parse(varargin{:}); + + s = p.Results; + + n = s.Size; + + params = otp.rangangermann.RangAngermannParameters; + params.Size = n; + params.ForcingU = @(t, x, y) (2*x + y).*cos(t) + 2*x.*sin(t) + y.*sin(t) ... + - exactu(t, x, y) + exactv(t, x, y); + params.ForcingV = @(t, x, y) exactu(t, x, y).^2 + exactv(t, x, y).^2; + params.BoundaryConditionsU = exactu; + params.BoundaryConditionsV = exactv; + + x = linspace(0, 1, n + 2); + y = linspace(0, 1, n + 2); + [xfull, yfull] = meshgrid(x, y); + + x = xfull(2:(end-1), 2:(end-1)); + y = yfull(2:(end-1), 2:(end-1)); + uv0 = [reshape(exactu(0, x, y), [], 1); reshape(exactv(0, x, y), [], 1)]; + + tspan = [0, 1]; + + obj = obj@otp.rangangermann.RangAngermannProblem(tspan, uv0, params); + end + end +end diff --git a/src/+otp/+rangangermann/RangAngermannParameters.m b/src/+otp/+rangangermann/RangAngermannParameters.m new file mode 100644 index 00000000..bdc68cfb --- /dev/null +++ b/src/+otp/+rangangermann/RangAngermannParameters.m @@ -0,0 +1,19 @@ +classdef RangAngermannParameters + % User-configurable parameters for the Rang-Angermann Problem + % + % See also otp.rangangermann.RangAngermannProblem + + properties + %Size is the dimension of the problem + Size %MATLAB ONLY: (1,1) {mustBeInteger} + %ForcingU is a forcing function + ForcingU %MATLAB ONLY: {mustBeA(ForcingU, {'numeric', 'function_handle'})} + %ForcingV is a forcing function + ForcingV %MATLAB ONLY: {mustBeA(ForcingU, {'numeric', 'function_handle'})} + %BoundaryConditionsU is a boundary + BoundaryConditionsU %MATLAB ONLY: {mustBeA(BoundaryConditionsU, {'numeric', 'function_handle'})} + %BoundaryConditionsV is a boundary + BoundaryConditionsV %MATLAB ONLY: {mustBeA(BoundaryConditionsV, {'numeric', 'function_handle'})} + + end +end diff --git a/src/+otp/+rangangermann/RangAngermannProblem.m b/src/+otp/+rangangermann/RangAngermannProblem.m new file mode 100644 index 00000000..43477d01 --- /dev/null +++ b/src/+otp/+rangangermann/RangAngermannProblem.m @@ -0,0 +1,57 @@ +classdef RangAngermannProblem < otp.Problem + %RangAngermannProblem + + methods + function obj = RangAngermannProblem(timeSpan, y0, parameters) + obj@otp.Problem('Rang-Angermann Problem', [], timeSpan, y0, parameters); + end + end + + methods (Access = protected) + function onSettingsChanged(obj) + n = obj.Parameters.Size; + forcingu = obj.Parameters.ForcingU; + forcingv = obj.Parameters.ForcingV; + FBCu = obj.Parameters.BoundaryConditionsU; + FBCv = obj.Parameters.BoundaryConditionsV; + + if obj.NumVars ~= n^2 + warning('OTP:inconsistentNumVars', ... + 'NumVars is %d, but there are %d grid points', ... + obj.NumVars, n^2); + end + + domainx = [0, 1]; + domainy = [0, 1]; + domain = [domainx; domainy]; + L = otp.utils.pde.laplacian([n+2 n+2], domain, [1, 1], 'DD'); + Dx = otp.utils.pde.Dd([n+2 n+2], domainx, 1, 2, 'D'); + Dy = otp.utils.pde.Dd([n+2 n+2], domainy, 2, 2, 'D'); + + xb = linspace(0, 1, n + 2); + yb = linspace(0, 1, n + 2); + + [x, y] = meshgrid(xb, yb); + + xsmall = x(2:(end - 1), 2:(end - 1)); + ysmall = y(2:(end - 1), 2:(end - 1)); + + hx = 1/(n + 3); + hy = 1/(n + 3); + domainxsmall = [hx, 1 - hx]; + domainysmall = [hy, 1 - hy]; + domainsmall = [domainxsmall; domainysmall]; + Lsmall = otp.utils.pde.laplacian([n, n], domainsmall, [1, 1], 'DD'); + Dxsmall = otp.utils.pde.Dd([n, n], domainxsmall, 1, 2, 'D'); + Dysmall = otp.utils.pde.Dd([n, n], domainysmall, 2, 2, 'D'); + + Md = [ones(n^2, 1); zeros(n^2, 1)]; + M = spdiags(Md, 0, 2*(n^2), 2*(n^2)); + + obj.RHS = otp.RHS(@(t, uv) otp.rangangermann.f(t, uv, L, Dx, Dy, x, y, forcingu, forcingv, FBCu, FBCv), ... + 'Jacobian', @(t, uv) otp.rangangermann.jacobian(t, uv, Lsmall, Dxsmall, Dysmall, xsmall, ysmall, forcingu, forcingv, FBCu, FBCv), ... + 'Mass', M); + + end + end +end diff --git a/src/+otp/+rangangermann/f.m b/src/+otp/+rangangermann/f.m new file mode 100644 index 00000000..b237876f --- /dev/null +++ b/src/+otp/+rangangermann/f.m @@ -0,0 +1,37 @@ +function duvdt = f(t, uv, L, Dx, Dy, x, y, forcingu, forcingv, FBCu, FBCv) + +usmall = uv(1:(end/2)); +vsmall = uv((end/2 + 1):end); + +n = sqrt(numel(usmall)); + +u = FBCu(t, x, y); +v = FBCv(t, x, y); + +u(2:(end - 1), 2:(end - 1)) = reshape(usmall, n, n); +v(2:(end - 1), 2:(end - 1)) = reshape(vsmall, n, n); + +u = reshape(u, [], 1); +v = reshape(v, [], 1); + +Lu = L*u; +Lv = L*v; +Dxu = Dx*u; +Dyu = Dy*u; + + +f1 = reshape(forcingu(t, x, y), [], 1); +f2 = reshape(forcingv(t, x, y), [], 1); + +dudt = f1 + Lu + Lv - reshape(x, [], 1).*Dxu - reshape(y, [], 1).*Dyu + u - v; +valg = f2 + Lu + Lv - u.^2 - v.^2; + +dudt = reshape(dudt, n + 2, n + 2); +dudt = reshape(dudt(2:(end - 1), 2:(end - 1)), [], 1); + +valg = reshape(valg, n + 2, n + 2); +valg = reshape(valg(2:(end - 1), 2:(end - 1)), [], 1); + +duvdt = [dudt; valg]; + +end diff --git a/src/+otp/+rangangermann/jacobian.m b/src/+otp/+rangangermann/jacobian.m new file mode 100644 index 00000000..829e85d3 --- /dev/null +++ b/src/+otp/+rangangermann/jacobian.m @@ -0,0 +1,19 @@ +function J = jacobian(~, uv, L, Dx, Dy, x, y, ~, ~, ~, ~) + +usmall = uv(1:(end/2)); +vsmall = uv((end/2 + 1):end); + +n = sqrt(numel(usmall)); + +ddudu = L ... + - spdiags(reshape(x, [], 1), 0, n^2, n^2)*Dx ... + - spdiags(reshape(y, [], 1), 0, n^2, n^2)*Dy ... + + speye(n^2, n^2); +ddudv = L - speye(n^2, n^2); +ddvdu = L - spdiags(reshape(2*usmall, [], 1), 0, n^2, n^2); +ddvdv = L - spdiags(reshape(2*vsmall, [], 1), 0, n^2, n^2); + + +J = [ddudu, ddudv; ddvdu, ddvdv]; + +end diff --git a/src/+otp/+rangangermann/mass.m b/src/+otp/+rangangermann/mass.m new file mode 100644 index 00000000..95452e94 --- /dev/null +++ b/src/+otp/+rangangermann/mass.m @@ -0,0 +1,5 @@ +function M = mass(n) + + + +end diff --git a/src/+otp/+rangangermann/test.m b/src/+otp/+rangangermann/test.m new file mode 100644 index 00000000..2fadd45c --- /dev/null +++ b/src/+otp/+rangangermann/test.m @@ -0,0 +1,84 @@ +n = 64; + +tend = 1; + +model = otp.rangangermann.presets.Canonical('Size', n); +sol = ode15s(model.RHS.F, [0, tend], model.Y0, ... + odeset('Mass', model.RHS.Mass, 'MassSingular','yes', 'RelTol', 1e-9)); +uv = sol.y(:, end); + + +xb = linspace(0, 1, n + 2); +yb = linspace(0, 1, n + 2); + +[x, y] = meshgrid(xb, yb); + +exactu = @(t, x, y) (2*x + y).*sin(t); +exactv = @(t, x, y) (x + 3*y).*cos(t); + + +for i = 1:numel(sol.x) + t = sol.x(i); + uv = sol.y(:, i); + usmall = uv(1:(end/2)); + vsmall = uv((end/2 + 1):end); + + u = exactu(t, x, y); + v = exactv(t, x, y); + + u(2:(end - 1), 2:(end - 1)) = reshape(usmall, n, n); + v(2:(end - 1), 2:(end - 1)) = reshape(vsmall, n, n); + + subplot(1, 2, 1) + imagesc(u); colorbar; + subplot(1, 2, 2) + imagesc(exactu(t, x, y)); colorbar; + + pause(0.01) + + drawnow; + + + +end + +%f = model.RHS.F; + +% Japprox = model.RHS.Jacobian(0, uv); +% +% J = zeros(2*n^2, 2*n^2); +% h = sqrt(eps); +% +% for i = 1:numel(uv) +% e = zeros(2*n^2, 1); +% e(i) = 1; +% +% J(:, i) = (f(0, uv + h*e) - f(0, uv - h*e))/(2*h); +% +% end +% +% imagesc(J - Japprox); colorbar + + +usmall = uv(1:(end/2)); +vsmall = uv((end/2 + 1):end); + +xb = linspace(0, 1, n + 2); +yb = linspace(0, 1, n + 2); + +[x, y] = meshgrid(xb, yb); + +xsmall = reshape(x(2:(end - 1), 2:(end - 1)), [], 1); +ysmall = reshape(y(2:(end - 1), 2:(end - 1)), [], 1); + +exactu = @(t, x, y) (2*x + y).*sin(t); +exactv = @(t, x, y) (x + 3*y).*cos(t); + +eu = exactu(tend, xsmall, ysmall); +ev = exactv(tend, xsmall, ysmall); + +uerr = norm(eu - usmall)/norm(eu) +verr = norm(ev - vsmall)/norm(ev) + + + From aced7105d74bc6592adeb9aded74942e8de23407 Mon Sep 17 00:00:00 2001 From: Andrey A Popov Date: Wed, 10 May 2023 17:30:06 -0500 Subject: [PATCH 2/3] Fix everything --- src/+otp/+rangangermann/+presets/Canonical.m | 2 +- .../+rangangermann/RangAngermannProblem.m | 40 +++++++++++-------- src/+otp/+rangangermann/f.m | 4 +- src/+otp/+rangangermann/test.m | 13 +++--- 4 files changed, 32 insertions(+), 27 deletions(-) diff --git a/src/+otp/+rangangermann/+presets/Canonical.m b/src/+otp/+rangangermann/+presets/Canonical.m index 11868e49..a95d2943 100644 --- a/src/+otp/+rangangermann/+presets/Canonical.m +++ b/src/+otp/+rangangermann/+presets/Canonical.m @@ -26,7 +26,7 @@ x = linspace(0, 1, n + 2); y = linspace(0, 1, n + 2); - [xfull, yfull] = meshgrid(x, y); + [yfull, xfull] = meshgrid(x, y); x = xfull(2:(end-1), 2:(end-1)); y = yfull(2:(end-1), 2:(end-1)); diff --git a/src/+otp/+rangangermann/RangAngermannProblem.m b/src/+otp/+rangangermann/RangAngermannProblem.m index 43477d01..11e3d141 100644 --- a/src/+otp/+rangangermann/RangAngermannProblem.m +++ b/src/+otp/+rangangermann/RangAngermannProblem.m @@ -20,30 +20,36 @@ function onSettingsChanged(obj) 'NumVars is %d, but there are %d grid points', ... obj.NumVars, n^2); end + + %full domain + nfull = n + 2; - domainx = [0, 1]; - domainy = [0, 1]; - domain = [domainx; domainy]; - L = otp.utils.pde.laplacian([n+2 n+2], domain, [1, 1], 'DD'); - Dx = otp.utils.pde.Dd([n+2 n+2], domainx, 1, 2, 'D'); - Dy = otp.utils.pde.Dd([n+2 n+2], domainy, 2, 2, 'D'); + hx = 1/(n + 3); + hy = 1/(n + 3); + Ddx = spdiags(repmat([-1 1 0]/(hx), nfull, 1), [-1, 0, 1], nfull, nfull); + Dx = kron(speye(nfull), Ddx); + + Ddy = spdiags(repmat([-1 1 0]/(hy), nfull, 1), [-1, 0, 1], nfull, nfull); + Dy = kron(Ddy, speye(nfull)); + + Ldx = spdiags(repmat([1 -2 1]/(hx^2), nfull, 1), [-1, 0, 1], nfull, nfull); + Ldy = spdiags(repmat([1 -2 1]/(hy^2), nfull, 1), [-1, 0, 1], nfull, nfull); - xb = linspace(0, 1, n + 2); - yb = linspace(0, 1, n + 2); - [x, y] = meshgrid(xb, yb); + L = kron(speye(nfull), Ldx) + kron(Ldy, speye(nfull)); + + xb = linspace(0, 1, nfull); + yb = linspace(0, 1, nfull); + + [y, x] = meshgrid(xb, yb); xsmall = x(2:(end - 1), 2:(end - 1)); ysmall = y(2:(end - 1), 2:(end - 1)); - hx = 1/(n + 3); - hy = 1/(n + 3); - domainxsmall = [hx, 1 - hx]; - domainysmall = [hy, 1 - hy]; - domainsmall = [domainxsmall; domainysmall]; - Lsmall = otp.utils.pde.laplacian([n, n], domainsmall, [1, 1], 'DD'); - Dxsmall = otp.utils.pde.Dd([n, n], domainxsmall, 1, 2, 'D'); - Dysmall = otp.utils.pde.Dd([n, n], domainysmall, 2, 2, 'D'); + Lsmall = kron(speye(n), Ldx(2:(end - 1), 2:(end - 1))) ... + + kron(Ldy(2:(end - 1), 2:(end - 1)), speye(n)); + Dxsmall = kron(speye(n), Ddx(2:(end - 1), 2:(end - 1))); + Dysmall = kron(Ddy(2:(end - 1), 2:(end - 1)), speye(n)); Md = [ones(n^2, 1); zeros(n^2, 1)]; M = spdiags(Md, 0, 2*(n^2), 2*(n^2)); diff --git a/src/+otp/+rangangermann/f.m b/src/+otp/+rangangermann/f.m index b237876f..9ad47a25 100644 --- a/src/+otp/+rangangermann/f.m +++ b/src/+otp/+rangangermann/f.m @@ -19,11 +19,13 @@ Dxu = Dx*u; Dyu = Dy*u; +xDxu = reshape(x, [], 1).*Dxu; +yDyu = reshape(y, [], 1).*Dyu; f1 = reshape(forcingu(t, x, y), [], 1); f2 = reshape(forcingv(t, x, y), [], 1); -dudt = f1 + Lu + Lv - reshape(x, [], 1).*Dxu - reshape(y, [], 1).*Dyu + u - v; +dudt = f1 + Lu + Lv - xDxu - yDyu + u - v; valg = f2 + Lu + Lv - u.^2 - v.^2; dudt = reshape(dudt, n + 2, n + 2); diff --git a/src/+otp/+rangangermann/test.m b/src/+otp/+rangangermann/test.m index 2fadd45c..ca179220 100644 --- a/src/+otp/+rangangermann/test.m +++ b/src/+otp/+rangangermann/test.m @@ -4,14 +4,14 @@ model = otp.rangangermann.presets.Canonical('Size', n); sol = ode15s(model.RHS.F, [0, tend], model.Y0, ... - odeset('Mass', model.RHS.Mass, 'MassSingular','yes', 'RelTol', 1e-9)); + odeset('Mass', model.RHS.Mass, 'MassSingular','yes', 'RelTol', 1e-3, 'Jacobian', model.RHS.Jacobian)); uv = sol.y(:, end); xb = linspace(0, 1, n + 2); yb = linspace(0, 1, n + 2); -[x, y] = meshgrid(xb, yb); +[y, x] = meshgrid(xb, yb); exactu = @(t, x, y) (2*x + y).*sin(t); exactv = @(t, x, y) (x + 3*y).*cos(t); @@ -42,7 +42,7 @@ end -%f = model.RHS.F; +f = model.RHS.F; % Japprox = model.RHS.Jacobian(0, uv); % @@ -58,16 +58,13 @@ % end % % imagesc(J - Japprox); colorbar +% +% return usmall = uv(1:(end/2)); vsmall = uv((end/2 + 1):end); -xb = linspace(0, 1, n + 2); -yb = linspace(0, 1, n + 2); - -[x, y] = meshgrid(xb, yb); - xsmall = reshape(x(2:(end - 1), 2:(end - 1)), [], 1); ysmall = reshape(y(2:(end - 1), 2:(end - 1)), [], 1); From f0bdbb89f80b18ff455c12e6805f4a6f515dc2c6 Mon Sep 17 00:00:00 2001 From: Andrey A Popov Date: Thu, 11 May 2023 11:32:06 -0500 Subject: [PATCH 3/3] Update to rang-angermann problem --- src/+otp/+rangangermann/+presets/Canonical.m | 7 +- .../+rangangermann/RangAngermannProblem.m | 61 +++++++++++--- src/+otp/+rangangermann/test.m | 81 ------------------- 3 files changed, 57 insertions(+), 92 deletions(-) delete mode 100644 src/+otp/+rangangermann/test.m diff --git a/src/+otp/+rangangermann/+presets/Canonical.m b/src/+otp/+rangangermann/+presets/Canonical.m index a95d2943..86c1176c 100644 --- a/src/+otp/+rangangermann/+presets/Canonical.m +++ b/src/+otp/+rangangermann/+presets/Canonical.m @@ -5,22 +5,25 @@ function obj = Canonical(varargin) p = inputParser; - p.addParameter('Size', 32); + p.addParameter('Size', 256); + % exact solution exactu = @(t, x, y) (2*x + y).*sin(t); exactv = @(t, x, y) (x + 3*y).*cos(t); p.parse(varargin{:}); s = p.Results; - n = s.Size; params = otp.rangangermann.RangAngermannParameters; params.Size = n; + % the forcing functions are based on the excact solution params.ForcingU = @(t, x, y) (2*x + y).*cos(t) + 2*x.*sin(t) + y.*sin(t) ... - exactu(t, x, y) + exactv(t, x, y); params.ForcingV = @(t, x, y) exactu(t, x, y).^2 + exactv(t, x, y).^2; + + % the boundary conditions are the exact solutions params.BoundaryConditionsU = exactu; params.BoundaryConditionsV = exactv; diff --git a/src/+otp/+rangangermann/RangAngermannProblem.m b/src/+otp/+rangangermann/RangAngermannProblem.m index 11e3d141..e7022a34 100644 --- a/src/+otp/+rangangermann/RangAngermannProblem.m +++ b/src/+otp/+rangangermann/RangAngermannProblem.m @@ -1,5 +1,11 @@ classdef RangAngermannProblem < otp.Problem - %RangAngermannProblem + % This is an Index-1 PDAE + % + % See + % Rang, J., & Angermann, L. (2005). + % New Rosenbrock W-methods of order 3 for partial differential algebraic equations of index 1. + % BIT Numerical Mathematics, 45, 761-787. + % methods function obj = RangAngermannProblem(timeSpan, y0, parameters) @@ -21,43 +27,80 @@ function onSettingsChanged(obj) obj.NumVars, n^2); end - %full domain + % First we construct the finite difference operators in the + % full domain including the boundary conditions nfull = n + 2; hx = 1/(n + 3); hy = 1/(n + 3); + % one dimensional derivative in the x direciton Ddx = spdiags(repmat([-1 1 0]/(hx), nfull, 1), [-1, 0, 1], nfull, nfull); + % two dimensional derivative in the x direction Dx = kron(speye(nfull), Ddx); + % one dimensional derivative in the y direciton Ddy = spdiags(repmat([-1 1 0]/(hy), nfull, 1), [-1, 0, 1], nfull, nfull); + % two dimensional derivative in the y direction Dy = kron(Ddy, speye(nfull)); + % both one dimensional Laplacians Ldx = spdiags(repmat([1 -2 1]/(hx^2), nfull, 1), [-1, 0, 1], nfull, nfull); Ldy = spdiags(repmat([1 -2 1]/(hy^2), nfull, 1), [-1, 0, 1], nfull, nfull); - + % two dimensional Laplacian L = kron(speye(nfull), Ldx) + kron(Ldy, speye(nfull)); + % construct the x-y grid xb = linspace(0, 1, nfull); yb = linspace(0, 1, nfull); [y, x] = meshgrid(xb, yb); - xsmall = x(2:(end - 1), 2:(end - 1)); - ysmall = y(2:(end - 1), 2:(end - 1)); - Lsmall = kron(speye(n), Ldx(2:(end - 1), 2:(end - 1))) ... + % internal point grid + xinternal = x(2:(end - 1), 2:(end - 1)); + yinternal = y(2:(end - 1), 2:(end - 1)); + + % construct the x derivative operator on the internal grid + Dxinternal = kron(speye(n), Ddx(2:(end - 1), 2:(end - 1))); + % construct the y derivative operator on the internal grid + Dyinternal = kron(Ddy(2:(end - 1), 2:(end - 1)), speye(n)); + % construct the Laplacian on the internal grid + Linternal = kron(speye(n), Ldx(2:(end - 1), 2:(end - 1))) ... + kron(Ldy(2:(end - 1), 2:(end - 1)), speye(n)); - Dxsmall = kron(speye(n), Ddx(2:(end - 1), 2:(end - 1))); - Dysmall = kron(Ddy(2:(end - 1), 2:(end - 1)), speye(n)); + % construct the mass matrix Md = [ones(n^2, 1); zeros(n^2, 1)]; M = spdiags(Md, 0, 2*(n^2), 2*(n^2)); + % construct the RHS obj.RHS = otp.RHS(@(t, uv) otp.rangangermann.f(t, uv, L, Dx, Dy, x, y, forcingu, forcingv, FBCu, FBCv), ... - 'Jacobian', @(t, uv) otp.rangangermann.jacobian(t, uv, Lsmall, Dxsmall, Dysmall, xsmall, ysmall, forcingu, forcingv, FBCu, FBCv), ... + 'Jacobian', @(t, uv) otp.rangangermann.jacobian(t, uv, Linternal, Dxinternal, Dyinternal, xinternal, yinternal, forcingu, forcingv, FBCu, FBCv), ... 'Mass', M); end + + function uv = internalSolveExactly(obj, t) + n = obj.Parameters.Size; + nfull = n + 2; + + % construct the x-y grid + xb = linspace(0, 1, nfull); + yb = linspace(0, 1, nfull); + + [y, x] = meshgrid(xb, yb); + + % internal point grid + xinternal = x(2:(end - 1), 2:(end - 1)); + yinternal = y(2:(end - 1), 2:(end - 1)); + + for i = length(t):-1:1 + u = obj.Parameters.BoundaryConditionsU(t(i), xinternal, yinternal); + u = reshape(u, [], 1); + v = obj.Parameters.BoundaryConditionsV(t(i), xinternal, yinternal); + v = reshape(v, [], 1); + uv(:, i) = [u; v]; + end + end end end diff --git a/src/+otp/+rangangermann/test.m b/src/+otp/+rangangermann/test.m deleted file mode 100644 index ca179220..00000000 --- a/src/+otp/+rangangermann/test.m +++ /dev/null @@ -1,81 +0,0 @@ -n = 64; - -tend = 1; - -model = otp.rangangermann.presets.Canonical('Size', n); -sol = ode15s(model.RHS.F, [0, tend], model.Y0, ... - odeset('Mass', model.RHS.Mass, 'MassSingular','yes', 'RelTol', 1e-3, 'Jacobian', model.RHS.Jacobian)); -uv = sol.y(:, end); - - -xb = linspace(0, 1, n + 2); -yb = linspace(0, 1, n + 2); - -[y, x] = meshgrid(xb, yb); - -exactu = @(t, x, y) (2*x + y).*sin(t); -exactv = @(t, x, y) (x + 3*y).*cos(t); - - -for i = 1:numel(sol.x) - t = sol.x(i); - uv = sol.y(:, i); - usmall = uv(1:(end/2)); - vsmall = uv((end/2 + 1):end); - - u = exactu(t, x, y); - v = exactv(t, x, y); - - u(2:(end - 1), 2:(end - 1)) = reshape(usmall, n, n); - v(2:(end - 1), 2:(end - 1)) = reshape(vsmall, n, n); - - subplot(1, 2, 1) - imagesc(u); colorbar; - subplot(1, 2, 2) - imagesc(exactu(t, x, y)); colorbar; - - pause(0.01) - - drawnow; - - - -end - -f = model.RHS.F; - -% Japprox = model.RHS.Jacobian(0, uv); -% -% J = zeros(2*n^2, 2*n^2); -% h = sqrt(eps); -% -% for i = 1:numel(uv) -% e = zeros(2*n^2, 1); -% e(i) = 1; -% -% J(:, i) = (f(0, uv + h*e) - f(0, uv - h*e))/(2*h); -% -% end -% -% imagesc(J - Japprox); colorbar -% -% return - - -usmall = uv(1:(end/2)); -vsmall = uv((end/2 + 1):end); - -xsmall = reshape(x(2:(end - 1), 2:(end - 1)), [], 1); -ysmall = reshape(y(2:(end - 1), 2:(end - 1)), [], 1); - -exactu = @(t, x, y) (2*x + y).*sin(t); -exactv = @(t, x, y) (x + 3*y).*cos(t); - -eu = exactu(tend, xsmall, ysmall); -ev = exactv(tend, xsmall, ysmall); - -uerr = norm(eu - usmall)/norm(eu) -verr = norm(ev - vsmall)/norm(ev) - - -