Mercurial > repos > public > sbplib

comparison +scheme/Utux2D.m @ 743:f4595f14d696 feature/utux2D

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Change schemes to work for special coefficients.
author Martin Almquist <malmquist@stanford.edu>
date Mon, 21 May 2018 23:11:34 -0700
parents 2d85f17a8aec
children 459eeb99130f
comparison
equal deleted inserted replaced
720:07f8311374c6 743:f4595f14d696
5 grid % Grid 5 grid % Grid
6 order % Order accuracy for the approximation 6 order % Order accuracy for the approximation
7 v0 % Initial data 7 v0 % Initial data
8 8
9 a % Wave speed a = [a1, a2]; 9 a % Wave speed a = [a1, a2];
10 % Can either be a constant vector or a cell array of function handles.
10 11
11 H % Discrete norm 12 H % Discrete norm
12 H_x, H_y % Norms in the x and y directions 13 H_x, H_y % Norms in the x and y directions
13 Hi, Hx, Hy, Hxi, Hyi % Kroneckered norms 14 Hi, Hx, Hy, Hxi, Hyi % Kroneckered norms
14 15
43 default_arg('interpolation_damping',{0,0}); 44 default_arg('interpolation_damping',{0,0});
44 default_arg('interpolation_type','AWW'); 45 default_arg('interpolation_type','AWW');
45 default_arg('coupling_type','upwind'); 46 default_arg('coupling_type','upwind');
46 default_arg('a',1/sqrt(2)*[1, 1]); 47 default_arg('a',1/sqrt(2)*[1, 1]);
47 default_arg('opSet',@sbp.D2Standard); 48 default_arg('opSet',@sbp.D2Standard);
49
48 assert(isa(g, 'grid.Cartesian')) 50 assert(isa(g, 'grid.Cartesian'))
51 if iscell(a)
52 a1 = grid.evalOn(g, a{1});
53 a2 = grid.evalOn(g, a{2});
54 a = {spdiag(a1), spdiag(a2)};
55 else
56 a = {a(1), a(2)};
57 end
49 58
50 m = g.size(); 59 m = g.size();
51 m_x = m(1); 60 m_x = m(1);
52 m_y = m(2); 61 m_y = m(2);
53 m_tot = g.N(); 62 m_tot = g.N();
94 obj.order = order; 103 obj.order = order;
95 obj.a = a; 104 obj.a = a;
96 obj.coupling_type = coupling_type; 105 obj.coupling_type = coupling_type;
97 obj.interpolation_type = interpolation_type; 106 obj.interpolation_type = interpolation_type;
98 obj.interpolation_damping = interpolation_damping; 107 obj.interpolation_damping = interpolation_damping;
99 obj.D = -(a(1)*obj.Dx + a(2)*obj.Dy); 108 obj.D = -(a{1}*obj.Dx + a{2}*obj.Dy);
100 109
101 end 110 end
102 % Closure functions return the opertors applied to the own domain to close the boundary 111 % Closure functions return the opertors applied to the own domain to close the boundary
103 % Penalty functions return the opertors to force the solution. In the case of an interface it returns the operator applied to the other doamin. 112 % Penalty functions return the opertors to force the solution. In the case of an interface it returns the operator applied to the other doamin.
104 % boundary is a string specifying the boundary e.g. 'l','r' or 'e','w','n','s'. 113 % boundary is a string specifying the boundary e.g. 'l','r' or 'e','w','n','s'.
110 default_arg('type','dirichlet'); 119 default_arg('type','dirichlet');
111 120
112 sigma = -1; % Scalar penalty parameter 121 sigma = -1; % Scalar penalty parameter
113 switch boundary 122 switch boundary
114 case {'w','W','west','West'} 123 case {'w','W','west','West'}
115 tau = sigma*obj.a(1)*obj.e_w*obj.H_y; 124 tau = sigma*obj.a{1}*obj.e_w*obj.H_y;
116 closure = obj.Hi*tau*obj.e_w'; 125 closure = obj.Hi*tau*obj.e_w';
117 126
118 case {'s','S','south','South'} 127 case {'s','S','south','South'}
119 tau = sigma*obj.a(2)*obj.e_s*obj.H_x; 128 tau = sigma*obj.a{2}*obj.e_s*obj.H_x;
120 closure = obj.Hi*tau*obj.e_s'; 129 closure = obj.Hi*tau*obj.e_s';
121 end 130 end
122 penalty = -obj.Hi*tau; 131 penalty = -obj.Hi*tau;
123 132
124 end 133 end
221 I_back_forth_ds = I_neighbour2local_ds*I_local2neighbour_ds; 230 I_back_forth_ds = I_neighbour2local_ds*I_local2neighbour_ds;
222 231
223 232
224 switch boundary 233 switch boundary
225 case {'w','W','west','West'} 234 case {'w','W','west','West'}
226 tau = sigma_ds*obj.a(1)*obj.e_w*obj.H_y; 235 tau = sigma_ds*obj.a{1}*obj.e_w*obj.H_y;
227 closure = obj.Hi*tau*obj.e_w'; 236 closure = obj.Hi*tau*obj.e_w';
228 penalty = -obj.Hi*tau*I_neighbour2local_ds*e_neighbour'; 237 penalty = -obj.Hi*tau*I_neighbour2local_ds*e_neighbour';
229 238
230 beta = int_damp_ds*obj.a(1)... 239 beta = int_damp_ds*obj.a{1}...
231 *obj.e_w*obj.H_y; 240 *obj.e_w*obj.H_y;
232 closure = closure + obj.Hi*beta*(I_back_forth_ds - I)*obj.e_w'; 241 closure = closure + obj.Hi*beta*(I_back_forth_ds - I)*obj.e_w';
233 case {'e','E','east','East'} 242 case {'e','E','east','East'}
234 tau = sigma_us*obj.a(1)*obj.e_e*obj.H_y; 243 tau = sigma_us*obj.a{1}*obj.e_e*obj.H_y;
235 closure = obj.Hi*tau*obj.e_e'; 244 closure = obj.Hi*tau*obj.e_e';
236 penalty = -obj.Hi*tau*I_neighbour2local_us*e_neighbour'; 245 penalty = -obj.Hi*tau*I_neighbour2local_us*e_neighbour';
237 246
238 beta = int_damp_us*obj.a(1)... 247 beta = int_damp_us*obj.a{1}...
239 *obj.e_e*obj.H_y; 248 *obj.e_e*obj.H_y;
240 closure = closure + obj.Hi*beta*(I_back_forth_us - I)*obj.e_e'; 249 closure = closure + obj.Hi*beta*(I_back_forth_us - I)*obj.e_e';
241 case {'s','S','south','South'} 250 case {'s','S','south','South'}
242 tau = sigma_ds*obj.a(2)*obj.e_s*obj.H_x; 251 tau = sigma_ds*obj.a{2}*obj.e_s*obj.H_x;
243 closure = obj.Hi*tau*obj.e_s'; 252 closure = obj.Hi*tau*obj.e_s';
244 penalty = -obj.Hi*tau*I_neighbour2local_ds*e_neighbour'; 253 penalty = -obj.Hi*tau*I_neighbour2local_ds*e_neighbour';
245 254
246 beta = int_damp_ds*obj.a(2)... 255 beta = int_damp_ds*obj.a{2}...
247 *obj.e_s*obj.H_x; 256 *obj.e_s*obj.H_x;
248 closure = closure + obj.Hi*beta*(I_back_forth_ds - I)*obj.e_s'; 257 closure = closure + obj.Hi*beta*(I_back_forth_ds - I)*obj.e_s';
249 case {'n','N','north','North'} 258 case {'n','N','north','North'}
250 tau = sigma_us*obj.a(2)*obj.e_n*obj.H_x; 259 tau = sigma_us*obj.a{2}*obj.e_n*obj.H_x;
251 closure = obj.Hi*tau*obj.e_n'; 260 closure = obj.Hi*tau*obj.e_n';
252 penalty = -obj.Hi*tau*I_neighbour2local_us*e_neighbour'; 261 penalty = -obj.Hi*tau*I_neighbour2local_us*e_neighbour';
253 262
254 beta = int_damp_us*obj.a(2)... 263 beta = int_damp_us*obj.a{2}...
255 *obj.e_n*obj.H_x; 264 *obj.e_n*obj.H_x;
256 closure = closure + obj.Hi*beta*(I_back_forth_us - I)*obj.e_n'; 265 closure = closure + obj.Hi*beta*(I_back_forth_us - I)*obj.e_n';
257 end 266 end
258 267
259 268