Mercurial > repos > public > sbplib
view +sbp/+implementations/d1_upwind_9.m @ 1337:bf2554f1825d feature/D2_boundary_opt
Add periodic D1 and D2 operators for orders 8,10,12
author | Vidar Stiernström <vidar.stiernstrom@it.uu.se> |
---|---|
date | Fri, 13 May 2022 13:28:10 +0200 |
parents | f7ac3cd6eeaa |
children |
line wrap: on
line source
function [H, HI, Dp, Dm, e_1, e_m] = d1_upwind_9(m,h) BP = 8; if(m<2*BP) error(['Operator requires at least ' num2str(2*BP) ' grid points']); end Hv=ones(m,1); Hv(1:8) = [0.1070017e7/0.3628800e7; 0.5537111e7/0.3628800e7; 0.103613e6/0.403200e6; 0.261115e6/0.145152e6; 0.298951e6/0.725760e6; 0.515677e6/0.403200e6; 0.3349879e7/0.3628800e7; 0.3662753e7/0.3628800e7;]; Hv(m-7:m)=rot90(Hv(1:8),2); Hv=Hv*h; H = spdiag(Hv,0); HI = spdiag(1./Hv,0); q_diags = [-4 -3 -2 -1 0 1 2 3 4 5]; q_stencil = [1/504 -1/42 +1/7 -2/3 -1/5 +1 -1/3 +2/21 -1/56 +1/630]; Qp = stripeMatrix(q_stencil, q_diags,m); Q_U =[ -0.5561e4/0.47263920e8 0.4186300102421e13/0.6193464076800e13 -0.377895002003e12/0.5806372572000e13 -0.16485548951749e14/0.111482353382400e15 -0.113245973003e12/0.3716078446080e13 0.355360297339e12/0.4645098057600e13 0.321012170669e12/0.55741176691200e14 -0.388397049437e12/0.26543417472000e14; -0.4178798062421e13/0.6193464076800e13 -0.493793e6/0.141791760e9 0.725405227507e12/0.2413037952000e13 0.3904159533697e13/0.9290196115200e13 0.2483046570341e13/0.13935294172800e14 -0.4336328670953e13/0.18580392230400e14 -0.1258688487061e13/0.37160784460800e14 0.12931584852209e14/0.278705883456000e15; 0.363359390003e12/0.5806372572000e13 -0.7539548734577e13/0.26543417472000e14 -0.69332623e8/0.2977626960e10 0.9994352248429e13/0.18580392230400e14 -0.8195655811631e13/0.18580392230400e14 0.7361486640463e13/0.61934640768000e14 0.5539855071347e13/0.92901961152000e14 -0.12898722943e11/0.422281641600e12; 0.16773595838149e14/0.111482353382400e15 -0.372477950627e12/0.844563283200e12 -0.8659050093229e13/0.18580392230400e14 -0.207799621e9/0.2977626960e10 0.1734921317461e13/0.2477385630720e13 0.2530020015841e13/0.18580392230400e14 0.441856623253e12/0.13935294172800e14 -0.115132773073e12/0.2654341747200e13; 0.108449122763e12/0.3716078446080e13 -0.2283566671541e13/0.13935294172800e14 0.6976424333231e13/0.18580392230400e14 -0.440819477447e12/0.825795210240e12 -0.55386253e8/0.425375280e9 0.2479572560009e13/0.3716078446080e13 -0.40258468963e11/0.120651897600e12 0.11808221047099e14/0.111482353382400e15; -0.32231128289e11/0.422281641600e12 0.4244793299753e13/0.18580392230400e14 -0.5173673584463e13/0.61934640768000e14 -0.4848139955041e13/0.18580392230400e14 -0.1506045711689e13/0.3716078446080e13 -0.526653889e9/0.2977626960e10 0.36411368691307e14/0.37160784460800e14 -0.825434105779e12/0.2903186286000e13; -0.316459841069e12/0.55741176691200e14 0.1277069729941e13/0.37160784460800e14 -0.6499182375347e13/0.92901961152000e14 0.355606625147e12/0.13935294172800e14 0.1519272420551e13/0.9290196115200e13 -0.2240079855137e13/0.3378253132800e13 -0.584765899e9/0.2977626960e10 0.2301241355533e13/0.2382101568000e13; 0.387779289437e12/0.26543417472000e14 -0.12908508708209e14/0.278705883456000e15 0.147710908133e12/0.4645098057600e13 0.534025841911e12/0.18580392230400e14 -0.4119981443899e13/0.111482353382400e15 0.279819152779e12/0.2903186286000e13 -0.1510324515533e13/0.2382101568000e13 -0.85017967e8/0.425375280e9; ]; Qp(1:8,1:8)=Q_U; Qp(m-7:m,m-7:m)=rot90(Q_U,2)'; %%% This is different from standard SBP Qm=-Qp'; e_1=sparse(m,1);e_1(1)=1; e_m=sparse(m,1);e_m(m)=1; Dp=HI*(Qp-1/2*(e_1*e_1')+1/2*(e_m*e_m')) ; Dm=HI*(Qm-1/2*(e_1*e_1')+1/2*(e_m*e_m')) ; end