Mercurial > repos > public > sbplib

changeset 1113:47e86b5270ad feature/timesteppers

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Change name of property k to dt in time.Timestepper
author Jonatan Werpers <jonatan@werpers.com>
date Wed, 10 Apr 2019 22:40:55 +0200
parents 835c8fa456ec
children f2988a63c3aa
files +time/Cdiff.m +time/CdiffImplicit.m +time/CdiffNonlin.m +time/CdiffTimeDep.m +time/Ode45.m +time/SBPInTime.m +time/SBPInTimeImplicitFormulation.m +time/SBPInTimeScaled.m +time/SBPInTimeSecondOrderForm.m +time/SBPInTimeSecondOrderFormImplicit.m +time/Timestepper.m
diffstat 11 files changed, 85 insertions(+), 85 deletions(-) [+]
line wrap: on
line diff
--- a/+time/Cdiff.m	Wed Apr 10 22:22:46 2019 +0200
+++ b/+time/Cdiff.m	Wed Apr 10 22:40:55 2019 +0200
@@ -3,7 +3,7 @@
         A, B, C
         AA, BB, CC
         G
-        k
+        dt
         t
         v, v_prev
         n
@@ -15,10 +15,10 @@
         %   A*v_tt + B*v_t + C*v = G(t)
         %   v(t0) = v0
         %   v_t(t0) = v0t
-        % starting at time t0 with timestep k
+        % starting at time t0 with timestep dt
         % Using
         % A*Dp*Dm*v_n + B*D0*v_n + C*v_n = G(t_n)
-        function obj = Cdiff(A, B, C, G, v0, v0t, k, t0)
+        function obj = Cdiff(A, B, C, G, v0, v0t, dt, t0)
             m = length(v0);
             default_arg('A', speye(m));
             default_arg('B', sparse(m,m));
@@ -31,14 +31,14 @@
             obj.G = G;
 
             % Rewrite as AA*v_(n+1) + BB*v_n + CC*v_(n-1) = G(t_n)
-            obj.AA =    A/k^2 + B/(2*k);
-            obj.BB = -2*A/k^2 + C;
-            obj.CC =    A/k^2 - B/(2*k);
+            obj.AA =    A/dt^2 + B/(2*dt);
+            obj.BB = -2*A/dt^2 + C;
+            obj.CC =    A/dt^2 - B/(2*dt);
 
-            obj.k = k;
+            obj.dt = dt;
             obj.v_prev = v0;
-            obj.v = v0 + k*v0t;
-            obj.t = t0+k;
+            obj.v = v0 + dt*v0t;
+            obj.t = t0+dt;
             obj.n = 1;
         end
 
@@ -48,14 +48,14 @@
         end
 
         function [vt,t] = getVt(obj)
-            vt = (obj.v-obj.v_prev)/obj.k; % Could be improved using u_tt = f(u))
+            vt = (obj.v-obj.v_prev)/obj.dt; % Could be improved using u_tt = f(u))
             t = obj.t;
         end
 
         function E = getEnergy(obj)
             v  = obj.v;
             vp = obj.v_prev;
-            vt = (obj.v - obj.v_prev)/obj.k;
+            vt = (obj.v - obj.v_prev)/obj.dt;
 
             E = vt'*obj.A*vt + v'*obj.C*vp;
         end
@@ -65,7 +65,7 @@
 
             obj.v_prev = obj.v;
             obj.v      = v_next;
-            obj.t = obj.t + obj.k;
+            obj.t = obj.t + obj.dt;
             obj.n = obj.n + 1;
         end
     end
--- a/+time/CdiffImplicit.m	Wed Apr 10 22:22:46 2019 +0200
+++ b/+time/CdiffImplicit.m	Wed Apr 10 22:40:55 2019 +0200
@@ -3,7 +3,7 @@
         A, B, C
         AA, BB, CC
         G
-        k
+        dt
         t
         v, v_prev
         n
@@ -20,7 +20,7 @@
         % starting at time t0 with timestep
         % Using
         % A*Dp*Dm*v_n + B*D0*v_n + C*I0*v_n = G(t_n)
