Mercurial > repos > public > sbplib
annotate +sbp/+implementations/d4_variable_2.m @ 1037:2d7ba44340d0 feature/burgers1d
Pass scheme specific parameters as cell array. This will enabale constructDiffOps to be more general. In addition, allow for schemes returning function handles as diffOps, which is currently how non-linear schemes such as Burgers1d are implemented.
author | Vidar Stiernström <vidar.stiernstrom@it.uu.se> |
---|---|
date | Fri, 18 Jan 2019 09:02:02 +0100 |
parents | 43d02533bea3 |
children |
rev | line source |
---|---|
261
6009f2712d13
Moved and renamned all implementations.
Martin Almquist <martin.almquist@it.uu.se>
parents:
diff
changeset
|
1 % Returns D2 as a function handle |
312 | 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) |
261
6009f2712d13
Moved and renamned all implementations.
Martin Almquist <martin.almquist@it.uu.se>
parents:
diff
changeset
|
3 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
6009f2712d13
Moved and renamned all implementations.
Martin Almquist <martin.almquist@it.uu.se>
parents:
diff
changeset
|
4 %%% 4:de ordn. SBP Finita differens %%% |
6009f2712d13
Moved and renamned all implementations.
Martin Almquist <martin.almquist@it.uu.se>
parents:
diff
changeset
|
5 %%% operatorer framtagna av Ken Mattsson %%% |
6009f2712d13
Moved and renamned all implementations.
Martin Almquist <martin.almquist@it.uu.se>
parents:
diff
changeset
|
6 %%% %%% |
6009f2712d13
Moved and renamned all implementations.
Martin Almquist <martin.almquist@it.uu.se>
parents:
diff
changeset
|
7 %%% 6 randpunkter, diagonal norm %%% |
6009f2712d13
Moved and renamned all implementations.
Martin Almquist <martin.almquist@it.uu.se>
parents:
diff
changeset
|
8 %%% %%% |
6009f2712d13
Moved and renamned all implementations.
Martin Almquist <martin.almquist@it.uu.se>
parents:
diff
changeset
|
9 %%% Datum: 2013-11-11 %%% |
6009f2712d13
Moved and renamned all implementations.
Martin Almquist <martin.almquist@it.uu.se>
parents:
diff
changeset
|
10 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
6009f2712d13
Moved and renamned all implementations.
Martin Almquist <martin.almquist@it.uu.se>
parents:
diff
changeset
|
11 |
341
43d02533bea3
Fixed error in boundary point check.
Jonatan Werpers <jonatan@werpers.com>
parents:
329
diff
changeset
|
12 BP = 2; |
313 | 13 if(m < 2*BP) |
14 error('Operator requires at least %d grid points', 2*BP); | |
266
bfa130b7abf6
Added error message for too few grid points to all implementation files.
Martin Almquist <martin.almquist@it.uu.se>
parents:
261
diff
changeset
|
15 end |
261
6009f2712d13
Moved and renamned all implementations.
Martin Almquist <martin.almquist@it.uu.se>
parents:
diff
changeset
|
16 |
313 | 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); | |
261
6009f2712d13
Moved and renamned all implementations.
Martin Almquist <martin.almquist@it.uu.se>
parents:
diff
changeset
|
24 |
313 | 25 % Boundary operators |
26 e_l = sparse(m,1); | |
27 e_l(1) = 1; | |
28 e_r = rot90(e_l, 2); | |
261
6009f2712d13
Moved and renamned all implementations.
Martin Almquist <martin.almquist@it.uu.se>
parents:
diff
changeset
|
29 |
313 | 30 d1_l = sparse(m,1); |
31 d1_l(1:3) = 1/h*[-3/2 2 -1/2]; | |
326
b19e142fcae1
Fixed bug in setting of boundary derivative.
Jonatan Werpers <jonatan@werpers.com>
parents:
314
diff
changeset
|
32 d1_r = -rot90(d1_l, 2); |
313 | 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]; | |
314
88584b0cfba1
Corrections and clean up order 4
Jonatan Werpers <jonatan@werpers.com>
parents:
313
diff
changeset
|
40 d3_r = -rot90(d3_l, 2); |
261
6009f2712d13
Moved and renamned all implementations.
Martin Almquist <martin.almquist@it.uu.se>
parents:
diff
changeset
|
41 |
6009f2712d13
Moved and renamned all implementations.
Martin Almquist <martin.almquist@it.uu.se>
parents:
diff
changeset
|
42 |
6009f2712d13
Moved and renamned all implementations.
Martin Almquist <martin.almquist@it.uu.se>
parents:
diff
changeset
|
43 % First derivative SBP operator, 1st order accurate at first 6 boundary points |
313 | 44 stencil = [-1/2, 0, 1/2]; |
45 diags = [-1 0 1]; | |
267
f7ac3cd6eeaa
Sparsified all implementation files, removed all matlab warnings, fixed small bugs on minimum grid points.
Martin Almquist <martin.almquist@it.uu.se>
parents:
266
diff
changeset
|
46 Q = stripeMatrix(stencil, diags, m); |
261
6009f2712d13
Moved and renamned all implementations.
Martin Almquist <martin.almquist@it.uu.se>
parents:
diff
changeset
|
47 |
314
88584b0cfba1
Corrections and clean up order 4
Jonatan Werpers <jonatan@werpers.com>
parents:
313
diff
changeset
|
48 D1 = HI*(Q - 1/2*e_l*e_l' + 1/2*e_r*e_r'); |
261
6009f2712d13
Moved and renamned all implementations.
Martin Almquist <martin.almquist@it.uu.se>
parents:
diff
changeset
|
49 |
6009f2712d13
Moved and renamned all implementations.
Martin Almquist <martin.almquist@it.uu.se>
parents:
diff
changeset
|
50 % Second derivative, 1st order accurate at first boundary points |
313 | 51 M = sparse(m,m); |
261
6009f2712d13
Moved and renamned all implementations.
Martin Almquist <martin.almquist@it.uu.se>
parents:
diff
changeset
|
52 |
6009f2712d13
Moved and renamned all implementations.
Martin Almquist <martin.almquist@it.uu.se>
parents:
diff
changeset
|
53 scheme_width = 3; |
6009f2712d13
Moved and renamned all implementations.
Martin Almquist <martin.almquist@it.uu.se>
parents:
diff
changeset
|
54 scheme_radius = (scheme_width-1)/2; |
6009f2712d13
Moved and renamned all implementations.
Martin Almquist <martin.almquist@it.uu.se>
parents:
diff
changeset
|
55 r = (1+scheme_radius):(m-scheme_radius); |
6009f2712d13
Moved and renamned all implementations.
Martin Almquist <martin.almquist@it.uu.se>
parents:
diff
changeset
|
56 |
6009f2712d13
Moved and renamned all implementations.
Martin Almquist <martin.almquist@it.uu.se>
parents:
diff
changeset
|
57 function D2 = D2_fun(c) |
6009f2712d13
Moved and renamned all implementations.
Martin Almquist <martin.almquist@it.uu.se>
parents:
diff
changeset
|
58 Mm1 = -c(r-1)/2 - c(r)/2; |
6009f2712d13
Moved and renamned all implementations.
Martin Almquist <martin.almquist@it.uu.se>
parents:
diff
changeset
|
59 M0 = c(r-1)/2 + c(r) + c(r+1)/2; |
6009f2712d13
Moved and renamned all implementations.
Martin Almquist <martin.almquist@it.uu.se>
parents:
diff
changeset
|
60 Mp1 = -c(r)/2 - c(r+1)/2; |
6009f2712d13
Moved and renamned all implementations.
Martin Almquist <martin.almquist@it.uu.se>
parents:
diff
changeset
|
61 |
6009f2712d13
Moved and renamned all implementations.
Martin Almquist <martin.almquist@it.uu.se>
parents:
diff
changeset
|
62 M(r,:) = spdiags([Mm1 M0 Mp1],0:2*scheme_radius,length(r),m); |
6009f2712d13
Moved and renamned all implementations.
Martin Almquist <martin.almquist@it.uu.se>
parents:
diff
changeset
|
63 |
313 | 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; | |
261
6009f2712d13
Moved and renamned all implementations.
Martin Almquist <martin.almquist@it.uu.se>
parents:
diff
changeset
|
67 |
329
bf801c3709be
Bug fixes in operators.
Jonatan Werpers <jonatan@werpers.com>
parents:
326
diff
changeset
|
68 D2 = HI*(-M - c(1)*e_l*d1_l' + c(m)*e_r*d1_r'); |
261
6009f2712d13
Moved and renamned all implementations.
Martin Almquist <martin.almquist@it.uu.se>
parents:
diff
changeset
|
69 end |
6009f2712d13
Moved and renamned all implementations.
Martin Almquist <martin.almquist@it.uu.se>
parents:
diff
changeset
|
70 D2 = @D2_fun; |
6009f2712d13
Moved and renamned all implementations.
Martin Almquist <martin.almquist@it.uu.se>
parents:
diff
changeset
|
71 |
313 | 72 % Fourth derivative, 0th order accurate at first 6 boundary points |
73 stencil = [1, -4, 6, -4, 1]; | |
74 diags = -2:2; | |
267
f7ac3cd6eeaa
Sparsified all implementation files, removed all matlab warnings, fixed small bugs on minimum grid points.
Martin Almquist <martin.almquist@it.uu.se>
parents:
266
diff
changeset
|
75 M4 = stripeMatrix(stencil, diags, m); |
261
6009f2712d13
Moved and renamned all implementations.
Martin Almquist <martin.almquist@it.uu.se>
parents:
diff
changeset
|
76 |
313 | 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; | |
312 | 82 ]; |
261
6009f2712d13
Moved and renamned all implementations.
Martin Almquist <martin.almquist@it.uu.se>
parents:
diff
changeset
|
83 |
313 | 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; | |
261
6009f2712d13
Moved and renamned all implementations.
Martin Almquist <martin.almquist@it.uu.se>
parents:
diff
changeset
|
87 |
314
88584b0cfba1
Corrections and clean up order 4
Jonatan Werpers <jonatan@werpers.com>
parents:
313
diff
changeset
|
88 D4=HI*(M4 - e_l*d3_l'+e_r*d3_r' + d1_l*d2_l'-d1_r*d2_r'); |
261
6009f2712d13
Moved and renamned all implementations.
Martin Almquist <martin.almquist@it.uu.se>
parents:
diff
changeset
|
89 end |