Mercurial > repos > public > sbplib

annotate +sbp/+implementations/d4_variable_2.m @ 774:66eb4a2bbb72 feature/grids

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Remove default scaling of the system. The scaling doens't seem to help actual solutions. One example that fails in the flexural code. With large timesteps the solutions seems to blow up. One particular example is profilePresentation on the tdb_presentation_figures branch with k = 0.0005
author Jonatan Werpers <jonatan@werpers.com>
date Wed, 18 Jul 2018 15:42:52 -0700
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