-        function obj = CdiffImplicit(A, B, C, G, v0, v0t, k, t0)
+        function obj = CdiffImplicit(A, B, C, G, v0, v0t, dt, t0)
             m = length(v0);
             default_arg('A', speye(m));
             default_arg('B', sparse(m,m));
@@ -33,9 +33,9 @@
             obj.G = G;
 
             % Rewrite as AA*v_(n+1) + BB*v_n + CC*v_(n-1) = G(t_n)
-            AA =    A/k^2 + B/(2*k) + C/2;
-            BB = -2*A/k^2;
-            CC =    A/k^2 - B/(2*k) + C/2;
+            AA =    A/dt^2 + B/(2*dt) + C/2;
+            BB = -2*A/dt^2;
+            CC =    A/dt^2 - B/(2*dt) + C/2;
 
             obj.AA = AA;
             obj.BB = BB;
@@ -43,7 +43,7 @@
 
             v_prev = v0;
             I = speye(m);
-            v = v0 + k*v0t;
+            v = v0 + dt*v0t;
 
             if ~issparse(A) || ~issparse(B) || ~issparse(C)
                 error('LU factorization with full pivoting only works for sparse matrices.')
@@ -56,8 +56,8 @@
             obj.p = p;
             obj.q = q;
 
-            obj.k = k;
-            obj.t = t0+k;
+            obj.dt = dt;
+            obj.t = t0+dt;
             obj.n = 1;
             obj.v = v;
             obj.v_prev = v_prev;
@@ -69,7 +69,7 @@
         end
 
         function [vt,t] = getVt(obj)
-            vt = (obj.v-obj.v_prev)/obj.k; % Could be improved using u_tt = f(u))
+            vt = (obj.v-obj.v_prev)/obj.dt; % Could be improved using u_tt = f(u))
             t = obj.t;
         end
 
@@ -77,7 +77,7 @@
         function E = getEnergy(obj)
             v  = obj.v;
             vp = obj.v_prev;
-            vt = (obj.v - obj.v_prev)/obj.k;
+            vt = (obj.v - obj.v_prev)/obj.dt;
 
             E = vt'*obj.A*vt + 1/2*(v'*obj.C*v + vp'*obj.C*vp);
         end
@@ -95,7 +95,7 @@
             obj.v(obj.q) = z;
 
             % Update time
-            obj.t = obj.t + obj.k;
+            obj.t = obj.t + obj.dt;
             obj.n = obj.n + 1;
         end
     end
--- a/+time/CdiffNonlin.m	Wed Apr 10 22:22:46 2019 +0200
+++ b/+time/CdiffNonlin.m	Wed Apr 10 22:40:55 2019 +0200
@@ -3,7 +3,7 @@
         D
         E
         S
-        k
+        dt
         t
         v
         v_prev
@@ -12,7 +12,7 @@
 
 
     methods
-        function obj = CdiffNonlin(D, E, S, k, t0,n0, v, v_prev)
+        function obj = CdiffNonlin(D, E, S, dt, t0,n0, v, v_prev)
             m = size(D(v),1);
             default_arg('E',0);
             default_arg('S',0);
@@ -33,7 +33,7 @@
             obj.D = D;
             obj.E = E;
             obj.S = S;
-            obj.k = k;
+            obj.dt = dt;
             obj.t = t0;
             obj.n = n0;
             obj.v = v;
@@ -46,7 +46,7 @@
         end
 
         function [vt,t] = getVt(obj)
-            vt = (obj.v-obj.v_prev)/obj.k; % Could be improved using u_tt = f(u))
+            vt = (obj.v-obj.v_prev)/obj.dt; % Could be improved using u_tt = f(u))
             t = obj.t;
         end
 
@@ -65,12 +65,12 @@
 
 
             %% Calculate matrices need for the timestep
-            % Before optimization:  A =  1/k^2 * I - 1/(2*k)*E;
-            k = obj.k;
+            % Before optimization:  A =  1/dt^2 * I - 1/(2*dt)*E;
+            dt = obj.dt;
 
