From 4029abc5de88f9c3186daea3912080956731823c Mon Sep 17 00:00:00 2001 From: Andrey A Popov Date: Mon, 22 May 2023 17:56:25 -0500 Subject: [PATCH 1/7] Initial commit of GFDM and KdV --- .../+kortewegdevries/+presets/Canonical.m | 62 +++++++++ .../KortewegDeVriesParameters.m | 27 ++++ .../+kortewegdevries/KortewegDeVriesProblem.m | 55 ++++++++ src/+otp/+kortewegdevries/f.m | 16 +++ src/+otp/+utils/+compatibility/pagemfun.m | 41 ++++++ src/+otp/+utils/+compatibility/pagemldivide.m | 21 ++++ src/+otp/+utils/+compatibility/pagemtimes.m | 38 +----- src/+otp/+utils/+pde/+gfdm/evalAB.m | 24 ++++ src/+otp/+utils/+pde/+gfdm/getAB.m | 118 ++++++++++++++++++ src/+otp/+utils/+rbf/gaussian.m | 5 + 10 files changed, 371 insertions(+), 36 deletions(-) create mode 100644 src/+otp/+kortewegdevries/+presets/Canonical.m create mode 100644 src/+otp/+kortewegdevries/KortewegDeVriesParameters.m create mode 100644 src/+otp/+kortewegdevries/KortewegDeVriesProblem.m create mode 100644 src/+otp/+kortewegdevries/f.m create mode 100644 src/+otp/+utils/+compatibility/pagemfun.m create mode 100644 src/+otp/+utils/+compatibility/pagemldivide.m create mode 100644 src/+otp/+utils/+pde/+gfdm/evalAB.m create mode 100644 src/+otp/+utils/+pde/+gfdm/getAB.m create mode 100644 src/+otp/+utils/+rbf/gaussian.m diff --git a/src/+otp/+kortewegdevries/+presets/Canonical.m b/src/+otp/+kortewegdevries/+presets/Canonical.m new file mode 100644 index 00000000..782c4976 --- /dev/null +++ b/src/+otp/+kortewegdevries/+presets/Canonical.m @@ -0,0 +1,62 @@ +classdef Canonical < otp.kortewegdevries.KortewegDeVriesProblem + methods + function obj = Canonical + N = 200; + + L = 10; + + mesh = linspace(-L, L, N + 1); + mesh = mesh(1:(end-1)); + meshBC = []; + order = 10; + + u0 = 6*(sech(mesh).^2); + + BC = @(t, x) 0*x; + + theta = 0; + nu = -1; + alpha = -3; + rho = 0; + + hfun = @(x, y) hfunfull(x, y, L); + + dist = @(x, y) abs(hfun(x, y)); + r = 0.2; + weightfun = @(mesh, meshi) otp.utils.rbf.gaussian(dist(mesh, meshi)/r); + + params = otp.kortewegdevries.KortewegDeVriesParameters; + params.GFDMOrder = order; + params.GFDMMesh = mesh; + params.GFDMMeshBC = meshBC; + params.GFDMHFun = hfun; + params.GFDMWeightFun = weightfun; + params.BCFun = BC; + params.Theta = theta; + params.Nu = nu; + params.Alpha = alpha; + params.Rho = rho; + + timespan = [0, 2]; + + obj = obj@otp.kortewegdevries.KortewegDeVriesProblem(timespan, u0, params); + + end + + + end + + +end + +function h = hfunfull(x, y, L) + +hs = [x - y; 2*L + x - y; -2*L + x - y]; + +[~, I] = min(abs(hs)); + +J = 1:numel(x); + +h = hs(sub2ind(size(hs), I, J)); + +end diff --git a/src/+otp/+kortewegdevries/KortewegDeVriesParameters.m b/src/+otp/+kortewegdevries/KortewegDeVriesParameters.m new file mode 100644 index 00000000..399edfad --- /dev/null +++ b/src/+otp/+kortewegdevries/KortewegDeVriesParameters.m @@ -0,0 +1,27 @@ +classdef KortewegDeVriesParameters + +properties + % + GFDMOrder + % + GFDMMesh + % + GFDMMeshBC + % + GFDMHFun + % + GFDMWeightFun + % + BCFun + % + Theta + % + Nu + % + Alpha + % + Rho +end + + +end \ No newline at end of file diff --git a/src/+otp/+kortewegdevries/KortewegDeVriesProblem.m b/src/+otp/+kortewegdevries/KortewegDeVriesProblem.m new file mode 100644 index 00000000..b28cfc2b --- /dev/null +++ b/src/+otp/+kortewegdevries/KortewegDeVriesProblem.m @@ -0,0 +1,55 @@ +classdef KortewegDeVriesProblem < otp.Problem + + methods + function obj = KortewegDeVriesProblem(timeSpan, y0, parameters) + + obj@otp.Problem('Korteweg de Vries', [], ... + timeSpan, y0, parameters); + + end + end + + properties (SetAccess = private) + DistanceFunction + end + + methods (Access = protected) + function onSettingsChanged(obj) + + mesh = obj.Parameters.GFDMMesh; + meshBC = obj.Parameters.GFDMMeshBC; + order = obj.Parameters.GFDMOrder; + weightfun = obj.Parameters.GFDMWeightFun; + hfun = obj.Parameters.GFDMHFun; + BC = obj.Parameters.BCFun; + theta = obj.Parameters.Theta; + alpha = obj.Parameters.Alpha; + nu = obj.Parameters.Nu; + rho = obj.Parameters.Rho; + + gridPts = size(mesh, 2); + + if obj.NumVars ~= gridPts + warning('OTP:inconsistentNumVars', ... + 'NumVars is %d, but there are %d grid points', ... + obj.NumVars, gridPts); + end + + %Re = obj.Parameters.ReynoldsNumber; + %Ro = obj.Parameters.RossbyNumber; + + [A, B] = otp.utils.pde.gfdm.getAB(mesh, meshBC, order, weightfun, hfun); + + obj.RHS = otp.RHS(@(t, u) ... + otp.kortewegdevries.f(t, u, BC, meshBC, A, B, theta, alpha, nu, rho)); + + %% Distance function + obj.DistanceFunction = @(t, y, i, j) otp.quasigeostrophic.distanceFunction(t, y, i, j, nx, ny); + end + + function sol = internalSolve(obj, varargin) + sol = internalSolve@otp.Problem(obj, 'Solver', otp.utils.Solver.Nonstiff, varargin{:}); + end + + end +end diff --git a/src/+otp/+kortewegdevries/f.m b/src/+otp/+kortewegdevries/f.m new file mode 100644 index 00000000..1b9d3a9d --- /dev/null +++ b/src/+otp/+kortewegdevries/f.m @@ -0,0 +1,16 @@ +function dudt = f(t, u, BC, meshBC, A, B, theta, alpha, nu, rho) + +uBC = BC(t, meshBC); + +dud = otp.utils.pde.gfdm.evalAB(A, B, u.', uBC); + +dudx = dud(1, :).'; +dudxxx = dud(3, :).'; + +du2d = otp.utils.pde.gfdm.evalAB(A, B, (u.^2).', uBC); + +du2dx = du2d(1, :).'; + +dudt = alpha*(theta*du2dx + 2*(1 - theta)*(u.*dudx)) + rho*dudx + nu*dudxxx; + +end diff --git a/src/+otp/+utils/+compatibility/pagemfun.m b/src/+otp/+utils/+compatibility/pagemfun.m new file mode 100644 index 00000000..f7c200ca --- /dev/null +++ b/src/+otp/+utils/+compatibility/pagemfun.m @@ -0,0 +1,41 @@ +function X = pagemfun(A, transpA, B, transpB, mfun) + +if ndims(A) > 3 || ndims(B) > 3 + error('OTP:notImplemented', 'Dimensions greater than 3 are not supported.'); +end + +% repmat A or B depending on the third dimension +N = max(size(A, 3), size(B, 3)); + +if size(A, 3) == 1 + repmat(A, 1, 1, N); +elseif size(A, 3) ~= N + error('OTP:dimensionMismatch', 'The dimensions of A and B do not match.'); +end + +if size(B, 3) == 1 + repmat(B, 1, 1, N); +elseif size(B, 3) ~= N + error('OTP:dimensionMismatch', 'The dimensions of A and B do not match.'); +end + +% transpose A and B as required +switch transpA + case 'transpose' + A = permute(A, [2, 1, 3]); + case 'ctranspose' + A = permute(conj(A), [2, 1, 3]); +end + +switch transpB + case 'transpose' + B = permute(B, [2, 1, 3]); + case 'ctranspose' + B = permute(conj(B), [2, 1, 3]); +end + +for i = N:-1:1 + X(:, :, i) = mfun(A(:, :, i), B(:, :, i)); +end + +end diff --git a/src/+otp/+utils/+compatibility/pagemldivide.m b/src/+otp/+utils/+compatibility/pagemldivide.m new file mode 100644 index 00000000..84910cea --- /dev/null +++ b/src/+otp/+utils/+compatibility/pagemldivide.m @@ -0,0 +1,21 @@ +function X = pagemldivide(A, B, C, D) +% OCTAVE FIX: This function is a naive implementation of the pagemldivide +% functionality for the purposes of supporting octave and legacy matlab +% installations + +if nargin == 2 + transpA = 'none'; + transpB = 'none'; +elseif nargin == 4 + transpA = B; + transpB = D; + B = C; +else + error('OTP:invalidNumberOfArguments', 'The number of arguments has to be 2 or 4.'); +end + +mfun = @mldivide; + +X = otp.utils.compatibility.pagemfun(A, transpA, B, transpB, mfun); + +end diff --git a/src/+otp/+utils/+compatibility/pagemtimes.m b/src/+otp/+utils/+compatibility/pagemtimes.m index b1515485..57a66383 100644 --- a/src/+otp/+utils/+compatibility/pagemtimes.m +++ b/src/+otp/+utils/+compatibility/pagemtimes.m @@ -14,42 +14,8 @@ error('OTP:invalidNumberOfArguments', 'The number of arguments has to be 2 or 4.'); end -if ndims(A) > 3 || ndims(B) > 3 - error('OTP:notImplemented', 'Dimensions greater than 3 are not supported.'); -end - -% repmat A or B depending on the third dimension -N = max(size(A, 3), size(B, 3)); - -if size(A, 3) == 1 - repmat(A, 1, 1, N); -elseif size(A, 3) ~= N - error('OTP:dimensionMismatch', 'The dimensions of A and B do not match.'); -end - -if size(B, 3) == 1 - repmat(B, 1, 1, N); -elseif size(B, 3) ~= N - error('OTP:dimensionMismatch', 'The dimensions of A and B do not match.'); -end +mfun = @mtimes; -% transpose A and B as required -switch transpA - case 'transpose' - A = permute(A, [2, 1, 3]); - case 'ctranspose' - A = permute(conj(A), [2, 1, 3]); -end - -switch transpB - case 'transpose' - B = permute(B, [2, 1, 3]); - case 'ctranspose' - B = permute(conj(B), [2, 1, 3]); -end - -for i = N:-1:1 - X(:, :, i) = A(:, :, i)*B(:, :, i); -end +X = otp.utils.compatibility.pagemfun(A, transpA, B, transpB, mfun); end diff --git a/src/+otp/+utils/+pde/+gfdm/evalAB.m b/src/+otp/+utils/+pde/+gfdm/evalAB.m new file mode 100644 index 00000000..e8000227 --- /dev/null +++ b/src/+otp/+utils/+pde/+gfdm/evalAB.m @@ -0,0 +1,24 @@ +function dfd = evalAB(A, B, f, fBC) + +nmesh = size(f, 2); +nmeshbc = size(fBC, 2); + +nmeshfull = nmesh + nmeshbc; +ndivterms = size(A, 1); + +ffull = [f, fBC]; + +fdiff = reshape(-f + ffull.', 1, nmeshfull, nmesh); +B = reshape(sum(B.*fdiff, 2), ndivterms, 1, nmesh); + +% OCTAVE FIX: find out if the function pagemldivide exists, and +% if it does not, replace it with a compatible function +if exist('pagemldivide', 'builtin') == 0 + pmld = @otp.utils.compatibility.pagemldivide; +else + pmld = @pagemldivide; +end + +dfd = reshape(pmld(A, B), ndivterms, []); + +end diff --git a/src/+otp/+utils/+pde/+gfdm/getAB.m b/src/+otp/+utils/+pde/+gfdm/getAB.m new file mode 100644 index 00000000..be85d40a --- /dev/null +++ b/src/+otp/+utils/+pde/+gfdm/getAB.m @@ -0,0 +1,118 @@ +function [A, B, alldivs] = getAB(mesh, meshBC, order, weightfun, hfun) + +[spatialdim, nmesh] = size(mesh); +[~, nmeshbc] = size(meshBC); +nmeshfull = nmesh + nmeshbc; +meshfull = [mesh, meshBC]; + +alldivs = struct('div', {}, 'mult', {}); + +for ord = 1:order + unique = getallpartial(spatialdim, ord); + + alldivs = [alldivs, unique]; +end + +%weightfun = @(d, r) exp(-0.5*(d/r).^2); + +ndivterms = numel(alldivs); + +W = sparse(nmeshfull, nmesh); + +%% build A and B matrix +A = zeros(ndivterms, ndivterms, nmesh); +B = zeros(ndivterms, nmeshfull, nmesh); +for mi = 1:nmesh +meshi = mesh(:, mi); + +h = hfun(meshfull, meshi); + +w = weightfun(meshfull, meshi); + +W(:, mi) = w.'; + +for i = 1:ndivterms + dvi = alldivs(i).div; + divmi = 1/prod(1:numel(dvi)); + divmi = divmi*(alldivs(i).mult); + + for j = 1:ndivterms + dvj = alldivs(j).div; + divmj = 1/prod(1:numel(dvj)); + divmj = divmj*(alldivs(j).mult); + + b = divmi*(w.^2); + + a = divmj*divmi*(w.^2); + for ii = 1:numel(dvi) + divc = dvi(ii); + a = a.*h(divc, :); + b = b.*h(divc, :); + end + B(i, :, mi) = b; + for jj = 1:numel(dvj) + divc = dvj(jj); + a = a.*h(divc, :); + end + A(i, j, mi) = sum(a); + + end +end + +end + +end + + +function unique = getallpartial(dim, div) + +alldivs = combn(dim, div); + +unique = struct; +unique.div = alldivs{1}; +unique.mult = 1; + + +for i = 2:numel(alldivs) + + wat = alldivs{i}; + + isthere = false; + % check if it is in unique + for j = 1:numel(unique) + wat2 = unique(j).div; + if wat == wat2 + % if it is increase multiplicity + isthere = true; + unique(j).mult = unique(j).mult + 1; + break + end + end + + if ~isthere + unique(end + 1) = struct('div', wat, 'mult', 1); + end + +end + +end + + +function strs = combn(objs, ways) + +strs = {}; + +if ways == 1 + for i = 1:objs + strs{end + 1} = i; + end +else + for i = 1:objs + ends = combn(objs, ways - 1); + for j = 1:numel(ends) + strs{end + 1} = sort([i, ends{j}]); + end + end +end + +end \ No newline at end of file diff --git a/src/+otp/+utils/+rbf/gaussian.m b/src/+otp/+utils/+rbf/gaussian.m new file mode 100644 index 00000000..df4108b9 --- /dev/null +++ b/src/+otp/+utils/+rbf/gaussian.m @@ -0,0 +1,5 @@ +function w = gaussian(k) + +w = exp(-0.5*(k.^2)); + +end From 7cef6b6e6f537a502f264cef2ed4592737f0c5e9 Mon Sep 17 00:00:00 2001 From: Andrey A Popov Date: Mon, 22 May 2023 18:03:25 -0500 Subject: [PATCH 2/7] A fix to make all tests pass --- src/+otp/+kortewegdevries/+presets/Canonical.m | 2 +- src/tests/validateallpresets.m | 3 +++ 2 files changed, 4 insertions(+), 1 deletion(-) diff --git a/src/+otp/+kortewegdevries/+presets/Canonical.m b/src/+otp/+kortewegdevries/+presets/Canonical.m index 782c4976..65478208 100644 --- a/src/+otp/+kortewegdevries/+presets/Canonical.m +++ b/src/+otp/+kortewegdevries/+presets/Canonical.m @@ -10,7 +10,7 @@ meshBC = []; order = 10; - u0 = 6*(sech(mesh).^2); + u0 = 6*(sech(mesh).^2).'; BC = @(t, x) 0*x; diff --git a/src/tests/validateallpresets.m b/src/tests/validateallpresets.m index f7be4a62..f2b57507 100644 --- a/src/tests/validateallpresets.m +++ b/src/tests/validateallpresets.m @@ -45,6 +45,9 @@ if strcmp(problemname, 'nbody') problem.TimeSpan = [0, 0.01]; end + if strcmp(problemname, 'kortewegdevries') + problem.TimeSpan = [0, 0.01]; + end problem.solve(); fprintf('PASS '); assert(true) From ec851fcd2478de874c53973f5f55eec4a3c9d5c2 Mon Sep 17 00:00:00 2001 From: Andrey A Popov Date: Mon, 22 May 2023 19:46:31 -0500 Subject: [PATCH 3/7] reduced order and mesh size in Canonical --- src/+otp/+kortewegdevries/+presets/Canonical.m | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/+otp/+kortewegdevries/+presets/Canonical.m b/src/+otp/+kortewegdevries/+presets/Canonical.m index 65478208..6a93e44c 100644 --- a/src/+otp/+kortewegdevries/+presets/Canonical.m +++ b/src/+otp/+kortewegdevries/+presets/Canonical.m @@ -1,14 +1,14 @@ classdef Canonical < otp.kortewegdevries.KortewegDeVriesProblem methods function obj = Canonical - N = 200; + N = 150; L = 10; mesh = linspace(-L, L, N + 1); mesh = mesh(1:(end-1)); meshBC = []; - order = 10; + order = 8; u0 = 6*(sech(mesh).^2).'; From 43a6b66d948907ea32aa4b1d0cd0a1b54d5a68e4 Mon Sep 17 00:00:00 2001 From: Andrey A Popov Date: Mon, 22 May 2023 19:55:12 -0500 Subject: [PATCH 4/7] further reduce order --- src/+otp/+kortewegdevries/+presets/Canonical.m | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/+otp/+kortewegdevries/+presets/Canonical.m b/src/+otp/+kortewegdevries/+presets/Canonical.m index 6a93e44c..e3b77057 100644 --- a/src/+otp/+kortewegdevries/+presets/Canonical.m +++ b/src/+otp/+kortewegdevries/+presets/Canonical.m @@ -8,7 +8,7 @@ mesh = linspace(-L, L, N + 1); mesh = mesh(1:(end-1)); meshBC = []; - order = 8; + order = 6; u0 = 6*(sech(mesh).^2).'; From 40aa990a77c00d91257aebf2e206be566fbec89e Mon Sep 17 00:00:00 2001 From: Andrey A Popov Date: Tue, 23 May 2023 16:42:41 -0500 Subject: [PATCH 5/7] Added invariant, increased order to preserve it, and added an alternate preset --- .../+presets/AscherMcLachlan.m | 62 +++++++++++++++++++ .../+kortewegdevries/+presets/Canonical.m | 4 +- .../+kortewegdevries/KortewegDeVriesProblem.m | 29 +++++++-- src/+otp/+kortewegdevries/invariant.m | 14 +++++ 4 files changed, 102 insertions(+), 7 deletions(-) create mode 100644 src/+otp/+kortewegdevries/+presets/AscherMcLachlan.m create mode 100644 src/+otp/+kortewegdevries/invariant.m diff --git a/src/+otp/+kortewegdevries/+presets/AscherMcLachlan.m b/src/+otp/+kortewegdevries/+presets/AscherMcLachlan.m new file mode 100644 index 00000000..5382df57 --- /dev/null +++ b/src/+otp/+kortewegdevries/+presets/AscherMcLachlan.m @@ -0,0 +1,62 @@ +classdef AscherMcLachlan < otp.kortewegdevries.KortewegDeVriesProblem + methods + function obj = AscherMcLachlan + N = 200; + + L = 1; + + mesh = linspace(-L, L, N + 1); + mesh = mesh(1:(end-1)); + meshBC = []; + order = 14; + + u0 = cospi(mesh).'; + + BC = @(t, x) 0*x; + + theta = 0; + nu = -(2/3)*(10^-3); + alpha = -3/8; + rho = -1/10; + + hfun = @(x, y) hfunfull(x, y, L); + + dist = @(x, y) abs(hfun(x, y)); + r = 2*2/N; + weightfun = @(mesh, meshi) otp.utils.rbf.gaussian(dist(mesh, meshi)/r); + + params = otp.kortewegdevries.KortewegDeVriesParameters; + params.GFDMOrder = order; + params.GFDMMesh = mesh; + params.GFDMMeshBC = meshBC; + params.GFDMHFun = hfun; + params.GFDMWeightFun = weightfun; + params.BCFun = BC; + params.Theta = theta; + params.Nu = nu; + params.Alpha = alpha; + params.Rho = rho; + + timespan = [0, 5]; + + obj = obj@otp.kortewegdevries.KortewegDeVriesProblem(timespan, u0, params); + + end + + + end + + +end + +function h = hfunfull(x, y, L) + +hs = [x - y; 2*L + x - y; -2*L + x - y]; + +[~, I] = min(abs(hs)); + +J = 1:numel(x); + +h = hs(sub2ind(size(hs), I, J)); + +end diff --git a/src/+otp/+kortewegdevries/+presets/Canonical.m b/src/+otp/+kortewegdevries/+presets/Canonical.m index e3b77057..3abb92b9 100644 --- a/src/+otp/+kortewegdevries/+presets/Canonical.m +++ b/src/+otp/+kortewegdevries/+presets/Canonical.m @@ -1,14 +1,14 @@ classdef Canonical < otp.kortewegdevries.KortewegDeVriesProblem methods function obj = Canonical - N = 150; + N = 250; L = 10; mesh = linspace(-L, L, N + 1); mesh = mesh(1:(end-1)); meshBC = []; - order = 6; + order = 12; u0 = 6*(sech(mesh).^2).'; diff --git a/src/+otp/+kortewegdevries/KortewegDeVriesProblem.m b/src/+otp/+kortewegdevries/KortewegDeVriesProblem.m index b28cfc2b..530950b4 100644 --- a/src/+otp/+kortewegdevries/KortewegDeVriesProblem.m +++ b/src/+otp/+kortewegdevries/KortewegDeVriesProblem.m @@ -11,6 +11,7 @@ properties (SetAccess = private) DistanceFunction + Invariant end methods (Access = protected) @@ -34,22 +35,40 @@ function onSettingsChanged(obj) 'NumVars is %d, but there are %d grid points', ... obj.NumVars, gridPts); end - - %Re = obj.Parameters.ReynoldsNumber; - %Ro = obj.Parameters.RossbyNumber; [A, B] = otp.utils.pde.gfdm.getAB(mesh, meshBC, order, weightfun, hfun); obj.RHS = otp.RHS(@(t, u) ... otp.kortewegdevries.f(t, u, BC, meshBC, A, B, theta, alpha, nu, rho)); - %% Distance function - obj.DistanceFunction = @(t, y, i, j) otp.quasigeostrophic.distanceFunction(t, y, i, j, nx, ny); + %% Distance function and Invariant + dist = @(x, y) abs(hfun(x, y)); + + obj.DistanceFunction = @(t, y, i, j) dist(mesh(i), mesh(j)); + + % find the volume of the full mesh + meshfull = [mesh, meshBC]; + volume = 0; + for i = 1:numel(mesh) + ds = sort(dist(meshfull, mesh(i)), "ascend"); + volume = volume + sum(ds(1:3))/2; + end + + obj.Invariant = @(t, u) ... + otp.kortewegdevries.invariant(t, u, BC, meshBC, A, B, theta, alpha, nu, rho, volume); + end function sol = internalSolve(obj, varargin) sol = internalSolve@otp.Problem(obj, 'Solver', otp.utils.Solver.Nonstiff, varargin{:}); end + + function mov = internalMovie(obj, t, y, varargin) + % BUG This movie is technically broken, and will not work for + % arbitrary meshes. + mov = otp.utils.movie.LineMovie('title', obj.Name, varargin{:}); + mov.record(t, y); + end end end diff --git a/src/+otp/+kortewegdevries/invariant.m b/src/+otp/+kortewegdevries/invariant.m new file mode 100644 index 00000000..19ca664b --- /dev/null +++ b/src/+otp/+kortewegdevries/invariant.m @@ -0,0 +1,14 @@ +function c = invariant(t, u, BC, meshBC, A, B, ~, alpha, nu, rho, volume) + +c = zeros(3, size(u, 2)); + +c(1, :) = volume*mean(u, 1); +c(2, :) = volume*mean(u.^2, 1); + +uBC = BC(t, meshBC); +dud = otp.utils.pde.gfdm.evalAB(A, B, u.', uBC); +dudx = dud(1, :).'; + +c(3, :) = volume*mean((alpha/3)*(u.^3) + (rho/2)*(u.^2) - (nu/2)*(dudx.^2), 1); + +end From f46645eececb57d08e4e986b32696bc9f212b149 Mon Sep 17 00:00:00 2001 From: Andrey A Popov Date: Tue, 23 May 2023 17:03:51 -0500 Subject: [PATCH 6/7] octave fix --- src/+otp/+kortewegdevries/+presets/AscherMcLachlan.m | 6 +----- 1 file changed, 1 insertion(+), 5 deletions(-) diff --git a/src/+otp/+kortewegdevries/+presets/AscherMcLachlan.m b/src/+otp/+kortewegdevries/+presets/AscherMcLachlan.m index 5382df57..af8fdf0b 100644 --- a/src/+otp/+kortewegdevries/+presets/AscherMcLachlan.m +++ b/src/+otp/+kortewegdevries/+presets/AscherMcLachlan.m @@ -10,7 +10,7 @@ meshBC = []; order = 14; - u0 = cospi(mesh).'; + u0 = cos(pi*mesh).'; BC = @(t, x) 0*x; @@ -42,11 +42,7 @@ obj = obj@otp.kortewegdevries.KortewegDeVriesProblem(timespan, u0, params); end - - end - - end function h = hfunfull(x, y, L) From f119e85a09b2d2e0713aaf805aafb3c93eb4f33d Mon Sep 17 00:00:00 2001 From: Andrey A Popov Date: Thu, 25 May 2023 16:49:12 -0500 Subject: [PATCH 7/7] massive speed increase for evaluating gfdm --- src/+otp/+utils/+pde/+gfdm/evalAB.m | 29 +++++++++++++++-------------- 1 file changed, 15 insertions(+), 14 deletions(-) diff --git a/src/+otp/+utils/+pde/+gfdm/evalAB.m b/src/+otp/+utils/+pde/+gfdm/evalAB.m index e8000227..812b7efb 100644 --- a/src/+otp/+utils/+pde/+gfdm/evalAB.m +++ b/src/+otp/+utils/+pde/+gfdm/evalAB.m @@ -1,24 +1,25 @@ function dfd = evalAB(A, B, f, fBC) -nmesh = size(f, 2); -nmeshbc = size(fBC, 2); - -nmeshfull = nmesh + nmeshbc; -ndivterms = size(A, 1); - -ffull = [f, fBC]; - -fdiff = reshape(-f + ffull.', 1, nmeshfull, nmesh); -B = reshape(sum(B.*fdiff, 2), ndivterms, 1, nmesh); - -% OCTAVE FIX: find out if the function pagemldivide exists, and -% if it does not, replace it with a compatible function +% OCTAVE FIX: find out if the functions pagemtimes and pagemldivide exist, and +% if they does not, replace them with a compatible function if exist('pagemldivide', 'builtin') == 0 pmld = @otp.utils.compatibility.pagemldivide; else pmld = @pagemldivide; end +if exist('pagemtimes', 'builtin') == 0 + pmt = @otp.utils.compatibility.pagemtimes; +else + pmt = @pagemtimes; +end + +ndivterms = size(A, 1); + +ffull = [f, fBC]; + +fdiff = permute(-f + ffull.', [3, 1, 2]); +Bfdiff = (pmt(B, 'none', fdiff, 'transpose')); -dfd = reshape(pmld(A, B), ndivterms, []); +dfd = reshape(pmld(A, Bfdiff), ndivterms, []); end