Mercurial > repos > public > sbplib

comparison +sbp/D1StaggeredUpwind.m @ 1261:a3f2c1781612 feature/poroelastic

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Fix bug with incorrect grids returned by D1StaggeredUpwind
author Martin Almquist <malmquist@stanford.edu>
date Wed, 29 Apr 2020 21:02:05 -0700
parents 99f92bfc1157
children
comparison
equal deleted inserted replaced
1260:7fbee3de0ab2 1261:a3f2c1781612
8 % D1_dual takes FROM primal grid TO dual grid 8 % D1_dual takes FROM primal grid TO dual grid
9 9
10 D1_primal % SBP operator approximating first derivative 10 D1_primal % SBP operator approximating first derivative
11 D1_dual % SBP operator approximating first derivative 11 D1_dual % SBP operator approximating first derivative
12 12
13 Dplus_primal % Upwind operator on primal grid 13 Dplus_primal % Upwind operator on primal grid
14 Dminus_primal % Upwind operator on primal grid 14 Dminus_primal % Upwind operator on primal grid
15 Dplus_dual % Upwind operator on dual grid 15 Dplus_dual % Upwind operator on dual grid
16 Dminus_dual % Upwind operator on dual grid 16 Dminus_dual % Upwind operator on dual grid
17 17
18 H_primal % Norm matrix 18 H_primal % Norm matrix
19 H_dual % Norm matrix 19 H_dual % Norm matrix
20 H_primalI % H^-1 20 H_primalI % H^-1
21 H_dualI % H^-1 21 H_dualI % H^-1
34 end 34 end
35 35
36 methods 36 methods
37 function obj = D1StaggeredUpwind(m,lim,order) 37 function obj = D1StaggeredUpwind(m,lim,order)
38 38
39 x_l = lim{1}; 39 xl = lim{1};
40 x_r = lim{2}; 40 xr = lim{2};
41 L = x_r-x_l; 41 L = xr-xl;
42 h = L/(m-1);
42 43
43 m_primal = m; 44 m_primal = m;
44 m_dual = m+1; 45 m_dual = m+1;
45 46
46 switch order 47 switch order
47 case 2 48 case 2
48 [obj.x_primal, obj.x_dual, obj.H_primal, obj.H_dual,... 49 [~, ~, obj.H_primal, obj.H_dual,...
49 obj.H_primalI, obj.H_dualI,... 50 obj.H_primalI, obj.H_dualI,...
50 obj.D1_primal, obj.D1_dual, obj.Dplus_primal, obj.Dminus_primal,... 51 obj.D1_primal, obj.D1_dual, obj.Dplus_primal, obj.Dminus_primal,...
51 obj.Dplus_dual, obj.Dminus_dual] = sbp.implementations.d1_staggered_upwind_2(m, L); 52 obj.Dplus_dual, obj.Dminus_dual] = sbp.implementations.d1_staggered_upwind_2(m, L);
52 case 4 53 case 4
53 [obj.x_primal, obj.x_dual, obj.H_primal, obj.H_dual,... 54 [~, ~, obj.H_primal, obj.H_dual,...
54 obj.H_primalI, obj.H_dualI,... 55 obj.H_primalI, obj.H_dualI,...
55 obj.D1_primal, obj.D1_dual, obj.Dplus_primal, obj.Dminus_primal,... 56 obj.D1_primal, obj.D1_dual, obj.Dplus_primal, obj.Dminus_primal,...
56 obj.Dplus_dual, obj.Dminus_dual] = sbp.implementations.d1_staggered_upwind_4(m, L); 57 obj.Dplus_dual, obj.Dminus_dual] = sbp.implementations.d1_staggered_upwind_4(m, L);
57 case 6 58 case 6
58 [obj.x_primal, obj.x_dual, obj.H_primal, obj.H_dual,... 59 [~, ~, obj.H_primal, obj.H_dual,...
59 obj.H_primalI, obj.H_dualI,... 60 obj.H_primalI, obj.H_dualI,...
60 obj.D1_primal, obj.D1_dual, obj.Dplus_primal, obj.Dminus_primal,... 61 obj.D1_primal, obj.D1_dual, obj.Dplus_primal, obj.Dminus_primal,...
61 obj.Dplus_dual, obj.Dminus_dual] = sbp.implementations.d1_staggered_upwind_6(m, L); 62 obj.Dplus_dual, obj.Dminus_dual] = sbp.implementations.d1_staggered_upwind_6(m, L);
62 case 8 63 case 8
63 [obj.x_primal, obj.x_dual, obj.H_primal, obj.H_dual,... 64 [~, ~, obj.H_primal, obj.H_dual,...
64 obj.H_primalI, obj.H_dualI,... 65 obj.H_primalI, obj.H_dualI,...
65 obj.D1_primal, obj.D1_dual, obj.Dplus_primal, obj.Dminus_primal,... 66 obj.D1_primal, obj.D1_dual, obj.Dplus_primal, obj.Dminus_primal,...
66 obj.Dplus_dual, obj.Dminus_dual] = sbp.implementations.d1_staggered_upwind_8(m, L); 67 obj.Dplus_dual, obj.Dminus_dual] = sbp.implementations.d1_staggered_upwind_8(m, L);
67 otherwise 68 otherwise
68 error('Invalid operator order %d.',order); 69 error('Invalid operator order %d.',order);
69 end 70 end
70 71
71 obj.m = m; 72 obj.m = m;
72 obj.m_primal = m_primal; 73 obj.m_primal = m_primal;
73 obj.m_dual = m_dual; 74 obj.m_dual = m_dual;
74 obj.h = L/(m-1); 75 obj.h = h;
76
77 obj.x_primal = linspace(xl, xr, m)';
78 obj.x_dual = [xl, linspace(xl+h/2, xr-h/2, m-1), xr]';
75 79
76 obj.e_primal_l = sparse(m_primal,1); 80 obj.e_primal_l = sparse(m_primal,1);
77 obj.e_primal_r = sparse(m_primal,1); 81 obj.e_primal_r = sparse(m_primal,1);
78 obj.e_primal_l(1) = 1; 82 obj.e_primal_l(1) = 1;
79 obj.e_primal_r(m_primal) = 1; 83 obj.e_primal_r(m_primal) = 1;