Mercurial > repos > public > sbplib
changeset 325:72468bc9b63f feature/beams
Renamed some operator implementations.
| author | Jonatan Werpers <jonatan@werpers.com> |
|---|---|
| date | Mon, 26 Sep 2016 09:55:16 +0200 |
| parents | c0cbffcf6513 |
| children | b19e142fcae1 |
| files | +sbp/+implementations/d4_lonely_4_min_boundary_points.m +sbp/+implementations/d4_lonely_6_2.m +sbp/+implementations/d4_lonely_6_3.m +sbp/+implementations/d4_lonely_6_min_boundary_points.m +sbp/+implementations/d4_lonely_8_higher_boundary_order.m +sbp/+implementations/d4_lonely_8_min_boundary_points.m +sbp/+implementations/d4_variable_4_min_boundary_points.m +sbp/+implementations/d4_variable_6_2.m +sbp/+implementations/d4_variable_6_3.m +sbp/+implementations/d4_variable_6_min_boundary_points.m +sbp/+implementations/d4_variable_8_higher_boundary_order.m +sbp/+implementations/d4_variable_8_min_boundary_points.m +sbp/D4Lonely.m |
| diffstat | 13 files changed, 455 insertions(+), 455 deletions(-) [+] |
line wrap: on
line diff
diff -r c0cbffcf6513 -r 72468bc9b63f +sbp/+implementations/d4_lonely_4_min_boundary_points.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+sbp/+implementations/d4_lonely_4_min_boundary_points.m Mon Sep 26 09:55:16 2016 +0200 @@ -0,0 +1,81 @@ +function [H, HI, D4, e_l, e_r, M4, d2_l, d2_r, d3_l, d3_r, d1_l, d1_r] = d4_variable_4_min_boundary_points(m,h) + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + %%% 4:de ordn. SBP Finita differens %%% + %%% operatorer framtagna av Mark Carpenter %%% + %%% %%% + %%% H (Normen) %%% + %%% D1=H^(-1)Q (approx f?rsta derivatan) %%% + %%% D2 (approx andra derivatan) %%% + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + %H?r med endast 4 randpunkter + + + BP = 4; + if(m<2*BP) + error(['Operator requires at least ' num2str(2*BP) ' grid points']); + end + + + % Norm + Hv = ones(m,1); + Hv(1:4) = [17/48 59/48 43/48 49/48]; + Hv(m-3:m) = rot90(Hv(1:4),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:4) = 1/h*[-11/6 3 -3/2 1/3]; + d1_r = -rot90(d1_l); + + d2_l = sparse(m,1); + d2_l(1:4) = 1/h^2*[2 -5 4 -1]; + d2_r = rot90(d2_l, 2); + + d3_l = sparse(m,1); + d3_l(1:4) = 1/h^3*[-1 3 -3 1]; + d3_r = -rot90(d3_l, 2); + + + % First derivative + stencil = [1/12 -2/3 0 2/3 -1/12]; + diags = [-1 0 1]; + + Q_U = [ + 0 0.59e2/0.96e2 -0.1e1/0.12e2 -0.1e1/0.32e2; + -0.59e2/0.96e2 0 0.59e2/0.96e2 0; + 0.1e1/0.12e2 -0.59e2/0.96e2 0 0.59e2/0.96e2; + 0.1e1/0.32e2 0 -0.59e2/0.96e2 0; + ]; + + Q = stripeMatrix(stencil, diags, m); + Q(1:4,1:4)=Q_U; + Q(m-3:m,m-3:m) = -rot90(Q_U, 2); + + D1 = HI*(Q - 1/2*e_l*e_l' + 1/2*e_r*e_r'); + + % Fourth derivative + stencil = [-1/6, 2, -13/2, 28/3, -13/2, 2, -1/6]; + diags = -3:3; + M4 = stripeMatrix(stencil, diags, m); + + M4_U=[ + 0.8e1/0.3e1 -0.37e2/0.6e1 0.13e2/0.3e1 -0.5e1/0.6e1; + -0.37e2/0.6e1 0.47e2/0.3e1 -13 0.11e2/0.3e1; + 0.13e2/0.3e1 -13 0.44e2/0.3e1 -0.47e2/0.6e1; + -0.5e1/0.6e1 0.11e2/0.3e1 -0.47e2/0.6e1 0.29e2/0.3e1; + ]; + + + M4(1:4,1:4) = M4_U; + M4(m-3:m,m-3:m) = rot90(M4_U, 2); + M4 = 1/h^3*M4; + + D4=HI*(M4 - e_l*d3_l'+e_r*d3_r' + d1_l*d2_l'-d1_r*d2_r'); +end \ No newline at end of file
diff -r c0cbffcf6513 -r 72468bc9b63f +sbp/+implementations/d4_lonely_6_2.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+sbp/+implementations/d4_lonely_6_2.m Mon Sep 26 09:55:16 2016 +0200 @@ -0,0 +1,76 @@ +function [H, HI, D4, e_l, e_r, M4, d2_l, d2_r, d3_l, d3_r, d1_l, d1_r] = d4_variable_6_2(m,h) + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + %%% 6:te ordn. SBP Finita differens %%% + %%% operatorer med diagonal norm %%% + %%% Extension to variable koeff %%% + %%% %%% + %%% H (Normen) %%% + %%% D1=H^(-1)Q (approx f?rsta derivatan) %%% + %%% D2 (approx andra derivatan) %%% + %%% D2=HI*(R+C*D*S %%% + %%% %%% + %%% R=-D1'*H*C*D1-RR %%% + %%% %%% + %%% RR ?r dissipation) %%% + %%% Dissipationen uppbyggd av D4: %%% + %%% DI=D4*B*H*D4 %%% + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + % H?r med 6 RP ist?llet f?r 8 f?r D4 operatorn, dock samma randderivator + % Denna ?r noggrannare, och har 2a ordningens randdslutning och b?r ge 6te + % ordningens konvergens. Hade dock ingen fri parameter att optimera + + BP = 6; + if(m<2*BP) + error(['Operator requires at least ' num2str(2*BP) ' grid points']); + end + + % Norm + Hv = ones(m,1); + Hv(1:6) = [0.181e3/0.576e3, 0.1343e4/0.960e3, 0.293e3/0.480e3, 0.1811e4/0.1440e4, 0.289e3/0.320e3, 0.65e2/0.64e2]; + Hv(m-5:m) = rot90(Hv(1:6),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:6) = [-0.137e3/0.60e2 5 -5 0.10e2/0.3e1 -0.5e1/0.4e1 0.1e1/0.5e1;]/h; + d1_r = -rot90(d1_l); + + d2_l = sparse(m,1); + d2_l(1:6) = [0.15e2/0.4e1 -0.77e2/0.6e1 0.107e3/0.6e1 -13 0.61e2/0.12e2 -0.5e1/0.6e1;]/h^2; + d2_r = rot90(d2_l, 2); + + d3_l = sparse(m,1); + d3_l(1:6) = [-0.17e2/0.4e1 0.71e2/0.4e1 -0.59e2/0.2e1 0.49e2/0.2e1 -0.41e2/0.4e1 0.7e1/0.4e1;]/h^3; + d3_r = -rot90(d3_l, 2); + + + % Fourth derivative, 1th order accurate at first 8 boundary points (still + % yield 5th order convergence if stable: for example u_tt = -u_xxxx + stencil = [7/240, -2/5, 169/60, -122/15, 91/8, -122/15, 169/60, -2/5, 7/240]; + diags = -4:4; + M4 = stripeMatrix(stencil, diags, m); + + M4_U = [ + 0.1009e4/0.192e3 -0.7657e4/0.480e3 0.9307e4/0.480e3 -0.509e3/0.40e2 0.4621e4/0.960e3 -0.25e2/0.32e2; + -0.7657e4/0.480e3 0.49513e5/0.960e3 -0.4007e4/0.60e2 0.21799e5/0.480e3 -0.8171e4/0.480e3 0.2657e4/0.960e3; + 0.9307e4/0.480e3 -0.4007e4/0.60e2 0.1399e4/0.15e2 -0.2721e4/0.40e2 0.12703e5/0.480e3 -0.521e3/0.120e3; + -0.509e3/0.40e2 0.21799e5/0.480e3 -0.2721e4/0.40e2 0.3349e4/0.60e2 -0.389e3/0.15e2 0.559e3/0.96e2; + 0.4621e4/0.960e3 -0.8171e4/0.480e3 0.12703e5/0.480e3 -0.389e3/0.15e2 0.17857e5/0.960e3 -0.1499e4/0.160e3; + -0.25e2/0.32e2 0.2657e4/0.960e3 -0.521e3/0.120e3 0.559e3/0.96e2 -0.1499e4/0.160e3 0.2225e4/0.192e3; + ]; + + + M4(1:6,1:6) = M4_U; + M4(m-5:m,m-5:m) = rot90(M4_U, 2); + M4 = 1/h^3*M4; + + D4=HI*(M4 - e_l*d3_l'+e_r*d3_r' + d1_l*d2_l'-d1_r*d2_r'); +end
diff -r c0cbffcf6513 -r 72468bc9b63f +sbp/+implementations/d4_lonely_6_3.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+sbp/+implementations/d4_lonely_6_3.m Mon Sep 26 09:55:16 2016 +0200 @@ -0,0 +1,70 @@ +function [H, HI, D4, e_l, e_r, M4, d2_l, d2_r, d3_l, d3_r, d1_l, d1_r] = d4_variable_6_3(m,h) + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + %%% 6:te ordn. SBP Finita differens %%% + %%% operatorer med diagonal norm %%% + %%% Extension to variable koeff %%% + %%% %%% + %%% H (Normen) %%% + %%% D1=H^(-1)Q (approx f?rsta derivatan) %%% + %%% D2 (approx andra derivatan) %%% + %%% D2=HI*(R+C*D*S %%% + %%% %%% + %%% R=-D1'*H*C*D1-RR %%% + %%% %%% + %%% RR ?r dissipation) %%% + %%% Dissipationen uppbyggd av D4: %%% + %%% DI=D4*B*H*D4 %%% + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + % H?r med 7 RP ist?llet f?r 8 f?r D4 operatorn, dock samma randderivator + % Denna ?r noggrannare, och har 2a ordningens randdslutning och b?r ge 6te + % ordningens konvergens. Hade 2 fria parametrar att optimera + + % Norm + Hv = ones(m,1); + Hv(1:7) = [0.414837907e9/0.1191965760e10, 0.475278367e9/0.397321920e9, 0.13872751e8/0.12416310e8, 0.346739027e9/0.595982880e9, 0.560227469e9/0.397321920e9, 0.322971631e9/0.397321920e9, 0.616122491e9/0.595982880e9]; + Hv(m-6:m) = rot90(Hv(1:7),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:6) = [-0.137e3/0.60e2 5 -5 0.10e2/0.3e1 -0.5e1/0.4e1 0.1e1/0.5e1;]/h; + d1_r = -rot90(d1_l); + + d2_l = sparse(m,1); + d2_l(1:6) = [0.15e2/0.4e1 -0.77e2/0.6e1 0.107e3/0.6e1 -13 0.61e2/0.12e2 -0.5e1/0.6e1;]/h^2; + d2_r = rot90(d2_l, 2); + + d3_l = sparse(m,1); + d3_l(1:6) = [-0.17e2/0.4e1 0.71e2/0.4e1 -0.59e2/0.2e1 0.49e2/0.2e1 -0.41e2/0.4e1 0.7e1/0.4e1;]/h^3; + d3_r = -rot90(d3_l, 2); + + + % Fourth derivative, 1th order accurate at first 8 boundary points + stencil = [7/240, -2/5, 169/60, -122/15, 91/8, -122/15, 169/60, -2/5, 7/240]; + diags = -4:4; + M4 = stripeMatrix(stencil, diags, m); + + M4_U = [ + 0.1399708478939e13/0.263487168000e12 -0.13482796013041e14/0.834376032000e12 0.344344095859e12/0.17565811200e11 -0.3166261424681e13/0.250312809600e12 0.1508605165681e13/0.333750412800e12 -0.486270829441e12/0.834376032000e12 -0.221976356359e12/0.5006256192000e13; + -0.13482796013041e14/0.834376032000e12 0.7260475818391e13/0.139062672000e12 -0.27224036353e11/0.406022400e9 0.1847477458951e13/0.41718801600e11 -0.848984558161e12/0.55625068800e11 0.247494925991e12/0.139062672000e12 0.165585445559e12/0.834376032000e12; + 0.344344095859e12/0.17565811200e11 -0.27224036353e11/0.406022400e9 0.2044938640393e13/0.22250027520e11 -0.1071086785417e13/0.16687520640e11 0.502199537033e12/0.22250027520e11 -0.143589154441e12/0.55625068800e11 -0.88181965559e11/0.333750412800e12; + -0.3166261424681e13/0.250312809600e12 0.1847477458951e13/0.41718801600e11 -0.1071086785417e13/0.16687520640e11 0.628860435593e12/0.12515640480e11 -0.73736245829e11/0.3337504128e10 0.195760572271e12/0.41718801600e11 -0.81156046361e11/0.250312809600e12; + 0.1508605165681e13/0.333750412800e12 -0.848984558161e12/0.55625068800e11 0.502199537033e12/0.22250027520e11 -0.73736245829e11/0.3337504128e10 0.76725285869e11/0.4450005504e10 -0.3912429433e10/0.406022400e9 0.53227370659e11/0.17565811200e11; + -0.486270829441e12/0.834376032000e12 0.247494925991e12/0.139062672000e12 -0.143589154441e12/0.55625068800e11 0.195760572271e12/0.41718801600e11 -0.3912429433e10/0.406022400e9 0.1699707221791e13/0.139062672000e12 -0.6959018412841e13/0.834376032000e12; + -0.221976356359e12/0.5006256192000e13 0.165585445559e12/0.834376032000e12 -0.88181965559e11/0.333750412800e12 -0.81156046361e11/0.250312809600e12 0.53227370659e11/0.17565811200e11 -0.6959018412841e13/0.834376032000e12 0.3012195053939e13/0.263487168000e12; + ]; + + M4(1:7,1:7) = M4_U; + M4(m-6:m,m-6:m) = rot90(M4_U, 2); + M4 = 1/h^3*M4; + + D4=HI*(M4 - e_l*d3_l'+e_r*d3_r' + d1_l*d2_l'-d1_r*d2_r'); +end
diff -r c0cbffcf6513 -r 72468bc9b63f +sbp/+implementations/d4_lonely_6_min_boundary_points.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+sbp/+implementations/d4_lonely_6_min_boundary_points.m Mon Sep 26 09:55:16 2016 +0200 @@ -0,0 +1,74 @@ +function [H, HI, D4, e_l, e_r, M4, d2_l, d2_r, d3_l, d3_r, d1_l, d1_r] = d4_variable_6_min_boundary_points(m,h) + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + %%% 6:te ordn. SBP Finita differens %%% + %%% operatorer med diagonal norm %%% + %%% Extension to variable koeff %%% + %%% %%% + %%% H (Normen) %%% + %%% D1=H^(-1)Q (approx f?rsta derivatan) %%% + %%% D2 (approx andra derivatan) %%% + %%% D2=HI*(R+C*D*S %%% + %%% %%% + %%% R=-D1'*H*C*D1-RR %%% + %%% %%% + %%% RR ?r dissipation) %%% + %%% Dissipationen uppbyggd av D4: %%% + %%% DI=D4*B*H*D4 %%% + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + % H?r med 6 RP ist?llet f?r 8 f?r D4 operatorn, dock samma randderivator + + BP = 6; + if(m<2*BP) + error(['Operator requires at least ' num2str(2*BP) ' grid points']); + end + + % Norm + Hv = ones(m,1); + Hv(1:6) = [13649/43200,12013/8640,2711/4320,5359/4320,7877/8640, 43801/43200]; + Hv(m-5:m) = rot90(Hv(1:6),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:5) = [-25/12, 4, -3, 4/3, -1/4]/h; + d1_r = -rot90(d1_l); + + d2_l = sparse(m,1); + d2_l(1:5) = [0.35e2/0.12e2 -0.26e2/0.3e1 0.19e2/0.2e1 -0.14e2/0.3e1 0.11e2/0.12e2;]/h^2; + d2_r = rot90(d2_l, 2); + + d3_l = sparse(m,1); + d3_l(1:5) = [-0.5e1/0.2e1 9 -12 7 -0.3e1/0.2e1;]/h^3; + d3_r = -rot90(d3_l, 2); + + + % Fourth derivative, 1th order accurate at first 8 boundary points (still + % yield 5th order convergence if stable: for example u_tt=-u_xxxx + + stencil = [7/240, -2/5, 169/60, -122/15, 91/8, -122/15, 169/60, -2/5, 7/240]; + diags = -4:4; + M4 = stripeMatrix(stencil, diags, m); + + M4_U=[ + 0.3504379e7/0.907200e6 -0.4613983e7/0.453600e6 0.4260437e7/0.453600e6 -0.418577e6/0.113400e6 0.524579e6/0.907200e6 0.535e3/0.18144e5; + -0.4613983e7/0.453600e6 0.5186159e7/0.181440e6 -0.81121e5/0.2835e4 0.218845e6/0.18144e5 -0.159169e6/0.90720e5 -0.94669e5/0.907200e6; + 0.4260437e7/0.453600e6 -0.81121e5/0.2835e4 0.147695e6/0.4536e4 -0.384457e6/0.22680e5 0.339653e6/0.90720e5 -0.18233e5/0.113400e6; + -0.418577e6/0.113400e6 0.218845e6/0.18144e5 -0.384457e6/0.22680e5 0.65207e5/0.4536e4 -0.22762e5/0.2835e4 0.1181753e7/0.453600e6; + 0.524579e6/0.907200e6 -0.159169e6/0.90720e5 0.339653e6/0.90720e5 -0.22762e5/0.2835e4 0.2006171e7/0.181440e6 -0.3647647e7/0.453600e6; + 0.535e3/0.18144e5 -0.94669e5/0.907200e6 -0.18233e5/0.113400e6 0.1181753e7/0.453600e6 -0.3647647e7/0.453600e6 0.10305271e8/0.907200e6; + ]; + + M4(1:6,1:6) = M4_U; + M4(m-5:m,m-5:m) = rot90(M4_U, 2); + M4 = 1/h^3*M4; + + D4=HI*(M4 - e_l*d3_l'+e_r*d3_r' + d1_l*d2_l'-d1_r*d2_r'); +end
diff -r c0cbffcf6513 -r 72468bc9b63f +sbp/+implementations/d4_lonely_8_higher_boundary_order.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+sbp/+implementations/d4_lonely_8_higher_boundary_order.m Mon Sep 26 09:55:16 2016 +0200 @@ -0,0 +1,75 @@ +function [H, HI, D4, e_l, e_r, M4, d2_l, d2_r, d3_l, d3_r, d1_l, d1_r] = d4_variable_8_higher_boundary_order(m,h) + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + %%% 8:te ordn. SBP Finita differens %%% + %%% operatorer med diagonal norm %%% + %%% %%% + %%% %%% + %%% H (Normen) %%% + %%% D1=H^(-1)Q (approx f?rsta derivatan) %%% + %%% D2 (approx andra derivatan) %%% + %%% D2=HI*(R+C*D*S %%% + %%% %%% + %%% R=-D1'*H*C*D1-RR %%% + %%% %%% + %%% RR ?r dissipation) %%% + %%% Dissipationen uppbyggd av D4: %%% + %%% DI=D4*B*H*D4 %%% + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + %This is 3rd order accurate at the boundary. Not same norm as D1 operator + + BP = 8; + if(m<2*BP) + error(['Operator requires at least ' num2str(2*BP) ' grid points']); + end + + % Norm + Hv = ones(m,1); + Hv(1:8) = [0.7488203e7/0.25401600e8, 0.5539027e7/0.3628800e7, 0.308923e6/0.1209600e7, 0.1307491e7/0.725760e6, 0.59407e5/0.145152e6, 0.1548947e7/0.1209600e7, 0.3347963e7/0.3628800e7, 0.25641187e8/0.25401600e8]; + Hv(m-7:m) = rot90(Hv(1:8),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:7) = [-0.49e2/0.20e2 6 -0.15e2/0.2e1 0.20e2/0.3e1 -0.15e2/0.4e1 0.6e1/0.5e1 -0.1e1/0.6e1]/h; + d1_r = -rot90(d1_l); + + d2_l = sparse(m,1); + d2_l(1:7) = [0.203e3/0.45e2 -0.87e2/0.5e1 0.117e3/0.4e1 -0.254e3/0.9e1 0.33e2/0.2e1 -0.27e2/0.5e1 0.137e3/0.180e3]/h^2; + d2_r = rot90(d2_l, 2); + + d3_l = sparse(m,1); + d3_l(1:7) = [-0.49e2/0.8e1 29 -0.461e3/0.8e1 62 -0.307e3/0.8e1 13 -0.15e2/0.8e1]/h^3; + d3_r = -rot90(d3_l, 2); + + + + % Fourth derivative, 1th order accurate at first 8 boundary points (still + % yield 5th order convergence if stable: for example u_tt = -u_xxxx + + stencil = [-0.41e2/0.7560e4, 0.1261e4/0.15120e5, -0.541e3/0.840e3, 0.4369e4/0.1260e4, -0.1669e4/0.180e3, 0.1529e4/0.120e3, -0.1669e4/0.180e3, 0.4369e4/0.1260e4, -0.541e3/0.840e3, 0.1261e4/0.15120e5,-0.41e2/0.7560e4]; + diags = -5:5; + + M4_U = [ + 0.1031569831e10/0.155675520e9 -0.32874237931e11/0.1452971520e10 0.3069551773e10/0.90810720e8 -0.658395212131e12/0.21794572800e11 0.31068454007e11/0.1816214400e10 -0.39244130657e11/0.7264857600e10 0.1857767503e10/0.2724321600e10 0.1009939e7/0.49420800e8; + -0.32874237931e11/0.1452971520e10 0.12799022387e11/0.155675520e9 -0.134456503627e12/0.1037836800e10 0.15366749479e11/0.129729600e9 -0.207640325549e12/0.3113510400e10 0.5396424073e10/0.259459200e9 -0.858079351e9/0.345945600e9 -0.19806607e8/0.170270100e9; + 0.3069551773e10/0.90810720e8 -0.134456503627e12/0.1037836800e10 0.6202056779e10/0.28828800e8 -0.210970327081e12/0.1037836800e10 0.2127730129e10/0.18532800e8 -0.4048692749e10/0.115315200e9 0.1025943959e10/0.259459200e9 0.71054663e8/0.290594304e9; + -0.658395212131e12/0.21794572800e11 0.15366749479e11/0.129729600e9 -0.210970327081e12/0.1037836800e10 0.31025293213e11/0.155675520e9 -0.1147729001e10/0.9884160e7 0.1178067773e10/0.32432400e8 -0.13487255581e11/0.3113510400e10 -0.231082547e9/0.1816214400e10; + 0.31068454007e11/0.1816214400e10 -0.207640325549e12/0.3113510400e10 0.2127730129e10/0.18532800e8 -0.1147729001e10/0.9884160e7 0.11524865123e11/0.155675520e9 -0.29754506009e11/0.1037836800e10 0.14231221e8/0.2316600e7 -0.15030629699e11/0.21794572800e11; + -0.39244130657e11/0.7264857600e10 0.5396424073e10/0.259459200e9 -0.4048692749e10/0.115315200e9 0.1178067773e10/0.32432400e8 -0.29754506009e11/0.1037836800e10 0.572247737e9/0.28828800e8 -0.11322059051e11/0.1037836800e10 0.3345834083e10/0.908107200e9; + 0.1857767503e10/0.2724321600e10 -0.858079351e9/0.345945600e9 0.1025943959e10/0.259459200e9 -0.13487255581e11/0.3113510400e10 0.14231221e8/0.2316600e7 -0.11322059051e11/0.1037836800e10 0.10478882597e11/0.778377600e9 -0.68446325191e11/0.7264857600e10; + 0.1009939e7/0.49420800e8 -0.19806607e8/0.170270100e9 0.71054663e8/0.290594304e9 -0.231082547e9/0.1816214400e10 -0.15030629699e11/0.21794572800e11 0.3345834083e10/0.908107200e9 -0.68446325191e11/0.7264857600e10 0.9944747557e10/0.778377600e9; + ]; + + M4(1:8,1:8) = M4_U; + M4(m-7:m,m-7:m) = rot90(M4_U, 2); + M4 = 1/h^3*M4; + + D4=HI*(M4 - e_l*d3_l'+e_r*d3_r' + d1_l*d2_l'-d1_r*d2_r'); +end
diff -r c0cbffcf6513 -r 72468bc9b63f +sbp/+implementations/d4_lonely_8_min_boundary_points.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+sbp/+implementations/d4_lonely_8_min_boundary_points.m Mon Sep 26 09:55:16 2016 +0200 @@ -0,0 +1,73 @@ +function [H, HI, D4, e_l, e_r, M4, d2_l, d2_r, d3_l, d3_r, d1_l, d1_r] = d4_variable_8_min_boundary_points(m,h) + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + %%% 8:te ordn. SBP Finita differens %%% + %%% operatorer med diagonal norm %%% + %%% %%% + %%% %%% + %%% H (Normen) %%% + %%% D1=H^(-1)Q (approx f?rsta derivatan) %%% + %%% D2 (approx andra derivatan) %%% + %%% D2=HI*(R+C*D*S %%% + %%% %%% + %%% R=-D1'*H*C*D1-RR %%% + %%% %%% + %%% RR ?r dissipation) %%% + %%% Dissipationen uppbyggd av D4: %%% + %%% DI=D4*B*H*D4 %%% + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + BP = 8; + if(m<2*BP) + error(['Operator requires at least ' num2str(2*BP) ' grid points']); + end + + % Norm + Hv = ones(m,1); + Hv(1:8) = [1498139/5080320, 1107307/725760, 20761/80640, 1304999/725760, 299527/725760, 103097/80640, 670091/725760, 5127739/5080320]; + Hv(m-7:m) = rot90(Hv(1:8),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:6) = [-0.137e3/0.60e2 5 -5 0.10e2/0.3e1 -0.5e1/0.4e1 0.1e1/0.5e1;]/h; + d1_r = -rot90(d1_l); + + d2_l = sparse(m,1); + d2_l(1:6) = [0.15e2/0.4e1 -0.77e2/0.6e1 0.107e3/0.6e1 -13 0.61e2/0.12e2 -0.5e1/0.6e1;]/h^2; + d2_r = rot90(d2_l, 2); + + d3_l = sparse(m,1); + d3_l(1:6) = [-0.17e2/0.4e1 0.71e2/0.4e1 -0.59e2/0.2e1 0.49e2/0.2e1 -0.41e2/0.4e1 0.7e1/0.4e1;]/h^3; + d3_r = -rot90(d3_l, 2); + + + % Fourth derivative, 1th order accurate at first 8 boundary points + + stencil = [-0.41e2/0.7560e4, 0.1261e4/0.15120e5,-0.541e3/0.840e3,0.4369e4/0.1260e4,-0.1669e4/0.180e3,0.1529e4/0.120e3,-0.1669e4/0.180e3,0.4369e4/0.1260e4,-0.541e3/0.840e3, 0.1261e4/0.15120e5,-0.41e2/0.7560e4]; + diags = -5:5; + M4 = stripeMatrix(stencil, diags, m); + + M4_U = [ + 0.151705142321e12/0.29189160000e11 -0.25643455801727e14/0.1634592960000e13 0.286417898677e12/0.15135120000e11 -0.4038072020317e13/0.326918592000e12 0.96455968907e11/0.20432412000e11 -0.151076916769e12/0.181621440000e12 0.14511526363e11/0.408648240000e12 -0.196663079e9/0.33359040000e11; + -0.25643455801727e14/0.1634592960000e13 0.735383382473e12/0.14594580000e11 -0.5035391734409e13/0.77837760000e11 0.20392440917e11/0.467026560e9 -0.109540902413e12/0.6671808000e10 0.2488686539e10/0.884520000e9 -0.2798067539e10/0.33359040000e11 0.6433463591e10/0.408648240000e12; + 0.286417898677e12/0.15135120000e11 -0.5035391734409e13/0.77837760000e11 0.145019791981e12/0.1621620000e10 -0.333577111061e12/0.5189184000e10 0.18928722391e11/0.778377600e9 -0.93081704557e11/0.25945920000e11 -0.372660319e9/0.3243240000e10 0.2861399869e10/0.544864320000e12; + -0.4038072020317e13/0.326918592000e12 0.20392440917e11/0.467026560e9 -0.333577111061e12/0.5189184000e10 0.59368471277e11/0.1167566400e10 -0.201168708569e12/0.9340531200e10 0.1492314487e10/0.432432000e9 0.1911896257e10/0.9340531200e10 0.24383341e8/0.2554051500e10; + 0.96455968907e11/0.20432412000e11 -0.109540902413e12/0.6671808000e10 0.18928722391e11/0.778377600e9 -0.201168708569e12/0.9340531200e10 0.1451230301e10/0.106142400e9 -0.103548247007e12/0.15567552000e11 0.27808437809e11/0.11675664000e11 -0.36870830713e11/0.65383718400e11; + -0.151076916769e12/0.181621440000e12 0.2488686539e10/0.884520000e9 -0.93081704557e11/0.25945920000e11 0.1492314487e10/0.432432000e9 -0.103548247007e12/0.15567552000e11 0.1229498243e10/0.115830000e9 -0.32222519717e11/0.3706560000e10 0.470092704233e12/0.136216080000e12; + 0.14511526363e11/0.408648240000e12 -0.2798067539e10/0.33359040000e11 -0.372660319e9/0.3243240000e10 0.1911896257e10/0.9340531200e10 0.27808437809e11/0.11675664000e11 -0.32222519717e11/0.3706560000e10 0.11547819313e11/0.912161250e9 -0.15187033999199e14/0.1634592960000e13; + -0.196663079e9/0.33359040000e11 0.6433463591e10/0.408648240000e12 0.2861399869e10/0.544864320000e12 0.24383341e8/0.2554051500e10 -0.36870830713e11/0.65383718400e11 0.470092704233e12/0.136216080000e12 -0.15187033999199e14/0.1634592960000e13 0.33832994693e11/0.2653560000e10; + ]; + + M4(1:8,1:8) = M4_U; + M4(m-7:m,m-7:m) = rot90(M4_U, 2); + M4 = 1/h^3*M4; + + D4=HI*(M4 - e_l*d3_l'+e_r*d3_r' + d1_l*d2_l'-d1_r*d2_r'); +end
diff -r c0cbffcf6513 -r 72468bc9b63f +sbp/+implementations/d4_variable_4_min_boundary_points.m --- a/+sbp/+implementations/d4_variable_4_min_boundary_points.m Mon Sep 26 09:51:45 2016 +0200 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,81 +0,0 @@ -function [H, HI, D4, e_l, e_r, M4, d2_l, d2_r, d3_l, d3_r, d1_l, d1_r] = d4_variable_4_min_boundary_points(m,h) - %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - %%% 4:de ordn. SBP Finita differens %%% - %%% operatorer framtagna av Mark Carpenter %%% - %%% %%% - %%% H (Normen) %%% - %%% D1=H^(-1)Q (approx f?rsta derivatan) %%% - %%% D2 (approx andra derivatan) %%% - %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - %H?r med endast 4 randpunkter - - - BP = 4; - if(m<2*BP) - error(['Operator requires at least ' num2str(2*BP) ' grid points']); - end - - - % Norm - Hv = ones(m,1); - Hv(1:4) = [17/48 59/48 43/48 49/48]; - Hv(m-3:m) = rot90(Hv(1:4),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:4) = 1/h*[-11/6 3 -3/2 1/3]; - d1_r = -rot90(d1_l); - - d2_l = sparse(m,1); - d2_l(1:4) = 1/h^2*[2 -5 4 -1]; - d2_r = rot90(d2_l, 2); - - d3_l = sparse(m,1); - d3_l(1:4) = 1/h^3*[-1 3 -3 1]; - d3_r = -rot90(d3_l, 2); - - - % First derivative - stencil = [1/12 -2/3 0 2/3 -1/12]; - diags = [-1 0 1]; - - Q_U = [ - 0 0.59e2/0.96e2 -0.1e1/0.12e2 -0.1e1/0.32e2; - -0.59e2/0.96e2 0 0.59e2/0.96e2 0; - 0.1e1/0.12e2 -0.59e2/0.96e2 0 0.59e2/0.96e2; - 0.1e1/0.32e2 0 -0.59e2/0.96e2 0; - ]; - - Q = stripeMatrix(stencil, diags, m); - Q(1:4,1:4)=Q_U; - Q(m-3:m,m-3:m) = -rot90(Q_U, 2); - - D1 = HI*(Q - 1/2*e_l*e_l' + 1/2*e_r*e_r'); - - % Fourth derivative - stencil = [-1/6, 2, -13/2, 28/3, -13/2, 2, -1/6]; - diags = -3:3; - M4 = stripeMatrix(stencil, diags, m); - - M4_U=[ - 0.8e1/0.3e1 -0.37e2/0.6e1 0.13e2/0.3e1 -0.5e1/0.6e1; - -0.37e2/0.6e1 0.47e2/0.3e1 -13 0.11e2/0.3e1; - 0.13e2/0.3e1 -13 0.44e2/0.3e1 -0.47e2/0.6e1; - -0.5e1/0.6e1 0.11e2/0.3e1 -0.47e2/0.6e1 0.29e2/0.3e1; - ]; - - - M4(1:4,1:4) = M4_U; - M4(m-3:m,m-3:m) = rot90(M4_U, 2); - M4 = 1/h^3*M4; - - D4=HI*(M4 - e_l*d3_l'+e_r*d3_r' + d1_l*d2_l'-d1_r*d2_r'); -end \ No newline at end of file
diff -r c0cbffcf6513 -r 72468bc9b63f +sbp/+implementations/d4_variable_6_2.m --- a/+sbp/+implementations/d4_variable_6_2.m Mon Sep 26 09:51:45 2016 +0200 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,76 +0,0 @@ -function [H, HI, D4, e_l, e_r, M4, d2_l, d2_r, d3_l, d3_r, d1_l, d1_r] = d4_variable_6_2(m,h) - %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - %%% 6:te ordn. SBP Finita differens %%% - %%% operatorer med diagonal norm %%% - %%% Extension to variable koeff %%% - %%% %%% - %%% H (Normen) %%% - %%% D1=H^(-1)Q (approx f?rsta derivatan) %%% - %%% D2 (approx andra derivatan) %%% - %%% D2=HI*(R+C*D*S %%% - %%% %%% - %%% R=-D1'*H*C*D1-RR %%% - %%% %%% - %%% RR ?r dissipation) %%% - %%% Dissipationen uppbyggd av D4: %%% - %%% DI=D4*B*H*D4 %%% - %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - - % H?r med 6 RP ist?llet f?r 8 f?r D4 operatorn, dock samma randderivator - % Denna ?r noggrannare, och har 2a ordningens randdslutning och b?r ge 6te - % ordningens konvergens. Hade dock ingen fri parameter att optimera - - BP = 6; - if(m<2*BP) - error(['Operator requires at least ' num2str(2*BP) ' grid points']); - end - - % Norm - Hv = ones(m,1); - Hv(1:6) = [0.181e3/0.576e3, 0.1343e4/0.960e3, 0.293e3/0.480e3, 0.1811e4/0.1440e4, 0.289e3/0.320e3, 0.65e2/0.64e2]; - Hv(m-5:m) = rot90(Hv(1:6),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:6) = [-0.137e3/0.60e2 5 -5 0.10e2/0.3e1 -0.5e1/0.4e1 0.1e1/0.5e1;]/h; - d1_r = -rot90(d1_l); - - d2_l = sparse(m,1); - d2_l(1:6) = [0.15e2/0.4e1 -0.77e2/0.6e1 0.107e3/0.6e1 -13 0.61e2/0.12e2 -0.5e1/0.6e1;]/h^2; - d2_r = rot90(d2_l, 2); - - d3_l = sparse(m,1); - d3_l(1:6) = [-0.17e2/0.4e1 0.71e2/0.4e1 -0.59e2/0.2e1 0.49e2/0.2e1 -0.41e2/0.4e1 0.7e1/0.4e1;]/h^3; - d3_r = -rot90(d3_l, 2); - - - % Fourth derivative, 1th order accurate at first 8 boundary points (still - % yield 5th order convergence if stable: for example u_tt = -u_xxxx - stencil = [7/240, -2/5, 169/60, -122/15, 91/8, -122/15, 169/60, -2/5, 7/240]; - diags = -4:4; - M4 = stripeMatrix(stencil, diags, m); - - M4_U = [ - 0.1009e4/0.192e3 -0.7657e4/0.480e3 0.9307e4/0.480e3 -0.509e3/0.40e2 0.4621e4/0.960e3 -0.25e2/0.32e2; - -0.7657e4/0.480e3 0.49513e5/0.960e3 -0.4007e4/0.60e2 0.21799e5/0.480e3 -0.8171e4/0.480e3 0.2657e4/0.960e3; - 0.9307e4/0.480e3 -0.4007e4/0.60e2 0.1399e4/0.15e2 -0.2721e4/0.40e2 0.12703e5/0.480e3 -0.521e3/0.120e3; - -0.509e3/0.40e2 0.21799e5/0.480e3 -0.2721e4/0.40e2 0.3349e4/0.60e2 -0.389e3/0.15e2 0.559e3/0.96e2; - 0.4621e4/0.960e3 -0.8171e4/0.480e3 0.12703e5/0.480e3 -0.389e3/0.15e2 0.17857e5/0.960e3 -0.1499e4/0.160e3; - -0.25e2/0.32e2 0.2657e4/0.960e3 -0.521e3/0.120e3 0.559e3/0.96e2 -0.1499e4/0.160e3 0.2225e4/0.192e3; - ]; - - - M4(1:6,1:6) = M4_U; - M4(m-5:m,m-5:m) = rot90(M4_U, 2); - M4 = 1/h^3*M4; - - D4=HI*(M4 - e_l*d3_l'+e_r*d3_r' + d1_l*d2_l'-d1_r*d2_r'); -end
diff -r c0cbffcf6513 -r 72468bc9b63f +sbp/+implementations/d4_variable_6_3.m --- a/+sbp/+implementations/d4_variable_6_3.m Mon Sep 26 09:51:45 2016 +0200 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,70 +0,0 @@ -function [H, HI, D4, e_l, e_r, M4, d2_l, d2_r, d3_l, d3_r, d1_l, d1_r] = d4_variable_6_3(m,h) - %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - %%% 6:te ordn. SBP Finita differens %%% - %%% operatorer med diagonal norm %%% - %%% Extension to variable koeff %%% - %%% %%% - %%% H (Normen) %%% - %%% D1=H^(-1)Q (approx f?rsta derivatan) %%% - %%% D2 (approx andra derivatan) %%% - %%% D2=HI*(R+C*D*S %%% - %%% %%% - %%% R=-D1'*H*C*D1-RR %%% - %%% %%% - %%% RR ?r dissipation) %%% - %%% Dissipationen uppbyggd av D4: %%% - %%% DI=D4*B*H*D4 %%% - %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - - % H?r med 7 RP ist?llet f?r 8 f?r D4 operatorn, dock samma randderivator - % Denna ?r noggrannare, och har 2a ordningens randdslutning och b?r ge 6te - % ordningens konvergens. Hade 2 fria parametrar att optimera - - % Norm - Hv = ones(m,1); - Hv(1:7) = [0.414837907e9/0.1191965760e10, 0.475278367e9/0.397321920e9, 0.13872751e8/0.12416310e8, 0.346739027e9/0.595982880e9, 0.560227469e9/0.397321920e9, 0.322971631e9/0.397321920e9, 0.616122491e9/0.595982880e9]; - Hv(m-6:m) = rot90(Hv(1:7),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:6) = [-0.137e3/0.60e2 5 -5 0.10e2/0.3e1 -0.5e1/0.4e1 0.1e1/0.5e1;]/h; - d1_r = -rot90(d1_l); - - d2_l = sparse(m,1); - d2_l(1:6) = [0.15e2/0.4e1 -0.77e2/0.6e1 0.107e3/0.6e1 -13 0.61e2/0.12e2 -0.5e1/0.6e1;]/h^2; - d2_r = rot90(d2_l, 2); - - d3_l = sparse(m,1); - d3_l(1:6) = [-0.17e2/0.4e1 0.71e2/0.4e1 -0.59e2/0.2e1 0.49e2/0.2e1 -0.41e2/0.4e1 0.7e1/0.4e1;]/h^3; - d3_r = -rot90(d3_l, 2); - - - % Fourth derivative, 1th order accurate at first 8 boundary points - stencil = [7/240, -2/5, 169/60, -122/15, 91/8, -122/15, 169/60, -2/5, 7/240]; - diags = -4:4; - M4 = stripeMatrix(stencil, diags, m); - - M4_U = [ - 0.1399708478939e13/0.263487168000e12 -0.13482796013041e14/0.834376032000e12 0.344344095859e12/0.17565811200e11 -0.3166261424681e13/0.250312809600e12 0.1508605165681e13/0.333750412800e12 -0.486270829441e12/0.834376032000e12 -0.221976356359e12/0.5006256192000e13; - -0.13482796013041e14/0.834376032000e12 0.7260475818391e13/0.139062672000e12 -0.27224036353e11/0.406022400e9 0.1847477458951e13/0.41718801600e11 -0.848984558161e12/0.55625068800e11 0.247494925991e12/0.139062672000e12 0.165585445559e12/0.834376032000e12; - 0.344344095859e12/0.17565811200e11 -0.27224036353e11/0.406022400e9 0.2044938640393e13/0.22250027520e11 -0.1071086785417e13/0.16687520640e11 0.502199537033e12/0.22250027520e11 -0.143589154441e12/0.55625068800e11 -0.88181965559e11/0.333750412800e12; - -0.3166261424681e13/0.250312809600e12 0.1847477458951e13/0.41718801600e11 -0.1071086785417e13/0.16687520640e11 0.628860435593e12/0.12515640480e11 -0.73736245829e11/0.3337504128e10 0.195760572271e12/0.41718801600e11 -0.81156046361e11/0.250312809600e12; - 0.1508605165681e13/0.333750412800e12 -0.848984558161e12/0.55625068800e11 0.502199537033e12/0.22250027520e11 -0.73736245829e11/0.3337504128e10 0.76725285869e11/0.4450005504e10 -0.3912429433e10/0.406022400e9 0.53227370659e11/0.17565811200e11; - -0.486270829441e12/0.834376032000e12 0.247494925991e12/0.139062672000e12 -0.143589154441e12/0.55625068800e11 0.195760572271e12/0.41718801600e11 -0.3912429433e10/0.406022400e9 0.1699707221791e13/0.139062672000e12 -0.6959018412841e13/0.834376032000e12; - -0.221976356359e12/0.5006256192000e13 0.165585445559e12/0.834376032000e12 -0.88181965559e11/0.333750412800e12 -0.81156046361e11/0.250312809600e12 0.53227370659e11/0.17565811200e11 -0.6959018412841e13/0.834376032000e12 0.3012195053939e13/0.263487168000e12; - ]; - - M4(1:7,1:7) = M4_U; - M4(m-6:m,m-6:m) = rot90(M4_U, 2); - M4 = 1/h^3*M4; - - D4=HI*(M4 - e_l*d3_l'+e_r*d3_r' + d1_l*d2_l'-d1_r*d2_r'); -end
diff -r c0cbffcf6513 -r 72468bc9b63f +sbp/+implementations/d4_variable_6_min_boundary_points.m --- a/+sbp/+implementations/d4_variable_6_min_boundary_points.m Mon Sep 26 09:51:45 2016 +0200 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,74 +0,0 @@ -function [H, HI, D4, e_l, e_r, M4, d2_l, d2_r, d3_l, d3_r, d1_l, d1_r] = d4_variable_6_min_boundary_points(m,h) - %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - %%% 6:te ordn. SBP Finita differens %%% - %%% operatorer med diagonal norm %%% - %%% Extension to variable koeff %%% - %%% %%% - %%% H (Normen) %%% - %%% D1=H^(-1)Q (approx f?rsta derivatan) %%% - %%% D2 (approx andra derivatan) %%% - %%% D2=HI*(R+C*D*S %%% - %%% %%% - %%% R=-D1'*H*C*D1-RR %%% - %%% %%% - %%% RR ?r dissipation) %%% - %%% Dissipationen uppbyggd av D4: %%% - %%% DI=D4*B*H*D4 %%% - %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - - % H?r med 6 RP ist?llet f?r 8 f?r D4 operatorn, dock samma randderivator - - BP = 6; - if(m<2*BP) - error(['Operator requires at least ' num2str(2*BP) ' grid points']); - end - - % Norm - Hv = ones(m,1); - Hv(1:6) = [13649/43200,12013/8640,2711/4320,5359/4320,7877/8640, 43801/43200]; - Hv(m-5:m) = rot90(Hv(1:6),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:5) = [-25/12, 4, -3, 4/3, -1/4]/h; - d1_r = -rot90(d1_l); - - d2_l = sparse(m,1); - d2_l(1:5) = [0.35e2/0.12e2 -0.26e2/0.3e1 0.19e2/0.2e1 -0.14e2/0.3e1 0.11e2/0.12e2;]/h^2; - d2_r = rot90(d2_l, 2); - - d3_l = sparse(m,1); - d3_l(1:5) = [-0.5e1/0.2e1 9 -12 7 -0.3e1/0.2e1;]/h^3; - d3_r = -rot90(d3_l, 2); - - - % Fourth derivative, 1th order accurate at first 8 boundary points (still - % yield 5th order convergence if stable: for example u_tt=-u_xxxx - - stencil = [7/240, -2/5, 169/60, -122/15, 91/8, -122/15, 169/60, -2/5, 7/240]; - diags = -4:4; - M4 = stripeMatrix(stencil, diags, m); - - M4_U=[ - 0.3504379e7/0.907200e6 -0.4613983e7/0.453600e6 0.4260437e7/0.453600e6 -0.418577e6/0.113400e6 0.524579e6/0.907200e6 0.535e3/0.18144e5; - -0.4613983e7/0.453600e6 0.5186159e7/0.181440e6 -0.81121e5/0.2835e4 0.218845e6/0.18144e5 -0.159169e6/0.90720e5 -0.94669e5/0.907200e6; - 0.4260437e7/0.453600e6 -0.81121e5/0.2835e4 0.147695e6/0.4536e4 -0.384457e6/0.22680e5 0.339653e6/0.90720e5 -0.18233e5/0.113400e6; - -0.418577e6/0.113400e6 0.218845e6/0.18144e5 -0.384457e6/0.22680e5 0.65207e5/0.4536e4 -0.22762e5/0.2835e4 0.1181753e7/0.453600e6; - 0.524579e6/0.907200e6 -0.159169e6/0.90720e5 0.339653e6/0.90720e5 -0.22762e5/0.2835e4 0.2006171e7/0.181440e6 -0.3647647e7/0.453600e6; - 0.535e3/0.18144e5 -0.94669e5/0.907200e6 -0.18233e5/0.113400e6 0.1181753e7/0.453600e6 -0.3647647e7/0.453600e6 0.10305271e8/0.907200e6; - ]; - - M4(1:6,1:6) = M4_U; - M4(m-5:m,m-5:m) = rot90(M4_U, 2); - M4 = 1/h^3*M4; - - D4=HI*(M4 - e_l*d3_l'+e_r*d3_r' + d1_l*d2_l'-d1_r*d2_r'); -end
diff -r c0cbffcf6513 -r 72468bc9b63f +sbp/+implementations/d4_variable_8_higher_boundary_order.m --- a/+sbp/+implementations/d4_variable_8_higher_boundary_order.m Mon Sep 26 09:51:45 2016 +0200 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,75 +0,0 @@ -function [H, HI, D4, e_l, e_r, M4, d2_l, d2_r, d3_l, d3_r, d1_l, d1_r] = d4_variable_8_higher_boundary_order(m,h) - %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - %%% 8:te ordn. SBP Finita differens %%% - %%% operatorer med diagonal norm %%% - %%% %%% - %%% %%% - %%% H (Normen) %%% - %%% D1=H^(-1)Q (approx f?rsta derivatan) %%% - %%% D2 (approx andra derivatan) %%% - %%% D2=HI*(R+C*D*S %%% - %%% %%% - %%% R=-D1'*H*C*D1-RR %%% - %%% %%% - %%% RR ?r dissipation) %%% - %%% Dissipationen uppbyggd av D4: %%% - %%% DI=D4*B*H*D4 %%% - %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - %This is 3rd order accurate at the boundary. Not same norm as D1 operator - - BP = 8; - if(m<2*BP) - error(['Operator requires at least ' num2str(2*BP) ' grid points']); - end - - % Norm - Hv = ones(m,1); - Hv(1:8) = [0.7488203e7/0.25401600e8, 0.5539027e7/0.3628800e7, 0.308923e6/0.1209600e7, 0.1307491e7/0.725760e6, 0.59407e5/0.145152e6, 0.1548947e7/0.1209600e7, 0.3347963e7/0.3628800e7, 0.25641187e8/0.25401600e8]; - Hv(m-7:m) = rot90(Hv(1:8),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:7) = [-0.49e2/0.20e2 6 -0.15e2/0.2e1 0.20e2/0.3e1 -0.15e2/0.4e1 0.6e1/0.5e1 -0.1e1/0.6e1]/h; - d1_r = -rot90(d1_l); - - d2_l = sparse(m,1); - d2_l(1:7) = [0.203e3/0.45e2 -0.87e2/0.5e1 0.117e3/0.4e1 -0.254e3/0.9e1 0.33e2/0.2e1 -0.27e2/0.5e1 0.137e3/0.180e3]/h^2; - d2_r = rot90(d2_l, 2); - - d3_l = sparse(m,1); - d3_l(1:7) = [-0.49e2/0.8e1 29 -0.461e3/0.8e1 62 -0.307e3/0.8e1 13 -0.15e2/0.8e1]/h^3; - d3_r = -rot90(d3_l, 2); - - - - % Fourth derivative, 1th order accurate at first 8 boundary points (still - % yield 5th order convergence if stable: for example u_tt = -u_xxxx - - stencil = [-0.41e2/0.7560e4, 0.1261e4/0.15120e5, -0.541e3/0.840e3, 0.4369e4/0.1260e4, -0.1669e4/0.180e3, 0.1529e4/0.120e3, -0.1669e4/0.180e3, 0.4369e4/0.1260e4, -0.541e3/0.840e3, 0.1261e4/0.15120e5,-0.41e2/0.7560e4]; - diags = -5:5; - - M4_U = [ - 0.1031569831e10/0.155675520e9 -0.32874237931e11/0.1452971520e10 0.3069551773e10/0.90810720e8 -0.658395212131e12/0.21794572800e11 0.31068454007e11/0.1816214400e10 -0.39244130657e11/0.7264857600e10 0.1857767503e10/0.2724321600e10 0.1009939e7/0.49420800e8; - -0.32874237931e11/0.1452971520e10 0.12799022387e11/0.155675520e9 -0.134456503627e12/0.1037836800e10 0.15366749479e11/0.129729600e9 -0.207640325549e12/0.3113510400e10 0.5396424073e10/0.259459200e9 -0.858079351e9/0.345945600e9 -0.19806607e8/0.170270100e9; - 0.3069551773e10/0.90810720e8 -0.134456503627e12/0.1037836800e10 0.6202056779e10/0.28828800e8 -0.210970327081e12/0.1037836800e10 0.2127730129e10/0.18532800e8 -0.4048692749e10/0.115315200e9 0.1025943959e10/0.259459200e9 0.71054663e8/0.290594304e9; - -0.658395212131e12/0.21794572800e11 0.15366749479e11/0.129729600e9 -0.210970327081e12/0.1037836800e10 0.31025293213e11/0.155675520e9 -0.1147729001e10/0.9884160e7 0.1178067773e10/0.32432400e8 -0.13487255581e11/0.3113510400e10 -0.231082547e9/0.1816214400e10; - 0.31068454007e11/0.1816214400e10 -0.207640325549e12/0.3113510400e10 0.2127730129e10/0.18532800e8 -0.1147729001e10/0.9884160e7 0.11524865123e11/0.155675520e9 -0.29754506009e11/0.1037836800e10 0.14231221e8/0.2316600e7 -0.15030629699e11/0.21794572800e11; - -0.39244130657e11/0.7264857600e10 0.5396424073e10/0.259459200e9 -0.4048692749e10/0.115315200e9 0.1178067773e10/0.32432400e8 -0.29754506009e11/0.1037836800e10 0.572247737e9/0.28828800e8 -0.11322059051e11/0.1037836800e10 0.3345834083e10/0.908107200e9; - 0.1857767503e10/0.2724321600e10 -0.858079351e9/0.345945600e9 0.1025943959e10/0.259459200e9 -0.13487255581e11/0.3113510400e10 0.14231221e8/0.2316600e7 -0.11322059051e11/0.1037836800e10 0.10478882597e11/0.778377600e9 -0.68446325191e11/0.7264857600e10; - 0.1009939e7/0.49420800e8 -0.19806607e8/0.170270100e9 0.71054663e8/0.290594304e9 -0.231082547e9/0.1816214400e10 -0.15030629699e11/0.21794572800e11 0.3345834083e10/0.908107200e9 -0.68446325191e11/0.7264857600e10 0.9944747557e10/0.778377600e9; - ]; - - M4(1:8,1:8) = M4_U; - M4(m-7:m,m-7:m) = rot90(M4_U, 2); - M4 = 1/h^3*M4; - - D4=HI*(M4 - e_l*d3_l'+e_r*d3_r' + d1_l*d2_l'-d1_r*d2_r'); -end
diff -r c0cbffcf6513 -r 72468bc9b63f +sbp/+implementations/d4_variable_8_min_boundary_points.m --- a/+sbp/+implementations/d4_variable_8_min_boundary_points.m Mon Sep 26 09:51:45 2016 +0200 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,73 +0,0 @@ -function [H, HI, D4, e_l, e_r, M4, d2_l, d2_r, d3_l, d3_r, d1_l, d1_r] = d4_variable_8_min_boundary_points(m,h) - %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - %%% 8:te ordn. SBP Finita differens %%% - %%% operatorer med diagonal norm %%% - %%% %%% - %%% %%% - %%% H (Normen) %%% - %%% D1=H^(-1)Q (approx f?rsta derivatan) %%% - %%% D2 (approx andra derivatan) %%% - %%% D2=HI*(R+C*D*S %%% - %%% %%% - %%% R=-D1'*H*C*D1-RR %%% - %%% %%% - %%% RR ?r dissipation) %%% - %%% Dissipationen uppbyggd av D4: %%% - %%% DI=D4*B*H*D4 %%% - %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - - BP = 8; - if(m<2*BP) - error(['Operator requires at least ' num2str(2*BP) ' grid points']); - end - - % Norm - Hv = ones(m,1); - Hv(1:8) = [1498139/5080320, 1107307/725760, 20761/80640, 1304999/725760, 299527/725760, 103097/80640, 670091/725760, 5127739/5080320]; - Hv(m-7:m) = rot90(Hv(1:8),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:6) = [-0.137e3/0.60e2 5 -5 0.10e2/0.3e1 -0.5e1/0.4e1 0.1e1/0.5e1;]/h; - d1_r = -rot90(d1_l); - - d2_l = sparse(m,1); - d2_l(1:6) = [0.15e2/0.4e1 -0.77e2/0.6e1 0.107e3/0.6e1 -13 0.61e2/0.12e2 -0.5e1/0.6e1;]/h^2; - d2_r = rot90(d2_l, 2); - - d3_l = sparse(m,1); - d3_l(1:6) = [-0.17e2/0.4e1 0.71e2/0.4e1 -0.59e2/0.2e1 0.49e2/0.2e1 -0.41e2/0.4e1 0.7e1/0.4e1;]/h^3; - d3_r = -rot90(d3_l, 2); - - - % Fourth derivative, 1th order accurate at first 8 boundary points - - stencil = [-0.41e2/0.7560e4, 0.1261e4/0.15120e5,-0.541e3/0.840e3,0.4369e4/0.1260e4,-0.1669e4/0.180e3,0.1529e4/0.120e3,-0.1669e4/0.180e3,0.4369e4/0.1260e4,-0.541e3/0.840e3, 0.1261e4/0.15120e5,-0.41e2/0.7560e4]; - diags = -5:5; - M4 = stripeMatrix(stencil, diags, m); - - M4_U = [ - 0.151705142321e12/0.29189160000e11 -0.25643455801727e14/0.1634592960000e13 0.286417898677e12/0.15135120000e11 -0.4038072020317e13/0.326918592000e12 0.96455968907e11/0.20432412000e11 -0.151076916769e12/0.181621440000e12 0.14511526363e11/0.408648240000e12 -0.196663079e9/0.33359040000e11; - -0.25643455801727e14/0.1634592960000e13 0.735383382473e12/0.14594580000e11 -0.5035391734409e13/0.77837760000e11 0.20392440917e11/0.467026560e9 -0.109540902413e12/0.6671808000e10 0.2488686539e10/0.884520000e9 -0.2798067539e10/0.33359040000e11 0.6433463591e10/0.408648240000e12; - 0.286417898677e12/0.15135120000e11 -0.5035391734409e13/0.77837760000e11 0.145019791981e12/0.1621620000e10 -0.333577111061e12/0.5189184000e10 0.18928722391e11/0.778377600e9 -0.93081704557e11/0.25945920000e11 -0.372660319e9/0.3243240000e10 0.2861399869e10/0.544864320000e12; - -0.4038072020317e13/0.326918592000e12 0.20392440917e11/0.467026560e9 -0.333577111061e12/0.5189184000e10 0.59368471277e11/0.1167566400e10 -0.201168708569e12/0.9340531200e10 0.1492314487e10/0.432432000e9 0.1911896257e10/0.9340531200e10 0.24383341e8/0.2554051500e10; - 0.96455968907e11/0.20432412000e11 -0.109540902413e12/0.6671808000e10 0.18928722391e11/0.778377600e9 -0.201168708569e12/0.9340531200e10 0.1451230301e10/0.106142400e9 -0.103548247007e12/0.15567552000e11 0.27808437809e11/0.11675664000e11 -0.36870830713e11/0.65383718400e11; - -0.151076916769e12/0.181621440000e12 0.2488686539e10/0.884520000e9 -0.93081704557e11/0.25945920000e11 0.1492314487e10/0.432432000e9 -0.103548247007e12/0.15567552000e11 0.1229498243e10/0.115830000e9 -0.32222519717e11/0.3706560000e10 0.470092704233e12/0.136216080000e12; - 0.14511526363e11/0.408648240000e12 -0.2798067539e10/0.33359040000e11 -0.372660319e9/0.3243240000e10 0.1911896257e10/0.9340531200e10 0.27808437809e11/0.11675664000e11 -0.32222519717e11/0.3706560000e10 0.11547819313e11/0.912161250e9 -0.15187033999199e14/0.1634592960000e13; - -0.196663079e9/0.33359040000e11 0.6433463591e10/0.408648240000e12 0.2861399869e10/0.544864320000e12 0.24383341e8/0.2554051500e10 -0.36870830713e11/0.65383718400e11 0.470092704233e12/0.136216080000e12 -0.15187033999199e14/0.1634592960000e13 0.33832994693e11/0.2653560000e10; - ]; - - M4(1:8,1:8) = M4_U; - M4(m-7:m,m-7:m) = rot90(M4_U, 2); - M4 = 1/h^3*M4; - - D4=HI*(M4 - e_l*d3_l'+e_r*d3_r' + d1_l*d2_l'-d1_r*d2_r'); -end
diff -r c0cbffcf6513 -r 72468bc9b63f +sbp/D4Lonely.m --- a/+sbp/D4Lonely.m Mon Sep 26 09:51:45 2016 +0200 +++ b/+sbp/D4Lonely.m Mon Sep 26 09:55:16 2016 +0200 @@ -34,7 +34,7 @@ switch opt case 'min_boundary_points' [H, HI, D4, e_l, e_r, M4, d2_l, d2_r, d3_l, d3_r, d1_l, d1_r] = ... - sbp.implementations.d4_variable_4_min_boundary_points(m, obj.h); + sbp.implementations.d4_lonely_4_min_boundary_points(m, obj.h); % obj.borrowing.N.S2 = 0.5055; % obj.borrowing.N.S3 = 0.9290; otherwise @@ -48,17 +48,17 @@ switch opt case '2' [H, HI, D4, e_l, e_r, M4, d2_l, d2_r, d3_l, d3_r, d1_l, d1_r] = ... - sbp.implementations.d4_variable_6_2(m, obj.h); + sbp.implementations.d4_lonely_6_2(m, obj.h); % obj.borrowing.N.S2 = 0.3259; % obj.borrowing.N.S3 = 0.1580; case '3' [H, HI, D4, e_l, e_r, M4, d2_l, d2_r, d3_l, d3_r, d1_l, d1_r] = ... - sbp.implementations.d4_variable_6_3(m, obj.h); + sbp.implementations.d4_lonely_6_3(m, obj.h); % obj.borrowing.N.S2 = 0.3259; % obj.borrowing.N.S3 = 0.1580; case 'min_boundary_points' [H, HI, D4, e_l, e_r, M4, d2_l, d2_r, d3_l, d3_r, d1_l, d1_r] = ... - sbp.implementations.d4_variable_6_min_boundary_points(m, obj.h); + sbp.implementations.d4_lonely_6_min_boundary_points(m, obj.h); % obj.borrowing.N.S2 = 0.3259; % obj.borrowing.N.S3 = 0.1580; otherwise @@ -72,12 +72,12 @@ switch opt case 'min_boundary_points' [H, HI, D4, e_l, e_r, M4, d2_l, d2_r, d3_l, d3_r, d1_l, d1_r] = ... - sbp.implementations.d4_variable_8_min_boundary_points(m, obj.h); + sbp.implementations.d4_lonely_8_min_boundary_points(m, obj.h); % obj.borrowing.N.S2 = 0.3259; % obj.borrowing.N.S3 = 0.1580; otherwise [H, HI, D4, e_l, e_r, M4, d2_l, d2_r, d3_l, d3_r, d1_l, d1_r] = ... - sbp.implementations.d4_variable_8_higher_boundary_order(m, obj.h); + sbp.implementations.d4_lonely_8_higher_boundary_order(m, obj.h); % obj.borrowing.N.S2 = 0.3259; % obj.borrowing.N.S3 = 0.1580; end