-            Aj = 1/k^2 * I(j,j) - 1/(2*k)*E(j,j);
-            B =  2/k^2 * I + D;
-            C = -1/k^2 * I - 1/(2*k)*E;
+            Aj = 1/dt^2 * I(j,j) - 1/(2*dt)*E(j,j);
+            B =  2/dt^2 * I + D;
+            C = -1/dt^2 * I - 1/(2*dt)*E;
 
             %% Take the timestep
             v = obj.v;
@@ -81,13 +81,13 @@
 
             % Before optimization:  obj.v = A\b;
 
-            obj.v(i) = k^2*b(i);
+            obj.v(i) = dt^2*b(i);
             obj.v(j) =  Aj\b(j);
 
             obj.v_prev = v;
 
             %% Update state of the timestepper
-            obj.t = obj.t + obj.k;
+            obj.t = obj.t + obj.dt;
             obj.n = obj.n + 1;
         end
     end
--- a/+time/CdiffTimeDep.m	Wed Apr 10 22:22:46 2019 +0200
+++ b/+time/CdiffTimeDep.m	Wed Apr 10 22:40:55 2019 +0200
@@ -3,7 +3,7 @@
         D
         E
         S
-        k
+        dt
         t
         v
         v_prev
@@ -15,8 +15,8 @@
         % Solves u_tt = Du + E(t)u_t + S(t)
         % D, E, S can either all be constants or all be function handles,
         % They can also be omitted by setting them equal to the empty matrix.
-        % CdiffTimeDep(D, E, S, k, t0, n0, v, v_prev)
-        function obj = CdiffTimeDep(D, E, S, k, t0, n0, v, v_prev)
+        % CdiffTimeDep(D, E, S, dt, t0, n0, v, v_prev)
+        function obj = CdiffTimeDep(D, E, S, dt, t0, n0, v, v_prev)
             m = length(v);
             default_arg('E', @(t)sparse(m,m));
             default_arg('S', @(t)sparse(m,1));
@@ -25,7 +25,7 @@
             obj.E = E;
             obj.S = S;
 
-            obj.k = k;
+            obj.dt = dt;
             obj.t = t0;
             obj.n = n0;
             obj.v = v;
@@ -38,13 +38,13 @@
         end
 
         function [vt,t] = getVt(obj)
-            vt = (obj.v-obj.v_prev)/obj.k; % Could be improved using u_tt = f(u))
+            vt = (obj.v-obj.v_prev)/obj.dt; % Could be improved using u_tt = f(u))
             t = obj.t;
         end
 
         function obj = step(obj)
-            [obj.v, obj.v_prev] = time.cdiff.cdiff(obj.v, obj.v_prev, obj.k, obj.D, obj.E(obj.t), obj.S(obj.t));
-            obj.t = obj.t + obj.k;
+            [obj.v, obj.v_prev] = time.cdiff.cdiff(obj.v, obj.v_prev, obj.dt, obj.D, obj.E(obj.t), obj.S(obj.t));
+            obj.t = obj.t + obj.dt;
             obj.n = obj.n + 1;
         end
     end
--- a/+time/Ode45.m	Wed Apr 10 22:22:46 2019 +0200
+++ b/+time/Ode45.m	Wed Apr 10 22:40:55 2019 +0200
@@ -1,7 +1,7 @@
 classdef Ode45 < time.Timestepper
     properties
         F
-        k
+        dt
         t
         w
         m
@@ -15,7 +15,7 @@
 
 
     methods
-        function obj = Ode45(D, E, S, k, t0, v0, v0t)
+        function obj = Ode45(D, E, S, dt, t0, v0, v0t)
             obj.D = D;
             obj.E = E;
             obj.S = S;
@@ -31,7 +31,7 @@
                 obj.C = [zeros(obj.m,1), S];
             end
 
-            obj.k = k;
+            obj.dt = dt;
             obj.t = t0;
             obj.w = [v0; v0t];
 
@@ -49,7 +49,7 @@
         end
 
         function obj = step(obj)
