Mercurial > repos > public > sbplib

annotate +sbp/+implementations/d4_lonely_6_2.m @ 1198:2924b3a9b921 feature/d2_compatible

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Add OpSet for fully compatible D2Variable, created from regular D2Variable by replacing d1 by first row of D1. Formal reduction by one order of accuracy at the boundary point.
author Martin Almquist <malmquist@stanford.edu>
date Fri, 16 Aug 2019 14:30:28 -0700
parents b19e142fcae1
children
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
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
9230c056a574 Fixed formatting.
Jonatan Werpers <jonatan@werpers.com>
parents: 310
diff changeset
2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
9230c056a574 Fixed formatting.
Jonatan Werpers <jonatan@werpers.com>
parents: 310
diff changeset
3 %%% 6:te ordn. SBP Finita differens %%%
9230c056a574 Fixed formatting.
Jonatan Werpers <jonatan@werpers.com>
parents: 310
diff changeset
4 %%% operatorer med diagonal norm %%%
9230c056a574 Fixed formatting.
Jonatan Werpers <jonatan@werpers.com>
parents: 310
diff changeset
5 %%% Extension to variable koeff %%%
9230c056a574 Fixed formatting.
Jonatan Werpers <jonatan@werpers.com>
parents: 310
diff changeset
6 %%% %%%
9230c056a574 Fixed formatting.
Jonatan Werpers <jonatan@werpers.com>
parents: 310
diff changeset
7 %%% H (Normen) %%%
9230c056a574 Fixed formatting.
Jonatan Werpers <jonatan@werpers.com>
parents: 310
diff changeset
8 %%% D1=H^(-1)Q (approx f?rsta derivatan) %%%
9230c056a574 Fixed formatting.
Jonatan Werpers <jonatan@werpers.com>
parents: 310
diff changeset
9 %%% D2 (approx andra derivatan) %%%
9230c056a574 Fixed formatting.
Jonatan Werpers <jonatan@werpers.com>
parents: 310
diff changeset
10 %%% D2=HI*(R+C*D*S %%%
9230c056a574 Fixed formatting.
Jonatan Werpers <jonatan@werpers.com>
parents: 310
diff changeset
11 %%% %%%
9230c056a574 Fixed formatting.
Jonatan Werpers <jonatan@werpers.com>
parents: 310
diff changeset
12 %%% R=-D1'*H*C*D1-RR %%%
9230c056a574 Fixed formatting.
Jonatan Werpers <jonatan@werpers.com>
parents: 310
diff changeset
13 %%% %%%
9230c056a574 Fixed formatting.
Jonatan Werpers <jonatan@werpers.com>
parents: 310
diff changeset
14 %%% RR ?r dissipation) %%%
9230c056a574 Fixed formatting.
Jonatan Werpers <jonatan@werpers.com>
parents: 310
diff changeset
15 %%% Dissipationen uppbyggd av D4: %%%
9230c056a574 Fixed formatting.
Jonatan Werpers <jonatan@werpers.com>
parents: 310
diff changeset
16 %%% DI=D4*B*H*D4 %%%
9230c056a574 Fixed formatting.
Jonatan Werpers <jonatan@werpers.com>
parents: 310
diff changeset
17 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
246
fe26791489e0 Added a bunch of new operators. Still non-functional.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
18
312
9230c056a574 Fixed formatting.
Jonatan Werpers <jonatan@werpers.com>
parents: 310
diff changeset
19 % H?r med 6 RP ist?llet f?r 8 f?r D4 operatorn, dock samma randderivator
9230c056a574 Fixed formatting.
Jonatan Werpers <jonatan@werpers.com>
parents: 310
diff changeset
20 % Denna ?r noggrannare, och har 2a ordningens randdslutning och b?r ge 6te
9230c056a574 Fixed formatting.
Jonatan Werpers <jonatan@werpers.com>
parents: 310
diff changeset
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
9230c056a574 Fixed formatting.
Jonatan Werpers <jonatan@werpers.com>
parents: 310
diff changeset
55 % Fourth derivative, 1th order accurate at first 8 boundary points (still
9230c056a574 Fixed formatting.
Jonatan Werpers <jonatan@werpers.com>
parents: 310
diff changeset
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
9230c056a574 Fixed formatting.
Jonatan Werpers <jonatan@werpers.com>
parents: 310
diff changeset
61 M4_U = [
9230c056a574 Fixed formatting.
Jonatan Werpers <jonatan@werpers.com>
parents: 310
diff changeset
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;
9230c056a574 Fixed formatting.
Jonatan Werpers <jonatan@werpers.com>
parents: 310
diff changeset
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;
9230c056a574 Fixed formatting.
Jonatan Werpers <jonatan@werpers.com>
parents: 310
diff changeset
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;
9230c056a574 Fixed formatting.
Jonatan Werpers <jonatan@werpers.com>
parents: 310
diff changeset
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;
9230c056a574 Fixed formatting.
Jonatan Werpers <jonatan@werpers.com>
parents: 310
diff changeset
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;
9230c056a574 Fixed formatting.
Jonatan Werpers <jonatan@werpers.com>
parents: 310
diff changeset
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;
9230c056a574 Fixed formatting.
Jonatan Werpers <jonatan@werpers.com>
parents: 310
diff changeset
68 ];
9230c056a574 Fixed formatting.
Jonatan Werpers <jonatan@werpers.com>
parents: 310
diff changeset
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
9230c056a574 Fixed formatting.
Jonatan Werpers <jonatan@werpers.com>
parents: 310
diff changeset
76 end