Mercurial > repos > public > sbplib

comparison +sbp/D1StaggeredUpwind.m @ 643:62147d8ce062 feature/d1_staggered

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Add staggered-upwind operators by Ken and Ossian.
author Martin Almquist <malmquist@stanford.edu>
date Mon, 13 Nov 2017 11:49:08 -0800
parents
children a3f2c1781612
comparison
equal deleted inserted replaced
642:8f0ee6ba6313 643:62147d8ce062
1 classdef D1StaggeredUpwind < sbp.OpSet
2 % Compatible staggered and upwind operators by Ken Mattsson and Ossian O'reilly
3 properties
4 % x_primal: "primal" grid with m points. Equidistant. Called Plus grid in Ossian's paper.
5 % x_dual: "dual" grid with m+1 points. Called Minus grid in Ossian's paper.
6
7 % D1_primal takes FROM dual grid TO primal grid
8 % D1_dual takes FROM primal grid TO dual grid
9
10 D1_primal % SBP operator approximating first derivative
11 D1_dual % SBP operator approximating first derivative
12
13 Dplus_primal % Upwind operator on primal grid
14 Dminus_primal % Upwind operator on primal grid
15 Dplus_dual % Upwind operator on dual grid
16 Dminus_dual % Upwind operator on dual grid
17
18 H_primal % Norm matrix
19 H_dual % Norm matrix
20 H_primalI % H^-1
21 H_dualI % H^-1
22 e_primal_l % Left boundary operator
23 e_dual_l % Left boundary operator
24 e_primal_r % Right boundary operator
25 e_dual_r % Right boundary operator
26 m % Number of grid points.
27 m_primal % Number of grid points.
28 m_dual % Number of grid points.
29 h % Step size
30 x_primal % grid
31 x_dual % grid
32 x
33 borrowing % Struct with borrowing limits for different norm matrices
34 end
35
36 methods
37 function obj = D1StaggeredUpwind(m,lim,order)
38
39 x_l = lim{1};
40 x_r = lim{2};
41 L = x_r-x_l;
42
43 m_primal = m;
44 m_dual = m+1;
45
46 switch order
47 case 2
48 [obj.x_primal, obj.x_dual, obj.H_primal, obj.H_dual,...
49 obj.H_primalI, obj.H_dualI,...
50 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 case 4
53 [obj.x_primal, obj.x_dual, obj.H_primal, obj.H_dual,...
54 obj.H_primalI, obj.H_dualI,...
55 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 case 6
58 [obj.x_primal, obj.x_dual, obj.H_primal, obj.H_dual,...
59 obj.H_primalI, obj.H_dualI,...
60 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 case 8
63 [obj.x_primal, obj.x_dual, obj.H_primal, obj.H_dual,...
64 obj.H_primalI, obj.H_dualI,...
65 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 otherwise
68 error('Invalid operator order %d.',order);
69 end
70
71 obj.m = m;
72 obj.m_primal = m_primal;
73 obj.m_dual = m_dual;
74 obj.h = L/(m-1);
75
76 obj.e_primal_l = sparse(m_primal,1);
77 obj.e_primal_r = sparse(m_primal,1);
78 obj.e_primal_l(1) = 1;
79 obj.e_primal_r(m_primal) = 1;
80
81 obj.e_dual_l = sparse(m_dual,1);
82 obj.e_dual_r = sparse(m_dual,1);
83 obj.e_dual_l(1) = 1;
84 obj.e_dual_r(m_dual) = 1;
85
86 obj.borrowing = [];
87 obj.x = [];
88
89 end
90
91 function str = string(obj)
92 str = [class(obj) '_' num2str(obj.order)];
93 end
94 end
95 end