diff --git a/src/+otp/+kortewegdevries/+presets/AscherMcLachlan.m b/src/+otp/+kortewegdevries/+presets/AscherMcLachlan.m new file mode 100644 index 00000000..af8fdf0b --- /dev/null +++ b/src/+otp/+kortewegdevries/+presets/AscherMcLachlan.m @@ -0,0 +1,58 @@ +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 = cos(pi*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 new file mode 100644 index 00000000..3abb92b9 --- /dev/null +++ b/src/+otp/+kortewegdevries/+presets/Canonical.m @@ -0,0 +1,62 @@ +classdef Canonical < otp.kortewegdevries.KortewegDeVriesProblem + methods + function obj = Canonical + N = 250; + + L = 10; + + mesh = linspace(-L, L, N + 1); + mesh = mesh(1:(end-1)); + meshBC = []; + order = 12; + + 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..530950b4 --- /dev/null +++ b/src/+otp/+kortewegdevries/KortewegDeVriesProblem.m @@ -0,0 +1,74 @@ +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 + Invariant + 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 + + [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 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/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/+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 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..812b7efb --- /dev/null +++ b/src/+otp/+utils/+pde/+gfdm/evalAB.m @@ -0,0 +1,25 @@ +function dfd = evalAB(A, B, f, fBC) + +% 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, Bfdiff), 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 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)