Mercurial > repos > public > sbplib

comparison +scheme/Heat2dCurvilinear.m @ 868:57760d7088ad bcSetupExperiment

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Merge with feature/grids
author Jonatan Werpers <jonatan@werpers.com>
date Fri, 27 Jul 2018 10:39:12 -0700
parents 08f3ffe63f48
children 78db023a7fe3
comparison
equal deleted inserted replaced
867:d634d4deb263 868:57760d7088ad
1 classdef Heat2dCurvilinear < scheme.Scheme
2
3 % Discretizes the Laplacian with variable coefficent, curvilinear,
4 % in the Heat equation way (i.e., the discretization matrix is not necessarily
5 % symmetric)
6 % u_t = div * (kappa * grad u )
7 % opSet should be cell array of opSets, one per dimension. This
8 % is useful if we have periodic BC in one direction.
9
10 properties
11 m % Number of points in each direction, possibly a vector
12 h % Grid spacing
13
14 grid
15 dim
16
17 order % Order of accuracy for the approximation
18
19 % Diagonal matrix for variable coefficients
20 KAPPA % Variable coefficient
21
22 D % Total operator
23 D1 % First derivatives
24
25 % Second derivatives
26 D2_kappa
27
28 H, Hi % Inner products
29 e_l, e_r
30 d1_l, d1_r % Normal derivatives at the boundary
31 alpha % Vector of borrowing constants
32
33 % Boundary inner products
34 H_boundary_l, H_boundary_r
35
36 % Metric coefficients
37 b % Cell matrix of size dim x dim
38 J, Ji
39 beta % Cell array of scale factors
40
41 % Numerical boundary flux operators
42 flux_l, flux_r
43
44 end
45
46 methods
47
48 function obj = Heat2dCurvilinear(g ,order, kappa_fun, opSet)
49 default_arg('opSet',{@sbp.D2Variable, @sbp.D2Variable});
50 default_arg('kappa_fun', @(x,y) 0*x+1);
51 dim = 2;
52
53 kappa = grid.evalOn(g, kappa_fun);
54 m = g.size();
55 m_tot = g.N();
56
57 % 1D operators
58 ops = cell(dim,1);
59 for i = 1:dim
60 ops{i} = opSet{i}(m(i), {0, 1}, order);
61 end
62
63 I = cell(dim,1);
64 D1 = cell(dim,1);
65 D2 = cell(dim,1);
66 H = cell(dim,1);
67 Hi = cell(dim,1);
68 e_l = cell(dim,1);
69 e_r = cell(dim,1);
70 d1_l = cell(dim,1);
71 d1_r = cell(dim,1);
72
73 for i = 1:dim
74 I{i} = speye(m(i));
75 D1{i} = ops{i}.D1;
76 D2{i} = ops{i}.D2;
77 H{i} = ops{i}.H;
78 Hi{i} = ops{i}.HI;
79 e_l{i} = ops{i}.e_l;
80 e_r{i} = ops{i}.e_r;
81 d1_l{i} = ops{i}.d1_l;
82 d1_r{i} = ops{i}.d1_r;
83 end
84
85 %====== Assemble full operators ========
86 KAPPA = spdiag(kappa);
87 obj.KAPPA = KAPPA;
88
89 % Allocate
90 obj.D1 = cell(dim,1);
91 obj.D2_kappa = cell(dim,1);
92 obj.e_l = cell(dim,1);
93 obj.e_r = cell(dim,1);
94 obj.d1_l = cell(dim,1);
95 obj.d1_r = cell(dim,1);
96
97 % D1
98 obj.D1{1} = kron(D1{1},I{2});
99 obj.D1{2} = kron(I{1},D1{2});
100
101 % -- Metric coefficients ----
102 coords = g.points();
103 x = coords(:,1);
104 y = coords(:,2);
105
106 % Use non-periodic difference operators for metric even if opSet is periodic.
107 xmax = max(ops{1}.x);
108 ymax = max(ops{2}.x);
109 opSetMetric{1} = sbp.D2Variable(m(1), {0, xmax}, order);
110 opSetMetric{2} = sbp.D2Variable(m(2), {0, ymax}, order);
111 D1Metric{1} = kron(opSetMetric{1}.D1, I{2});
112 D1Metric{2} = kron(I{1}, opSetMetric{2}.D1);
113
114 x_xi = D1Metric{1}*x;
115 x_eta = D1Metric{2}*x;
116 y_xi = D1Metric{1}*y;
117 y_eta = D1Metric{2}*y;
118
119 J = x_xi.*y_eta - x_eta.*y_xi;
120
121 b = cell(dim,dim);
122 b{1,1} = y_eta./J;
123 b{1,2} = -x_eta./J;
124 b{2,1} = -y_xi./J;
125 b{2,2} = x_xi./J;
126
127 % Scale factors for boundary integrals
128 beta = cell(dim,1);
129 beta{1} = sqrt(x_eta.^2 + y_eta.^2);
130 beta{2} = sqrt(x_xi.^2 + y_xi.^2);
131
132 J = spdiag(J);
133 Ji = inv(J);
134 for i = 1:dim
135 beta{i} = spdiag(beta{i});
136 for j = 1:dim
137 b{i,j} = spdiag(b{i,j});
138 end
139 end
140 obj.J = J;
141 obj.Ji = Ji;
142 obj.b = b;
143 obj.beta = beta;
144 %----------------------------
145
146 % Boundary operators
147 obj.e_l{1} = kron(e_l{1},I{2});
148 obj.e_l{2} = kron(I{1},e_l{2});
149 obj.e_r{1} = kron(e_r{1},I{2});
150 obj.e_r{2} = kron(I{1},e_r{2});
151
152 obj.d1_l{1} = kron(d1_l{1},I{2});
153 obj.d1_l{2} = kron(I{1},d1_l{2});
154 obj.d1_r{1} = kron(d1_r{1},I{2});
155 obj.d1_r{2} = kron(I{1},d1_r{2});
156
157 % D2 coefficients
158 kappa_coeff = cell(dim,dim);
159 for j = 1:dim
160 obj.D2_kappa{j} = sparse(m_tot,m_tot);
161 kappa_coeff{j} = sparse(m_tot,1);
162 for i = 1:dim
163 kappa_coeff{j} = kappa_coeff{j} + b{i,j}*J*b{i,j}*kappa;
164 end
165 end
166 ind = grid.funcToMatrix(g, 1:m_tot);
167
168 % x-dir
169 j = 1;
170 for col = 1:m(2)
171 D_kappa = D2{1}(kappa_coeff{j}(ind(:,col)));
172
173 p = ind(:,col);
174 obj.D2_kappa{j}(p,p) = D_kappa;
175 end
176
177 % y-dir
178 j = 2;
179 for row = 1:m(1)
180 D_kappa = D2{2}(kappa_coeff{j}(ind(row,:)));
181
182 p = ind(row,:);
183 obj.D2_kappa{j}(p,p) = D_kappa;
184 end
185
186 % Quadratures
187 obj.H = kron(H{1},H{2});
188 obj.Hi = inv(obj.H);
189 obj.H_boundary_l = cell(dim,1);
190 obj.H_boundary_l{1} = obj.e_l{1}'*beta{1}*obj.e_l{1}*H{2};
191 obj.H_boundary_l{2} = obj.e_l{2}'*beta{2}*obj.e_l{2}*H{1};
192 obj.H_boundary_r = cell(dim,1);
193 obj.H_boundary_r{1} = obj.e_r{1}'*beta{1}*obj.e_r{1}*H{2};
194 obj.H_boundary_r{2} = obj.e_r{2}'*beta{2}*obj.e_r{2}*H{1};
195
196 %=== Differentiation matrix D (without SAT) ===
197 D2_kappa = obj.D2_kappa;
198 D1 = obj.D1;
199 D = sparse(m_tot,m_tot);
200
201 d = @kroneckerDelta; % Kronecker delta
202 db = @(i,j) 1-d(i,j); % Logical not of Kronecker delta
203
204 % 2nd derivatives
205 for j = 1:dim
206 D = D + Ji*D2_kappa{j};
207 end
208
209 % Mixed terms
210 for i = 1:dim
211 for j = 1:dim
212 for k = 1:dim
213 D = D + db(i,j)*Ji*D1{j}*b{i,j}*J*KAPPA*b{i,k}*D1{k};
214 end
215 end
216 end
217 obj.D = D;
218 %=========================================%
219
220 % Normal flux operators for BC.
221 flux_l = cell(dim,1);
222 flux_r = cell(dim,1);
223
224 d1_l = obj.d1_l;
225 d1_r = obj.d1_r;
226 e_l = obj.e_l;
227 e_r = obj.e_r;
228
229 % Loop over boundaries
230 for j = 1:dim
231 flux_l{j} = sparse(m_tot,m_tot);
232 flux_r{j} = sparse(m_tot,m_tot);
233
234 % Loop over dummy index
235 for i = 1:dim
236 % Loop over dummy index
237 for k = 1:dim
238 flux_l{j} = flux_l{j} ...
239 - beta{j}\b{i,j}*J*KAPPA*b{i,k}*( d(j,k)*e_l{k}*d1_l{k}' + db(j,k)*D1{k} );
240
241 flux_r{j} = flux_r{j} ...
242 + beta{j}\b{i,j}*J*KAPPA*b{i,k}*( d(j,k)*e_r{k}*d1_r{k}' + db(j,k)*D1{k} );
243 end
244
245 end
246 end
247 obj.flux_l = flux_l;
248 obj.flux_r = flux_r;
249
250 % Misc.
251 obj.m = m;
252 obj.h = g.scaling();
253 obj.order = order;
254 obj.grid = g;
255 obj.dim = dim;
256 obj.alpha = [ops{1}.borrowing.M.d1, ops{2}.borrowing.M.d1];
257
258 end
259
260
261 % Closure functions return the operators applied to the own domain to close the boundary
262 % Penalty functions return the operators to force the solution. In the case of an interface it returns the operator applied to the other doamin.
263 % boundary is a string specifying the boundary e.g. 'l','r' or 'e','w','n','s'.
264 % type is a string specifying the type of boundary condition.
265 % data is a function returning the data that should be applied at the boundary.
266 % neighbour_scheme is an instance of Scheme that should be interfaced to.
267 % neighbour_boundary is a string specifying which boundary to interface to.
268 function [closure, penalty] = boundary_condition(obj, boundary, type, symmetric, tuning)
269 default_arg('type','Neumann');
270 default_arg('symmetric', false);
271 default_arg('tuning',1.2);
272
273 % j is the coordinate direction of the boundary
274 % nj: outward unit normal component.
275 % nj = -1 for west, south, bottom boundaries
276 % nj = 1 for east, north, top boundaries
277 [j, nj] = obj.get_boundary_number(boundary);
278 switch nj
279 case 1
280 e = obj.e_r{j};
281 flux = obj.flux_r{j};
282 H_gamma = obj.H_boundary_r{j};
283 case -1
284 e = obj.e_l{j};
285 flux = obj.flux_l{j};
286 H_gamma = obj.H_boundary_l{j};
287 end
288
289 Hi = obj.Hi;
290 Ji = obj.Ji;
291 KAPPA = obj.KAPPA;
292 kappa_gamma = e'*KAPPA*e;
293 h = obj.h(j);
294 alpha = h*obj.alpha(j);
295
296 switch type
297
298 % Dirichlet boundary condition
299 case {'D','d','dirichlet','Dirichlet'}
300
301 if ~symmetric
302 closure = -Ji*Hi*flux'*e*H_gamma*(e' );
303 penalty = Ji*Hi*flux'*e*H_gamma;
304 else
305 closure = Ji*Hi*flux'*e*H_gamma*(e' )...
306 -tuning*2/alpha*Ji*Hi*e*kappa_gamma*H_gamma*(e' ) ;
307 penalty = -Ji*Hi*flux'*e*H_gamma ...
308 +tuning*2/alpha*Ji*Hi*e*kappa_gamma*H_gamma;
309 end
310
311 % Normal flux boundary condition
312 case {'N','n','neumann','Neumann'}
313 closure = -Ji*Hi*e*H_gamma*(e'*flux );
314 penalty = Ji*Hi*e*H_gamma;
315
316 % Unknown boundary condition
317 otherwise
318 error('No such boundary condition: type = %s',type);
319 end
320 end
321
322 function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary)
323 % u denotes the solution in the own domain
324 % v denotes the solution in the neighbour domain
325 error('Interface not implemented');
326 end
327
328 % Returns the coordinate number and outward normal component for the boundary specified by the string boundary.
329 function [j, nj] = get_boundary_number(obj, boundary)
330
331 switch boundary
332 case {'w','W','west','West', 'e', 'E', 'east', 'East'}
333 j = 1;
334 case {'s','S','south','South', 'n', 'N', 'north', 'North'}
335 j = 2;
336 otherwise
337 error('No such boundary: boundary = %s',boundary);
338 end
339
340 switch boundary
341 case {'w','W','west','West','s','S','south','South'}
342 nj = -1;
343 case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'}
344 nj = 1;
345 end
346 end
347
348 % Returns the coordinate number and outward normal component for the boundary specified by the string boundary.
349 function [return_op] = get_boundary_operator(obj, op, boundary)
350
351 switch boundary
352 case {'w','W','west','West', 'e', 'E', 'east', 'East'}
353 j = 1;
354 case {'s','S','south','South', 'n', 'N', 'north', 'North'}
355 j = 2;
356 otherwise
357 error('No such boundary: boundary = %s',boundary);
358 end
359
360 switch op
361 case 'e'
362 switch boundary
363 case {'w','W','west','West','s','S','south','South'}
364 return_op = obj.e_l{j};
365 case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'}
366 return_op = obj.e_r{j};
367 end
368 case 'd'
369 switch boundary
370 case {'w','W','west','West','s','S','south','South'}
371 return_op = obj.d1_l{j};
372 case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'}
373 return_op = obj.d1_r{j};
374 end
375 otherwise
376 error(['No such operator: operatr = ' op]);
377 end
378
379 end
380
381 function N = size(obj)
382 N = prod(obj.m);
383 end
384 end
385 end