From b849181deb4fdb92309ea07cd63e4b3ed5d5d4a5 Mon Sep 17 00:00:00 2001 From: Andrey A Popov Date: Fri, 4 Sep 2020 21:02:13 -0400 Subject: [PATCH 01/15] Initial Commit of KS Simple 1D Chaotic test problem that is significantly better understood than Lorenz '96 --- .../+kuramotosivashinsky/+presets/Canonical.m | 29 +++++++++++++++++ .../KuramotoSivashinskyProblem.m | 32 +++++++++++++++++++ src/+otp/+kuramotosivashinsky/f.m | 5 +++ src/+otp/+kuramotosivashinsky/jac.m | 5 +++ 4 files changed, 71 insertions(+) create mode 100644 src/+otp/+kuramotosivashinsky/+presets/Canonical.m create mode 100644 src/+otp/+kuramotosivashinsky/KuramotoSivashinskyProblem.m create mode 100644 src/+otp/+kuramotosivashinsky/f.m create mode 100644 src/+otp/+kuramotosivashinsky/jac.m diff --git a/src/+otp/+kuramotosivashinsky/+presets/Canonical.m b/src/+otp/+kuramotosivashinsky/+presets/Canonical.m new file mode 100644 index 00000000..430df49b --- /dev/null +++ b/src/+otp/+kuramotosivashinsky/+presets/Canonical.m @@ -0,0 +1,29 @@ +classdef Canonical < otp.kuramotosivashinsky.KuramotoSivashinskyProblem + + methods + function obj = Canonical(varargin) + + p = inputParser; + addParameter(p, 'Size', 64, @isscalar); + addParameter(p, 'L', 25, @isscalar); + + parse(p, varargin{:}); + + s = p.Results; + + n = s.Size; + + + params.n = n; + params.l = s.L; + + x = linspace(0, 10*pi, n + 1); + + u0 = 4*cos(x(1:end-1)).'; + + tspan = [0, 100]; + + obj = obj@otp.kuramotosivashinsky.KuramotoSivashinskyProblem(tspan, u0, params); + end + end +end diff --git a/src/+otp/+kuramotosivashinsky/KuramotoSivashinskyProblem.m b/src/+otp/+kuramotosivashinsky/KuramotoSivashinskyProblem.m new file mode 100644 index 00000000..2e24150e --- /dev/null +++ b/src/+otp/+kuramotosivashinsky/KuramotoSivashinskyProblem.m @@ -0,0 +1,32 @@ +classdef KuramotoSivashinskyProblem < otp.Problem + + methods + function obj = KuramotoSivashinskyProblem(timeSpan, y0, parameters) + obj@otp.Problem('Kuramoto Sivashinsky Problem', [], timeSpan, y0, parameters); + end + end + + methods (Access = protected) + function onSettingsChanged(obj) + n = obj.Parameters.n; + l = obj.Parameters.l; + + domain = [-l, l]; + D = otp.utils.pde.D(n, domain, 'C'); + L = otp.utils.pde.laplacian(n, domain, 1, 'C'); + + obj.Rhs = otp.Rhs(@(t, u) otp.kuramotosivashinsky.f(t, u, D, L), ... + otp.Rhs.FieldNames.Jacobian, @(t, u) otp.kuramotosivashinsky.jac(t, u, D, L)); + + end + + function validateNewState(obj, newTimeSpan, newY0, newParameters) + validateNewState@otp.Problem(obj, ... + newTimeSpan, newY0, newParameters) + + otp.utils.StructParser(newParameters) ... + .checkField('n', 'scalar', 'integer', 'finite', 'positive') ... + .checkField('l', 'scalar', 'integer', 'finite', 'positive'); + end + end +end diff --git a/src/+otp/+kuramotosivashinsky/f.m b/src/+otp/+kuramotosivashinsky/f.m new file mode 100644 index 00000000..15f328d2 --- /dev/null +++ b/src/+otp/+kuramotosivashinsky/f.m @@ -0,0 +1,5 @@ +function du = f(~, u, D, L) + +du = - L*(L*u + u) - u.*(D*u); + +end diff --git a/src/+otp/+kuramotosivashinsky/jac.m b/src/+otp/+kuramotosivashinsky/jac.m new file mode 100644 index 00000000..aaf9b104 --- /dev/null +++ b/src/+otp/+kuramotosivashinsky/jac.m @@ -0,0 +1,5 @@ +function j = jac(~, u, D, L) + +j = -L*L - L - spdiags(u, 0, numel(u), numel(u))*D - spdiags(D*u, 0, numel(u), numel(u)); + +end From c7beadc9e2cde6220682b4be9215f06ef0ab65e1 Mon Sep 17 00:00:00 2001 From: Andrey A Popov Date: Tue, 8 Sep 2020 20:35:57 -0400 Subject: [PATCH 02/15] Completely changed the problem to be spectral --- .../+kuramotosivashinsky/+presets/Canonical.m | 20 ++++---- .../KuramotoSivashinskyProblem.m | 46 +++++++++++++++---- src/+otp/+kuramotosivashinsky/f.m | 6 ++- src/+otp/+kuramotosivashinsky/jac.m | 4 +- src/+otp/+kuramotosivashinsky/jvp.m | 5 ++ 5 files changed, 59 insertions(+), 22 deletions(-) create mode 100644 src/+otp/+kuramotosivashinsky/jvp.m diff --git a/src/+otp/+kuramotosivashinsky/+presets/Canonical.m b/src/+otp/+kuramotosivashinsky/+presets/Canonical.m index 430df49b..8f54369c 100644 --- a/src/+otp/+kuramotosivashinsky/+presets/Canonical.m +++ b/src/+otp/+kuramotosivashinsky/+presets/Canonical.m @@ -4,26 +4,30 @@ function obj = Canonical(varargin) p = inputParser; - addParameter(p, 'Size', 64, @isscalar); - addParameter(p, 'L', 25, @isscalar); + addParameter(p, 'Size', 128); + addParameter(p, 'L', 32*pi); parse(p, varargin{:}); s = p.Results; - n = s.Size; + N = s.Size; + L = s.L; + params.L = L; - params.n = n; - params.l = s.L; + h=L/N; - x = linspace(0, 10*pi, n + 1); + x=h*(1:N).'; - u0 = 4*cos(x(1:end-1)).'; + u0 = cos(x/16).*(1+sin(x/16)); + + u0hat = fft(u0); tspan = [0, 100]; - obj = obj@otp.kuramotosivashinsky.KuramotoSivashinskyProblem(tspan, u0, params); + obj = obj@otp.kuramotosivashinsky.KuramotoSivashinskyProblem(tspan, u0hat, params); + end end end diff --git a/src/+otp/+kuramotosivashinsky/KuramotoSivashinskyProblem.m b/src/+otp/+kuramotosivashinsky/KuramotoSivashinskyProblem.m index 2e24150e..521d2f54 100644 --- a/src/+otp/+kuramotosivashinsky/KuramotoSivashinskyProblem.m +++ b/src/+otp/+kuramotosivashinsky/KuramotoSivashinskyProblem.m @@ -2,31 +2,57 @@ methods function obj = KuramotoSivashinskyProblem(timeSpan, y0, parameters) - obj@otp.Problem('Kuramoto Sivashinsky Problem', [], timeSpan, y0, parameters); + obj@otp.Problem('Kuramoto-Sivashinsky', [], timeSpan, y0, parameters); end end + methods + function soly = solution2real(soly) + + if isstruct(soly) + soly.y = real(ifft(soly.y)); + else + soly = real(ifft(soly.')).'; + end + + end + + end + methods (Access = protected) function onSettingsChanged(obj) - n = obj.Parameters.n; - l = obj.Parameters.l; - domain = [-l, l]; - D = otp.utils.pde.D(n, domain, 'C'); - L = otp.utils.pde.laplacian(n, domain, 1, 'C'); + L = obj.Parameters.L; - obj.Rhs = otp.Rhs(@(t, u) otp.kuramotosivashinsky.f(t, u, D, L), ... - otp.Rhs.FieldNames.Jacobian, @(t, u) otp.kuramotosivashinsky.jac(t, u, D, L)); + N = numel(obj.Y0); + + div = L/(2*pi); + + k = (1i*[0:(N/2 - 1), 0, (-N/2 + 1):-1].'/div); + k2 = k.^2; + k4 = k.^4; + + Dfft = dftmtx(N); + Difft = conj(Dfft)/N; + + obj.Rhs = otp.Rhs(@(t, u) otp.kuramotosivashinsky.f(t, u, k, k2, k4), ... + otp.Rhs.FieldNames.Jacobian, ... + @(t, u) otp.kuramotosivashinsky.jac(t,u, k, k2, k4, Dfft, Difft), ... + otp.Rhs.FieldNames.JacobianVectorProduct, ... + @(t, u, v) otp.kuramotosivashinsky.jvp(t,u, k, k2, k4, v)); end function validateNewState(obj, newTimeSpan, newY0, newParameters) + validateNewState@otp.Problem(obj, ... newTimeSpan, newY0, newParameters) + assert(mod(numel(newY0), 2) == 0); + otp.utils.StructParser(newParameters) ... - .checkField('n', 'scalar', 'integer', 'finite', 'positive') ... - .checkField('l', 'scalar', 'integer', 'finite', 'positive'); + .checkField('L', 'scalar', 'finite', 'positive'); + end end end diff --git a/src/+otp/+kuramotosivashinsky/f.m b/src/+otp/+kuramotosivashinsky/f.m index 15f328d2..1c58fb1c 100644 --- a/src/+otp/+kuramotosivashinsky/f.m +++ b/src/+otp/+kuramotosivashinsky/f.m @@ -1,5 +1,7 @@ -function du = f(~, u, D, L) +function ut = f(~, u, k, k2, k4) -du = - L*(L*u + u) - u.*(D*u); +u2 = fft(real(ifft(u)).^2); + +ut = - k2.*u - k4.*u -(k/2).*u2; end diff --git a/src/+otp/+kuramotosivashinsky/jac.m b/src/+otp/+kuramotosivashinsky/jac.m index aaf9b104..b75bcacf 100644 --- a/src/+otp/+kuramotosivashinsky/jac.m +++ b/src/+otp/+kuramotosivashinsky/jac.m @@ -1,5 +1,5 @@ -function j = jac(~, u, D, L) +function j = jac(~, u, k, k2, k4, Dfft, Difft) -j = -L*L - L - spdiags(u, 0, numel(u), numel(u))*D - spdiags(D*u, 0, numel(u), numel(u)); +j = -diag(k2) - diag(k4) - diag(k)*Dfft*diag(ifft(u))*Difft; end diff --git a/src/+otp/+kuramotosivashinsky/jvp.m b/src/+otp/+kuramotosivashinsky/jvp.m new file mode 100644 index 00000000..f340ed33 --- /dev/null +++ b/src/+otp/+kuramotosivashinsky/jvp.m @@ -0,0 +1,5 @@ +function jv = jvp(~, u, k, k2, k4, v) + +jv = -k2.*v - k4.*v - k.*fft(real(ifft(u)).*real(ifft(v))); + +end From e2f3c6a008dffe5d185e010c1577e6ecb22f7255 Mon Sep 17 00:00:00 2001 From: Andrey A Popov Date: Tue, 8 Sep 2020 22:44:32 -0400 Subject: [PATCH 03/15] Made jacobian more optimal; other Note: the jacobian is dense for all non-trivial inputs. for the initial condition it might seem like the jacobian is sparse, but this is a red herring. The jacobian is dense, and the way in which it is computed in this code is as efficient as I dare make it. --- src/+otp/+kuramotosivashinsky/KuramotoSivashinskyProblem.m | 7 ++----- src/+otp/+kuramotosivashinsky/jac.m | 4 ++-- 2 files changed, 4 insertions(+), 7 deletions(-) diff --git a/src/+otp/+kuramotosivashinsky/KuramotoSivashinskyProblem.m b/src/+otp/+kuramotosivashinsky/KuramotoSivashinskyProblem.m index 521d2f54..bed9637d 100644 --- a/src/+otp/+kuramotosivashinsky/KuramotoSivashinskyProblem.m +++ b/src/+otp/+kuramotosivashinsky/KuramotoSivashinskyProblem.m @@ -24,20 +24,17 @@ function onSettingsChanged(obj) L = obj.Parameters.L; - N = numel(obj.Y0); + N = obj.NumVars; div = L/(2*pi); k = (1i*[0:(N/2 - 1), 0, (-N/2 + 1):-1].'/div); k2 = k.^2; k4 = k.^4; - - Dfft = dftmtx(N); - Difft = conj(Dfft)/N; obj.Rhs = otp.Rhs(@(t, u) otp.kuramotosivashinsky.f(t, u, k, k2, k4), ... otp.Rhs.FieldNames.Jacobian, ... - @(t, u) otp.kuramotosivashinsky.jac(t,u, k, k2, k4, Dfft, Difft), ... + @(t, u) otp.kuramotosivashinsky.jac(t,u, k, k2, k4), ... otp.Rhs.FieldNames.JacobianVectorProduct, ... @(t, u, v) otp.kuramotosivashinsky.jvp(t,u, k, k2, k4, v)); diff --git a/src/+otp/+kuramotosivashinsky/jac.m b/src/+otp/+kuramotosivashinsky/jac.m index b75bcacf..8be8bc71 100644 --- a/src/+otp/+kuramotosivashinsky/jac.m +++ b/src/+otp/+kuramotosivashinsky/jac.m @@ -1,5 +1,5 @@ -function j = jac(~, u, k, k2, k4, Dfft, Difft) +function j = jac(~, u, k, k2, k4) -j = -diag(k2) - diag(k4) - diag(k)*Dfft*diag(ifft(u))*Difft; +j = -diag(k2) - diag(k4) - diag(k)*ifft(fft(diag(ifft(u))).').'; end From dd2f38af2d1167efa171fd723fdd3603ee5a0cbf Mon Sep 17 00:00:00 2001 From: Andrey A Popov Date: Wed, 9 Sep 2020 11:08:48 -0400 Subject: [PATCH 04/15] Added Jacobian Adjoint (conjugate as well) vector product. --- src/+otp/+kuramotosivashinsky/KuramotoSivashinskyProblem.m | 4 +++- src/+otp/+kuramotosivashinsky/javp.m | 5 +++++ src/+otp/+kuramotosivashinsky/jvp.m | 2 +- 3 files changed, 9 insertions(+), 2 deletions(-) create mode 100644 src/+otp/+kuramotosivashinsky/javp.m diff --git a/src/+otp/+kuramotosivashinsky/KuramotoSivashinskyProblem.m b/src/+otp/+kuramotosivashinsky/KuramotoSivashinskyProblem.m index bed9637d..61c00aee 100644 --- a/src/+otp/+kuramotosivashinsky/KuramotoSivashinskyProblem.m +++ b/src/+otp/+kuramotosivashinsky/KuramotoSivashinskyProblem.m @@ -36,7 +36,9 @@ function onSettingsChanged(obj) otp.Rhs.FieldNames.Jacobian, ... @(t, u) otp.kuramotosivashinsky.jac(t,u, k, k2, k4), ... otp.Rhs.FieldNames.JacobianVectorProduct, ... - @(t, u, v) otp.kuramotosivashinsky.jvp(t,u, k, k2, k4, v)); + @(t, u, v) otp.kuramotosivashinsky.jvp(t,u, k, k2, k4, v), ... + otp.Rhs.FieldNames.JacobianAdjointVectorProduct, ... + @(t, u, v) otp.kuramotosivashinsky.javp(t,u, k, k2, k4, v)); end diff --git a/src/+otp/+kuramotosivashinsky/javp.m b/src/+otp/+kuramotosivashinsky/javp.m new file mode 100644 index 00000000..05577830 --- /dev/null +++ b/src/+otp/+kuramotosivashinsky/javp.m @@ -0,0 +1,5 @@ +function jv = javp(~, u, k, k2, k4, v) + +jv = -k2.*v - k4.*v - fft(real(ifft(u)).*ifft(conj(k).*v)); + +end diff --git a/src/+otp/+kuramotosivashinsky/jvp.m b/src/+otp/+kuramotosivashinsky/jvp.m index f340ed33..da28a3d8 100644 --- a/src/+otp/+kuramotosivashinsky/jvp.m +++ b/src/+otp/+kuramotosivashinsky/jvp.m @@ -1,5 +1,5 @@ function jv = jvp(~, u, k, k2, k4, v) -jv = -k2.*v - k4.*v - k.*fft(real(ifft(u)).*real(ifft(v))); +jv = -k2.*v - k4.*v - k.*fft(real(ifft(u)).*ifft(v)); end From 09e0966e088f80f0ed27c1d183a9fed276a00dec Mon Sep 17 00:00:00 2001 From: Andrey A Popov Date: Wed, 9 Sep 2020 13:05:21 -0400 Subject: [PATCH 05/15] added citation and fixed the tspan --- .../+kuramotosivashinsky/+presets/Canonical.m | 20 ++++++++++++++++--- 1 file changed, 17 insertions(+), 3 deletions(-) diff --git a/src/+otp/+kuramotosivashinsky/+presets/Canonical.m b/src/+otp/+kuramotosivashinsky/+presets/Canonical.m index 8f54369c..a56b80a8 100644 --- a/src/+otp/+kuramotosivashinsky/+presets/Canonical.m +++ b/src/+otp/+kuramotosivashinsky/+presets/Canonical.m @@ -1,3 +1,17 @@ +% The Kuramoto-Sivashinky equation is a chaotic problem. +% +% In this particular discretization, we are applying a spectral +% method, therefore the boundary conditions will be chosen to be as cyclic +% on the domain [0, L]. Note that this is different from another typical +% domain of [-L, L]. The larger the L, the more interesting the problem is +% but the more points are required to do a good discretization. The current +% canonical implementation with the size, L, and is used in +% +% Kassam, Aly-Khan, and Lloyd N. Trefethen. +% "Fourth-order time-stepping for stiff PDEs." +% SIAM Journal on Scientific Computing 26, no. 4 (2005): 1214-1233. +% + classdef Canonical < otp.kuramotosivashinsky.KuramotoSivashinskyProblem methods @@ -16,15 +30,15 @@ params.L = L; - h=L/N; + h = L/N; - x=h*(1:N).'; + x = h*(1:N).'; u0 = cos(x/16).*(1+sin(x/16)); u0hat = fft(u0); - tspan = [0, 100]; + tspan = [0, 150]; obj = obj@otp.kuramotosivashinsky.KuramotoSivashinskyProblem(tspan, u0hat, params); From 7e63fe29a33be0e486e2136879e334757d45b671 Mon Sep 17 00:00:00 2001 From: Andrey A Popov Date: Wed, 9 Sep 2020 15:56:25 -0400 Subject: [PATCH 06/15] applied fixes --- src/+otp/+kuramotosivashinsky/+presets/Canonical.m | 4 +++- .../+kuramotosivashinsky/KuramotoSivashinskyProblem.m | 8 +++++--- src/+otp/+kuramotosivashinsky/javp.m | 2 +- src/+otp/+kuramotosivashinsky/jvp.m | 2 +- 4 files changed, 10 insertions(+), 6 deletions(-) diff --git a/src/+otp/+kuramotosivashinsky/+presets/Canonical.m b/src/+otp/+kuramotosivashinsky/+presets/Canonical.m index a56b80a8..6e993ccb 100644 --- a/src/+otp/+kuramotosivashinsky/+presets/Canonical.m +++ b/src/+otp/+kuramotosivashinsky/+presets/Canonical.m @@ -32,7 +32,9 @@ h = L/N; - x = h*(1:N).'; + % exclude the left boundary point as it is identical to the + % right boundary point + x = linspace(h, L, N).'; u0 = cos(x/16).*(1+sin(x/16)); diff --git a/src/+otp/+kuramotosivashinsky/KuramotoSivashinskyProblem.m b/src/+otp/+kuramotosivashinsky/KuramotoSivashinskyProblem.m index 61c00aee..a18d8a45 100644 --- a/src/+otp/+kuramotosivashinsky/KuramotoSivashinskyProblem.m +++ b/src/+otp/+kuramotosivashinsky/KuramotoSivashinskyProblem.m @@ -36,9 +36,9 @@ function onSettingsChanged(obj) otp.Rhs.FieldNames.Jacobian, ... @(t, u) otp.kuramotosivashinsky.jac(t,u, k, k2, k4), ... otp.Rhs.FieldNames.JacobianVectorProduct, ... - @(t, u, v) otp.kuramotosivashinsky.jvp(t,u, k, k2, k4, v), ... + @(t, u, v) otp.kuramotosivashinsky.jvp(t, u, v, k, k2, k4), ... otp.Rhs.FieldNames.JacobianAdjointVectorProduct, ... - @(t, u, v) otp.kuramotosivashinsky.javp(t,u, k, k2, k4, v)); + @(t, u, v) otp.kuramotosivashinsky.javp(t, u, v, k, k2, k4)); end @@ -47,7 +47,9 @@ function validateNewState(obj, newTimeSpan, newY0, newParameters) validateNewState@otp.Problem(obj, ... newTimeSpan, newY0, newParameters) - assert(mod(numel(newY0), 2) == 0); + if mod(numel(newY0), 2) ~= 0 + error('The problem size has to be an even integer.'); + end otp.utils.StructParser(newParameters) ... .checkField('L', 'scalar', 'finite', 'positive'); diff --git a/src/+otp/+kuramotosivashinsky/javp.m b/src/+otp/+kuramotosivashinsky/javp.m index 05577830..e26d9b10 100644 --- a/src/+otp/+kuramotosivashinsky/javp.m +++ b/src/+otp/+kuramotosivashinsky/javp.m @@ -1,4 +1,4 @@ -function jv = javp(~, u, k, k2, k4, v) +function jv = javp(~, u, v, k, k2, k4) jv = -k2.*v - k4.*v - fft(real(ifft(u)).*ifft(conj(k).*v)); diff --git a/src/+otp/+kuramotosivashinsky/jvp.m b/src/+otp/+kuramotosivashinsky/jvp.m index da28a3d8..c208d029 100644 --- a/src/+otp/+kuramotosivashinsky/jvp.m +++ b/src/+otp/+kuramotosivashinsky/jvp.m @@ -1,4 +1,4 @@ -function jv = jvp(~, u, k, k2, k4, v) +function jv = jvp(~, u, v, k, k2, k4) jv = -k2.*v - k4.*v - k.*fft(real(ifft(u)).*ifft(v)); From 91d439c8d5691993f606cb279b34400cecf7a923 Mon Sep 17 00:00:00 2001 From: Andrey A Popov Date: Thu, 8 Jul 2021 15:59:01 -0400 Subject: [PATCH 07/15] Steven's changes all implemented. --- .../+kuramotosivashinsky/KuramotoSivashinskyProblem.m | 11 ++++++----- src/+otp/+kuramotosivashinsky/f.m | 4 ++-- src/+otp/+kuramotosivashinsky/jac.m | 4 ++-- src/+otp/+kuramotosivashinsky/javp.m | 4 ++-- src/+otp/+kuramotosivashinsky/jvp.m | 4 ++-- 5 files changed, 14 insertions(+), 13 deletions(-) diff --git a/src/+otp/+kuramotosivashinsky/KuramotoSivashinskyProblem.m b/src/+otp/+kuramotosivashinsky/KuramotoSivashinskyProblem.m index a18d8a45..3c1d8686 100644 --- a/src/+otp/+kuramotosivashinsky/KuramotoSivashinskyProblem.m +++ b/src/+otp/+kuramotosivashinsky/KuramotoSivashinskyProblem.m @@ -7,7 +7,7 @@ end methods - function soly = solution2real(soly) + function soly = convert2grid(soly) if isstruct(soly) soly.y = real(ifft(soly.y)); @@ -31,14 +31,15 @@ function onSettingsChanged(obj) k = (1i*[0:(N/2 - 1), 0, (-N/2 + 1):-1].'/div); k2 = k.^2; k4 = k.^4; + k24 = k2 + k4; - obj.Rhs = otp.Rhs(@(t, u) otp.kuramotosivashinsky.f(t, u, k, k2, k4), ... + obj.Rhs = otp.Rhs(@(t, u) otp.kuramotosivashinsky.f(t, u, k, k24), ... otp.Rhs.FieldNames.Jacobian, ... - @(t, u) otp.kuramotosivashinsky.jac(t,u, k, k2, k4), ... + @(t, u) otp.kuramotosivashinsky.jac(t,u, k, k24), ... otp.Rhs.FieldNames.JacobianVectorProduct, ... - @(t, u, v) otp.kuramotosivashinsky.jvp(t, u, v, k, k2, k4), ... + @(t, u, v) otp.kuramotosivashinsky.jvp(t, u, v, k, k24), ... otp.Rhs.FieldNames.JacobianAdjointVectorProduct, ... - @(t, u, v) otp.kuramotosivashinsky.javp(t, u, v, k, k2, k4)); + @(t, u, v) otp.kuramotosivashinsky.javp(t, u, v, k, k24)); end diff --git a/src/+otp/+kuramotosivashinsky/f.m b/src/+otp/+kuramotosivashinsky/f.m index 1c58fb1c..bcd01ce3 100644 --- a/src/+otp/+kuramotosivashinsky/f.m +++ b/src/+otp/+kuramotosivashinsky/f.m @@ -1,7 +1,7 @@ -function ut = f(~, u, k, k2, k4) +function ut = f(~, u, k, k24) u2 = fft(real(ifft(u)).^2); -ut = - k2.*u - k4.*u -(k/2).*u2; +ut = -k24.*u - (k/2).*u2; end diff --git a/src/+otp/+kuramotosivashinsky/jac.m b/src/+otp/+kuramotosivashinsky/jac.m index 8be8bc71..df64dada 100644 --- a/src/+otp/+kuramotosivashinsky/jac.m +++ b/src/+otp/+kuramotosivashinsky/jac.m @@ -1,5 +1,5 @@ -function j = jac(~, u, k, k2, k4) +function j = jac(~, u, k, k24) -j = -diag(k2) - diag(k4) - diag(k)*ifft(fft(diag(ifft(u))).').'; +j = -diag(k24) - k.*ifft(fft(diag(ifft(u))).').'; end diff --git a/src/+otp/+kuramotosivashinsky/javp.m b/src/+otp/+kuramotosivashinsky/javp.m index e26d9b10..141f7574 100644 --- a/src/+otp/+kuramotosivashinsky/javp.m +++ b/src/+otp/+kuramotosivashinsky/javp.m @@ -1,5 +1,5 @@ -function jv = javp(~, u, v, k, k2, k4) +function jv = javp(~, u, v, k, k24) -jv = -k2.*v - k4.*v - fft(real(ifft(u)).*ifft(conj(k).*v)); +jv = -k24.*v - fft(real(ifft(u)).*ifft(conj(k).*v)); end diff --git a/src/+otp/+kuramotosivashinsky/jvp.m b/src/+otp/+kuramotosivashinsky/jvp.m index c208d029..9e01ff41 100644 --- a/src/+otp/+kuramotosivashinsky/jvp.m +++ b/src/+otp/+kuramotosivashinsky/jvp.m @@ -1,5 +1,5 @@ -function jv = jvp(~, u, v, k, k2, k4) +function jv = jvp(~, u, v, k, k24) -jv = -k2.*v - k4.*v - k.*fft(real(ifft(u)).*ifft(v)); +jv = -k24.*v - k.*fft(real(ifft(u)).*ifft(v)); end From 9a8a63c3e8967c7dfb8e536d34456bff8d2636da Mon Sep 17 00:00:00 2001 From: Andrey A Popov Date: Thu, 8 Jul 2021 18:00:11 -0400 Subject: [PATCH 08/15] more of steven's fixes --- .../+kuramotosivashinsky/KuramotoSivashinskyProblem.m | 8 ++------ src/+otp/+kuramotosivashinsky/jac.m | 2 +- 2 files changed, 3 insertions(+), 7 deletions(-) diff --git a/src/+otp/+kuramotosivashinsky/KuramotoSivashinskyProblem.m b/src/+otp/+kuramotosivashinsky/KuramotoSivashinskyProblem.m index 3c1d8686..e6d4f2a6 100644 --- a/src/+otp/+kuramotosivashinsky/KuramotoSivashinskyProblem.m +++ b/src/+otp/+kuramotosivashinsky/KuramotoSivashinskyProblem.m @@ -7,13 +7,9 @@ end methods - function soly = convert2grid(soly) + function soly = convert2grid(~, soly) - if isstruct(soly) - soly.y = real(ifft(soly.y)); - else - soly = real(ifft(soly.')).'; - end + soly = real(ifft(soly)); end diff --git a/src/+otp/+kuramotosivashinsky/jac.m b/src/+otp/+kuramotosivashinsky/jac.m index df64dada..08b6a47e 100644 --- a/src/+otp/+kuramotosivashinsky/jac.m +++ b/src/+otp/+kuramotosivashinsky/jac.m @@ -1,5 +1,5 @@ function j = jac(~, u, k, k24) -j = -diag(k24) - k.*ifft(fft(diag(ifft(u))).').'; +j = -diag(k24) - k.*ifft(fft(diag(real(ifft(u)))).').'; end From 3eb99cd86b4b9c17d3c5b406de6da3c3f08c892b Mon Sep 17 00:00:00 2001 From: Andrey A Popov Date: Thu, 8 Jul 2021 18:01:11 -0400 Subject: [PATCH 09/15] added a comment for future idiot me --- src/+otp/+kuramotosivashinsky/KuramotoSivashinskyProblem.m | 1 + 1 file changed, 1 insertion(+) diff --git a/src/+otp/+kuramotosivashinsky/KuramotoSivashinskyProblem.m b/src/+otp/+kuramotosivashinsky/KuramotoSivashinskyProblem.m index e6d4f2a6..ea45f664 100644 --- a/src/+otp/+kuramotosivashinsky/KuramotoSivashinskyProblem.m +++ b/src/+otp/+kuramotosivashinsky/KuramotoSivashinskyProblem.m @@ -24,6 +24,7 @@ function onSettingsChanged(obj) div = L/(2*pi); + % note that k already has "i" in it k = (1i*[0:(N/2 - 1), 0, (-N/2 + 1):-1].'/div); k2 = k.^2; k4 = k.^4; From c1bf87ef9bc7e97e3ee1450d992738bd0b772564 Mon Sep 17 00:00:00 2001 From: Andrey A Popov Date: Thu, 8 Jul 2021 19:03:27 -0400 Subject: [PATCH 10/15] Change initial condition to always be periodic. --- src/+otp/+kuramotosivashinsky/+presets/Canonical.m | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/+otp/+kuramotosivashinsky/+presets/Canonical.m b/src/+otp/+kuramotosivashinsky/+presets/Canonical.m index 6e993ccb..e84a6147 100644 --- a/src/+otp/+kuramotosivashinsky/+presets/Canonical.m +++ b/src/+otp/+kuramotosivashinsky/+presets/Canonical.m @@ -36,7 +36,7 @@ % right boundary point x = linspace(h, L, N).'; - u0 = cos(x/16).*(1+sin(x/16)); + u0 = cospi(2*x/L).*(1+sinpi(2*x/L)); u0hat = fft(u0); From 1b53abc4b0b19779860af8829665d19835a95477 Mon Sep 17 00:00:00 2001 From: Andrey A Popov Date: Mon, 12 Jul 2021 10:53:17 -0400 Subject: [PATCH 11/15] Added the preset from Sol ODES2 --- .../+presets/HairerWanner.m | 40 +++++++++++++++++++ 1 file changed, 40 insertions(+) create mode 100644 src/+otp/+kuramotosivashinsky/+presets/HairerWanner.m diff --git a/src/+otp/+kuramotosivashinsky/+presets/HairerWanner.m b/src/+otp/+kuramotosivashinsky/+presets/HairerWanner.m new file mode 100644 index 00000000..414a1e30 --- /dev/null +++ b/src/+otp/+kuramotosivashinsky/+presets/HairerWanner.m @@ -0,0 +1,40 @@ +classdef HairerWanner < otp.kuramotosivashinsky.KuramotoSivashinskyProblem + + methods + function obj = HairerWanner(varargin) + + p = inputParser; + addParameter(p, 'Size', 128); + addParameter(p, 'L', 80*pi); + + parse(p, varargin{:}); + + s = p.Results; + + N = s.Size; + L = s.L; + + params.L = L; + + h = L/N; + + % exclude the left boundary point as it is identical to the + % right boundary point + x = linspace(h, L, N).'; + + eta1 = min(x/L, 0.1 - x/L); + eta2 = 20*(x/L - 0.2).*(0.3 - x/L); + eta3 = min(x/L - 0.6, 0.7 - x/L); + eta4 = min(x/L - 0.9, 1 - x/L); + + u0 = 16*max(0, max(eta1, max(eta2, max(eta3, eta4)))); + + u0hat = fft(u0); + + tspan = [0, 150]; + + obj = obj@otp.kuramotosivashinsky.KuramotoSivashinskyProblem(tspan, u0hat, params); + + end + end +end From 33b0987156a495577ee06b3b16fb37b800832e31 Mon Sep 17 00:00:00 2001 From: Andrey A Popov Date: Mon, 12 Jul 2021 11:56:51 -0400 Subject: [PATCH 12/15] Minor fix to N --- src/+otp/+kuramotosivashinsky/+presets/HairerWanner.m | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/+otp/+kuramotosivashinsky/+presets/HairerWanner.m b/src/+otp/+kuramotosivashinsky/+presets/HairerWanner.m index 414a1e30..7b08608b 100644 --- a/src/+otp/+kuramotosivashinsky/+presets/HairerWanner.m +++ b/src/+otp/+kuramotosivashinsky/+presets/HairerWanner.m @@ -4,7 +4,7 @@ function obj = HairerWanner(varargin) p = inputParser; - addParameter(p, 'Size', 128); + addParameter(p, 'Size', 1024); addParameter(p, 'L', 80*pi); parse(p, varargin{:}); From 09f0c921ebd8f1d23fe639e07822827c521c711f Mon Sep 17 00:00:00 2001 From: Andrey A Popov Date: Mon, 12 Jul 2021 12:13:41 -0400 Subject: [PATCH 13/15] tspan fix --- src/+otp/+kuramotosivashinsky/+presets/HairerWanner.m | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/+otp/+kuramotosivashinsky/+presets/HairerWanner.m b/src/+otp/+kuramotosivashinsky/+presets/HairerWanner.m index 7b08608b..93824d14 100644 --- a/src/+otp/+kuramotosivashinsky/+presets/HairerWanner.m +++ b/src/+otp/+kuramotosivashinsky/+presets/HairerWanner.m @@ -31,7 +31,7 @@ u0hat = fft(u0); - tspan = [0, 150]; + tspan = [0, 100]; obj = obj@otp.kuramotosivashinsky.KuramotoSivashinskyProblem(tspan, u0hat, params); From 8c0e7a3463b0412fe71b33b3d8e3c65f1727176d Mon Sep 17 00:00:00 2001 From: Steven Roberts Date: Wed, 28 Jul 2021 20:35:48 -0400 Subject: [PATCH 14/15] Improvements for Jacobian and JAVP --- .../+kuramotosivashinsky/+presets/Canonical.m | 14 ++++------- .../+presets/HairerWanner.m | 12 ++++------ .../KuramotoSivashinskyProblem.m | 24 +++++++------------ src/+otp/+kuramotosivashinsky/f.m | 6 ++--- src/+otp/+kuramotosivashinsky/jac.m | 5 ++-- src/+otp/+kuramotosivashinsky/javp.m | 4 ++-- src/+otp/+kuramotosivashinsky/jvp.m | 4 ++-- 7 files changed, 27 insertions(+), 42 deletions(-) diff --git a/src/+otp/+kuramotosivashinsky/+presets/Canonical.m b/src/+otp/+kuramotosivashinsky/+presets/Canonical.m index e84a6147..7800a0b6 100644 --- a/src/+otp/+kuramotosivashinsky/+presets/Canonical.m +++ b/src/+otp/+kuramotosivashinsky/+presets/Canonical.m @@ -7,8 +7,8 @@ % but the more points are required to do a good discretization. The current % canonical implementation with the size, L, and is used in % -% Kassam, Aly-Khan, and Lloyd N. Trefethen. -% "Fourth-order time-stepping for stiff PDEs." +% Kassam, Aly-Khan, and Lloyd N. Trefethen. +% "Fourth-order time-stepping for stiff PDEs." % SIAM Journal on Scientific Computing 26, no. 4 (2005): 1214-1233. % @@ -20,26 +20,22 @@ p = inputParser; addParameter(p, 'Size', 128); addParameter(p, 'L', 32*pi); - parse(p, varargin{:}); s = p.Results; N = s.Size; L = s.L; - params.L = L; - h = L/N; - % exclude the left boundary point as it is identical to the % right boundary point - x = linspace(h, L, N).'; + x = linspace(L / N, L, N).'; - u0 = cospi(2*x/L).*(1+sinpi(2*x/L)); + u0 = cospi(2 * x / L) .* (1 + sinpi(2 * x / L)); u0hat = fft(u0); - + tspan = [0, 150]; obj = obj@otp.kuramotosivashinsky.KuramotoSivashinskyProblem(tspan, u0hat, params); diff --git a/src/+otp/+kuramotosivashinsky/+presets/HairerWanner.m b/src/+otp/+kuramotosivashinsky/+presets/HairerWanner.m index 93824d14..04030db2 100644 --- a/src/+otp/+kuramotosivashinsky/+presets/HairerWanner.m +++ b/src/+otp/+kuramotosivashinsky/+presets/HairerWanner.m @@ -5,8 +5,8 @@ p = inputParser; addParameter(p, 'Size', 1024); - addParameter(p, 'L', 80*pi); - + addParameter(p, 'L', 80 * pi); + parse(p, varargin{:}); s = p.Results; @@ -16,21 +16,19 @@ params.L = L; - h = L/N; - % exclude the left boundary point as it is identical to the % right boundary point - x = linspace(h, L, N).'; + x = linspace(L/N, L, N).'; eta1 = min(x/L, 0.1 - x/L); eta2 = 20*(x/L - 0.2).*(0.3 - x/L); eta3 = min(x/L - 0.6, 0.7 - x/L); eta4 = min(x/L - 0.9, 1 - x/L); - u0 = 16*max(0, max(eta1, max(eta2, max(eta3, eta4)))); + u0 = 16*max(0, max([eta1, eta2, eta3, eta4], [], 2)); u0hat = fft(u0); - + tspan = [0, 100]; obj = obj@otp.kuramotosivashinsky.KuramotoSivashinskyProblem(tspan, u0hat, params); diff --git a/src/+otp/+kuramotosivashinsky/KuramotoSivashinskyProblem.m b/src/+otp/+kuramotosivashinsky/KuramotoSivashinskyProblem.m index ea45f664..e9c6ccdb 100644 --- a/src/+otp/+kuramotosivashinsky/KuramotoSivashinskyProblem.m +++ b/src/+otp/+kuramotosivashinsky/KuramotoSivashinskyProblem.m @@ -9,7 +9,7 @@ methods function soly = convert2grid(~, soly) - soly = real(ifft(soly)); + soly = abs(ifft(soly)); end @@ -17,27 +17,19 @@ methods (Access = protected) function onSettingsChanged(obj) - - L = obj.Parameters.L; - N = obj.NumVars; - div = L/(2*pi); + k = 2 * pi * [0:(N/2 - 1), 0, (-N/2 + 1):-1].' / obj.Parameters.L; + ik = 1i * k; + k24 = k.^2 - k.^4; - % note that k already has "i" in it - k = (1i*[0:(N/2 - 1), 0, (-N/2 + 1):-1].'/div); - k2 = k.^2; - k4 = k.^4; - k24 = k2 + k4; - - obj.Rhs = otp.Rhs(@(t, u) otp.kuramotosivashinsky.f(t, u, k, k24), ... + obj.Rhs = otp.Rhs(@(t, u) otp.kuramotosivashinsky.f(t, u, ik, k24), ... otp.Rhs.FieldNames.Jacobian, ... - @(t, u) otp.kuramotosivashinsky.jac(t,u, k, k24), ... + @(t, u) otp.kuramotosivashinsky.jac(t,u, ik, k24), ... otp.Rhs.FieldNames.JacobianVectorProduct, ... - @(t, u, v) otp.kuramotosivashinsky.jvp(t, u, v, k, k24), ... + @(t, u, v) otp.kuramotosivashinsky.jvp(t, u, v, ik, k24), ... otp.Rhs.FieldNames.JacobianAdjointVectorProduct, ... - @(t, u, v) otp.kuramotosivashinsky.javp(t, u, v, k, k24)); - + @(t, u, v) otp.kuramotosivashinsky.javp(t, u, v, ik, k24)); end function validateNewState(obj, newTimeSpan, newY0, newParameters) diff --git a/src/+otp/+kuramotosivashinsky/f.m b/src/+otp/+kuramotosivashinsky/f.m index bcd01ce3..2aeb2721 100644 --- a/src/+otp/+kuramotosivashinsky/f.m +++ b/src/+otp/+kuramotosivashinsky/f.m @@ -1,7 +1,5 @@ -function ut = f(~, u, k, k24) +function ut = f(~, u, ik, k24) -u2 = fft(real(ifft(u)).^2); - -ut = -k24.*u - (k/2).*u2; +ut = k24 .* u - (ik/2) .* fft(ifft(u).^2); end diff --git a/src/+otp/+kuramotosivashinsky/jac.m b/src/+otp/+kuramotosivashinsky/jac.m index 08b6a47e..0a0532c2 100644 --- a/src/+otp/+kuramotosivashinsky/jac.m +++ b/src/+otp/+kuramotosivashinsky/jac.m @@ -1,5 +1,6 @@ -function j = jac(~, u, k, k24) +function j = jac(~, u, ik, k24) -j = -diag(k24) - k.*ifft(fft(diag(real(ifft(u)))).').'; +circulantU = toeplitz(u, [u(1); flipud(u(2:end))]); +j = diag(k24) - ik / length(u) .* circulantU; end diff --git a/src/+otp/+kuramotosivashinsky/javp.m b/src/+otp/+kuramotosivashinsky/javp.m index 141f7574..62145033 100644 --- a/src/+otp/+kuramotosivashinsky/javp.m +++ b/src/+otp/+kuramotosivashinsky/javp.m @@ -1,5 +1,5 @@ -function jv = javp(~, u, v, k, k24) +function jv = javp(~, u, v, ik, k24) -jv = -k24.*v - fft(real(ifft(u)).*ifft(conj(k).*v)); +jv = k24 .* v + fft(conj(ifft(u)) .* ifft(ik .* v)); end diff --git a/src/+otp/+kuramotosivashinsky/jvp.m b/src/+otp/+kuramotosivashinsky/jvp.m index 9e01ff41..995bbac3 100644 --- a/src/+otp/+kuramotosivashinsky/jvp.m +++ b/src/+otp/+kuramotosivashinsky/jvp.m @@ -1,5 +1,5 @@ -function jv = jvp(~, u, v, k, k24) +function jv = jvp(~, u, v, ik, k24) -jv = -k24.*v - k.*fft(real(ifft(u)).*ifft(v)); +jv = k24 .* v - ik .* fft(ifft(u) .* ifft(v)); end From f9af33c5177731e7b0ac669582843811424cfea6 Mon Sep 17 00:00:00 2001 From: Andrey A Popov Date: Mon, 8 May 2023 13:22:34 -0500 Subject: [PATCH 15/15] Updated for new format, and validated all derivatives. --- .../+kuramotosivashinsky/+presets/Canonical.m | 28 +++++++++---------- .../KuramotoSivashinskyParameters.m | 7 +++++ .../KuramotoSivashinskyProblem.m | 25 +++-------------- 3 files changed, 25 insertions(+), 35 deletions(-) create mode 100644 src/+otp/+kuramotosivashinsky/KuramotoSivashinskyParameters.m diff --git a/src/+otp/+kuramotosivashinsky/+presets/Canonical.m b/src/+otp/+kuramotosivashinsky/+presets/Canonical.m index 7800a0b6..1117e61f 100644 --- a/src/+otp/+kuramotosivashinsky/+presets/Canonical.m +++ b/src/+otp/+kuramotosivashinsky/+presets/Canonical.m @@ -1,18 +1,17 @@ -% The Kuramoto-Sivashinky equation is a chaotic problem. -% -% In this particular discretization, we are applying a spectral -% method, therefore the boundary conditions will be chosen to be as cyclic -% on the domain [0, L]. Note that this is different from another typical -% domain of [-L, L]. The larger the L, the more interesting the problem is -% but the more points are required to do a good discretization. The current -% canonical implementation with the size, L, and is used in -% -% Kassam, Aly-Khan, and Lloyd N. Trefethen. -% "Fourth-order time-stepping for stiff PDEs." -% SIAM Journal on Scientific Computing 26, no. 4 (2005): 1214-1233. -% - classdef Canonical < otp.kuramotosivashinsky.KuramotoSivashinskyProblem + % The Kuramoto-Sivashinky equation is a chaotic problem. + % + % In this particular discretization, we are applying a spectral + % method, therefore the boundary conditions will be chosen to be as cyclic + % on the domain [0, L]. Note that this is different from another typical + % domain of [-L, L]. The larger the L, the more interesting the problem is + % but the more points are required to do a good discretization. The current + % canonical implementation with the size, L, and is used in + % + % Kassam, Aly-Khan, and Lloyd N. Trefethen. + % "Fourth-order time-stepping for stiff PDEs." + % SIAM Journal on Scientific Computing 26, no. 4 (2005): 1214-1233. + % methods function obj = Canonical(varargin) @@ -26,6 +25,7 @@ N = s.Size; L = s.L; + params = otp.kuramotosivashinsky.KuramotoSivashinskyParameters; params.L = L; % exclude the left boundary point as it is identical to the diff --git a/src/+otp/+kuramotosivashinsky/KuramotoSivashinskyParameters.m b/src/+otp/+kuramotosivashinsky/KuramotoSivashinskyParameters.m new file mode 100644 index 00000000..357f0140 --- /dev/null +++ b/src/+otp/+kuramotosivashinsky/KuramotoSivashinskyParameters.m @@ -0,0 +1,7 @@ +classdef KuramotoSivashinskyParameters + % Parameters for the Kuramoto Sivashinsky problem + properties + % Represents the length-scale of the problem + L + end +end diff --git a/src/+otp/+kuramotosivashinsky/KuramotoSivashinskyProblem.m b/src/+otp/+kuramotosivashinsky/KuramotoSivashinskyProblem.m index e9c6ccdb..39d9c53d 100644 --- a/src/+otp/+kuramotosivashinsky/KuramotoSivashinskyProblem.m +++ b/src/+otp/+kuramotosivashinsky/KuramotoSivashinskyProblem.m @@ -8,11 +8,8 @@ methods function soly = convert2grid(~, soly) - soly = abs(ifft(soly)); - end - end methods (Access = protected) @@ -23,27 +20,13 @@ function onSettingsChanged(obj) ik = 1i * k; k24 = k.^2 - k.^4; - obj.Rhs = otp.Rhs(@(t, u) otp.kuramotosivashinsky.f(t, u, ik, k24), ... - otp.Rhs.FieldNames.Jacobian, ... + obj.RHS = otp.RHS(@(t, u) otp.kuramotosivashinsky.f(t, u, ik, k24), ... + 'Jacobian', ... @(t, u) otp.kuramotosivashinsky.jac(t,u, ik, k24), ... - otp.Rhs.FieldNames.JacobianVectorProduct, ... + 'JacobianVectorProduct', ... @(t, u, v) otp.kuramotosivashinsky.jvp(t, u, v, ik, k24), ... - otp.Rhs.FieldNames.JacobianAdjointVectorProduct, ... + 'JacobianAdjointVectorProduct', ... @(t, u, v) otp.kuramotosivashinsky.javp(t, u, v, ik, k24)); end - - function validateNewState(obj, newTimeSpan, newY0, newParameters) - - validateNewState@otp.Problem(obj, ... - newTimeSpan, newY0, newParameters) - - if mod(numel(newY0), 2) ~= 0 - error('The problem size has to be an even integer.'); - end - - otp.utils.StructParser(newParameters) ... - .checkField('L', 'scalar', 'finite', 'positive'); - - end end end