Mercurial > repos > public > sbplib
changeset 398:fccd746d8573 feature/SBPinTimeClarity
Refactor for clarity
author | Jonatan Werpers <jonatan@werpers.com> |
---|---|
date | Thu, 02 Feb 2017 09:57:54 +0100 |
parents | 4fd5bfe5d0bb |
children | c92d2f8319c2 |
files | +time/SBPInTime.m |
diffstat | 1 files changed, 42 insertions(+), 73 deletions(-) [+] |
line wrap: on
line diff
diff -r 4fd5bfe5d0bb -r fccd746d8573 +time/SBPInTime.m --- a/+time/SBPInTime.m Thu Feb 02 09:47:21 2017 +0100 +++ b/+time/SBPInTime.m Thu Feb 02 09:57:54 2017 +0100 @@ -1,89 +1,83 @@ classdef SBPInTime < time.Timestepper % The SBP in time method. % Implemented for v_t = A*v + f(t) - % k_local -- time-step - % Nblock -- number of points in each block - % nodes -- points such that t_n + nodes are the points in block n. - % Each "step" takes one block step and thus advances - % k = k_local*(Nblock-1) in time. - % M -- matrix used in every solve. - % [L,U,P,Q] = lu(M); + % + % Each "step" takes one block step and thus advances + % k = k_local*(blockSize-1) in time. properties - M - L - U - P - Q + M % System matrix + L,U,P,Q % LU factorization of M A Et_r penalty f - k_local - k + k_local % step size within a block + k % Time size of a block k/(blockSize-1) = k_local t v m n - Nblock + blockSize % number of points in each block order nodes end methods - function obj = SBPInTime(A, f, k, order, Nblock, t0, v0, TYPE) - default_arg('TYPE','equidistant'); - default_arg('Nblock',time.SBPInTime.smallestBlockSize(order,TYPE)); - + function obj = SBPInTime(A, f, k, t0, v0, TYPE, order, blockSize) + default_arg('TYPE','minimal'); + default_arg('order', 8); + default_arg('blockSize',time.SBPInTime.smallestBlockSize(order,TYPE)); + obj.A = A; obj.f = f; - obj.k_local = k; - obj.k = k*(Nblock-1); - obj.Nblock = Nblock; + obj.k_local = k/(blockSize-1); + obj.k = k; + obj.blockSize = blockSize; obj.t = t0; obj.m = length(v0); obj.n = 0; - + %==== Build the time discretization matrix =====% switch TYPE case 'equidistant' - ops = sbp.D2Standard(Nblock,{0,obj.k},order); + ops = sbp.D2Standard(blockSize,{0,obj.k},order); case 'optimal' - ops = sbp.D1Nonequidistant(Nblock,{0,obj.k},order); + ops = sbp.D1Nonequidistant(blockSize,{0,obj.k},order); case 'minimal' - ops = sbp.D1Nonequidistant(Nblock,{0,obj.k},order,'minimal'); + ops = sbp.D1Nonequidistant(blockSize,{0,obj.k},order,'minimal'); end - + D1 = ops.D1; HI = ops.HI; e_l = ops.e_l; e_r = ops.e_r; obj.nodes = ops.x; - + Ix = speye(size(A)); - It = speye(Nblock,Nblock); - - obj.Et_r = kron(e_r,Ix); - + It = speye(blockSize,blockSize); + + obj.Et_r = kron(e_r,Ix); + % Time derivative + penalty tau = 1; - Mt = D1 + tau*HI*(e_l*e_l'); - + Mt = D1 + tau*HI*(e_l*e_l'); + % penalty to impose "data" penalty = tau*HI*e_l; obj.penalty = kron(penalty,Ix); - + Mx = kron(It,A); - Mt = kron(Mt,Ix); + Mt = kron(Mt,Ix); obj.M = Mt - Mx; %==============================================% - + % LU factorization [obj.L,obj.U,obj.P,obj.Q] = lu(obj.M); - + % Pretend that the initial condition is the last level % of a previous step. obj.v = obj.Et_r * v0; - + end function [v,t] = getV(obj) @@ -93,39 +87,20 @@ function obj = step(obj) obj.v = time.sbp.sbpintime(obj.v, obj.t, obj.nodes,... - obj.penalty, obj.f, obj.Nblock,... + obj.penalty, obj.f, obj.blockSize,... obj.Et_r,... obj.L, obj.U, obj.P, obj.Q); obj.t = obj.t + obj.k; - obj.n = obj.n + obj.Nblock-1; + obj.n = obj.n + 1; end end - - + methods(Static) - - % - function [k,numberOfBlocks] = alignedTimeStep(k,Tend,Nblock) - - % input k is the desired time-step - % Nblock is the number of points per block. - - % Make sure that we reach the final time by advancing - % an integer number of blocks - kblock = (Nblock-1)*k; - numberOfBlocks = ceil(Tend/kblock); - kblock = Tend/(numberOfBlocks); - - % Corrected time step - k = kblock/(Nblock-1); - - end - function N = smallestBlockSize(order,TYPE) default_arg('TYPE','equidistant') - + switch TYPE - + case 'equidistant' switch order case 2 @@ -143,9 +118,9 @@ otherwise error('Operator does not exist'); end - + case 'optimal' - + switch order case 4 N = 8; @@ -160,9 +135,9 @@ otherwise error('Operator does not exist'); end - + case 'minimal' - + switch order case 4 N = 6; @@ -177,13 +152,7 @@ otherwise error('Operator does not exist'); end - end - end - end - - - end