-            [t,w] = ode45(@(t,w)(obj.F(w,t)),[obj.t obj.t+obj.k],obj.w);
+            [t,w] = ode45(@(t,w)(obj.F(w,t)),[obj.t obj.t+obj.dt],obj.w);
 
             obj.t = t(end);
             obj.w = w(end,:)';
@@ -59,8 +59,8 @@
 
 
     methods (Static)
-        function k = getTimeStep(lambda)
-            k = rk4.get_rk4_time_step(lambda);
+        function dt = getTimeStep(lambda)
+            dt = rk4.get_rk4_time_step(lambda);
         end
     end
 
--- a/+time/SBPInTime.m	Wed Apr 10 22:22:46 2019 +0200
+++ b/+time/SBPInTime.m	Wed Apr 10 22:40:55 2019 +0200
@@ -3,7 +3,7 @@
     % Implemented for v_t = A*v + f(t)
     %
     % Each "step" takes one block step and thus advances
-    % k = k_local*(blockSize-1) in time.
+    % dt = k_local*(blockSize-1) in time.
     properties
         M     % System matrix
         L,U,P,Q % LU factorization of M
@@ -12,7 +12,7 @@
         penalty
         f
         k_local % step size within a block
-        k % Time size of a block  k/(blockSize-1) = k_local
+        dt % Time size of a block  dt/(blockSize-1) = k_local
         t
         v
         m
@@ -23,7 +23,7 @@
     end
 
     methods
-        function obj = SBPInTime(A, f, k, t0, v0, TYPE, order, blockSize)
+        function obj = SBPInTime(A, f, dt, t0, v0, TYPE, order, blockSize)
 
             default_arg('TYPE','gauss');
 
@@ -37,8 +37,8 @@
 
             obj.A = A;
             obj.f = f;
-            obj.k_local = k/(blockSize-1);
-            obj.k = k;
+            obj.k_local = dt/(blockSize-1);
+            obj.dt = dt;
             obj.blockSize = blockSize;
             obj.t = t0;
             obj.m = length(v0);
@@ -47,13 +47,13 @@
             %==== Build the time discretization matrix =====%
             switch TYPE
                 case 'equidistant'
-                    ops = sbp.D2Standard(blockSize,{0,obj.k},order);
+                    ops = sbp.D2Standard(blockSize,{0,obj.dt},order);
                 case 'optimal'
-                    ops = sbp.D1Nonequidistant(blockSize,{0,obj.k},order);
+                    ops = sbp.D1Nonequidistant(blockSize,{0,obj.dt},order);
                 case 'minimal'
-                    ops = sbp.D1Nonequidistant(blockSize,{0,obj.k},order,'minimal');
+                    ops = sbp.D1Nonequidistant(blockSize,{0,obj.dt},order,'minimal');
                 case 'gauss'
-                    ops = sbp.D1Gauss(blockSize,{0,obj.k});
+                    ops = sbp.D1Gauss(blockSize,{0,obj.dt});
             end
 
             D1 = ops.D1;
@@ -99,7 +99,7 @@
                               obj.penalty, obj.f, obj.blockSize,...
                               obj.Et_r,...
                               obj.L, obj.U, obj.P, obj.Q);
-            obj.t = obj.t + obj.k;
+            obj.t = obj.t + obj.dt;
             obj.n = obj.n + 1;
         end
     end
--- a/+time/SBPInTimeImplicitFormulation.m	Wed Apr 10 22:22:46 2019 +0200
+++ b/+time/SBPInTimeImplicitFormulation.m	Wed Apr 10 22:40:55 2019 +0200
@@ -5,7 +5,7 @@
         A,B
         f
 
-        k % total time step.
+        dt % total time step.
 
         blockSize % number of points in each block
         N % Number of components
@@ -24,7 +24,7 @@
     end
 
     methods
-        function obj = SBPInTimeImplicitFormulation(A, B, f, k, t0, v0, TYPE, order, blockSize)
+        function obj = SBPInTimeImplicitFormulation(A, B, f, dt, t0, v0, TYPE, order, blockSize)
 
             default_arg('TYPE','gauss');
             default_arg('f',[]);
@@ -46,7 +46,7 @@
                 obj.f = @(t)sparse(length(v0),1);
             end
 
