Mercurial > repos > public > sbplib

comparison +scheme/Heat2dVariable.m @ 812:6b83dcb46f54 feature/grids

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Merge with feature/poroelastic
author Martin Almquist <malmquist@stanford.edu>
date Fri, 27 Jul 2018 10:31:51 -0700
parents e9e15d64f9f9
children 21394c78c72e
comparison
equal deleted inserted replaced
798:e76321b89c1e 812:6b83dcb46f54
1 classdef Heat2dVariable < scheme.Scheme
2
3 % Discretizes the Laplacian with variable coefficent,
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 H_boundary % Boundary inner products
34
35 end
36
37 methods
38
39 function obj = Heat2dVariable(g ,order, kappa_fun, opSet)
40 default_arg('opSet',{@sbp.D2Variable, @sbp.D2Variable});
41 default_arg('kappa_fun', @(x,y) 0*x+1);
42 dim = 2;
43
44 assert(isa(g, 'grid.Cartesian'))
45
46 kappa = grid.evalOn(g, kappa_fun);
47 m = g.size();
48 m_tot = g.N();
49
50 h = g.scaling();
51 lim = g.lim;
52
53 % 1D operators
54 ops = cell(dim,1);
55 for i = 1:dim
56 ops{i} = opSet{i}(m(i), lim{i}, order);
57 end
58
59 I = cell(dim,1);
60 D1 = cell(dim,1);
61 D2 = cell(dim,1);
62 H = cell(dim,1);
63 Hi = cell(dim,1);
64 e_l = cell(dim,1);
65 e_r = cell(dim,1);
66 d1_l = cell(dim,1);
67 d1_r = cell(dim,1);
68
69 for i = 1:dim
70 I{i} = speye(m(i));
71 D1{i} = ops{i}.D1;
72 D2{i} = ops{i}.D2;
73 H{i} = ops{i}.H;
74 Hi{i} = ops{i}.HI;
75 e_l{i} = ops{i}.e_l;
76 e_r{i} = ops{i}.e_r;
77 d1_l{i} = ops{i}.d1_l;
78 d1_r{i} = ops{i}.d1_r;
79 end
80
81 %====== Assemble full operators ========
82 KAPPA = spdiag(kappa);
83 obj.KAPPA = KAPPA;
84
85 obj.D1 = cell(dim,1);
86 obj.D2_kappa = cell(dim,1);
87 obj.e_l = cell(dim,1);
88 obj.e_r = cell(dim,1);
89 obj.d1_l = cell(dim,1);
90 obj.d1_r = cell(dim,1);
91
92 % D1
93 obj.D1{1} = kron(D1{1},I{2});
94 obj.D1{2} = kron(I{1},D1{2});
95
96 % Boundary operators
97 obj.e_l{1} = kron(e_l{1},I{2});
98 obj.e_l{2} = kron(I{1},e_l{2});
99 obj.e_r{1} = kron(e_r{1},I{2});
100 obj.e_r{2} = kron(I{1},e_r{2});
101
102 obj.d1_l{1} = kron(d1_l{1},I{2});
103 obj.d1_l{2} = kron(I{1},d1_l{2});
104 obj.d1_r{1} = kron(d1_r{1},I{2});
105 obj.d1_r{2} = kron(I{1},d1_r{2});
106
107 % D2
108 for i = 1:dim
109 obj.D2_kappa{i} = sparse(m_tot);
110 end
111 ind = grid.funcToMatrix(g, 1:m_tot);
112
113 for i = 1:m(2)
114 D_kappa = D2{1}(kappa(ind(:,i)));
115 p = ind(:,i);
116 obj.D2_kappa{1}(p,p) = D_kappa;
117 end
118
119 for i = 1:m(1)
120 D_kappa = D2{2}(kappa(ind(i,:)));
121 p = ind(i,:);
122 obj.D2_kappa{2}(p,p) = D_kappa;
123 end
124
125 % Quadratures
126 obj.H = kron(H{1},H{2});
127 obj.Hi = inv(obj.H);
128 obj.H_boundary = cell(dim,1);
129 obj.H_boundary{1} = H{2};
130 obj.H_boundary{2} = H{1};
131
132 % Differentiation matrix D (without SAT)
133 D2_kappa = obj.D2_kappa;
134 D1 = obj.D1;
135 D = sparse(m_tot,m_tot);
136 for i = 1:dim
137 D = D + D2_kappa{i};
138 end
139 obj.D = D;
140 %=========================================%
141
142 % Misc.
143 obj.m = m;
144 obj.h = h;
145 obj.order = order;
146 obj.grid = g;
147 obj.dim = dim;
148 obj.alpha = [ops{1}.borrowing.M.d1, ops{2}.borrowing.M.d1];
149
150 end
151
152
153 % Closure functions return the operators applied to the own domain to close the boundary
154 % Penalty functions return the operators to force the solution. In the case of an interface it returns the operator applied to the other doamin.
155 % boundary is a string specifying the boundary e.g. 'l','r' or 'e','w','n','s'.
156 % type is a string specifying the type of boundary condition.
157 % data is a function returning the data that should be applied at the boundary.
158 % neighbour_scheme is an instance of Scheme that should be interfaced to.
159 % neighbour_boundary is a string specifying which boundary to interface to.
160 function [closure, penalty] = boundary_condition(obj, boundary, type, symmetric, tuning)
161 default_arg('type','Neumann');
162 default_arg('symmetric', false);
163 default_arg('tuning',1.2);
164
165 % j is the coordinate direction of the boundary
166 % nj: outward unit normal component.
167 % nj = -1 for west, south, bottom boundaries
168 % nj = 1 for east, north, top boundaries
169 [j, nj] = obj.get_boundary_number(boundary);
170 switch nj
171 case 1
172 e = obj.e_r;
173 d = obj.d1_r;
174 case -1
175 e = obj.e_l;
176 d = obj.d1_l;
177 end
178
179 Hi = obj.Hi;
180 H_gamma = obj.H_boundary{j};
181 KAPPA = obj.KAPPA;
182 kappa_gamma = e{j}'*KAPPA*e{j};
183 h = obj.h(j);
184 alpha = h*obj.alpha(j);
185
186 switch type
187
188 % Dirichlet boundary condition
189 case {'D','d','dirichlet','Dirichlet'}
190
191 if ~symmetric
192 closure = -nj*Hi*d{j}*kappa_gamma*H_gamma*(e{j}' );
193 penalty = nj*Hi*d{j}*kappa_gamma*H_gamma;
194 else
195 closure = nj*Hi*d{j}*kappa_gamma*H_gamma*(e{j}' )...
196 -tuning*2/alpha*Hi*e{j}*kappa_gamma*H_gamma*(e{j}' ) ;
197 penalty = -nj*Hi*d{j}*kappa_gamma*H_gamma ...
198 +tuning*2/alpha*Hi*e{j}*kappa_gamma*H_gamma;
199 end
200
201 % Free boundary condition
202 case {'N','n','neumann','Neumann'}
203 closure = -nj*Hi*e{j}*kappa_gamma*H_gamma*(d{j}' );
204 penalty = Hi*e{j}*kappa_gamma*H_gamma;
205 % penalty is for normal derivative and not for derivative, hence the sign.
206
207 % Unknown boundary condition
208 otherwise
209 error('No such boundary condition: type = %s',type);
210 end
211 end
212
213 function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary)
214 % u denotes the solution in the own domain
215 % v denotes the solution in the neighbour domain
216 error('Interface not implemented');
217 end
218
219 % Returns the coordinate number and outward normal component for the boundary specified by the string boundary.
220 function [j, nj] = get_boundary_number(obj, boundary)
221
222 switch boundary
223 case {'w','W','west','West', 'e', 'E', 'east', 'East'}
224 j = 1;
225 case {'s','S','south','South', 'n', 'N', 'north', 'North'}
226 j = 2;
227 otherwise
228 error('No such boundary: boundary = %s',boundary);
229 end
230
231 switch boundary
232 case {'w','W','west','West','s','S','south','South'}
233 nj = -1;
234 case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'}
235 nj = 1;
236 end
237 end
238
239 % Returns the coordinate number and outward normal component for the boundary specified by the string boundary.
240 function [return_op] = get_boundary_operator(obj, op, boundary)
241
242 switch boundary
243 case {'w','W','west','West', 'e', 'E', 'east', 'East'}
244 j = 1;
245 case {'s','S','south','South', 'n', 'N', 'north', 'North'}
246 j = 2;
247 otherwise
248 error('No such boundary: boundary = %s',boundary);
249 end
250
251 switch op
252 case 'e'
253 switch boundary
254 case {'w','W','west','West','s','S','south','South'}
255 return_op = obj.e_l{j};
256 case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'}
257 return_op = obj.e_r{j};
258 end
259 case 'd'
260 switch boundary
261 case {'w','W','west','West','s','S','south','South'}
262 return_op = obj.d1_l{j};
263 case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'}
264 return_op = obj.d1_r{j};
265 end
266 otherwise
267 error(['No such operator: operatr = ' op]);
268 end
269
270 end
271
272 function N = size(obj)
273 N = prod(obj.m);
274 end
275 end
276 end