Mercurial > repos > public > sbplib

comparison +sbp/+implementations/d4_variable_2.m @ 886:8894e9c49e40 feature/timesteppers

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Merge with default for latest changes
author Vidar Stiernström <vidar.stiernstrom@it.uu.se>
date Thu, 15 Nov 2018 16:36:21 -0800
parents 43d02533bea3
children
comparison
equal deleted inserted replaced
816:b5e5b195da1e 886:8894e9c49e40
1 % Returns D2 as a function handle
2 function [H, HI, D1, D2, D4, e_l, e_r, M4, d2_l, d2_r, d3_l, d3_r, d1_l, d1_r] = d4_variable_2(m,h)
3 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4 %%% 4:de ordn. SBP Finita differens %%%
5 %%% operatorer framtagna av Ken Mattsson %%%
6 %%% %%%
7 %%% 6 randpunkter, diagonal norm %%%
8 %%% %%%
9 %%% Datum: 2013-11-11 %%%
10 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
11
12 BP = 2;
13 if(m < 2*BP)
14 error('Operator requires at least %d grid points', 2*BP);
15 end
16
17 % Norm
18 Hv = ones(m,1);
19 Hv(1) = 1/2;
20 Hv(m) = 1/2;
21 Hv = h*Hv;
22 H = spdiag(Hv, 0);
23 HI = spdiag(1./Hv, 0);
24
25 % Boundary operators
26 e_l = sparse(m,1);
27 e_l(1) = 1;
28 e_r = rot90(e_l, 2);
29
30 d1_l = sparse(m,1);
31 d1_l(1:3) = 1/h*[-3/2 2 -1/2];
32 d1_r = -rot90(d1_l, 2);
33
34 d2_l = sparse(m,1);
35 d2_l(1:3) = 1/h^2*[1 -2 1];
36 d2_r = rot90(d2_l, 2);
37
38 d3_l = sparse(m,1);
39 d3_l(1:4) = 1/h^3*[-1 3 -3 1];
40 d3_r = -rot90(d3_l, 2);
41
42
43 % First derivative SBP operator, 1st order accurate at first 6 boundary points
44 stencil = [-1/2, 0, 1/2];
45 diags = [-1 0 1];
46 Q = stripeMatrix(stencil, diags, m);
47
48 D1 = HI*(Q - 1/2*e_l*e_l' + 1/2*e_r*e_r');
49
50 % Second derivative, 1st order accurate at first boundary points
51 M = sparse(m,m);
52
53 scheme_width = 3;
54 scheme_radius = (scheme_width-1)/2;
55 r = (1+scheme_radius):(m-scheme_radius);
56
57 function D2 = D2_fun(c)
58 Mm1 = -c(r-1)/2 - c(r)/2;
59 M0 = c(r-1)/2 + c(r) + c(r+1)/2;
60 Mp1 = -c(r)/2 - c(r+1)/2;
61
62 M(r,:) = spdiags([Mm1 M0 Mp1],0:2*scheme_radius,length(r),m);
63
64 M(1:2,1:2) = [c(1)/2 + c(2)/2 -c(1)/2 - c(2)/2; -c(1)/2 - c(2)/2 c(1)/2 + c(2) + c(3)/2;];
65 M(m-1:m,m-1:m) = [c(m-2)/2 + c(m-1) + c(m)/2 -c(m-1)/2 - c(m)/2; -c(m-1)/2 - c(m)/2 c(m-1)/2 + c(m)/2;];
66 M = 1/h*M;
67
68 D2 = HI*(-M - c(1)*e_l*d1_l' + c(m)*e_r*d1_r');
69 end
70 D2 = @D2_fun;
71
72 % Fourth derivative, 0th order accurate at first 6 boundary points
73 stencil = [1, -4, 6, -4, 1];
74 diags = -2:2;
75 M4 = stripeMatrix(stencil, diags, m);
76
77 M4_U = [
78 0.13e2/0.10e2 -0.12e2/0.5e1 0.9e1/0.10e2 0.1e1/0.5e1;
79 -0.12e2/0.5e1 0.26e2/0.5e1 -0.16e2/0.5e1 0.2e1/0.5e1;
80 0.9e1/0.10e2 -0.16e2/0.5e1 0.47e2/0.10e2 -0.17e2/0.5e1;
81 0.1e1/0.5e1 0.2e1/0.5e1 -0.17e2/0.5e1 0.29e2/0.5e1;
82 ];
83
84 M4(1:4,1:4) = M4_U;
85 M4(m-3:m,m-3:m) = rot90(M4_U, 2);
86 M4 = 1/h^3*M4;
87
88 D4=HI*(M4 - e_l*d3_l'+e_r*d3_r' + d1_l*d2_l'-d1_r*d2_r');
89 end