Mercurial > repos > public > sbplib

changeset 979:7a5e770974ed feature/timesteppers

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Merge with default
author Jonatan Werpers <jonatan@werpers.com>
date Mon, 07 Jan 2019 16:26:00 +0100
parents 34b3d092a4d0 (diff) 1a30dbe99c7c (current diff)
children b89379fb0814
files +sbp/+implementations/intOpAWW_orders_2to2_ratio2to1.m +sbp/+implementations/intOpAWW_orders_2to2_ratio_2to1_accC2F1_accF2C2.m +sbp/+implementations/intOpAWW_orders_2to2_ratio_2to1_accC2F2_accF2C1.m +sbp/+implementations/intOpAWW_orders_4to4_ratio2to1.m +sbp/+implementations/intOpAWW_orders_4to4_ratio_2to1_accC2F2_accF2C3.m +sbp/+implementations/intOpAWW_orders_4to4_ratio_2to1_accC2F3_accF2C2.m +sbp/+implementations/intOpAWW_orders_6to6_ratio2to1.m +sbp/+implementations/intOpAWW_orders_6to6_ratio_2to1_accC2F3_accF2C4.m +sbp/+implementations/intOpAWW_orders_6to6_ratio_2to1_accC2F4_accF2C3.m +sbp/+implementations/intOpAWW_orders_8to8_ratio2to1.m +sbp/+implementations/intOpAWW_orders_8to8_ratio_2to1_accC2F4_accF2C5.m +sbp/+implementations/intOpAWW_orders_8to8_ratio_2to1_accC2F5_accF2C4.m +sbp/InterpAWW.m +sbp/InterpMC.m +scheme/Wave.m +time/+cdiff/cdiff.m +time/Cdiff.m
diffstat 16 files changed, 286 insertions(+), 167 deletions(-) [+]
line wrap: on
line diff
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/+time/+rk/butcherTableau.m	Mon Jan 07 16:26:00 2019 +0100
@@ -0,0 +1,51 @@
+% Returns the coefficients used in a RK method as defined by a Butcher Tableau.
+%
+% @param method - string specifying which Runge-Kutta method to be used.
+% @return s - number of stages
+% @return a - coefficents for intermediate stages
+% @return b - weights for summing stages
+% @return c - time step coefficents for intermediate stages
+function [s,a,b,c] = butcherTableau(method)
+switch method
+    case "tvd-3"
+        % TVD (Total Variational Diminishing)
+        s = 3;
+        a = zeros(s,s-1);
+        a(2,1) = 1;
+        a(3,1) = 1/4; a(3,2) = 1/4;
+        b = [1/6, 1/6, 2/3];
+        c = [0 1 1/2];
+    case "rk4"
+        % Standard RK4
+        s = 4;
+        a = zeros(s,s-1);
+        a(2,1) = 1/2;
+        a(3,1) = 0; a(3,2) = 1/2;
+        a(4,1) = 0; a(4,2) = 0; a(4,3) = 1;
+        b = [1/6 1/3 1/3 1/6];
+        c = [0, 1/2, 1/2, 1];
+    case "rk4-3/8"
+        % 3/8 RK4 (Kuttas method). Lower truncation error, more flops.
+        % Irreducible, unlike standard rk4.
+        s = 4;
+        a = zeros(s,s-1);
+        a(2,1) = 1/3;
+        a(3,1) = -1/3; a(3,2) = 1;
+        a(4,1) = 1; a(4,2) = -1; a(4,3) = 1;
+        b = [1/8 3/8 3/8 1/8];
+        c = [0, 1/3, 2/3, 1];
+    case "rk6"
+        % Runge-Kutta 6 from Alshina07
+        s = 7;
+        a = zeros(s,s-1);
+        a(2,1) = 4/7;
+        a(3,1) = 115/112; a(3,2) = -5/16;
+        a(4,1) = 589/630; a(4,2) = 5/18; a(4,3) = -16/45;
+        a(5,1) = 229/1200 - 29/6000*sqrt(5); a(5,2) = 119/240 - 187/1200*sqrt(5); a(5,3) = -14/75 + 34/375*sqrt(5); a(5,4) = -3/100*sqrt(5);
+        a(6,1) = 71/2400 - 587/12000*sqrt(5); a(6,2) = 187/480 - 391/2400*sqrt(5); a(6,3) = -38/75 + 26/375*sqrt(5); a(6,4) = 27/80 - 3/400*sqrt(5); a(6,5) = (1+sqrt(5))/4;
+        a(7,1) = -49/480 + 43/160*sqrt(5); a(7,2) = -425/96 + 51/32*sqrt(5); a(7,3) = 52/15 - 4/5*sqrt(5); a(7,4) = -27/16 + 3/16*sqrt(5); a(7,5) = 5/4 - 3/4*sqrt(5); a(7,6) = 5/2 - 1/2*sqrt(5);
+        b = [1/12 0 0 0 5/12 5/12 1/12];
+        c = [0, 4/7, 5/7, 6/7, (5-sqrt(5))/10, (5+sqrt(5))/10, 1];
+    otherwise
+        error('Runge-Kutta method %s is not implemented', method)
+end
\ No newline at end of file
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/+time/+rk/get_rk4_time_step.m	Mon Jan 07 16:26:00 2019 +0100
@@ -0,0 +1,21 @@
+% Calculate the size of the largest time step given the largest evalue for a operator with pure imaginary e.values.
+function k = get_rk4_time_step(lambda,l_type)
+    default_arg('l_type','complex')
+
+    rad = abs(lambda);
+    if strcmp(l_type,'real')
+        % Real eigenvalue
+        % kl > -2.7852
+        k = 2.7852/rad;
+
+    elseif strcmp(l_type,'imag')
+        % Imaginary eigenvalue
+        % |kl| < 2.8284
+        k = 2.8284/rad;
+    elseif strcmp(l_type,'complex')
+        % |kl| < 2.5
+        k = 2.5/rad;
+    else
+        error('l_type must be one of ''real'',''imag'' or ''complex''.')
+    end
+end
\ No newline at end of file
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/+time/+rk/rk4_stability.m	Mon Jan 07 16:26:00 2019 +0100
@@ -0,0 +1,58 @@
+function rk_stability()
+    ruku4 = @(z)(abs(1 + z +(1/2)*z.^2 + (1/6)*z.^3 + (1/24)*z.^4));
+    circ  = @(z)(abs(z));
+
+
+    % contour(X,Y,z)
+    ax = [-4 2 -3 3];
+    % hold on
+    fcontour(ruku4,[1,1],[-3, 0.6],[-3.2, 3.2])
+    hold on
+    r = 2.6;
+    fcontour(circ,[r,r],[-3, 0.6],[-3.2, 3.2],'r')
+    hold off
+    % contour(X,Y,z,[1,1],'b')
+    axis(ax)
+    title('4th order Runge-Kutta stability region')
+    xlabel('Re')
+    ylabel('Im')
+    axis equal
+    grid on
+    box on
+    hold off
+    % surf(X,Y,z)
+
+
+    rk4roots()
+end
+
+function fcontour(f,levels,x_lim,y_lim,opt)
+    default_arg('opt','b')
+    x = linspace(x_lim(1),x_lim(2));
+    y = linspace(y_lim(1),y_lim(2));
+    [X,Y] = meshgrid(x,y);
+    mu = X+ 1i*Y;
+
+    z = f(mu);
+
+    contour(X,Y,z,levels,opt)
+
+end
+
+
+function rk4roots()
+    ruku4 = @(z)(abs(1 + z +(1/2)*z.^2 + (1/6)*z.^3 + (1/24)*z.^4));
+    % Roots for real evalues:
+    F = @(x)(abs(ruku4(x))-1);
+    real_x = fzero(F,-3);
+
+    % Roots for imaginary evalues:
+    F = @(x)(abs(ruku4(1i*x))-1);
+    imag_x1 = fzero(F,-3);
+    imag_x2 = fzero(F,3);
+
+
+    fprintf('Real x = %f\n',real_x)
+    fprintf('Imag x = %f\n',imag_x1)
+    fprintf('Imag x = %f\n',imag_x2)
+end
\ No newline at end of file
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/+time/+rk/rungekutta.m	Mon Jan 07 16:26:00 2019 +0100
@@ -0,0 +1,23 @@
+% Takes one time step of size dt using the rungekutta method
+% starting from @arg v and where the function F(v,t) gives the
+% time derivatives. coeffs is a struct holding the RK coefficients
+% for the specific method.
+% Also returns the stage approximations (V) and stage rates (K).
+function [v, V, K] = rungekutta(v, t , dt, F, coeffs)
+    % Compute the intermediate stages k
+    K = zeros(length(v), coeffs.s);
+    V = zeros(length(v), coeffs.s);
+    for i = 1:coeffs.s
+        u = v;
+        for j = 1:i-1
+            u = u + dt*coeffs.a(i,j)*K(:,j);
+        end
+        V(:,i) = u;
+        K(:,i) = F(u,t+coeffs.c(i)*dt);
+    end
+    % Compute the updated solution as a linear combination
+    % of the intermediate stages.
+    for i = 1:coeffs.s
+        v = v + dt*coeffs.b(i)*k(:,i);
+    end
+end
\ No newline at end of file
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/+time/+rk/rungekuttaDiscreteData.m	Mon Jan 07 16:26:00 2019 +0100
@@ -0,0 +1,30 @@
+% Takes one time step of size dt using the rungekutta method
+% starting from @arg v.
+%
+% discreteData contains (a part of) the forcing function, already
+% evaluated on the space-time grid.
+%
+% ODE: dv/dt = F(v,t) + discreteData(:, nt), where nt denotes the current time-point.
+%
+% coeffs is a struct holding the RK coefficients
+% for the specific method.
+% Also returns the stage approximations (V) and stage rates (K).
+function [v, V, K] = rungekuttaDiscreteData(v, t , dt, F, coeffs, discreteData, n)
+    % Compute the intermediate stages k
+    K = zeros(length(v), coeffs.s);
+    V = zeros(length(v), coeffs.s);
+    for i = 1:coeffs.s
+        u = v;
+        for j = 1:i-1
+            u = u + dt*coeffs.a(i,j)*K(:,j);
+        end
+        V(:,i) = u;
+        K(:,i) = F(u,t+coeffs.c(i)*dt);
+        K(:,i) = K(:,i) + discreteData(:, n*coeffs.s + i);
+    end
+    % Compute the updated solution as a linear combination
+    % of the intermediate stages.
+    for i = 1:coeffs.s
+        v = v + dt*coeffs.b(i)*k(:,i);
+    end
+end
\ No newline at end of file
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/+time/+rk/rungekutta_4.m	Mon Jan 07 16:26:00 2019 +0100
@@ -0,0 +1,10 @@
+% Takes one time step of size k using the rungekutta method
+% starting from v_0 and where the function F(v,t) gives the
+% time derivatives.
+function v = rungekutta_4(v, t , k, F)
+    k1 = F(v         ,t      );
+    k2 = F(v+0.5*k*k1,t+0.5*k);
+    k3 = F(v+0.5*k*k2,t+0.5*k);
+    k4 = F(v+    k*k3,t+    k);
+    v = v + (1/6)*(k1+2*(k2+k3)+k4)*k;
+end
\ No newline at end of file
--- a/+time/+rk4/get_rk4_time_step.m	Mon Jan 07 16:15:49 2019 +0100
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,21 +0,0 @@
-% Calculate the size of the largest time step given the largest evalue for a operator with pure imaginary e.values.
-function k = get_rk4_time_step(lambda,l_type)
-    default_arg('l_type','complex')
-
-    rad = abs(lambda);
-    if strcmp(l_type,'real')
-        % Real eigenvalue
-        % kl > -2.7852
-        k = 2.7852/rad;
-
-    elseif strcmp(l_type,'imag')
-        % Imaginary eigenvalue
-        % |kl| < 2.8284
-        k = 2.8284/rad;
-    elseif strcmp(l_type,'complex')
-        % |kl| < 2.5
-        k = 2.5/rad;
-    else
-        error('l_type must be one of ''real'',''imag'' or ''complex''.')
-    end
-end
\ No newline at end of file
--- a/+time/+rk4/rk4_stability.m	Mon Jan 07 16:15:49 2019 +0100
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,58 +0,0 @@
-function rk_stability()
-    ruku4 = @(z)(abs(1 + z +(1/2)*z.^2 + (1/6)*z.^3 + (1/24)*z.^4));
-    circ  = @(z)(abs(z));
-
-
-    % contour(X,Y,z)
-    ax = [-4 2 -3 3];
-    % hold on
-    fcontour(ruku4,[1,1],[-3, 0.6],[-3.2, 3.2])
-    hold on
-    r = 2.6;
-    fcontour(circ,[r,r],[-3, 0.6],[-3.2, 3.2],'r')
-    hold off
-    % contour(X,Y,z,[1,1],'b')
-    axis(ax)
-    title('4th order Runge-Kutta stability region')
-    xlabel('Re')
-    ylabel('Im')
-    axis equal
-    grid on
-    box on
-    hold off
-    % surf(X,Y,z)
-
-
-    rk4roots()
-end
-
-function fcontour(f,levels,x_lim,y_lim,opt)
-    default_arg('opt','b')
-    x = linspace(x_lim(1),x_lim(2));
-    y = linspace(y_lim(1),y_lim(2));
-    [X,Y] = meshgrid(x,y);
-    mu = X+ 1i*Y;
-
-    z = f(mu);
-
-    contour(X,Y,z,levels,opt)
-
-end
-
-
-function rk4roots()
-    ruku4 = @(z)(abs(1 + z +(1/2)*z.^2 + (1/6)*z.^3 + (1/24)*z.^4));
-    % Roots for real evalues:
-    F = @(x)(abs(ruku4(x))-1);
-    real_x = fzero(F,-3);
-
-    % Roots for imaginary evalues:
-    F = @(x)(abs(ruku4(1i*x))-1);
-    imag_x1 = fzero(F,-3);
-    imag_x2 = fzero(F,3);
-
-
-    fprintf('Real x = %f\n',real_x)
-    fprintf('Imag x = %f\n',imag_x1)
-    fprintf('Imag x = %f\n',imag_x2)
-end
\ No newline at end of file
--- a/+time/+rk4/rungekutta_4.m	Mon Jan 07 16:15:49 2019 +0100
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,10 +0,0 @@
-% Takes one time step of size k using the rungekutta method
-% starting from v_0 and where the function F(v,t) gives the
-% time derivatives.
-function v = rungekutta_4(v, t , k, F)
-    k1 = F(v         ,t      );
-    k2 = F(v+0.5*k*k1,t+0.5*k);
-    k3 = F(v+0.5*k*k2,t+0.5*k);
-    k4 = F(v+    k*k3,t+    k);
-    v = v + (1/6)*(k1+2*(k2+k3)+k4)*k;
-end
\ No newline at end of file
--- a/+time/+rk4/rungekutta_6.m	Mon Jan 07 16:15:49 2019 +0100
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,31 +0,0 @@
-% Takes one time step of size k using the rungekutta method
-% starting from v_0 and where the function F(v,t) gives the
-% time derivatives.
-function v = rungekutta_6(v, t , k, F)
-    s = 7
-    k = zeros(length(v),s)
-    a = zeros(7,6);
-    c = [0, 4/7, 5/7, 6/7, (5-sqrt(5))/10, (5+sqrt(5))/10, 1];
-    b = [1/12, 0, 0, 0, 5/12, 5/12, 1/12];
-    a = [
-        0,                           0,                          0,                       0,                     0,                 0;
-        4/7,                         0,                          0,                       0,                     0,                 0;
-        115/112,                     -5/16,                      0,                       0,                     0,                 0;
-        589/630,                     5/18,                       -16/45,                  0,                     0,                 0;
-        229/1200 - 29/6000*sqrt(5),  119/240 - 187/1200*sqrt(5), -14/75 + 34/375*sqrt(5), -3/100*sqrt(5),        0,                 0;
-        71/2400 - 587/12000*sqrt(5), 187/480 - 391/2400*sqrt(5), -38/75 + 26/375*sqrt(5), 27/80 - 3/400*sqrt(5), (1+sqrt(5))/4,     0;
-        -49/480 + 43/160*sqrt(5),    -425/96 + 51/32*sqrt(5),    52/15 - 4/5*sqrt(5),     -27/16 + 3/16*sqrt(5), 5/4 - 3/4*sqrt(5), 5/2 - 1/2*sqrt(5);
-    ]
-
-    for i = 1:s
-        u = v
-        for j = 1: i-1
-            u = u + h*a(i,j) * k(:,j)
-        end
-        k(:,i) = F(t+c(i)*k,u)
-    end
-
-    for i = 1:s
-        v = v + k*b(i)*k(:,i)
-    end
-end
--- a/+time/Rk4SecondOrderNonlin.m	Mon Jan 07 16:15:49 2019 +0100
+++ b/+time/Rk4SecondOrderNonlin.m	Mon Jan 07 16:26:00 2019 +0100
@@ -61,7 +61,7 @@
         end
 
         function obj = step(obj)