-            obj.k = k;
+            obj.dt = dt;
             obj.blockSize = blockSize;
             obj.N = length(v0);
 
@@ -56,13 +56,13 @@
             %==== Build the time discretization matrix =====%
             switch TYPE
                 case 'equidistant'
-                    ops = sbp.D2Standard(blockSize,{0,obj.k},order);
+                    ops = sbp.D2Standard(blockSize,{0,obj.dt},order);
                 case 'optimal'
-                    ops = sbp.D1Nonequidistant(blockSize,{0,obj.k},order);
+                    ops = sbp.D1Nonequidistant(blockSize,{0,obj.dt},order);
                 case 'minimal'
-                    ops = sbp.D1Nonequidistant(blockSize,{0,obj.k},order,'minimal');
+                    ops = sbp.D1Nonequidistant(blockSize,{0,obj.dt},order,'minimal');
                 case 'gauss'
-                    ops = sbp.D1Gauss(blockSize,{0,obj.k});
+                    ops = sbp.D1Gauss(blockSize,{0,obj.dt});
             end
 
             I = speye(size(A));
@@ -112,7 +112,7 @@
 
             obj.v = obj.e_T'*w;
 
-            obj.t = obj.t + obj.k;
+            obj.t = obj.t + obj.dt;
             obj.n = obj.n + 1;
         end
     end
--- a/+time/SBPInTimeScaled.m	Wed Apr 10 22:22:46 2019 +0200
+++ b/+time/SBPInTimeScaled.m	Wed Apr 10 22:40:55 2019 +0200
@@ -7,7 +7,7 @@
         A,B
         f
 
-        k % total time step.
+        dt % total time step.
 
         blockSize % number of points in each block
         N % Number of components
@@ -29,7 +29,7 @@
     end
 
     methods
-        function obj = SBPInTimeScaled(A, B, f, k, t0, v0, scaling, TYPE, order, blockSize)
+        function obj = SBPInTimeScaled(A, B, f, dt, t0, v0, scaling, TYPE, order, blockSize)
             default_arg('TYPE','gauss');
             default_arg('f',[]);
 
@@ -51,7 +51,7 @@
                 obj.f = @(t)sparse(length(v0),1);
             end
 
-            obj.k = k;
+            obj.dt = dt;
             obj.blockSize = blockSize;
             obj.N = length(v0);
 
@@ -61,13 +61,13 @@
             %==== Build the time discretization matrix =====%
             switch TYPE
                 case 'equidistant'
-                    ops = sbp.D2Standard(blockSize,{0,obj.k},order);
+                    ops = sbp.D2Standard(blockSize,{0,obj.dt},order);
                 case 'optimal'
-                    ops = sbp.D1Nonequidistant(blockSize,{0,obj.k},order);
+                    ops = sbp.D1Nonequidistant(blockSize,{0,obj.dt},order);
                 case 'minimal'
-                    ops = sbp.D1Nonequidistant(blockSize,{0,obj.k},order,'minimal');
+                    ops = sbp.D1Nonequidistant(blockSize,{0,obj.dt},order,'minimal');
                 case 'gauss'
-                    ops = sbp.D1Gauss(blockSize,{0,obj.k});
+                    ops = sbp.D1Gauss(blockSize,{0,obj.dt});
             end
 
             I = speye(size(A));
@@ -121,7 +121,7 @@
 
             obj.vtilde = obj.e_T'*w;
 
-            obj.t = obj.t + obj.k;
+            obj.t = obj.t + obj.dt;
             obj.n = obj.n + 1;
         end
     end
--- a/+time/SBPInTimeSecondOrderForm.m	Wed Apr 10 22:22:46 2019 +0200
+++ b/+time/SBPInTimeSecondOrderForm.m	Wed Apr 10 22:40:55 2019 +0200
@@ -5,7 +5,7 @@
 
         n
         t
-        k
+        dt
 
         firstOrderTimeStepper
     end
@@ -14,7 +14,7 @@
         % Solves u_tt = Au + Bu_t + C
         % A, B can either both be constants or both be function handles,
         % They can also be omitted by setting them equal to the empty matrix.
