Mercurial > repos > public > sbplib

comparison +scheme/Heat2dCurvilinear.m @ 741:5a9acf282b34 feature/poroelastic

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