Mercurial > repos > public > sbplib

comparison +scheme/LaplaceCurvilinear.m @ 450:8d455e49364f feature/grids

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Copy Wave2dCurve to new scheme LaplaceCurvilinear
author Jonatan Werpers <jonatan@werpers.com>
date Tue, 09 May 2017 14:58:27 +0200
parents +scheme/Wave2dCurve.m@0707a7192bc3
children 56eb7c088bd4
comparison
equal deleted inserted replaced
449:0707a7192bc3 450:8d455e49364f
1 classdef LaplaceCurvilinear < scheme.Scheme
2 properties
3 m % Number of points in each direction, possibly a vector
4 h % Grid spacing
5
6 grid
7
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 lambda
29
30 Dx, Dy % Physical derivatives
31
32 x_u
33 x_v
34 y_u
35 y_v
36 end
37
38 methods
39 function obj = LaplaceCurvilinear(g ,order, c, opSet)
40 default_arg('opSet',@sbp.D2Variable);
41 default_arg('c', 1);
42
43 assert(isa(g, 'grid.Curvilinear'))
44
45 m = g.size();
46 m_u = m(1);
47 m_v = m(2);
48 m_tot = g.N();
49
50 h = g.scaling();
51 h_u = h(1);
52 h_v = h(2);
53
54 % Operators
55 ops_u = opSet(m_u, {0, 1}, order);
56 ops_v = opSet(m_v, {0, 1}, order);
57
58 I_u = speye(m_u);
59 I_v = speye(m_v);
60
61 D1_u = ops_u.D1;
62 D2_u = ops_u.D2;
63 H_u = ops_u.H;
64 Hi_u = ops_u.HI;
65 e_l_u = ops_u.e_l;
66 e_r_u = ops_u.e_r;
67 d1_l_u = ops_u.d1_l;
68 d1_r_u = ops_u.d1_r;
69
70 D1_v = ops_v.D1;
71 D2_v = ops_v.D2;
72 H_v = ops_v.H;
73 Hi_v = ops_v.HI;
74 e_l_v = ops_v.e_l;
75 e_r_v = ops_v.e_r;
76 d1_l_v = ops_v.d1_l;
77 d1_r_v = ops_v.d1_r;
78
79 Du = kr(D1_u,I_v);
80 Dv = kr(I_u,D1_v);
81
82 % Metric derivatives
83 coords = g.points();
84 x = coords(:,1);
85 y = coords(:,2);
86
87 x_u = Du*x;
88 x_v = Dv*x;
89 y_u = Du*y;
90 y_v = Dv*y;
91
92 J = x_u.*y_v - x_v.*y_u;
93 a11 = 1./J .* (x_v.^2 + y_v.^2);
94 a12 = -1./J .* (x_u.*x_v + y_u.*y_v);
95 a22 = 1./J .* (x_u.^2 + y_u.^2);
96 lambda = 1/2 * (a11 + a22 - sqrt((a11-a22).^2 + 4*a12.^2));
97
98 % Assemble full operators
99 L_12 = spdiags(a12, 0, m_tot, m_tot);
100 Duv = Du*L_12*Dv;
101 Dvu = Dv*L_12*Du;
102
103 Duu = sparse(m_tot);
104 Dvv = sparse(m_tot);
105 ind = grid.funcToMatrix(g, 1:m_tot);
106
107 for i = 1:m_v
108 D = D2_u(a11(ind(:,i)));
109 p = ind(:,i);
110 Duu(p,p) = D;
111 end
112
113 for i = 1:m_u
114 D = D2_v(a22(ind(i,:)));
115 p = ind(i,:);
116 Dvv(p,p) = D;
117 end
118
119 obj.H = kr(H_u,H_v);
120 obj.Hi = kr(Hi_u,Hi_v);
121 obj.Hu = kr(H_u,I_v);
122 obj.Hv = kr(I_u,H_v);
123 obj.Hiu = kr(Hi_u,I_v);
124 obj.Hiv = kr(I_u,Hi_v);
125
126 obj.e_w = kr(e_l_u,I_v);
127 obj.e_e = kr(e_r_u,I_v);
128 obj.e_s = kr(I_u,e_l_v);
129 obj.e_n = kr(I_u,e_r_v);
130 obj.du_w = kr(d1_l_u,I_v);
131 obj.dv_w = (obj.e_w'*Dv)';
132 obj.du_e = kr(d1_r_u,I_v);
133 obj.dv_e = (obj.e_e'*Dv)';
134 obj.du_s = (obj.e_s'*Du)';
135 obj.dv_s = kr(I_u,d1_l_v);
136 obj.du_n = (obj.e_n'*Du)';
137 obj.dv_n = kr(I_u,d1_r_v);
138
139 obj.x_u = x_u;
140 obj.x_v = x_v;
141 obj.y_u = y_u;
142 obj.y_v = y_v;
143
144 obj.m = m;
145 obj.h = [h_u h_v];
146 obj.order = order;
147 obj.grid = g;
148
149 obj.c = c;
150 obj.J = spdiags(J, 0, m_tot, m_tot);
151 obj.Ji = spdiags(1./J, 0, m_tot, m_tot);
152 obj.a11 = a11;
153 obj.a12 = a12;
154 obj.a22 = a22;
155 obj.D = obj.Ji*c^2*(Duu + Duv + Dvu + Dvv);
156 obj.lambda = lambda;
157
158 obj.Dx = spdiag( y_v./J)*Du + spdiag(-y_u./J)*Dv;
159 obj.Dy = spdiag(-x_v./J)*Du + spdiag( x_u./J)*Dv;
160
161 obj.gamm_u = h_u*ops_u.borrowing.M.d1;
162 obj.gamm_v = h_v*ops_v.borrowing.M.d1;
163 end
164
165
166 % Closure functions return the opertors applied to the own doamin to close the boundary
167 % Penalty functions return the opertors to force the solution. In the case of an interface it returns the operator applied to the other doamin.
168 % boundary is a string specifying the boundary e.g. 'l','r' or 'e','w','n','s'.
169 % type is a string specifying the type of boundary condition if there are several.
170 % data is a function returning the data that should be applied at the boundary.
171 % neighbour_scheme is an instance of Scheme that should be interfaced to.
172 % neighbour_boundary is a string specifying which boundary to interface to.
173 function [closure, penalty] = boundary_condition(obj, boundary, type, parameter)
174 default_arg('type','neumann');
175 default_arg('parameter', []);
176
177 [e, d_n, d_t, coeff_n, coeff_t, s, gamm, halfnorm_inv , ~, ~, ~, scale_factor] = obj.get_boundary_ops(boundary);
178 switch type
179 % Dirichlet boundary condition
180 case {'D','d','dirichlet'}
181 % v denotes the solution in the neighbour domain
182 tuning = 1.2;
183 % tuning = 20.2;
184 [e, d_n, d_t, coeff_n, coeff_t, s, gamm, halfnorm_inv_n, halfnorm_inv_t, halfnorm_t] = obj.get_boundary_ops(boundary);
185
186 a_n = spdiag(coeff_n);
187 a_t = spdiag(coeff_t);
188
189 F = (s * a_n * d_n' + s * a_t*d_t')';
190
191 u = obj;
192
193 b1 = gamm*u.lambda./u.a11.^2;
194 b2 = gamm*u.lambda./u.a22.^2;
195
196 tau = -1./b1 - 1./b2;
197 tau = tuning * spdiag(tau);
198 sig1 = 1;
199
200 penalty_parameter_1 = halfnorm_inv_n*(tau + sig1*halfnorm_inv_t*F*e'*halfnorm_t)*e;
201
202 closure = obj.Ji*obj.c^2 * penalty_parameter_1*e';
203 penalty = -obj.Ji*obj.c^2 * penalty_parameter_1;
204
205
206 % Neumann boundary condition
207 case {'N','n','neumann'}
208 c = obj.c;
209
210 a_n = spdiags(coeff_n,0,length(coeff_n),length(coeff_n));
211 a_t = spdiags(coeff_t,0,length(coeff_t),length(coeff_t));
212 d = (a_n * d_n' + a_t*d_t')';
213
214 tau1 = -s;
215 tau2 = 0;
216 tau = c.^2 * obj.Ji*(tau1*e + tau2*d);
217
218 closure = halfnorm_inv*tau*d';
219 penalty = -halfnorm_inv*tau;
220
221 % Characteristic boundary condition
222 case {'characteristic', 'char', 'c'}
223 default_arg('parameter', 1);
224 beta = parameter;
225 c = obj.c;
226
227 a_n = spdiags(coeff_n,0,length(coeff_n),length(coeff_n));
228 a_t = spdiags(coeff_t,0,length(coeff_t),length(coeff_t));
229 d = s*(a_n * d_n' + a_t*d_t')'; % outward facing normal derivative
230
231 tau = -c.^2 * 1/beta*obj.Ji*e;
232
233 closure{1} = halfnorm_inv*tau/c*spdiag(scale_factor)*e';
234 closure{2} = halfnorm_inv*tau*beta*d';
235 penalty = -halfnorm_inv*tau;
236
237 % Unknown, boundary condition
238 otherwise
239 error('No such boundary condition: type = %s',type);
240 end
241 end
242
243 function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary)
244 % u denotes the solution in the own domain
245 % v denotes the solution in the neighbour domain
246 tuning = 1.2;
247 % tuning = 20.2;
248 [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, I_u] = obj.get_boundary_ops(boundary);
249 [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, I_v] = neighbour_scheme.get_boundary_ops(neighbour_boundary);
250
251 a_n_u = spdiag(coeff_n_u);
252 a_t_u = spdiag(coeff_t_u);
253 a_n_v = spdiag(coeff_n_v);
254 a_t_v = spdiag(coeff_t_v);
255
256 F_u = (s_u * a_n_u * d_n_u' + s_u * a_t_u*d_t_u')';
257 F_v = (s_v * a_n_v * d_n_v' + s_v * a_t_v*d_t_v')';
258
259 u = obj;
260 v = neighbour_scheme;
261
262 b1_u = gamm_u*u.lambda(I_u)./u.a11(I_u).^2;
263 b2_u = gamm_u*u.lambda(I_u)./u.a22(I_u).^2;
264 b1_v = gamm_v*v.lambda(I_v)./v.a11(I_v).^2;
265 b2_v = gamm_v*v.lambda(I_v)./v.a22(I_v).^2;
266
267 tau = -1./(4*b1_u) -1./(4*b1_v) -1./(4*b2_u) -1./(4*b2_v);
268 tau = tuning * spdiag(tau);
269 sig1 = 1/2;
270 sig2 = -1/2;
271
272 penalty_parameter_1 = halfnorm_inv_u_n*(e_u*tau + sig1*halfnorm_inv_u_t*F_u*e_u'*halfnorm_u_t*e_u);
273 penalty_parameter_2 = halfnorm_inv_u_n * sig2 * e_u;
274
275
276 closure = obj.Ji*obj.c^2 * ( penalty_parameter_1*e_u' + penalty_parameter_2*F_u');
277 penalty = obj.Ji*obj.c^2 * (-penalty_parameter_1*e_v' + penalty_parameter_2*F_v');
278 end
279
280 % Ruturns the boundary ops and sign for the boundary specified by the string boundary.
281 % The right boundary is considered the positive boundary
282 %
283 % I -- the indecies of the boundary points in the grid matrix
284 function [e, d_n, d_t, coeff_n, coeff_t, s, gamm, halfnorm_inv_n, halfnorm_inv_t, halfnorm_t, I, scale_factor] = get_boundary_ops(obj, boundary)
285
286 % gridMatrix = zeros(obj.m(2),obj.m(1));
287 % gridMatrix(:) = 1:numel(gridMatrix);
288
289 ind = grid.funcToMatrix(obj.grid, 1:prod(obj.m));
290
291 switch boundary
292 case 'w'
293 e = obj.e_w;
294 d_n = obj.du_w;
295 d_t = obj.dv_w;
296 s = -1;
297
298 I = ind(1,:);
299 coeff_n = obj.a11(I);
300 coeff_t = obj.a12(I);
301 scale_factor = sqrt(obj.x_v(I).^2 + obj.y_v(I).^2);
302 case 'e'
303 e = obj.e_e;
304 d_n = obj.du_e;
305 d_t = obj.dv_e;
306 s = 1;
307
308 I = ind(end,:);
309 coeff_n = obj.a11(I);
310 coeff_t = obj.a12(I);
311 scale_factor = sqrt(obj.x_v(I).^2 + obj.y_v(I).^2);
312 case 's'
313 e = obj.e_s;
314 d_n = obj.dv_s;
315 d_t = obj.du_s;
316 s = -1;
317
318 I = ind(:,1)';
319 coeff_n = obj.a22(I);
320 coeff_t = obj.a12(I);
321 scale_factor = sqrt(obj.x_u(I).^2 + obj.y_u(I).^2);
322 case 'n'
323 e = obj.e_n;
324 d_n = obj.dv_n;
325 d_t = obj.du_n;
326 s = 1;
327
328 I = ind(:,end)';
329 coeff_n = obj.a22(I);
330 coeff_t = obj.a12(I);
331 scale_factor = sqrt(obj.x_u(I).^2 + obj.y_u(I).^2);
332 otherwise
333 error('No such boundary: boundary = %s',boundary);
334 end
335
336 switch boundary
337 case {'w','e'}
338 halfnorm_inv_n = obj.Hiu;
339 halfnorm_inv_t = obj.Hiv;
340 halfnorm_t = obj.Hv;
341 gamm = obj.gamm_u;
342 case {'s','n'}
343 halfnorm_inv_n = obj.Hiv;
344 halfnorm_inv_t = obj.Hiu;
345 halfnorm_t = obj.Hu;
346 gamm = obj.gamm_v;
347 end
348 end
349
350 function N = size(obj)
351 N = prod(obj.m);
352 end
353
354
355 end
356 end