-            obj.w = time.rk4.rungekutta_4(obj.w, obj.t, obj.k, obj.F);
+            obj.w = time.rk.rungekutta_4(obj.w, obj.t, obj.k, obj.F);
             obj.t = obj.t + obj.k;
             obj.n = obj.n + 1;
         end
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/+time/Rungekutta.m	Mon Jan 07 16:26:00 2019 +0100
@@ -0,0 +1,87 @@
+classdef Rungekutta < time.Timestepper
+    properties
+        F       % RHS of the ODE
+        dt      % Time step
+        t       % Time point
+        v       % Solution vector
+        n       % Time level
+        scheme  % The scheme used for the time stepping, e.g rk4, rk6 etc.
+        coeffs  % Butcher tableau coefficients
+        V       % All stage approximations in most recent time step
+        K       % All stage rates in most recent time step
+    end
+
+
+    methods
+        % Timesteps v_t = F(v,t), using the specified RK method from t = t0 with
+        % timestep dt and initial conditions v = v0
+        function obj = Rungekutta(F, dt, t0, v0, method, discreteData)
+            default_arg('method',"rk4");
+            default_arg('discreteData', []);
+            obj.F = F;
+            obj.dt = dt;
+            obj.t = t0;
+            obj.v = v0;
+            obj.n = 0;
+
+            % Extract the coefficients for the specified method
+            % used for the RK updates from the Butcher tableua.
+            [s,a,b,c] = time.rk.butcherTableau(method);
+            obj.coeffs = struct('s',s,'a',a,'b',b,'c',c);
+
+            if isempty(discreteData)
+                % TODO: method "rk4" is also implemented in the butcher tableau, but the rungekutta_4.m implementation
+                % might be slightly more efficient. Need to do some profiling before deciding whether or not to keep it.
+                if (method == "rk4")
+                    obj.scheme = @(v,t,n) time.rk.rungekutta_4(v ,t, dt, F);
+                else
+                    obj.scheme = @(v,t,n) time.rk.rungekutta(v, t, dt, F, obj.coeffs);
+                end
+            else
+                obj.scheme = @(v,t,n) time.rk.rungekuttaDiscreteData(v, t, dt, F, obj.coeffs, discreteData, n);
+            end
+        end
+
+        % v: Current solution
+        % t: Current time
+        % V: All stage approximations in most recent time step
+        % K: All stage rates in most recent time step
+        % T: Time points (corresponding to V and K) in most recent time step
+        function [v,t,V,T,K] = getV(obj)
+            v = obj.v;
+            t = obj.t;
+            V = obj.V;
+            K = obj.K;
+            T = obj.t + obj.dt*obj.coeffs.b;
+        end
+
+        function obj = step(obj)
+            [obj.v, obj.V, obj.K] = obj.scheme(obj.v, obj.t, obj.n);
+            obj.t = obj.t + obj.dt;
+            obj.n = obj.n + 1;
+        end
+
+        % Returns a vector of time points, including substage points,
+        % in the time interval [t0, tEnd].
+        % The time-step obj.dt is assumed to be aligned with [t0, tEnd] already.
+        function tvec = timePoints(obj, t0, tEnd)
+            N = round( (tEnd-t0)/obj.dt );
+            tvec = zeros(N*obj.s, 1);
+            s = obj.coeffs.s;
+            c = obj.coeffs.c;
+            for i = 1:N
+                ind = (i-1)*s+1 : i*s;
+                tvec(ind) = ((i-1) + c')*obj.dt;
+            end
+        end
+
+        % Returns a vector of quadrature weights corresponding to grid points
+        % in time interval [t0, tEnd], substage points included.
+        % The time-step obj.dt is assumed to be aligned with [t0, tEnd] already.
+        function weights = quadWeights(obj, t0, tEnd)
+            N = round( (tEnd-t0)/obj.dt );
+            b = obj.coeffs.b;
+            weights = repmat(b', N, 1);
+        end
+    end
+end
\ No newline at end of file
--- a/+time/Rungekutta4.m	Mon Jan 07 16:15:49 2019 +0100
+++ b/+time/Rungekutta4.m	Mon Jan 07 16:26:00 2019 +0100
@@ -39,7 +39,7 @@
         end
 
         function obj = step(obj)
-            obj.v = time.rk4.rungekutta_4(obj.v, obj.t, obj.k, obj.F);
+            obj.v = time.rk.rungekutta_4(obj.v, obj.t, obj.k, obj.F);
             obj.t = obj.t + obj.k;
             obj.n = obj.n + 1;
         end
--- a/+time/Rungekutta4SecondOrder.m	Mon Jan 07 16:15:49 2019 +0100
+++ b/+time/Rungekutta4SecondOrder.m	Mon Jan 07 16:26:00 2019 +0100
@@ -99,7 +99,7 @@
         end
 
         function obj = step(obj)
-            obj.w = time.rk4.rungekutta_4(obj.w, obj.t, obj.k, obj.F);
+            obj.w = time.rk.rungekutta_4(obj.w, obj.t, obj.k, obj.F);
             obj.t = obj.t + obj.k;
             obj.n = obj.n + 1;
         end
--- a/+time/Rungekutta4proper.m	Mon Jan 07 16:15:49 2019 +0100
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,42 +0,0 @@
-classdef Rungekutta4proper < time.Timestepper
-    properties
-        F
-        k
-        t
-        v
-        m
-        n
-    end
-
-
-    methods
-        % Timesteps v_t = F(v,t), using RK4 fromt t = t0 with timestep k and initial conditions v = v0
-        function obj = Rungekutta4proper(F, k, t0, v0)
-            obj.F = F;
-            obj.k = k;
-            obj.t = t0;
-            obj.v = v0;
-            obj.m = length(v0);
-            obj.n = 0;
-        end
-
-        function [v,t] = getV(obj)
-            v = obj.v;
-            t = obj.t;
-        end
-
-        function obj = step(obj)
-            obj.v = time.rk4.rungekutta_4(obj.v, obj.t, obj.k, obj.F);
-            obj.t = obj.t + obj.k;
-            obj.n = obj.n + 1;
-        end
-    end
-
-
-    methods (Static)
-        function k = getTimeStep(lambda)
-            k = rk4.get_rk4_time_step(lambda);
-        end
-    end
-
-end
\ No newline at end of file
--- a/+time/Timestepper.m	Mon Jan 07 16:15:49 2019 +0100
+++ b/+time/Timestepper.m	Mon Jan 07 16:26:00 2019 +0100
@@ -6,8 +6,9 @@
     end
 
     methods (Abstract)
-         [v,t] = getV(obj)
-         obj = step(obj)
+        % Returns the solution vector v at timestep t.
+        [v,t] = getV(obj)
+        obj = step(obj)
     end