Mercurial > repos > public > sbplib
changeset 375:55c14f7321a3 feature/hypsyst
Meged with default
| author | Ylva Rydin <ylva.rydin@telia.com> |
|---|---|
| date | Thu, 26 Jan 2017 14:03:20 +0100 |
| parents | 0fd6561964b0 (current diff) f1289f3b86b8 (diff) |
| children | b0e95c81943e |
| files | |
| diffstat | 6 files changed, 299 insertions(+), 0 deletions(-) [+] |
line wrap: on
line diff
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+sbp/+implementations/d2_variable_2.m Thu Jan 26 14:03:20 2017 +0100 @@ -0,0 +1,59 @@ +function [H, HI, D1, D2, e_l, e_r, d1_l, d1_r] = d2_variable_2(m,h) + + BP = 1; + if(m<2*BP) + error(['Operator requires at least ' num2str(2*BP) ' grid points']); + end + + % Norm + Hv = ones(m,1); + Hv(1) = 1/2; + Hv(m:m) = 1/2; + Hv = h*Hv; + H = spdiag(Hv, 0); + HI = spdiag(1./Hv, 0); + + + % Boundary operators + e_l = sparse(m,1); + e_l(1) = 1; + e_r = rot90(e_l, 2); + + d1_l = sparse(m,1); + d1_l(1:3) = 1/h*[-3/2 2 -1/2]; + d1_r = -rot90(d1_l, 2); + + % D1 operator + diags = -1:1; + stencil = [-1/2 0 1/2]; + D1 = stripeMatrix(stencil, diags, m); + + D1(1,1)=-1;D1(1,2)=1;D1(m,m-1)=-1;D1(m,m)=1; + D1(m,m-1)=-1;D1(m,m)=1; + D1=D1/h; + %Q=H*D1 + 1/2*(e_1*e_1') - 1/2*(e_m*e_m'); + + + M=sparse(m,m); + + scheme_width = 3; + scheme_radius = (scheme_width-1)/2; + r = (1+scheme_radius):(m-scheme_radius); + + function D2 = D2_fun(c) + + Mm1 = -c(r-1)/2 - c(r)/2; + M0 = c(r-1)/2 + c(r) + c(r+1)/2; + Mp1 = -c(r)/2 - c(r+1)/2; + + M(r,:) = spdiags([Mm1 M0 Mp1],0:2*scheme_radius,length(r),m); + + + M(1:2,1:2)=[c(1)/2 + c(2)/2 -c(1)/2 - c(2)/2; -c(1)/2 - c(2)/2 c(1)/2 + c(2) + c(3)/2;]; + M(m-1:m,m-1:m)=[c(m-2)/2 + c(m-1) + c(m)/2 -c(m-1)/2 - c(m)/2; -c(m-1)/2 - c(m)/2 c(m-1)/2 + c(m)/2;]; + M=M/h; + + D2=HI*(-M-c(1)*e_l*d1_l'+c(m)*e_r*d1_r'); + end + D2 = @D2_fun; +end \ No newline at end of file
--- a/+sbp/D2Variable.m Thu Jan 26 13:57:03 2017 +0100 +++ b/+sbp/D2Variable.m Thu Jan 26 14:03:20 2017 +0100 @@ -31,6 +31,13 @@ obj.e_r, obj.d1_l, obj.d1_r] = ... sbp.implementations.d2_variable_4(m,obj.h); obj.borrowing.M.S = 0.2505765857; + case 2 + [obj.H, obj.HI, obj.D1, obj.D2, obj.e_l,... + obj.e_r, obj.d1_l, obj.d1_r] = ... + sbp.implementations.d2_variable_2(m,obj.h); + obj.borrowing.M.S = 0.3636363636; + % Borrowing const taken from Virta 2014 + otherwise error('Invalid operator order %d.',order); end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+time/+sbp/sbpintime.m Thu Jan 26 14:03:20 2017 +0100 @@ -0,0 +1,28 @@ +% Takes one time step of size k using the sbp in time method +% starting from v_0 and where +% M v_t = E_r'*v + f(t) +function v = sbpintime(v, t, tvec, penalty, f, Nblock, E_r,... + L, U, P, Q) + +% Pick out last time step +v_r = E_r'*v; + +% Form RHS +RHS = penalty*v_r; + +% Build vector of f-values and add to RHS +if(~isempty(f)) + fvec = []; + for i = 1:Nblock + fvec = [fvec;f(t+tvec(i))]; + end + RHS = RHS + fvec; +end + +% Solve system +%v = M\RHS; +RHS = P*RHS; +RHS = L\RHS; +v = U\RHS; +v = Q*v; +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+time/SBPInTime.m Thu Jan 26 14:03:20 2017 +0100 @@ -0,0 +1,189 @@ +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); + properties + M + L + U + P + Q + A + Et_r + penalty + f + k_local + k + t + v + m + n + Nblock + 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)); + + obj.A = A; + obj.f = f; + obj.k_local = k; + obj.k = k*(Nblock-1); + obj.Nblock = Nblock; + 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); + case 'optimal' + ops = sbp.D1Nonequidistant(Nblock,{0,obj.k},order); + case 'minimal' + ops = sbp.D1Nonequidistant(Nblock,{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); + + % Time derivative + penalty + tau = 1; + 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); + 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) + v = obj.Et_r' * obj.v; + t = obj.t; + end + + function obj = step(obj) + obj.v = time.sbp.sbpintime(obj.v, obj.t, obj.nodes,... + obj.penalty, obj.f, obj.Nblock,... + obj.Et_r,... + obj.L, obj.U, obj.P, obj.Q); + obj.t = obj.t + obj.k; + obj.n = obj.n + obj.Nblock-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 + N = 2; + case 4 + N = 8; + case 6 + N = 12; + case 8 + N = 16; + case 10 + N = 20; + case 12 + N = 24; + otherwise + error('Operator does not exist'); + end + + case 'optimal' + + switch order + case 4 + N = 8; + case 6 + N = 12; + case 8 + N = 16; + case 10 + N = 20; + case 12 + N = 24; + otherwise + error('Operator does not exist'); + end + + case 'minimal' + + switch order + case 4 + N = 6; + case 6 + N = 10; + case 8 + N = 12; + case 10 + N = 16; + case 12 + N = 20; + otherwise + error('Operator does not exist'); + end + + end + + end + + end + + + +end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scaleAxisLim.m Thu Jan 26 14:03:20 2017 +0100 @@ -0,0 +1,13 @@ +function scaleAxisLim(ah, s) + x = ah.XLim; + y = ah.YLim; + + xl = x(2) - x(1); + yl = y(2) - y(1); + + dx = xl*(s-1); + dy = yl*(s-1); + + ah.XLim = x + [-dx/2 dx/2]; + ah.YLim = y + [-dy/2 dy/2]; +end \ No newline at end of file