Mercurial > repos > public > sbplib

comparison +scheme/Wave2dCurve.m @ 0:48b6fb693025

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Initial commit.
author Jonatan Werpers <jonatan@werpers.com>
date Thu, 17 Sep 2015 10:12:50 +0200
parents
children 5f6b0b6a012b
comparison
equal deleted inserted replaced
-1:000000000000 0:48b6fb693025
1 classdef Wave2dCurve < scheme.Scheme
2 properties
3 m % Number of points in each direction, possibly a vector
4 h % Grid spacing
5 u,v % Grid
6 x,y % Values of x and y for each grid point
7 X,Y % Grid point locations as matrices
8 order % Order accuracy for the approximation
9
10 D % non-stabalized scheme operator
11 M % Derivative norm
12 c
13 J, Ji
14 a11, a12, a22
15
16 H % Discrete norm
17 Hi
18 H_u, H_v % Norms in the x and y directions
19 Hu,Hv % Kroneckerd norms. 1'*Hx*v corresponds to integration in the x dir.
20 Hi_u, Hi_v
21 Hiu, Hiv
22 e_w, e_e, e_s, e_n
23 du_w, dv_w
24 du_e, dv_e
25 du_s, dv_s
26 du_n, dv_n
27 gamm_u, gamm_v
28 end
29
30 methods
31 function obj = Wave2dCurve(m,ti,order,c,opSet)
32 default_arg('opSet',@sbp.Variable);
33
34 if length(m) == 1
35 m = [m m];
36 end
37
38 m_u = m(1);
39 m_v = m(2);
40 m_tot = m_u*m_v;
41
42 [u, h_u] = util.get_grid(0, 1, m_u);
43 [v, h_v] = util.get_grid(0, 1, m_v);
44
45
46 % Operators
47 ops_u = opSet(m_u,h_u,order);
48 ops_v = opSet(m_v,h_v,order);
49
50 I_u = speye(m_u);
51 I_v = speye(m_v);
52
53 D1_u = sparse(ops_u.derivatives.D1);
54 D2_u = ops_u.derivatives.D2;
55 H_u = sparse(ops_u.norms.H);
56 Hi_u = sparse(ops_u.norms.HI);
57 % M_u = sparse(ops_u.norms.M);
58 e_l_u = sparse(ops_u.boundary.e_1);
59 e_r_u = sparse(ops_u.boundary.e_m);
60 d1_l_u = sparse(ops_u.boundary.S_1);
61 d1_r_u = sparse(ops_u.boundary.S_m);
62
63 D1_v = sparse(ops_v.derivatives.D1);
64 D2_v = ops_v.derivatives.D2;
65 H_v = sparse(ops_v.norms.H);
66 Hi_v = sparse(ops_v.norms.HI);
67 % M_v = sparse(ops_v.norms.M);
68 e_l_v = sparse(ops_v.boundary.e_1);
69 e_r_v = sparse(ops_v.boundary.e_m);
70 d1_l_v = sparse(ops_v.boundary.S_1);
71 d1_r_v = sparse(ops_v.boundary.S_m);
72
73
74 % Metric derivatives
75 [X,Y] = ti.map(u,v);
76
77 [x_u,x_v] = gridDerivatives(X,D1_u,D1_v);
78 [y_u,y_v] = gridDerivatives(Y,D1_u,D1_v);
79
80
81
82 J = x_u.*y_v - x_v.*y_u;
83 a11 = 1./J .* (x_v.^2 + y_v.^2); %% GÖR SOM MATRISER
84 a12 = -1./J .* (x_u.*x_v + y_u.*y_v);
85 a22 = 1./J .* (x_u.^2 + y_u.^2);
86 lambda = 1/2 * (a11 + a22 - sqrt((a11-a22).^2 + 4*a12.^2));
87
88 dof_order = reshape(1:m_u*m_v,m_v,m_u);
89
90 Duu = sparse(m_tot);
91 Dvv = sparse(m_tot);
92
93 for i = 1:m_v
94 D = D2_u(a11(i,:));
95 p = dof_order(i,:);
96 Duu(p,p) = D;
97 end
98
99 for i = 1:m_u
100 D = D2_v(a22(:,i));
101 p = dof_order(:,i);
102 Dvv(p,p) = D;
103 end
104
105 L_12 = spdiags(a12(:),0,m_tot,m_tot);
106 Du = kr(D1_u,I_v);
107 Dv = kr(I_u,D1_v);
108
109 Duv = Du*L_12*Dv;
110 Dvu = Dv*L_12*Du;
111
112
113
114 obj.H = kr(H_u,H_v);
115 obj.Hi = kr(Hi_u,Hi_v);
116 obj.Hu = kr(H_u,I_v);
117 obj.Hv = kr(I_u,H_v);
118 obj.Hiu = kr(Hi_u,I_v);
119 obj.Hiv = kr(I_u,Hi_v);
120
121 % obj.M = kr(M_u,H_v)+kr(H_u,M_v);
122 obj.e_w = kr(e_l_u,I_v);
123 obj.e_e = kr(e_r_u,I_v);
124 obj.e_s = kr(I_u,e_l_v);
125 obj.e_n = kr(I_u,e_r_v);
126 obj.du_w = kr(d1_l_u,I_v);
127 obj.dv_w = (obj.e_w'*Dv)';
128 obj.du_e = kr(d1_r_u,I_v);
129 obj.dv_e = (obj.e_e'*Dv)';
130 obj.du_s = (obj.e_s'*Du)';
131 obj.dv_s = kr(I_u,d1_l_v);
132 obj.du_n = (obj.e_n'*Du)';
133 obj.dv_n = kr(I_u,d1_r_v);
134
135 obj.m = m;
136 obj.h = [h_u h_v];
137 obj.order = order;
138
139
140 obj.c = c;
141 obj.J = spdiags(J(:),0,m_tot,m_tot);
142 obj.Ji = spdiags(1./J(:),0,m_tot,m_tot);
143 obj.a11 = a11;
144 obj.a12 = a12;
145 obj.a22 = a22;
146 obj.D = obj.Ji*c^2*(Duu + Duv + Dvu + Dvv);
147 obj.u = u;
148 obj.v = v;
149 obj.X = X;
150 obj.Y = Y;
151 obj.x = X(:);
152 obj.y = Y(:);
153
154 obj.gamm_u = h_u*ops_u.borrowing.M.S;
155 obj.gamm_v = h_v*ops_v.borrowing.M.S;
156 end
157
158
159 % Closure functions return the opertors applied to the own doamin to close the boundary
160 % Penalty functions return the opertors to force the solution. In the case of an interface it returns the operator applied to the other doamin.
161 % boundary is a string specifying the boundary e.g. 'l','r' or 'e','w','n','s'.
162 % type is a string specifying the type of boundary condition if there are several.
163 % data is a function returning the data that should be applied at the boundary.
164 % neighbour_scheme is an instance of Scheme that should be interfaced to.
165 % neighbour_boundary is a string specifying which boundary to interface to.
166 function [closure, penalty] = boundary_condition(obj,boundary,type,data)
167 default_arg('type','neumann');
168 default_arg('data',0);
169
170 [e, d_n, d_t, coeff_n, coeff_t, s, gamm, halfnorm_inv] = obj.get_boundary_ops(boundary);
171
172 switch type
173 % Dirichlet boundary condition
174 case {'D','d','dirichlet'}
175 error('not implemented')
176 alpha = obj.alpha;
177
178 % tau1 < -alpha^2/gamma
179 tuning = 1.1;
180 tau1 = -tuning*alpha/gamm;
181 tau2 = s*alpha;
182
183 p = tau1*e + tau2*d;
184
185 closure = halfnorm_inv*p*e';
186
187 pp = halfnorm_inv*p;
188 switch class(data)
189 case 'double'
190 penalty = pp*data;
191 case 'function_handle'
192 penalty = @(t)pp*data(t);
193 otherwise
194 error('Weird data argument!')
195 end
196
197
198 % Neumann boundary condition
199 case {'N','n','neumann'}
200 c = obj.c;
201
202
203 a_n = spdiags(coeff_n,0,length(coeff_n),length(coeff_n));
204 a_t = spdiags(coeff_t,0,length(coeff_t),length(coeff_t));
205 d = (a_n * d_n' + a_t*d_t')';
206
207 tau1 = -s;
208 tau2 = 0;
209 tau = c.^2 * obj.Ji*(tau1*e + tau2*d);
210
211 closure = halfnorm_inv*tau*d';
212
213 pp = halfnorm_inv*tau;
214 switch class(data)
215 case 'double'
216 penalty = pp*data;
217 case 'function_handle'
218 penalty = @(t)pp*data(t);
219 otherwise
220 error('Weird data argument!')
221 end
222
223 % Unknown, boundary condition
224 otherwise
225 error('No such boundary condition: type = %s',type);
226 end
227 end
228
229 function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary)
230 % u denotes the solution in the own domain
231 % v denotes the solution in the neighbour domain
232 [e_u, d_n_u, d_t_u, coeff_n_u, coeff_t_u, s_u, gamm_u, halfnorm_inv_u_n, halfnorm_inv_u_t, halfnorm_u_t] = obj.get_boundary_ops(boundary);
233 [e_v, d_n_v, d_t_v, coeff_n_v, coeff_t_v, s_v, gamm_v, halfnorm_inv_v_n, halfnorm_inv_v_t, halfnorm_v_t] = neighbour_scheme.get_boundary_ops(boundary);
234
235 F_u = s_u * a_n_u * d_n_u' + s_u * a_t_u*d_t_u';
236 F_v = s_v * a_n_v * d_n_v' + s_v * a_t_v*d_t_v';
237
238 tau = 111111111111111111111111111111;
239 sig1 = 1/2;
240 sig2 = -1/2;
241
242 penalty_parameter_1 = s_u*halfnorm_inv_u_n*(tau + sig1*halfnorm_inv_u_t*F_u'*halfnorm_u_t)*e_u;
243 penalty_parameter_2 = halfnorm_inv_u_n * sig2 * e_u;
244
245 tuning = 1.2;
246
247 alpha_u = obj.alpha;
248 alpha_v = neighbour_scheme.alpha;
249
250 % tau1 < -(alpha_u/gamm_u + alpha_v/gamm_v)
251
252 tau1 = -(alpha_u/gamm_u + alpha_v/gamm_v) * tuning;
253 tau2 = s_u*1/2*alpha_u;
254 sig1 = s_u*(-1/2);
255 sig2 = 0;
256
257 tau = tau1*e_u + tau2*d_u;
258 sig = sig1*e_u + sig2*d_u;
259
260 closure = halfnorm_inv*( tau*e_u' + sig*alpha_u*d_u');
261 penalty = halfnorm_inv*(-tau*e_v' - sig*alpha_v*d_v');
262 end
263
264 % Ruturns the boundary ops and sign for the boundary specified by the string boundary.
265 % The right boundary is considered the positive boundary
266 function [e, d_n, d_t, coeff_n, coeff_t, s, gamm, halfnorm_inv_n, halfnorm_inv_t, halfnorm_t] = get_boundary_ops(obj,boundary)
267 switch boundary
268 case 'w'
269 e = obj.e_w;
270 d_n = obj.du_w;
271 d_t = obj.dv_w;
272 s = -1;
273
274 coeff_n = obj.a11(:,1);
275 coeff_t = obj.a12(:,1);
276 case 'e'
277 e = obj.e_e;
278 d_n = obj.du_e;
279 d_t = obj.dv_e;
280 s = 1;
281
282 coeff_n = obj.a11(:,end);
283 coeff_t = obj.a12(:,end);
284 case 's'
285 e = obj.e_s;
286 d_n = obj.dv_s;
287 d_t = obj.du_s;
288 s = -1;
289
290 coeff_n = obj.a22(1,:)';
291 coeff_t = obj.a12(1,:)';
292 case 'n'
293 e = obj.e_n;
294 d_n = obj.dv_n;
295 d_t = obj.du_n;
296 s = 1;
297
298 coeff_n = obj.a22(end,:)';
299 coeff_t = obj.a12(end,:)';
300 otherwise
301 error('No such boundary: boundary = %s',boundary);
302 end
303
304 switch boundary
305 case {'w','e'}
306 halfnorm_inv_n = obj.Hiu;
307 halfnorm_inv_t = obj.Hiv;
308 halfnorm_t = obj.Hv;
309 gamm = obj.gamm_u;
310 case {'s','n'}
311 halfnorm_inv_n = obj.Hiv;
312 halfnorm_inv_t = obj.Hiu;
313 halfnorm_t = obj.Hu;
314 gamm = obj.gamm_v;
315 end
316 end
317
318 function N = size(obj)
319 N = prod(obj.m);
320 end
321
322 end
323
324 methods(Static)
325 % Calculates the matrcis need for the inteface coupling between boundary bound_u of scheme schm_u
326 % and bound_v of scheme schm_v.
327 % [uu, uv, vv, vu] = inteface_couplong(A,'r',B,'l')
328 function [uu, uv, vv, vu] = interface_coupling(schm_u,bound_u,schm_v,bound_v)
329 [uu,uv] = schm_u.interface(bound_u,schm_v,bound_v);
330 [vv,vu] = schm_v.interface(bound_v,schm_u,bound_u);
331 end
332 end
333 end