Mercurial > repos > public > sbplib

comparison +scheme/Heat2dVariable.m @ 943:21394c78c72e feature/utux2D

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Merge with default
author Martin Almquist <malmquist@stanford.edu>
date Tue, 04 Dec 2018 15:24:36 -0800
parents b9c98661ff5d e9e15d64f9f9
children 78db023a7fe3
comparison
equal deleted inserted replaced
942:35701c85e356 943:21394c78c72e
1 classdef Heat2dVariable < scheme.Scheme 1 classdef Heat2dVariable < scheme.Scheme
2 2
3 % Discretizes the Laplacian with variable coefficent, 3 % Discretizes the Laplacian with variable coefficent,
4 % In the Heat equation way (i.e., the discretization matrix is not necessarily 4 % In the Heat equation way (i.e., the discretization matrix is not necessarily
5 % symmetric) 5 % symmetric)
6 % u_t = div * (kappa * grad u ) 6 % u_t = div * (kappa * grad u )
7 % opSet should be cell array of opSets, one per dimension. This 7 % opSet should be cell array of opSets, one per dimension. This
8 % is useful if we have periodic BC in one direction. 8 % is useful if we have periodic BC in one direction.
9 9
10 properties 10 properties
11 m % Number of points in each direction, possibly a vector 11 m % Number of points in each direction, possibly a vector
26 D2_kappa 26 D2_kappa
27 27
28 H, Hi % Inner products 28 H, Hi % Inner products
29 e_l, e_r 29 e_l, e_r
30 d1_l, d1_r % Normal derivatives at the boundary 30 d1_l, d1_r % Normal derivatives at the boundary
31 31 alpha % Vector of borrowing constants
32
32 H_boundary % Boundary inner products 33 H_boundary % Boundary inner products
33 34
34 end 35 end
35 36
36 methods 37 methods
142 obj.m = m; 143 obj.m = m;
143 obj.h = h; 144 obj.h = h;
144 obj.order = order; 145 obj.order = order;
145 obj.grid = g; 146 obj.grid = g;
146 obj.dim = dim; 147 obj.dim = dim;
148 obj.alpha = [ops{1}.borrowing.M.d1, ops{2}.borrowing.M.d1];
147 149
148 end 150 end
149 151
150 152
151 % Closure functions return the operators applied to the own domain to close the boundary 153 % Closure functions return the operators applied to the own domain to close the boundary
153 % boundary is a string specifying the boundary e.g. 'l','r' or 'e','w','n','s'. 155 % boundary is a string specifying the boundary e.g. 'l','r' or 'e','w','n','s'.
154 % type is a string specifying the type of boundary condition. 156 % type is a string specifying the type of boundary condition.
155 % data is a function returning the data that should be applied at the boundary. 157 % data is a function returning the data that should be applied at the boundary.
156 % neighbour_scheme is an instance of Scheme that should be interfaced to. 158 % neighbour_scheme is an instance of Scheme that should be interfaced to.
157 % neighbour_boundary is a string specifying which boundary to interface to. 159 % neighbour_boundary is a string specifying which boundary to interface to.
158 function [closure, penalty] = boundary_condition(obj, boundary, type, parameter) 160 function [closure, penalty] = boundary_condition(obj, boundary, type, symmetric, tuning)
159 default_arg('type','Neumann'); 161 default_arg('type','Neumann');
160 default_arg('parameter', []); 162 default_arg('symmetric', false);
163 default_arg('tuning',1.2);
161 164
162 % j is the coordinate direction of the boundary 165 % j is the coordinate direction of the boundary
163 % nj: outward unit normal component. 166 % nj: outward unit normal component.
164 % nj = -1 for west, south, bottom boundaries 167 % nj = -1 for west, south, bottom boundaries
165 % nj = 1 for east, north, top boundaries 168 % nj = 1 for east, north, top boundaries
166 [j, nj] = obj.get_boundary_number(boundary); 169 [j, nj] = obj.get_boundary_number(boundary);
167 switch nj 170 switch nj
168 case 1 171 case 1
174 end 177 end
175 178
176 Hi = obj.Hi; 179 Hi = obj.Hi;
177 H_gamma = obj.H_boundary{j}; 180 H_gamma = obj.H_boundary{j};
178 KAPPA = obj.KAPPA; 181 KAPPA = obj.KAPPA;
179 kappa_gamma = e{j}'*KAPPA*e{j}; 182 kappa_gamma = e{j}'*KAPPA*e{j};
183 h = obj.h(j);
184 alpha = h*obj.alpha(j);
180 185
181 switch type 186 switch type
182 187
183 % Dirichlet boundary condition 188 % Dirichlet boundary condition
184 case {'D','d','dirichlet','Dirichlet'} 189 case {'D','d','dirichlet','Dirichlet'}
185 closure = -nj*Hi*d{j}*kappa_gamma*H_gamma*(e{j}' ); 190
191 if ~symmetric
192 closure = -nj*Hi*d{j}*kappa_gamma*H_gamma*(e{j}' );
186 penalty = nj*Hi*d{j}*kappa_gamma*H_gamma; 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
187 200
188 % Free boundary condition 201 % Free boundary condition
189 case {'N','n','neumann','Neumann'} 202 case {'N','n','neumann','Neumann'}
190 closure = -nj*Hi*e{j}*kappa_gamma*H_gamma*(d{j}' ); 203 closure = -nj*Hi*e{j}*kappa_gamma*H_gamma*(d{j}' );
191 penalty = nj*Hi*e{j}*kappa_gamma*H_gamma; 204 penalty = Hi*e{j}*kappa_gamma*H_gamma;
205 % penalty is for normal derivative and not for derivative, hence the sign.
192 206
193 % Unknown boundary condition 207 % Unknown boundary condition
194 otherwise 208 otherwise
195 error('No such boundary condition: type = %s',type); 209 error('No such boundary condition: type = %s',type);
196 end 210 end
197 end 211 end
198 212
199 function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary,opts) 213 function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary)
200 % u denotes the solution in the own domain 214 % u denotes the solution in the own domain
201 % v denotes the solution in the neighbour domain 215 % v denotes the solution in the neighbour domain
202 error('Interface not implemented'); 216 error('Interface not implemented');
203 end 217 end
204 218