Mercurial > repos > public > sbplib

annotate +scheme/Staggered1DAcoustics.m @ 682:ec0ac76006e2 feature/d1_staggered

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Modify BC in Staggered 1D acoustics to ensure that only physically meaningful changes of variables are used. See hypsyst varcoeff sec 3.
author Martin Almquist <malmquist@stanford.edu>
date Wed, 07 Feb 2018 20:07:34 -0800
parents 993aac771efd
children
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
682
ec0ac76006e2 Modify BC in Staggered 1D acoustics to ensure that only physically meaningful changes of variables are used. See hypsyst varcoeff sec 3.
Martin Almquist <malmquist@stanford.edu>
parents: 671
diff changeset
121 % Need a transformation T such that w = T^{−1}*u is a
ec0ac76006e2 Modify BC in Staggered 1D acoustics to ensure that only physically meaningful changes of variables are used. See hypsyst varcoeff sec 3.
Martin Almquist <malmquist@stanford.edu>
parents: 671
diff changeset
122 % meaningful change of variables and T^T*B*T is block-diagonal,
ec0ac76006e2 Modify BC in Staggered 1D acoustics to ensure that only physically meaningful changes of variables are used. See hypsyst varcoeff sec 3.
Martin Almquist <malmquist@stanford.edu>
parents: 671
diff changeset
123 % For linear acoustics, T = T_C meets these criteria.
ec0ac76006e2 Modify BC in Staggered 1D acoustics to ensure that only physically meaningful changes of variables are used. See hypsyst varcoeff sec 3.
Martin Almquist <malmquist@stanford.edu>
parents: 671
diff changeset
124
671
993aac771efd Correct characteristic BC in staggered acoustics 1D constant coeff.
Martin Almquist <malmquist@stanford.edu>
parents: 650
diff changeset
125 A = obj.A;
993aac771efd Correct characteristic BC in staggered acoustics 1D constant coeff.
Martin Almquist <malmquist@stanford.edu>
parents: 650
diff changeset
126 B = obj.B;
993aac771efd Correct characteristic BC in staggered acoustics 1D constant coeff.
Martin Almquist <malmquist@stanford.edu>
parents: 650
diff changeset
127 C = inv(A)*B;
993aac771efd Correct characteristic BC in staggered acoustics 1D constant coeff.
Martin Almquist <malmquist@stanford.edu>
parents: 650
diff changeset
128
682
ec0ac76006e2 Modify BC in Staggered 1D acoustics to ensure that only physically meaningful changes of variables are used. See hypsyst varcoeff sec 3.
Martin Almquist <malmquist@stanford.edu>
parents: 671
diff changeset
129 % Diagonalize C and use T_C to diagonalize B
ec0ac76006e2 Modify BC in Staggered 1D acoustics to ensure that only physically meaningful changes of variables are used. See hypsyst varcoeff sec 3.
Martin Almquist <malmquist@stanford.edu>
parents: 671
diff changeset
130 [T, ~] = eig(C);
ec0ac76006e2 Modify BC in Staggered 1D acoustics to ensure that only physically meaningful changes of variables are used. See hypsyst varcoeff sec 3.
Martin Almquist <malmquist@stanford.edu>
parents: 671
diff changeset
131 Lambda = T'*B*T;
642
8f0ee6ba6313 Add scheme Staggered1DAcoustics
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
132 lambda = diag(Lambda);
8f0ee6ba6313 Add scheme Staggered1DAcoustics
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
133
682
ec0ac76006e2 Modify BC in Staggered 1D acoustics to ensure that only physically meaningful changes of variables are used. See hypsyst varcoeff sec 3.
Martin Almquist <malmquist@stanford.edu>
parents: 671
diff changeset
134
671
993aac771efd Correct characteristic BC in staggered acoustics 1D constant coeff.
Martin Almquist <malmquist@stanford.edu>
parents: 650
diff changeset
135 % Identify in- and outgoing variables
642
8f0ee6ba6313 Add scheme Staggered1DAcoustics
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
136 Iplus = lambda > 0;
8f0ee6ba6313 Add scheme Staggered1DAcoustics
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
137 Iminus = lambda < 0;
8f0ee6ba6313 Add scheme Staggered1DAcoustics
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
138
8f0ee6ba6313 Add scheme Staggered1DAcoustics
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
139 switch boundary
8f0ee6ba6313 Add scheme Staggered1DAcoustics
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
140 case 'l'
8f0ee6ba6313 Add scheme Staggered1DAcoustics
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
141 Iin = Iplus;
8f0ee6ba6313 Add scheme Staggered1DAcoustics
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
142 case 'r'
8f0ee6ba6313 Add scheme Staggered1DAcoustics
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
143 Iin = Iminus;
8f0ee6ba6313 Add scheme Staggered1DAcoustics
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
144 end
8f0ee6ba6313 Add scheme Staggered1DAcoustics
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
145 Tin = T(:,Iin);
8f0ee6ba6313 Add scheme Staggered1DAcoustics
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
146
650
8e55298657b9 Add characteristic BC
Martin Almquist <malmquist@stanford.edu>
parents: 644
diff changeset
147 switch type
8e55298657b9 Add characteristic BC
Martin Almquist <malmquist@stanford.edu>
parents: 644
diff changeset
148 case 'p'
8e55298657b9 Add characteristic BC
Martin Almquist <malmquist@stanford.edu>
parents: 644
diff changeset
149 L = [1, 0];
8e55298657b9 Add characteristic BC
Martin Almquist <malmquist@stanford.edu>
parents: 644
diff changeset
150 case 'v'
8e55298657b9 Add characteristic BC
Martin Almquist <malmquist@stanford.edu>
parents: 644
diff changeset
151 L = [0, 1];
8e55298657b9 Add characteristic BC
Martin Almquist <malmquist@stanford.edu>
parents: 644
diff changeset
152 case 'characteristic'
671
993aac771efd Correct characteristic BC in staggered acoustics 1D constant coeff.
Martin Almquist <malmquist@stanford.edu>
parents: 650
diff changeset
153 % Diagonalize C
993aac771efd Correct characteristic BC in staggered acoustics 1D constant coeff.
Martin Almquist <malmquist@stanford.edu>
parents: 650
diff changeset
154 [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
155 lambda_C = diag(Lambda_C);
993aac771efd Correct characteristic BC in staggered acoustics 1D constant coeff.
Martin Almquist <malmquist@stanford.edu>
parents: 650
diff changeset
156 % 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
157 Iplus = lambda_C > 0;
993aac771efd Correct characteristic BC in staggered acoustics 1D constant coeff.
Martin Almquist <malmquist@stanford.edu>
parents: 650
diff changeset
158 Iminus = lambda_C < 0;
993aac771efd Correct characteristic BC in staggered acoustics 1D constant coeff.
Martin Almquist <malmquist@stanford.edu>
parents: 650
diff changeset
159
993aac771efd Correct characteristic BC in staggered acoustics 1D constant coeff.
Martin Almquist <malmquist@stanford.edu>
parents: 650
diff changeset
160 switch boundary
993aac771efd Correct characteristic BC in staggered acoustics 1D constant coeff.
Martin Almquist <malmquist@stanford.edu>
parents: 650
diff changeset
161 case 'l'
993aac771efd Correct characteristic BC in staggered acoustics 1D constant coeff.
Martin Almquist <malmquist@stanford.edu>
parents: 650
diff changeset
162 Iin_C = Iplus;
993aac771efd Correct characteristic BC in staggered acoustics 1D constant coeff.
Martin Almquist <malmquist@stanford.edu>
parents: 650
diff changeset
163 case 'r'
993aac771efd Correct characteristic BC in staggered acoustics 1D constant coeff.
Martin Almquist <malmquist@stanford.edu>
parents: 650
diff changeset
164 Iin_C = Iminus;
993aac771efd Correct characteristic BC in staggered acoustics 1D constant coeff.
Martin Almquist <malmquist@stanford.edu>
parents: 650
diff changeset
165 end
993aac771efd Correct characteristic BC in staggered acoustics 1D constant coeff.
Martin Almquist <malmquist@stanford.edu>
parents: 650
diff changeset
166 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
167 L = T_C_inv(Iin_C,:);
650
8e55298657b9 Add characteristic BC
Martin Almquist <malmquist@stanford.edu>
parents: 644
diff changeset
168 otherwise
8e55298657b9 Add characteristic BC
Martin Almquist <malmquist@stanford.edu>
parents: 644
diff changeset
169 error('Boundary condition not implemented.');
8e55298657b9 Add characteristic BC
Martin Almquist <malmquist@stanford.edu>
parents: 644
diff changeset
170 end
8e55298657b9 Add characteristic BC
Martin Almquist <malmquist@stanford.edu>
parents: 644
diff changeset
171
642
8f0ee6ba6313 Add scheme Staggered1DAcoustics
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
172 % Penalty parameters
8f0ee6ba6313 Add scheme Staggered1DAcoustics
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
173 sigma = [0; 0];
8f0ee6ba6313 Add scheme Staggered1DAcoustics
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
174 sigma(Iin) = lambda(Iin);
671
993aac771efd Correct characteristic BC in staggered acoustics 1D constant coeff.
Martin Almquist <malmquist@stanford.edu>
parents: 650
diff changeset
175
993aac771efd Correct characteristic BC in staggered acoustics 1D constant coeff.
Martin Almquist <malmquist@stanford.edu>
parents: 650
diff changeset
176 % Sparsify
993aac771efd Correct characteristic BC in staggered acoustics 1D constant coeff.
Martin Almquist <malmquist@stanford.edu>
parents: 650
diff changeset
177 A = sparse(A);
993aac771efd Correct characteristic BC in staggered acoustics 1D constant coeff.
Martin Almquist <malmquist@stanford.edu>
parents: 650
diff changeset
178 T = sparse(T);
993aac771efd Correct characteristic BC in staggered acoustics 1D constant coeff.
Martin Almquist <malmquist@stanford.edu>
parents: 650
diff changeset
179 sigma = sparse(sigma);
993aac771efd Correct characteristic BC in staggered acoustics 1D constant coeff.
Martin Almquist <malmquist@stanford.edu>
parents: 650
diff changeset
180 L = sparse(L);
993aac771efd Correct characteristic BC in staggered acoustics 1D constant coeff.
Martin Almquist <malmquist@stanford.edu>
parents: 650
diff changeset
181 Tin = sparse(Tin);
993aac771efd Correct characteristic BC in staggered acoustics 1D constant coeff.
Martin Almquist <malmquist@stanford.edu>
parents: 650
diff changeset
182
642
8f0ee6ba6313 Add scheme Staggered1DAcoustics
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
183 switch boundary
8f0ee6ba6313 Add scheme Staggered1DAcoustics
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
184 case 'l'
682
ec0ac76006e2 Modify BC in Staggered 1D acoustics to ensure that only physically meaningful changes of variables are used. See hypsyst varcoeff sec 3.
Martin Almquist <malmquist@stanford.edu>
parents: 671
diff changeset
185 tau = -1*obj.e_l * inv(A) * inv(T)' * sigma * inv(L*Tin);
642
8f0ee6ba6313 Add scheme Staggered1DAcoustics
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
186 closure = obj.Hi*tau*L*obj.e_l';
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 case 'r'
682
ec0ac76006e2 Modify BC in Staggered 1D acoustics to ensure that only physically meaningful changes of variables are used. See hypsyst varcoeff sec 3.
Martin Almquist <malmquist@stanford.edu>
parents: 671
diff changeset
189 tau = 1*obj.e_r * inv(A) * inv(T)' * sigma * inv(L*Tin);
642
8f0ee6ba6313 Add scheme Staggered1DAcoustics
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
190 closure = obj.Hi*tau*L*obj.e_r';
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 end
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 penalty = -obj.Hi*tau;
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 end
8f0ee6ba6313 Add scheme Staggered1DAcoustics
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
197
8f0ee6ba6313 Add scheme Staggered1DAcoustics
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
198 function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary)
8f0ee6ba6313 Add scheme Staggered1DAcoustics
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
199
8f0ee6ba6313 Add scheme Staggered1DAcoustics
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
200 error('Staggered1DAcoustics, interface not implemented');
8f0ee6ba6313 Add scheme Staggered1DAcoustics
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
201
8f0ee6ba6313 Add scheme Staggered1DAcoustics
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
202 switch boundary
8f0ee6ba6313 Add scheme Staggered1DAcoustics
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
203 % Upwind coupling
8f0ee6ba6313 Add scheme Staggered1DAcoustics
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
204 case {'l','left'}
8f0ee6ba6313 Add scheme Staggered1DAcoustics
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
205 tau = -1*obj.e_l;
8f0ee6ba6313 Add scheme Staggered1DAcoustics
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
206 closure = obj.Hi*tau*obj.e_l';
8f0ee6ba6313 Add scheme Staggered1DAcoustics
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
207 penalty = -obj.Hi*tau*neighbour_scheme.e_r';
8f0ee6ba6313 Add scheme Staggered1DAcoustics
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
208 case {'r','right'}
8f0ee6ba6313 Add scheme Staggered1DAcoustics
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
209 tau = 0*obj.e_r;
8f0ee6ba6313 Add scheme Staggered1DAcoustics
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
210 closure = obj.Hi*tau*obj.e_r';
8f0ee6ba6313 Add scheme Staggered1DAcoustics
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
211 penalty = -obj.Hi*tau*neighbour_scheme.e_l';
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 function N = size(obj)
8f0ee6ba6313 Add scheme Staggered1DAcoustics
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
217 N = obj.m;
8f0ee6ba6313 Add scheme Staggered1DAcoustics
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
218 end
8f0ee6ba6313 Add scheme Staggered1DAcoustics
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
219
8f0ee6ba6313 Add scheme Staggered1DAcoustics
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
220 end
8f0ee6ba6313 Add scheme Staggered1DAcoustics
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
221
8f0ee6ba6313 Add scheme Staggered1DAcoustics
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
222 methods(Static)
8f0ee6ba6313 Add scheme Staggered1DAcoustics
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
223 % 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
224 % and bound_v of scheme schm_v.
8f0ee6ba6313 Add scheme Staggered1DAcoustics
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
225 % [uu, uv, vv, vu] = inteface_coupling(A,'r',B,'l')
8f0ee6ba6313 Add scheme Staggered1DAcoustics
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
226 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
227 [uu,uv] = schm_u.interface(bound_u,schm_v,bound_v);
8f0ee6ba6313 Add scheme Staggered1DAcoustics
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
228 [vv,vu] = schm_v.interface(bound_v,schm_u,bound_u);
8f0ee6ba6313 Add scheme Staggered1DAcoustics
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
229 end
8f0ee6ba6313 Add scheme Staggered1DAcoustics
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
230 end
8f0ee6ba6313 Add scheme Staggered1DAcoustics
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
231 end