Mercurial > repos > public > sbplib
annotate +sbp/+implementations/d4_lonely_8_higher_boundary_order.m @ 329:bf801c3709be feature/beams
Bug fixes in operators.
author | Jonatan Werpers <jonatan@werpers.com> |
---|---|
date | Wed, 28 Sep 2016 16:31:06 +0200 |
parents | b19e142fcae1 |
children |
rev | line source |
---|---|
324
c0cbffcf6513
Created new operator class for lonely D4 operators. Removed some output parameters of implementations.
Jonatan Werpers <jonatan@werpers.com>
parents:
323
diff
changeset
|
1 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) |
312 | 2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
3 %%% 8:te ordn. SBP Finita differens %%% | |
4 %%% operatorer med diagonal norm %%% | |
5 %%% %%% | |
6 %%% %%% | |
7 %%% H (Normen) %%% | |
8 %%% D1=H^(-1)Q (approx f?rsta derivatan) %%% | |
9 %%% D2 (approx andra derivatan) %%% | |
10 %%% D2=HI*(R+C*D*S %%% | |
11 %%% %%% | |
12 %%% R=-D1'*H*C*D1-RR %%% | |
13 %%% %%% | |
14 %%% RR ?r dissipation) %%% | |
15 %%% Dissipationen uppbyggd av D4: %%% | |
16 %%% DI=D4*B*H*D4 %%% | |
17 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% | |
18 %This is 3rd order accurate at the boundary. Not same norm as D1 operator | |
246
fe26791489e0
Added a bunch of new operators. Still non-functional.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff
changeset
|
19 |
323
7579c2abbf9f
Cleaning d4_variable_8_higher_boundary_order.m
Jonatan Werpers <jonatan@werpers.com>
parents:
318
diff
changeset
|
20 BP = 8; |
7579c2abbf9f
Cleaning d4_variable_8_higher_boundary_order.m
Jonatan Werpers <jonatan@werpers.com>
parents:
318
diff
changeset
|
21 if(m<2*BP) |
7579c2abbf9f
Cleaning d4_variable_8_higher_boundary_order.m
Jonatan Werpers <jonatan@werpers.com>
parents:
318
diff
changeset
|
22 error(['Operator requires at least ' num2str(2*BP) ' grid points']); |
7579c2abbf9f
Cleaning d4_variable_8_higher_boundary_order.m
Jonatan Werpers <jonatan@werpers.com>
parents:
318
diff
changeset
|
23 end |
312 | 24 |
323
7579c2abbf9f
Cleaning d4_variable_8_higher_boundary_order.m
Jonatan Werpers <jonatan@werpers.com>
parents:
318
diff
changeset
|
25 % Norm |
7579c2abbf9f
Cleaning d4_variable_8_higher_boundary_order.m
Jonatan Werpers <jonatan@werpers.com>
parents:
318
diff
changeset
|
26 Hv = ones(m,1); |
7579c2abbf9f
Cleaning d4_variable_8_higher_boundary_order.m
Jonatan Werpers <jonatan@werpers.com>
parents:
318
diff
changeset
|
27 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]; |
7579c2abbf9f
Cleaning d4_variable_8_higher_boundary_order.m
Jonatan Werpers <jonatan@werpers.com>
parents:
318
diff
changeset
|
28 Hv(m-7:m) = rot90(Hv(1:8),2); |
7579c2abbf9f
Cleaning d4_variable_8_higher_boundary_order.m
Jonatan Werpers <jonatan@werpers.com>
parents:
318
diff
changeset
|
29 Hv = h*Hv; |
7579c2abbf9f
Cleaning d4_variable_8_higher_boundary_order.m
Jonatan Werpers <jonatan@werpers.com>
parents:
318
diff
changeset
|
30 H = spdiag(Hv, 0); |
7579c2abbf9f
Cleaning d4_variable_8_higher_boundary_order.m
Jonatan Werpers <jonatan@werpers.com>
parents:
318
diff
changeset
|
31 HI = spdiag(1./Hv, 0); |
312 | 32 |
33 | |
323
7579c2abbf9f
Cleaning d4_variable_8_higher_boundary_order.m
Jonatan Werpers <jonatan@werpers.com>
parents:
318
diff
changeset
|
34 % Boundary operators |
7579c2abbf9f
Cleaning d4_variable_8_higher_boundary_order.m
Jonatan Werpers <jonatan@werpers.com>
parents:
318
diff
changeset
|
35 e_l = sparse(m,1); |
7579c2abbf9f
Cleaning d4_variable_8_higher_boundary_order.m
Jonatan Werpers <jonatan@werpers.com>
parents:
318
diff
changeset
|
36 e_l(1) = 1; |
7579c2abbf9f
Cleaning d4_variable_8_higher_boundary_order.m
Jonatan Werpers <jonatan@werpers.com>
parents:
318
diff
changeset
|
37 e_r = rot90(e_l, 2); |
7579c2abbf9f
Cleaning d4_variable_8_higher_boundary_order.m
Jonatan Werpers <jonatan@werpers.com>
parents:
318
diff
changeset
|
38 |
7579c2abbf9f
Cleaning d4_variable_8_higher_boundary_order.m
Jonatan Werpers <jonatan@werpers.com>
parents:
318
diff
changeset
|
39 d1_l = sparse(m,1); |
7579c2abbf9f
Cleaning d4_variable_8_higher_boundary_order.m
Jonatan Werpers <jonatan@werpers.com>
parents:
318
diff
changeset
|
40 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; |
326
b19e142fcae1
Fixed bug in setting of boundary derivative.
Jonatan Werpers <jonatan@werpers.com>
parents:
325
diff
changeset
|
41 d1_r = -rot90(d1_l, 2); |
246
fe26791489e0
Added a bunch of new operators. Still non-functional.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff
changeset
|
42 |
323
7579c2abbf9f
Cleaning d4_variable_8_higher_boundary_order.m
Jonatan Werpers <jonatan@werpers.com>
parents:
318
diff
changeset
|
43 d2_l = sparse(m,1); |
7579c2abbf9f
Cleaning d4_variable_8_higher_boundary_order.m
Jonatan Werpers <jonatan@werpers.com>
parents:
318
diff
changeset
|
44 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; |
7579c2abbf9f
Cleaning d4_variable_8_higher_boundary_order.m
Jonatan Werpers <jonatan@werpers.com>
parents:
318
diff
changeset
|
45 d2_r = rot90(d2_l, 2); |
316
203afa156f59
Collected boundary operators.
Jonatan Werpers <jonatan@werpers.com>
parents:
312
diff
changeset
|
46 |
323
7579c2abbf9f
Cleaning d4_variable_8_higher_boundary_order.m
Jonatan Werpers <jonatan@werpers.com>
parents:
318
diff
changeset
|
47 d3_l = sparse(m,1); |
7579c2abbf9f
Cleaning d4_variable_8_higher_boundary_order.m
Jonatan Werpers <jonatan@werpers.com>
parents:
318
diff
changeset
|
48 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; |
7579c2abbf9f
Cleaning d4_variable_8_higher_boundary_order.m
Jonatan Werpers <jonatan@werpers.com>
parents:
318
diff
changeset
|
49 d3_r = -rot90(d3_l, 2); |
316
203afa156f59
Collected boundary operators.
Jonatan Werpers <jonatan@werpers.com>
parents:
312
diff
changeset
|
50 |
323
7579c2abbf9f
Cleaning d4_variable_8_higher_boundary_order.m
Jonatan Werpers <jonatan@werpers.com>
parents:
318
diff
changeset
|
51 |
246
fe26791489e0
Added a bunch of new operators. Still non-functional.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff
changeset
|
52 |
312 | 53 % Fourth derivative, 1th order accurate at first 8 boundary points (still |
54 % yield 5th order convergence if stable: for example u_tt = -u_xxxx | |
246
fe26791489e0
Added a bunch of new operators. Still non-functional.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff
changeset
|
55 |
323
7579c2abbf9f
Cleaning d4_variable_8_higher_boundary_order.m
Jonatan Werpers <jonatan@werpers.com>
parents:
318
diff
changeset
|
56 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]; |
7579c2abbf9f
Cleaning d4_variable_8_higher_boundary_order.m
Jonatan Werpers <jonatan@werpers.com>
parents:
318
diff
changeset
|
57 diags = -5:5; |
329
bf801c3709be
Bug fixes in operators.
Jonatan Werpers <jonatan@werpers.com>
parents:
326
diff
changeset
|
58 M4 = stripeMatrix(stencil, diags, m); |
246
fe26791489e0
Added a bunch of new operators. Still non-functional.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff
changeset
|
59 |
312 | 60 M4_U = [ |
61 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; | |
62 -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; | |
63 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; | |
64 -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; | |
65 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; | |
66 -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; | |
67 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; | |
68 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; | |
69 ]; | |
246
fe26791489e0
Added a bunch of new operators. Still non-functional.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff
changeset
|
70 |
312 | 71 M4(1:8,1:8) = M4_U; |
323
7579c2abbf9f
Cleaning d4_variable_8_higher_boundary_order.m
Jonatan Werpers <jonatan@werpers.com>
parents:
318
diff
changeset
|
72 M4(m-7:m,m-7:m) = rot90(M4_U, 2); |
7579c2abbf9f
Cleaning d4_variable_8_higher_boundary_order.m
Jonatan Werpers <jonatan@werpers.com>
parents:
318
diff
changeset
|
73 M4 = 1/h^3*M4; |
312 | 74 |
323
7579c2abbf9f
Cleaning d4_variable_8_higher_boundary_order.m
Jonatan Werpers <jonatan@werpers.com>
parents:
318
diff
changeset
|
75 D4=HI*(M4 - e_l*d3_l'+e_r*d3_r' + d1_l*d2_l'-d1_r*d2_r'); |
312 | 76 end |