Mercurial > repos > public > sbplib
changeset 404:d6d27fdc342a
Merged in feature/SBPinTimeClarity (pull request #5)
Feature/sbpintimeclarity
Approved-by: Martin Almquist
| author | Martin Almquist <martin.almquist@it.uu.se> |
|---|---|
| date | Thu, 02 Feb 2017 14:46:36 +0000 |
| parents | 18525f1bb941 (current diff) 0ae72c67b03f (diff) |
| children | 4d9d8064e58b 16a76db8c279 |
| files | |
| diffstat | 4 files changed, 123 insertions(+), 98 deletions(-) [+] |
line wrap: on
line diff
diff -r 18525f1bb941 -r d6d27fdc342a +sbp/D1Nonequidistant.m --- a/+sbp/D1Nonequidistant.m Thu Jan 26 13:07:51 2017 +0000 +++ b/+sbp/D1Nonequidistant.m Thu Feb 02 14:46:36 2017 +0000 @@ -11,23 +11,23 @@ x % grid borrowing % Struct with borrowing limits for different norm matrices end - + methods function obj = D1Nonequidistant(m,lim,order,option) - + default_arg('option','Accurate'); % 'Accurate' operators are optimized for accuracy % 'Minimal' operators have the smallest possible boundary % closure - + x_l = lim{1}; x_r = lim{2}; L = x_r-x_l; - + switch option - + case {'Accurate','accurate','A'} - + if order == 4 [obj.D1,obj.H,obj.x,obj.h] = ... sbp.implementations.d1_noneq_4(m,L); @@ -46,9 +46,9 @@ else error('Invalid operator order %d.',order); end - + case {'Minimal','minimal','M'} - + if order == 4 [obj.D1,obj.H,obj.x,obj.h] = ... sbp.implementations.d1_noneq_minimal_4(m,L); @@ -67,28 +67,20 @@ else error('Invalid operator order %d.',order); end - + end - + obj.x = obj.x + x_l; - + obj.e_l = sparse(m,1); obj.e_r = sparse(m,1); obj.e_l(1) = 1; obj.e_r(m) = 1; - + obj.HI = inv(obj.H); obj.Q = obj.H*obj.D1 - obj.e_r*obj.e_r' + obj.e_l*obj.e_l'; - + obj.borrowing = []; - end end - - end - - - - -
diff -r 18525f1bb941 -r d6d27fdc342a +sbp/D2Standard.m --- a/+sbp/D2Standard.m Thu Jan 26 13:07:51 2017 +0000 +++ b/+sbp/D2Standard.m Thu Feb 02 14:46:36 2017 +0000 @@ -14,7 +14,7 @@ h % Step size x % grid borrowing % Struct with borrowing limits for different norm matrices - + end methods @@ -63,11 +63,8 @@ end obj.m = m; - end end - - end
diff -r 18525f1bb941 -r d6d27fdc342a +time/SBPInTime.m --- a/+time/SBPInTime.m Thu Jan 26 13:07:51 2017 +0000 +++ b/+time/SBPInTime.m Thu Feb 02 14:46:36 2017 +0000 @@ -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
diff -r 18525f1bb941 -r d6d27fdc342a +time/SBPInTimeSecondOrderForm.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+time/SBPInTimeSecondOrderForm.m Thu Feb 02 14:46:36 2017 +0000 @@ -0,0 +1,67 @@ +classdef SBPInTimeSecondOrderForm < time.Timestepper + properties + A,B,C + M, f + + n + t + k + + firstOrderTimeStepper + end + + methods + % 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) + default_arg('TYPE', []); + default_arg('order', []); + default_arg('blockSize',[]); + + m = length(v0); + + default_arg('A', sparse(m, m)); + default_arg('B', sparse(m, m)); + default_arg('C', sparse(m, 1)); + + I = speye(m); + O = sparse(m,m); + + obj.M = [ + O, I; + A, B; + ]; + obj.f = @(t)[ + sparse(m,1); + C; + ]; + + w0 = [v0; v0t]; + + obj.k = k; + obj.t = t0; + obj.n = 0; + + obj.firstOrderTimeStepper = time.SBPInTime(obj.M, obj.f, obj.k, obj.t, w0, TYPE, order, blockSize); + end + + function [v,t] = getV(obj) + w = obj.firstOrderTimeStepper.getV(); + v = w(1:end/2); + t = obj.t; + end + + function [vt,t] = getVt(obj) + w = obj.firstOrderTimeStepper.getV(); + vt = w(end/2+1:end); + t = obj.t; + end + + function obj = step(obj) + obj.firstOrderTimeStepper.step(); + obj.t = obj.firstOrderTimeStepper.t; + obj.n = obj.firstOrderTimeStepper.n; + end + end +end