Mercurial > repos > public > sbplib

annotate +scheme/Staggered1DAcoustics.m @ 671:993aac771efd feature/d1_staggered

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Correct characteristic BC in staggered acoustics 1D constant coeff.
author Martin Almquist <malmquist@stanford.edu>
date Wed, 27 Dec 2017 17:18:51 +0100
parents 8e55298657b9
children ec0ac76006e2
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
642
8f0ee6ba6313 Add scheme Staggered1DAcoustics
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
1 classdef Staggered1DAcoustics < scheme.Scheme
8f0ee6ba6313 Add scheme Staggered1DAcoustics
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
2 properties
8f0ee6ba6313 Add scheme Staggered1DAcoustics
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
3 m % Number of points of primal grid in each direction, possibly a vector
8f0ee6ba6313 Add scheme Staggered1DAcoustics
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
4 h % Grid spacing
8f0ee6ba6313 Add scheme Staggered1DAcoustics
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
5
8f0ee6ba6313 Add scheme Staggered1DAcoustics
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
6 % Grids
8f0ee6ba6313 Add scheme Staggered1DAcoustics
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
7 grid % Total grid object
8f0ee6ba6313 Add scheme Staggered1DAcoustics
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
8 grid_primal
8f0ee6ba6313 Add scheme Staggered1DAcoustics
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
9 grid_dual
8f0ee6ba6313 Add scheme Staggered1DAcoustics
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
10
8f0ee6ba6313 Add scheme Staggered1DAcoustics
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
11 order % Order accuracy for the approximation
8f0ee6ba6313 Add scheme Staggered1DAcoustics
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
12
8f0ee6ba6313 Add scheme Staggered1DAcoustics
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
13 H % Combined norm
8f0ee6ba6313 Add scheme Staggered1DAcoustics
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
14 Hi % Inverse
8f0ee6ba6313 Add scheme Staggered1DAcoustics
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
15 D % Semi-discrete approximation matrix
8f0ee6ba6313 Add scheme Staggered1DAcoustics
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
16
8f0ee6ba6313 Add scheme Staggered1DAcoustics
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
17 D1_primal
8f0ee6ba6313 Add scheme Staggered1DAcoustics
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
18 D1_dual
8f0ee6ba6313 Add scheme Staggered1DAcoustics
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
19
8f0ee6ba6313 Add scheme Staggered1DAcoustics
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
20 % Pick out left or right boundary
8f0ee6ba6313 Add scheme Staggered1DAcoustics
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
21 e_l
8f0ee6ba6313 Add scheme Staggered1DAcoustics
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
22 e_r
8f0ee6ba6313 Add scheme Staggered1DAcoustics
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
23
8f0ee6ba6313 Add scheme Staggered1DAcoustics
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
24 % Initial data
8f0ee6ba6313 Add scheme Staggered1DAcoustics
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
25 v0
8f0ee6ba6313 Add scheme Staggered1DAcoustics
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
26
8f0ee6ba6313 Add scheme Staggered1DAcoustics
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
27 % Pick out primal or dual component
8f0ee6ba6313 Add scheme Staggered1DAcoustics
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
28 e_primal
8f0ee6ba6313 Add scheme Staggered1DAcoustics
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
29 e_dual
8f0ee6ba6313 Add scheme Staggered1DAcoustics
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
30
8f0ee6ba6313 Add scheme Staggered1DAcoustics
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
31 % System matrices
8f0ee6ba6313 Add scheme Staggered1DAcoustics
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
32 A
8f0ee6ba6313 Add scheme Staggered1DAcoustics
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
33 B
8f0ee6ba6313 Add scheme Staggered1DAcoustics
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
34 end
8f0ee6ba6313 Add scheme Staggered1DAcoustics
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
35
8f0ee6ba6313 Add scheme Staggered1DAcoustics
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
36
8f0ee6ba6313 Add scheme Staggered1DAcoustics
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
37 methods
8f0ee6ba6313 Add scheme Staggered1DAcoustics
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
38 % Scheme for A*u_t + B u_x = 0,
8f0ee6ba6313 Add scheme Staggered1DAcoustics
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
39 % u = [p, v];
8f0ee6ba6313 Add scheme Staggered1DAcoustics
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
40 % A: Diagonal and A > 0,
8f0ee6ba6313 Add scheme Staggered1DAcoustics
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
41 % A = [a1, 0;
8f0ee6ba6313 Add scheme Staggered1DAcoustics
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
42 % 0, a2]
8f0ee6ba6313 Add scheme Staggered1DAcoustics
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
43 %
8f0ee6ba6313 Add scheme Staggered1DAcoustics
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
44 % B = B^T and with diagonal entries = 0.
8f0ee6ba6313 Add scheme Staggered1DAcoustics
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
45 % B = [0 b
8f0ee6ba6313 Add scheme Staggered1DAcoustics
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
46 % b 0]
8f0ee6ba6313 Add scheme Staggered1DAcoustics
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
47 % Here we store p on the primal grid and v on the dual
8f0ee6ba6313 Add scheme Staggered1DAcoustics
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
48 function obj = Staggered1DAcoustics(g, order, A, B)
8f0ee6ba6313 Add scheme Staggered1DAcoustics
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
49 default_arg('A',[1, 0; 0, 1]);
8f0ee6ba6313 Add scheme Staggered1DAcoustics
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
50 default_arg('B',[0, 1; 1, 0]);
8f0ee6ba6313 Add scheme Staggered1DAcoustics
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
51
8f0ee6ba6313 Add scheme Staggered1DAcoustics
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
52 obj.order = order;
8f0ee6ba6313 Add scheme Staggered1DAcoustics
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
53 obj.A = A;
8f0ee6ba6313 Add scheme Staggered1DAcoustics
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
54 obj.B = B;
8f0ee6ba6313 Add scheme Staggered1DAcoustics
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
55
8f0ee6ba6313 Add scheme Staggered1DAcoustics
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
56 % Grids
8f0ee6ba6313 Add scheme Staggered1DAcoustics
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
57 obj.m = g.size();
8f0ee6ba6313 Add scheme Staggered1DAcoustics
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
58 xl = g.getBoundary('l');
8f0ee6ba6313 Add scheme Staggered1DAcoustics
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
59 xr = g.getBoundary('r');
8f0ee6ba6313 Add scheme Staggered1DAcoustics
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
60 xlim = {xl, xr};
8f0ee6ba6313 Add scheme Staggered1DAcoustics
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
61
8f0ee6ba6313 Add scheme Staggered1DAcoustics
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
62 obj.grid = g;
8f0ee6ba6313 Add scheme Staggered1DAcoustics
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
63 obj.grid_primal = g.grids{1};
8f0ee6ba6313 Add scheme Staggered1DAcoustics
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
64 obj.grid_dual = g.grids{2};
8f0ee6ba6313 Add scheme Staggered1DAcoustics
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
65
8f0ee6ba6313 Add scheme Staggered1DAcoustics
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
66 % Get operators
644
dc2918fb104d Switch to staggered-upwind staggered operators in Staggered1DAcoustics scheme.
Martin Almquist <malmquist@stanford.edu>
parents: 642
diff changeset
67 ops = sbp.D1StaggeredUpwind(obj.m, xlim, order);
642
8f0ee6ba6313 Add scheme Staggered1DAcoustics
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
68 obj.h = ops.h;
8f0ee6ba6313 Add scheme Staggered1DAcoustics
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
69
8f0ee6ba6313 Add scheme Staggered1DAcoustics
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
70 % Build combined operators
8f0ee6ba6313 Add scheme Staggered1DAcoustics
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
71 H_primal = ops.H_primal;
8f0ee6ba6313 Add scheme Staggered1DAcoustics
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
72 H_dual = ops.H_dual;
8f0ee6ba6313 Add scheme Staggered1DAcoustics
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
73 obj.H = blockmatrix.toMatrix( {H_primal, []; [], H_dual } );
8f0ee6ba6313 Add scheme Staggered1DAcoustics
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
74 obj.Hi = inv(obj.H);
8f0ee6ba6313 Add scheme Staggered1DAcoustics
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
75
8f0ee6ba6313 Add scheme Staggered1DAcoustics
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
76 D1_primal = ops.D1_primal;
8f0ee6ba6313 Add scheme Staggered1DAcoustics
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
77 D1_dual = ops.D1_dual;
8f0ee6ba6313 Add scheme Staggered1DAcoustics
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
78 D = {[], -1/A(1,1)*B(1,2)*D1_primal;...
8f0ee6ba6313 Add scheme Staggered1DAcoustics
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
79 -1/A(2,2)*B(2,1)*D1_dual, []};
8f0ee6ba6313 Add scheme Staggered1DAcoustics
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
80 obj.D = blockmatrix.toMatrix(D);
8f0ee6ba6313 Add scheme Staggered1DAcoustics
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
81 obj.D1_primal = D1_primal;
8f0ee6ba6313 Add scheme Staggered1DAcoustics
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
82 obj.D1_dual = D1_dual;
8f0ee6ba6313 Add scheme Staggered1DAcoustics
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
83
8f0ee6ba6313 Add scheme Staggered1DAcoustics
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
84 % Combined boundary operators
8f0ee6ba6313 Add scheme Staggered1DAcoustics
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
85 e_primal_l = ops.e_primal_l;
8f0ee6ba6313 Add scheme Staggered1DAcoustics
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
86 e_primal_r = ops.e_primal_r;
8f0ee6ba6313 Add scheme Staggered1DAcoustics
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
87 e_dual_l = ops.e_dual_l;
8f0ee6ba6313 Add scheme Staggered1DAcoustics
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
88 e_dual_r = ops.e_dual_r;
8f0ee6ba6313 Add scheme Staggered1DAcoustics
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
89 e_l = {e_primal_l, [];...
8f0ee6ba6313 Add scheme Staggered1DAcoustics
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
90 [] , e_dual_l};
8f0ee6ba6313 Add scheme Staggered1DAcoustics
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
91 e_r = {e_primal_r, [];...
8f0ee6ba6313 Add scheme Staggered1DAcoustics
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
92 [] , e_dual_r};
8f0ee6ba6313 Add scheme Staggered1DAcoustics
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
93 obj.e_l = blockmatrix.toMatrix(e_l);
8f0ee6ba6313 Add scheme Staggered1DAcoustics
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
94 obj.e_r = blockmatrix.toMatrix(e_r);
8f0ee6ba6313 Add scheme Staggered1DAcoustics
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
95
8f0ee6ba6313 Add scheme Staggered1DAcoustics
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
96 % Pick out first or second component of solution
8f0ee6ba6313 Add scheme Staggered1DAcoustics
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
97 N_primal = obj.grid_primal.N();
8f0ee6ba6313 Add scheme Staggered1DAcoustics
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
98 N_dual = obj.grid_dual.N();
8f0ee6ba6313 Add scheme Staggered1DAcoustics
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
99 obj.e_primal = [speye(N_primal, N_primal); sparse(N_dual, N_primal)];
8f0ee6ba6313 Add scheme Staggered1DAcoustics
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
100 obj.e_dual = [sparse(N_primal, N_dual); speye(N_dual, N_dual)];
8f0ee6ba6313 Add scheme Staggered1DAcoustics
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
101
8f0ee6ba6313 Add scheme Staggered1DAcoustics
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
102
8f0ee6ba6313 Add scheme Staggered1DAcoustics
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
103 end
8f0ee6ba6313 Add scheme Staggered1DAcoustics
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
104 % Closure functions return the operators applied to the own domain to close the boundary
8f0ee6ba6313 Add scheme Staggered1DAcoustics
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
105 % Penalty functions return the operators to force the solution. In the case of an interface it returns the operator applied to the other domain.
8f0ee6ba6313 Add scheme Staggered1DAcoustics
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
106 % boundary is a string specifying the boundary e.g. 'l','r' or 'e','w','n','s'.
8f0ee6ba6313 Add scheme Staggered1DAcoustics
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
107 % type is a string specifying the type of boundary condition if there are several.
8f0ee6ba6313 Add scheme Staggered1DAcoustics
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
108 % neighbour_scheme is an instance of Scheme that should be interfaced to.
8f0ee6ba6313 Add scheme Staggered1DAcoustics
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
109 % neighbour_boundary is a string specifying which boundary to interface to.
8f0ee6ba6313 Add scheme Staggered1DAcoustics
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
110 function [closure, penalty] = boundary_condition(obj, boundary, type)
8f0ee6ba6313 Add scheme Staggered1DAcoustics
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
111
8f0ee6ba6313 Add scheme Staggered1DAcoustics
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
112 default_arg('type','p');
8f0ee6ba6313 Add scheme Staggered1DAcoustics
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
113
8f0ee6ba6313 Add scheme Staggered1DAcoustics
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
114 % type = 'p' => boundary condition for p
8f0ee6ba6313 Add scheme Staggered1DAcoustics
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
115 % type = 'v' => boundary condition for v
671
993aac771efd Correct characteristic BC in staggered acoustics 1D constant coeff.
Martin Almquist <malmquist@stanford.edu>
parents: 650
diff changeset
116 % type = 'characteristic' => incoming characteristc = data
642
8f0ee6ba6313 Add scheme Staggered1DAcoustics
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
117 % No other types implemented yet
8f0ee6ba6313 Add scheme Staggered1DAcoustics
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
118
8f0ee6ba6313 Add scheme Staggered1DAcoustics
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
119 % BC on the form Lu - g = 0;
8f0ee6ba6313 Add scheme Staggered1DAcoustics
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
120
671
993aac771efd Correct characteristic BC in staggered acoustics 1D constant coeff.
Martin Almquist <malmquist@stanford.edu>
parents: 650
diff changeset
121 A = obj.A;
993aac771efd Correct characteristic BC in staggered acoustics 1D constant coeff.
Martin Almquist <malmquist@stanford.edu>
parents: 650
diff changeset
122 B = obj.B;
993aac771efd Correct characteristic BC in staggered acoustics 1D constant coeff.
Martin Almquist <malmquist@stanford.edu>
parents: 650
diff changeset
123 C = inv(A)*B;
993aac771efd Correct characteristic BC in staggered acoustics 1D constant coeff.
Martin Almquist <malmquist@stanford.edu>
parents: 650
diff changeset
124
642
8f0ee6ba6313 Add scheme Staggered1DAcoustics
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
125 % Diagonalize B
8f0ee6ba6313 Add scheme Staggered1DAcoustics
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
126 [T, Lambda] = eig(B);
8f0ee6ba6313 Add scheme Staggered1DAcoustics
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
127 lambda = diag(Lambda);
8f0ee6ba6313 Add scheme Staggered1DAcoustics
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
128
671
993aac771efd Correct characteristic BC in staggered acoustics 1D constant coeff.
Martin Almquist <malmquist@stanford.edu>
parents: 650
diff changeset
129 % Identify in- and outgoing variables
642
8f0ee6ba6313 Add scheme Staggered1DAcoustics
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
130 Iplus = lambda > 0;
8f0ee6ba6313 Add scheme Staggered1DAcoustics
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
131 Iminus = lambda < 0;
8f0ee6ba6313 Add scheme Staggered1DAcoustics
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
132
8f0ee6ba6313 Add scheme Staggered1DAcoustics
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
133 switch boundary
8f0ee6ba6313 Add scheme Staggered1DAcoustics
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
134 case 'l'
8f0ee6ba6313 Add scheme Staggered1DAcoustics
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
135 Iin = Iplus;
8f0ee6ba6313 Add scheme Staggered1DAcoustics
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
136 case 'r'
8f0ee6ba6313 Add scheme Staggered1DAcoustics
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
137 Iin = Iminus;
8f0ee6ba6313 Add scheme Staggered1DAcoustics
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
138 end
8f0ee6ba6313 Add scheme Staggered1DAcoustics
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
139 Tin = T(:,Iin);
8f0ee6ba6313 Add scheme Staggered1DAcoustics
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
140
650
8e55298657b9 Add characteristic BC
Martin Almquist <malmquist@stanford.edu>
parents: 644
diff changeset
141 switch type
8e55298657b9 Add characteristic BC
Martin Almquist <malmquist@stanford.edu>
parents: 644
diff changeset
142 case 'p'
8e55298657b9 Add characteristic BC
Martin Almquist <malmquist@stanford.edu>
parents: 644
diff changeset
143 L = [1, 0];
8e55298657b9 Add characteristic BC
Martin Almquist <malmquist@stanford.edu>
parents: 644
diff changeset
144 case 'v'
8e55298657b9 Add characteristic BC
Martin Almquist <malmquist@stanford.edu>
parents: 644
diff changeset
145 L = [0, 1];
8e55298657b9 Add characteristic BC
Martin Almquist <malmquist@stanford.edu>
parents: 644
diff changeset
146 case 'characteristic'
671
993aac771efd Correct characteristic BC in staggered acoustics 1D constant coeff.
Martin Almquist <malmquist@stanford.edu>
parents: 650
diff changeset
147 % Diagonalize C
993aac771efd Correct characteristic BC in staggered acoustics 1D constant coeff.
Martin Almquist <malmquist@stanford.edu>
parents: 650
diff changeset
148 [T_C, Lambda_C] = eig(C);
993aac771efd Correct characteristic BC in staggered acoustics 1D constant coeff.
Martin Almquist <malmquist@stanford.edu>
parents: 650
diff changeset
149 lambda_C = diag(Lambda_C);
993aac771efd Correct characteristic BC in staggered acoustics 1D constant coeff.
Martin Almquist <malmquist@stanford.edu>
parents: 650
diff changeset
150 % Identify in- and outgoing characteristic variables
993aac771efd Correct characteristic BC in staggered acoustics 1D constant coeff.
Martin Almquist <malmquist@stanford.edu>
parents: 650
diff changeset
151 Iplus = lambda_C > 0;
993aac771efd Correct characteristic BC in staggered acoustics 1D constant coeff.
Martin Almquist <malmquist@stanford.edu>
parents: 650
diff changeset
152 Iminus = lambda_C < 0;
993aac771efd Correct characteristic BC in staggered acoustics 1D constant coeff.
Martin Almquist <malmquist@stanford.edu>
parents: 650
diff changeset
153
993aac771efd Correct characteristic BC in staggered acoustics 1D constant coeff.
Martin Almquist <malmquist@stanford.edu>
parents: 650
diff changeset
154 switch boundary
993aac771efd Correct characteristic BC in staggered acoustics 1D constant coeff.
Martin Almquist <malmquist@stanford.edu>
parents: 650
diff changeset
155 case 'l'
993aac771efd Correct characteristic BC in staggered acoustics 1D constant coeff.
Martin Almquist <malmquist@stanford.edu>
parents: 650
diff changeset
156 Iin_C = Iplus;
993aac771efd Correct characteristic BC in staggered acoustics 1D constant coeff.
Martin Almquist <malmquist@stanford.edu>
parents: 650
diff changeset
157 case 'r'
993aac771efd Correct characteristic BC in staggered acoustics 1D constant coeff.
Martin Almquist <malmquist@stanford.edu>
parents: 650
diff changeset
158 Iin_C = Iminus;
993aac771efd Correct characteristic BC in staggered acoustics 1D constant coeff.
Martin Almquist <malmquist@stanford.edu>
parents: 650
diff changeset
159 end
993aac771efd Correct characteristic BC in staggered acoustics 1D constant coeff.
Martin Almquist <malmquist@stanford.edu>
parents: 650
diff changeset
160 T_C_inv = inv(T_C);
993aac771efd Correct characteristic BC in staggered acoustics 1D constant coeff.
Martin Almquist <malmquist@stanford.edu>
parents: 650
diff changeset
161 L = T_C_inv(Iin_C,:);
650
8e55298657b9 Add characteristic BC
Martin Almquist <malmquist@stanford.edu>
parents: 644
diff changeset
162 otherwise
8e55298657b9 Add characteristic BC
Martin Almquist <malmquist@stanford.edu>
parents: 644
diff changeset
163 error('Boundary condition not implemented.');
8e55298657b9 Add characteristic BC
Martin Almquist <malmquist@stanford.edu>
parents: 644
diff changeset
164 end
8e55298657b9 Add characteristic BC
Martin Almquist <malmquist@stanford.edu>
parents: 644
diff changeset
165
642
8f0ee6ba6313 Add scheme Staggered1DAcoustics
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
166 % Penalty parameters
8f0ee6ba6313 Add scheme Staggered1DAcoustics
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
167 sigma = [0; 0];
8f0ee6ba6313 Add scheme Staggered1DAcoustics
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
168 sigma(Iin) = lambda(Iin);
671
993aac771efd Correct characteristic BC in staggered acoustics 1D constant coeff.
Martin Almquist <malmquist@stanford.edu>
parents: 650
diff changeset
169
993aac771efd Correct characteristic BC in staggered acoustics 1D constant coeff.
Martin Almquist <malmquist@stanford.edu>
parents: 650
diff changeset
170 % Sparsify
993aac771efd Correct characteristic BC in staggered acoustics 1D constant coeff.
Martin Almquist <malmquist@stanford.edu>
parents: 650
diff changeset
171 A = sparse(A);
993aac771efd Correct characteristic BC in staggered acoustics 1D constant coeff.
Martin Almquist <malmquist@stanford.edu>
parents: 650
diff changeset
172 T = sparse(T);
993aac771efd Correct characteristic BC in staggered acoustics 1D constant coeff.
Martin Almquist <malmquist@stanford.edu>
parents: 650
diff changeset
173 sigma = sparse(sigma);
993aac771efd Correct characteristic BC in staggered acoustics 1D constant coeff.
Martin Almquist <malmquist@stanford.edu>
parents: 650
diff changeset
174 L = sparse(L);
993aac771efd Correct characteristic BC in staggered acoustics 1D constant coeff.
Martin Almquist <malmquist@stanford.edu>
parents: 650
diff changeset
175 Tin = sparse(Tin);
993aac771efd Correct characteristic BC in staggered acoustics 1D constant coeff.
Martin Almquist <malmquist@stanford.edu>
parents: 650
diff changeset
176
642
8f0ee6ba6313 Add scheme Staggered1DAcoustics
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
177 switch boundary
8f0ee6ba6313 Add scheme Staggered1DAcoustics
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
178 case 'l'
8f0ee6ba6313 Add scheme Staggered1DAcoustics
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
179 tau = -1*obj.e_l * inv(A) * T * sigma * inv(L*Tin);
8f0ee6ba6313 Add scheme Staggered1DAcoustics
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
180 closure = obj.Hi*tau*L*obj.e_l';
8f0ee6ba6313 Add scheme Staggered1DAcoustics
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
181
8f0ee6ba6313 Add scheme Staggered1DAcoustics
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
182 case 'r'
8f0ee6ba6313 Add scheme Staggered1DAcoustics
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
183 tau = 1*obj.e_r * inv(A) * T * sigma * inv(L*Tin);
8f0ee6ba6313 Add scheme Staggered1DAcoustics
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
184 closure = obj.Hi*tau*L*obj.e_r';
8f0ee6ba6313 Add scheme Staggered1DAcoustics
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
185
8f0ee6ba6313 Add scheme Staggered1DAcoustics
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
186 end
8f0ee6ba6313 Add scheme Staggered1DAcoustics
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
187
8f0ee6ba6313 Add scheme Staggered1DAcoustics
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
188 penalty = -obj.Hi*tau;
8f0ee6ba6313 Add scheme Staggered1DAcoustics
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
189
8f0ee6ba6313 Add scheme Staggered1DAcoustics
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
190 end
8f0ee6ba6313 Add scheme Staggered1DAcoustics
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
191
8f0ee6ba6313 Add scheme Staggered1DAcoustics
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
192 function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary)
8f0ee6ba6313 Add scheme Staggered1DAcoustics
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
193
8f0ee6ba6313 Add scheme Staggered1DAcoustics
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
194 error('Staggered1DAcoustics, interface not implemented');
8f0ee6ba6313 Add scheme Staggered1DAcoustics
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
195
8f0ee6ba6313 Add scheme Staggered1DAcoustics
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
196 switch boundary
8f0ee6ba6313 Add scheme Staggered1DAcoustics
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
197 % Upwind coupling
8f0ee6ba6313 Add scheme Staggered1DAcoustics
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
198 case {'l','left'}
8f0ee6ba6313 Add scheme Staggered1DAcoustics
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
199 tau = -1*obj.e_l;
8f0ee6ba6313 Add scheme Staggered1DAcoustics
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
200 closure = obj.Hi*tau*obj.e_l';
8f0ee6ba6313 Add scheme Staggered1DAcoustics
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
201 penalty = -obj.Hi*tau*neighbour_scheme.e_r';
8f0ee6ba6313 Add scheme Staggered1DAcoustics
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
202 case {'r','right'}
8f0ee6ba6313 Add scheme Staggered1DAcoustics
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
203 tau = 0*obj.e_r;
8f0ee6ba6313 Add scheme Staggered1DAcoustics
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
204 closure = obj.Hi*tau*obj.e_r';
8f0ee6ba6313 Add scheme Staggered1DAcoustics
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
205 penalty = -obj.Hi*tau*neighbour_scheme.e_l';
8f0ee6ba6313 Add scheme Staggered1DAcoustics
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
206 end
8f0ee6ba6313 Add scheme Staggered1DAcoustics
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
207
8f0ee6ba6313 Add scheme Staggered1DAcoustics
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
208 end
8f0ee6ba6313 Add scheme Staggered1DAcoustics
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
209
8f0ee6ba6313 Add scheme Staggered1DAcoustics
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
210 function N = size(obj)
8f0ee6ba6313 Add scheme Staggered1DAcoustics
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
211 N = obj.m;
8f0ee6ba6313 Add scheme Staggered1DAcoustics
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
212 end
8f0ee6ba6313 Add scheme Staggered1DAcoustics
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
213
8f0ee6ba6313 Add scheme Staggered1DAcoustics
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
214 end
8f0ee6ba6313 Add scheme Staggered1DAcoustics
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
215
8f0ee6ba6313 Add scheme Staggered1DAcoustics
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
216 methods(Static)
8f0ee6ba6313 Add scheme Staggered1DAcoustics
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
217 % Calculates the matrices needed for the inteface coupling between boundary bound_u of scheme schm_u
8f0ee6ba6313 Add scheme Staggered1DAcoustics
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
218 % and bound_v of scheme schm_v.
8f0ee6ba6313 Add scheme Staggered1DAcoustics
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
219 % [uu, uv, vv, vu] = inteface_coupling(A,'r',B,'l')
8f0ee6ba6313 Add scheme Staggered1DAcoustics
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
220 function [uu, uv, vv, vu] = interface_coupling(schm_u,bound_u,schm_v,bound_v)
8f0ee6ba6313 Add scheme Staggered1DAcoustics
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
221 [uu,uv] = schm_u.interface(bound_u,schm_v,bound_v);
8f0ee6ba6313 Add scheme Staggered1DAcoustics
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
222 [vv,vu] = schm_v.interface(bound_v,schm_u,bound_u);
8f0ee6ba6313 Add scheme Staggered1DAcoustics
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
223 end
8f0ee6ba6313 Add scheme Staggered1DAcoustics
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
224 end
8f0ee6ba6313 Add scheme Staggered1DAcoustics
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
225 end