-        function obj = SBPInTimeSecondOrderForm(A, B, C, k, t0, v0, v0t, TYPE, order, blockSize)
+        function obj = SBPInTimeSecondOrderForm(A, B, C, dt, t0, v0, v0t, TYPE, order, blockSize)
             default_arg('TYPE', []);
             default_arg('order', []);
             default_arg('blockSize',[]);
@@ -39,11 +39,11 @@
 
             w0 = [v0; v0t];
 
-            obj.k = k;
+            obj.dt = dt;
             obj.t = t0;
             obj.n = 0;
 
-            obj.firstOrderTimeStepper = time.SBPInTime(obj.M, obj.f, obj.k, obj.t, w0, TYPE, order, blockSize);
+            obj.firstOrderTimeStepper = time.SBPInTime(obj.M, obj.f, obj.dt, obj.t, w0, TYPE, order, blockSize);
         end
 
         function [v,t] = getV(obj)
--- a/+time/SBPInTimeSecondOrderFormImplicit.m	Wed Apr 10 22:22:46 2019 +0200
+++ b/+time/SBPInTimeSecondOrderFormImplicit.m	Wed Apr 10 22:40:55 2019 +0200
@@ -5,7 +5,7 @@
 
         n
         t
-        k
+        dt
 
         firstOrderTimeStepper
     end
@@ -14,7 +14,7 @@
         % Solves A*u_tt + B*u_t + C*u = f(t)
         % A, B can either both be constants or both be function handles,
         % They can also be omitted by setting them equal to the empty matrix.
-        function obj = SBPInTimeSecondOrderFormImplicit(A, B, C, f, k, t0, v0, v0t, do_scaling, TYPE, order, blockSize)
+        function obj = SBPInTimeSecondOrderFormImplicit(A, B, C, f, dt, t0, v0, v0t, do_scaling, TYPE, order, blockSize)
             default_arg('f', []);
             default_arg('TYPE', []);
             default_arg('order', []);
@@ -53,15 +53,15 @@
 
             w0 = [v0; v0t];
 
-            obj.k = k;
+            obj.dt = dt;
             obj.t = t0;
             obj.n = 0;
 
             if do_scaling
                 scaling = [ones(m,1); sqrt(diag(C))];
-                obj.firstOrderTimeStepper = time.SBPInTimeScaled(obj.AA, obj.BB, obj.ff, obj.k, obj.t, w0, scaling, TYPE, order, blockSize);
+                obj.firstOrderTimeStepper = time.SBPInTimeScaled(obj.AA, obj.BB, obj.ff, obj.dt, obj.t, w0, scaling, TYPE, order, blockSize);
             else
-                obj.firstOrderTimeStepper = time.SBPInTimeImplicitFormulation(obj.AA, obj.BB, obj.ff, obj.k, obj.t, w0, TYPE, order, blockSize);
+                obj.firstOrderTimeStepper = time.SBPInTimeImplicitFormulation(obj.AA, obj.BB, obj.ff, obj.dt, obj.t, w0, TYPE, order, blockSize);
             end
         end
 
--- a/+time/Timestepper.m	Wed Apr 10 22:22:46 2019 +0200
+++ b/+time/Timestepper.m	Wed Apr 10 22:40:55 2019 +0200
@@ -1,7 +1,7 @@
 classdef Timestepper < handle
     properties (Abstract)
         t
-        k % TBD: should this be a method instead?
+        dt % TBD: should this be a method instead?
         n
     end
 
@@ -81,7 +81,7 @@
         end
 
         function evolve_without_progress(obj, tend)
-            while obj.t < tend - obj.k/2
+            while obj.t < tend - obj.dt/2
                 obj.step();
             end
         end
@@ -93,7 +93,7 @@
             steps_since_update = 0;
             last_update = tic();
             s = util.replace_string('','   %d %%',0);
-            while obj.t < tend - obj.k/2
+            while obj.t < tend - obj.dt/2
                 obj.step();
 
                 steps_since_update = steps_since_update + 1;