Mercurial > repos > public > sbplib
annotate +sbp/+implementations/d4_lonely_6_2.m @ 1031:2ef20d00b386 feature/advectionRV
For easier comparison, return both the first order and residual viscosity when evaluating the residual. Add the first order and residual viscosity to the state of the RungekuttaRV time steppers
author | Vidar Stiernström <vidar.stiernstrom@it.uu.se> |
---|---|
date | Thu, 17 Jan 2019 10:25:06 +0100 |
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:
321
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_6_2(m,h) |
312 | 2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
3 %%% 6:te ordn. SBP Finita differens %%% | |
4 %%% operatorer med diagonal norm %%% | |
5 %%% Extension to variable koeff %%% | |
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 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% | |
246
fe26791489e0
Added a bunch of new operators. Still non-functional.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff
changeset
|
18 |
312 | 19 % H?r med 6 RP ist?llet f?r 8 f?r D4 operatorn, dock samma randderivator |
20 % Denna ?r noggrannare, och har 2a ordningens randdslutning och b?r ge 6te | |
21 % ordningens konvergens. Hade dock ingen fri parameter att optimera | |
246
fe26791489e0
Added a bunch of new operators. Still non-functional.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff
changeset
|
22 |
321
5c9e5ba1c1ab
Clean up of d4_variable_6_2.m
Jonatan Werpers <jonatan@werpers.com>
parents:
318
diff
changeset
|
23 BP = 6; |
5c9e5ba1c1ab
Clean up of d4_variable_6_2.m
Jonatan Werpers <jonatan@werpers.com>
parents:
318
diff
changeset
|
24 if(m<2*BP) |
5c9e5ba1c1ab
Clean up of d4_variable_6_2.m
Jonatan Werpers <jonatan@werpers.com>
parents:
318
diff
changeset
|
25 error(['Operator requires at least ' num2str(2*BP) ' grid points']); |
5c9e5ba1c1ab
Clean up of d4_variable_6_2.m
Jonatan Werpers <jonatan@werpers.com>
parents:
318
diff
changeset
|
26 end |
246
fe26791489e0
Added a bunch of new operators. Still non-functional.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff
changeset
|
27 |
321
5c9e5ba1c1ab
Clean up of d4_variable_6_2.m
Jonatan Werpers <jonatan@werpers.com>
parents:
318
diff
changeset
|
28 % Norm |
5c9e5ba1c1ab
Clean up of d4_variable_6_2.m
Jonatan Werpers <jonatan@werpers.com>
parents:
318
diff
changeset
|
29 Hv = ones(m,1); |
5c9e5ba1c1ab
Clean up of d4_variable_6_2.m
Jonatan Werpers <jonatan@werpers.com>
parents:
318
diff
changeset
|
30 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]; |
5c9e5ba1c1ab
Clean up of d4_variable_6_2.m
Jonatan Werpers <jonatan@werpers.com>
parents:
318
diff
changeset
|
31 Hv(m-5:m) = rot90(Hv(1:6),2); |
5c9e5ba1c1ab
Clean up of d4_variable_6_2.m
Jonatan Werpers <jonatan@werpers.com>
parents:
318
diff
changeset
|
32 Hv = h*Hv; |
5c9e5ba1c1ab
Clean up of d4_variable_6_2.m
Jonatan Werpers <jonatan@werpers.com>
parents:
318
diff
changeset
|
33 H = spdiag(Hv, 0); |
5c9e5ba1c1ab
Clean up of d4_variable_6_2.m
Jonatan Werpers <jonatan@werpers.com>
parents:
318
diff
changeset
|
34 HI = spdiag(1./Hv, 0); |
246
fe26791489e0
Added a bunch of new operators. Still non-functional.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff
changeset
|
35 |
fe26791489e0
Added a bunch of new operators. Still non-functional.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff
changeset
|
36 |
321
5c9e5ba1c1ab
Clean up of d4_variable_6_2.m
Jonatan Werpers <jonatan@werpers.com>
parents:
318
diff
changeset
|
37 % Boundary operators |
5c9e5ba1c1ab
Clean up of d4_variable_6_2.m
Jonatan Werpers <jonatan@werpers.com>
parents:
318
diff
changeset
|
38 e_l = sparse(m,1); |
5c9e5ba1c1ab
Clean up of d4_variable_6_2.m
Jonatan Werpers <jonatan@werpers.com>
parents:
318
diff
changeset
|
39 e_l(1) = 1; |
5c9e5ba1c1ab
Clean up of d4_variable_6_2.m
Jonatan Werpers <jonatan@werpers.com>
parents:
318
diff
changeset
|
40 e_r = rot90(e_l, 2); |
246
fe26791489e0
Added a bunch of new operators. Still non-functional.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff
changeset
|
41 |
321
5c9e5ba1c1ab
Clean up of d4_variable_6_2.m
Jonatan Werpers <jonatan@werpers.com>
parents:
318
diff
changeset
|
42 d1_l = sparse(m,1); |
5c9e5ba1c1ab
Clean up of d4_variable_6_2.m
Jonatan Werpers <jonatan@werpers.com>
parents:
318
diff
changeset
|
43 d1_l(1:6) = [-0.137e3/0.60e2 5 -5 0.10e2/0.3e1 -0.5e1/0.4e1 0.1e1/0.5e1;]/h; |
326
b19e142fcae1
Fixed bug in setting of boundary derivative.
Jonatan Werpers <jonatan@werpers.com>
parents:
325
diff
changeset
|
44 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
|
45 |
321
5c9e5ba1c1ab
Clean up of d4_variable_6_2.m
Jonatan Werpers <jonatan@werpers.com>
parents:
318
diff
changeset
|
46 d2_l = sparse(m,1); |
5c9e5ba1c1ab
Clean up of d4_variable_6_2.m
Jonatan Werpers <jonatan@werpers.com>
parents:
318
diff
changeset
|
47 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; |
5c9e5ba1c1ab
Clean up of d4_variable_6_2.m
Jonatan Werpers <jonatan@werpers.com>
parents:
318
diff
changeset
|
48 d2_r = rot90(d2_l, 2); |
246
fe26791489e0
Added a bunch of new operators. Still non-functional.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff
changeset
|
49 |
321
5c9e5ba1c1ab
Clean up of d4_variable_6_2.m
Jonatan Werpers <jonatan@werpers.com>
parents:
318
diff
changeset
|
50 d3_l = sparse(m,1); |
5c9e5ba1c1ab
Clean up of d4_variable_6_2.m
Jonatan Werpers <jonatan@werpers.com>
parents:
318
diff
changeset
|
51 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; |
5c9e5ba1c1ab
Clean up of d4_variable_6_2.m
Jonatan Werpers <jonatan@werpers.com>
parents:
318
diff
changeset
|
52 d3_r = -rot90(d3_l, 2); |
246
fe26791489e0
Added a bunch of new operators. Still non-functional.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff
changeset
|
53 |
fe26791489e0
Added a bunch of new operators. Still non-functional.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff
changeset
|
54 |
312 | 55 % Fourth derivative, 1th order accurate at first 8 boundary points (still |
56 % yield 5th order convergence if stable: for example u_tt = -u_xxxx | |
321
5c9e5ba1c1ab
Clean up of d4_variable_6_2.m
Jonatan Werpers <jonatan@werpers.com>
parents:
318
diff
changeset
|
57 stencil = [7/240, -2/5, 169/60, -122/15, 91/8, -122/15, 169/60, -2/5, 7/240]; |
5c9e5ba1c1ab
Clean up of d4_variable_6_2.m
Jonatan Werpers <jonatan@werpers.com>
parents:
318
diff
changeset
|
58 diags = -4:4; |
5c9e5ba1c1ab
Clean up of d4_variable_6_2.m
Jonatan Werpers <jonatan@werpers.com>
parents:
318
diff
changeset
|
59 M4 = stripeMatrix(stencil, diags, m); |
246
fe26791489e0
Added a bunch of new operators. Still non-functional.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff
changeset
|
60 |
312 | 61 M4_U = [ |
62 0.1009e4/0.192e3 -0.7657e4/0.480e3 0.9307e4/0.480e3 -0.509e3/0.40e2 0.4621e4/0.960e3 -0.25e2/0.32e2; | |
63 -0.7657e4/0.480e3 0.49513e5/0.960e3 -0.4007e4/0.60e2 0.21799e5/0.480e3 -0.8171e4/0.480e3 0.2657e4/0.960e3; | |
64 0.9307e4/0.480e3 -0.4007e4/0.60e2 0.1399e4/0.15e2 -0.2721e4/0.40e2 0.12703e5/0.480e3 -0.521e3/0.120e3; | |
65 -0.509e3/0.40e2 0.21799e5/0.480e3 -0.2721e4/0.40e2 0.3349e4/0.60e2 -0.389e3/0.15e2 0.559e3/0.96e2; | |
66 0.4621e4/0.960e3 -0.8171e4/0.480e3 0.12703e5/0.480e3 -0.389e3/0.15e2 0.17857e5/0.960e3 -0.1499e4/0.160e3; | |
67 -0.25e2/0.32e2 0.2657e4/0.960e3 -0.521e3/0.120e3 0.559e3/0.96e2 -0.1499e4/0.160e3 0.2225e4/0.192e3; | |
68 ]; | |
69 | |
246
fe26791489e0
Added a bunch of new operators. Still non-functional.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff
changeset
|
70 |
321
5c9e5ba1c1ab
Clean up of d4_variable_6_2.m
Jonatan Werpers <jonatan@werpers.com>
parents:
318
diff
changeset
|
71 M4(1:6,1:6) = M4_U; |
5c9e5ba1c1ab
Clean up of d4_variable_6_2.m
Jonatan Werpers <jonatan@werpers.com>
parents:
318
diff
changeset
|
72 M4(m-5:m,m-5:m) = rot90(M4_U, 2); |
5c9e5ba1c1ab
Clean up of d4_variable_6_2.m
Jonatan Werpers <jonatan@werpers.com>
parents:
318
diff
changeset
|
73 M4 = 1/h^3*M4; |
246
fe26791489e0
Added a bunch of new operators. Still non-functional.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff
changeset
|
74 |
321
5c9e5ba1c1ab
Clean up of d4_variable_6_2.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 |