Mercurial > repos > public > sbplib

annotate +sbp/+implementations/d4_variable_2.m @ 1347:ac54767ae1fb feature/poroelastic tip

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Add interface, not fully compatible.
author Martin Almquist <martin.almquist@it.uu.se>
date Tue, 30 Apr 2024 14:58:35 +0200
parents 43d02533bea3
children
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
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
9230c056a574 Fixed formatting.
Jonatan Werpers <jonatan@werpers.com>
parents: 311
diff changeset
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
52b4cdf27633 Cleaning order 2.
Jonatan Werpers <jonatan@werpers.com>
parents: 312
diff changeset
13 if(m < 2*BP)
52b4cdf27633 Cleaning order 2.
Jonatan Werpers <jonatan@werpers.com>
parents: 312
diff changeset
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
52b4cdf27633 Cleaning order 2.
Jonatan Werpers <jonatan@werpers.com>
parents: 312
diff changeset
17 % Norm
52b4cdf27633 Cleaning order 2.
Jonatan Werpers <jonatan@werpers.com>
parents: 312
diff changeset
18 Hv = ones(m,1);
52b4cdf27633 Cleaning order 2.
Jonatan Werpers <jonatan@werpers.com>
parents: 312
diff changeset
19 Hv(1) = 1/2;
52b4cdf27633 Cleaning order 2.
Jonatan Werpers <jonatan@werpers.com>
parents: 312
diff changeset
20 Hv(m) = 1/2;
52b4cdf27633 Cleaning order 2.
Jonatan Werpers <jonatan@werpers.com>
parents: 312
diff changeset
21 Hv = h*Hv;
52b4cdf27633 Cleaning order 2.
Jonatan Werpers <jonatan@werpers.com>
parents: 312
diff changeset
22 H = spdiag(Hv, 0);
52b4cdf27633 Cleaning order 2.
Jonatan Werpers <jonatan@werpers.com>
parents: 312
diff changeset
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
52b4cdf27633 Cleaning order 2.
Jonatan Werpers <jonatan@werpers.com>
parents: 312
diff changeset
25 % Boundary operators
52b4cdf27633 Cleaning order 2.
Jonatan Werpers <jonatan@werpers.com>
parents: 312
diff changeset
26 e_l = sparse(m,1);
52b4cdf27633 Cleaning order 2.
Jonatan Werpers <jonatan@werpers.com>
parents: 312
diff changeset
27 e_l(1) = 1;
52b4cdf27633 Cleaning order 2.
Jonatan Werpers <jonatan@werpers.com>
parents: 312
diff changeset
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
52b4cdf27633 Cleaning order 2.
Jonatan Werpers <jonatan@werpers.com>
parents: 312
diff changeset
30 d1_l = sparse(m,1);
52b4cdf27633 Cleaning order 2.
Jonatan Werpers <jonatan@werpers.com>
parents: 312
diff changeset
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
52b4cdf27633 Cleaning order 2.
Jonatan Werpers <jonatan@werpers.com>
parents: 312
diff changeset
33
52b4cdf27633 Cleaning order 2.
Jonatan Werpers <jonatan@werpers.com>
parents: 312
diff changeset
34 d2_l = sparse(m,1);
52b4cdf27633 Cleaning order 2.
Jonatan Werpers <jonatan@werpers.com>
parents: 312
diff changeset
35 d2_l(1:3) = 1/h^2*[1 -2 1];
52b4cdf27633 Cleaning order 2.
Jonatan Werpers <jonatan@werpers.com>
parents: 312
diff changeset
36 d2_r = rot90(d2_l, 2);
52b4cdf27633 Cleaning order 2.
Jonatan Werpers <jonatan@werpers.com>
parents: 312
diff changeset
37
52b4cdf27633 Cleaning order 2.
Jonatan Werpers <jonatan@werpers.com>
parents: 312
diff changeset
38 d3_l = sparse(m,1);
52b4cdf27633 Cleaning order 2.
Jonatan Werpers <jonatan@werpers.com>
parents: 312
diff changeset
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
52b4cdf27633 Cleaning order 2.
Jonatan Werpers <jonatan@werpers.com>
parents: 312
diff changeset
44 stencil = [-1/2, 0, 1/2];
52b4cdf27633 Cleaning order 2.
Jonatan Werpers <jonatan@werpers.com>
parents: 312
diff changeset
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
52b4cdf27633 Cleaning order 2.
Jonatan Werpers <jonatan@werpers.com>
parents: 312
diff changeset
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
52b4cdf27633 Cleaning order 2.
Jonatan Werpers <jonatan@werpers.com>
parents: 312
diff changeset
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;];
52b4cdf27633 Cleaning order 2.
Jonatan Werpers <jonatan@werpers.com>
parents: 312
diff changeset
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;];
52b4cdf27633 Cleaning order 2.
Jonatan Werpers <jonatan@werpers.com>
parents: 312
diff changeset
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
52b4cdf27633 Cleaning order 2.
Jonatan Werpers <jonatan@werpers.com>
parents: 312
diff changeset
72 % Fourth derivative, 0th order accurate at first 6 boundary points
52b4cdf27633 Cleaning order 2.
Jonatan Werpers <jonatan@werpers.com>
parents: 312
diff changeset
73 stencil = [1, -4, 6, -4, 1];
52b4cdf27633 Cleaning order 2.
Jonatan Werpers <jonatan@werpers.com>
parents: 312
diff changeset
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
52b4cdf27633 Cleaning order 2.
Jonatan Werpers <jonatan@werpers.com>
parents: 312
diff changeset
77 M4_U = [
52b4cdf27633 Cleaning order 2.
Jonatan Werpers <jonatan@werpers.com>
parents: 312
diff changeset
78 0.13e2/0.10e2 -0.12e2/0.5e1 0.9e1/0.10e2 0.1e1/0.5e1;
52b4cdf27633 Cleaning order 2.
Jonatan Werpers <jonatan@werpers.com>
parents: 312
diff changeset
79 -0.12e2/0.5e1 0.26e2/0.5e1 -0.16e2/0.5e1 0.2e1/0.5e1;
52b4cdf27633 Cleaning order 2.
Jonatan Werpers <jonatan@werpers.com>
parents: 312
diff changeset
80 0.9e1/0.10e2 -0.16e2/0.5e1 0.47e2/0.10e2 -0.17e2/0.5e1;
52b4cdf27633 Cleaning order 2.
Jonatan Werpers <jonatan@werpers.com>
parents: 312
diff changeset
81 0.1e1/0.5e1 0.2e1/0.5e1 -0.17e2/0.5e1 0.29e2/0.5e1;
312
9230c056a574 Fixed formatting.
Jonatan Werpers <jonatan@werpers.com>
parents: 311
diff changeset
82 ];
261
6009f2712d13 Moved and renamned all implementations.
Martin Almquist <martin.almquist@it.uu.se>
parents:
diff changeset
83
313
52b4cdf27633 Cleaning order 2.
Jonatan Werpers <jonatan@werpers.com>
parents: 312
diff changeset
84 M4(1:4,1:4) = M4_U;
52b4cdf27633 Cleaning order 2.
Jonatan Werpers <jonatan@werpers.com>
parents: 312
diff changeset
85 M4(m-3:m,m-3:m) = rot90(M4_U, 2);
52b4cdf27633 Cleaning order 2.
Jonatan Werpers <jonatan@werpers.com>
parents: 312
diff changeset
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