Mercurial > repos > public > sbplib
changeset 421:b6a5dc423990
Merged in feature/better_map (pull request #8)
Feature/better map
| author | Jonatan Werpers <jonatan.werpers@it.uu.se> |
|---|---|
| date | Tue, 07 Feb 2017 14:38:51 +0000 |
| parents | 1adfc42bbc1a (diff) 8bd2e36f1f8b (current diff) |
| children | 38173ea263ed a2cb0d4f4a02 5f4540e13f9b 0bc37a25ed88 |
| files | |
| diffstat | 7 files changed, 249 insertions(+), 100 deletions(-) [+] |
line wrap: on
line diff
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+sbp/+implementations/d1_gauss_4.m Tue Feb 07 14:38:51 2017 +0000 @@ -0,0 +1,63 @@ +function [D1,H,x,h,e_l,e_r] = d1_gauss_4(L) + +% L: Domain length +% N: Number of grid points +if(nargin < 2) + L = 1; +end + +N = 4; + +% Quadrature nodes on interval [-1, 1] +x = [ -0.8611363115940526; -0.3399810435848563; 0.3399810435848563; 0.8611363115940526]; + +% Shift nodes to [0,L] +x = (x+1)/2*L; + +% Boundary extrapolation operators +e_l = [1.5267881254572668; -0.8136324494869273; 0.4007615203116504; -0.1139171962819899]; +e_r = flipud(e_l); +e_l = sparse(e_l); +e_r = sparse(e_r); + +%%%% Compute approximate h %%%%%%%%%% +h = L/(N-1); +%%%%%%%%%%%%%%%%%%%%%%%%% + +%%%% Norm matrix on [-1,1] %%%%%%%% +P = sparse(N,N); +P(1,1) = 0.3478548451374539; +P(2,2) = 0.6521451548625461; +P(3,3) = 0.6521451548625461; +P(4,4) = 0.3478548451374539; +%%%%%%%%%%%%%%%%%%%%%%%%% + +%%%% Norm matrix on [0,L] %%%%%%%% +H = P*L/2; +%%%%%%%%%%%%%%%%%%%%%%%%% + +%%%% D1 on [-1,1] %%%%%%%% +D1 = sparse(N,N); +D1(1,1) = -3.3320002363522817; +D1(1,2) = 4.8601544156851962; +D1(1,3) = -2.1087823484951789; +D1(1,4) = 0.5806281691622644; + +D1(2,1) = -0.7575576147992339; +D1(2,2) = -0.3844143922232086; +D1(2,3) = 1.4706702312807167; +D1(2,4) = -0.3286982242582743; + +D1(3,1) = 0.3286982242582743; +D1(3,2) = -1.4706702312807167; +D1(3,3) = 0.3844143922232086; +D1(3,4) = 0.7575576147992339; + +D1(4,1) = -0.5806281691622644; +D1(4,2) = 2.1087823484951789; +D1(4,3) = -4.8601544156851962; +D1(4,4) = 3.3320002363522817; +%%%%%%%%%%%%%%%%%%%%%%%%% + +% D1 on [0,L] +D1 = D1*2/L; \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+sbp/D1Gauss.m Tue Feb 07 14:38:51 2017 +0000 @@ -0,0 +1,41 @@ +classdef D1Gauss < sbp.OpSet + % Diagonal-norm SBP operators based on the Gauss quadrature formula + % with m nodes, which is of degree 2m-1. Hence, The operator D1 is + % accurate of order m. + properties + D1 % SBP operator approximating first derivative + H % Norm matrix + HI % H^-1 + Q % Skew-symmetric matrix + e_l % Left boundary operator + e_r % Right boundary operator + m % Number of grid points. + h % Step size + x % grid + borrowing % Struct with borrowing limits for different norm matrices + end + + methods + function obj = D1Gauss(m,lim) + + x_l = lim{1}; + x_r = lim{2}; + L = x_r-x_l; + + switch m + case 4 + [obj.D1,obj.H,obj.x,obj.h,obj.e_l,obj.e_r] = ... + sbp.implementations.d1_gauss_4(L); + otherwise + error('Invalid number of points: %d.', m); + end + + + obj.x = obj.x + x_l; + 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
--- a/+sbp/D1Nonequidistant.m Tue Feb 07 14:38:51 2017 +0000 +++ b/+sbp/D1Nonequidistant.m Tue Feb 07 14:38:51 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 - - - - -
--- a/+sbp/D2Standard.m Tue Feb 07 14:38:51 2017 +0000 +++ b/+sbp/D2Standard.m Tue Feb 07 14:38:51 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
--- a/+time/SBPInTime.m Tue Feb 07 14:38:51 2017 +0000 +++ b/+time/SBPInTime.m Tue Feb 07 14:38:51 2017 +0000 @@ -1,89 +1,92 @@ 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','gauss'); + + if(strcmp(TYPE,'gauss')) + default_arg('order',4) + default_arg('blockSize',4) + else + default_arg('order', 8); + default_arg('blockSize',time.SBPInTime.smallestBlockSize(order,TYPE)); + end + 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'); + case 'gauss' + ops = sbp.D1Gauss(blockSize,{0,obj.k}); 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; - + obj.v = 1/(e_r'*e_r) * obj.Et_r * v0; + end function [v,t] = getV(obj) @@ -93,39 +96,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); + function N = smallestBlockSize(order,TYPE) + default_arg('TYPE','gauss') - % 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 +127,9 @@ otherwise error('Operator does not exist'); end - + case 'optimal' - + switch order case 4 N = 8; @@ -160,9 +144,9 @@ otherwise error('Operator does not exist'); end - + case 'minimal' - + switch order case 4 N = 6; @@ -177,13 +161,9 @@ otherwise error('Operator does not exist'); end - + case 'gauss' + N = 4; end - end - end - - - end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+time/SBPInTimeSecondOrderForm.m Tue Feb 07 14:38:51 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