Mercurial > repos > public > sbplib

comparison +scheme/Utux2D.m @ 914:50cafc4b9e40 feature/utux2D

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Make default_field overwrite if field exists but is empty. Refactor interface method in Utux2d.
author Martin Almquist <malmquist@stanford.edu>
date Sun, 25 Nov 2018 21:09:22 -0800
parents b9c98661ff5d
children 35701c85e356
comparison
equal deleted inserted replaced
913:95cd70f4b07d 914:50cafc4b9e40
18 18
19 % Boundary operators 19 % Boundary operators
20 e_w, e_e, e_s, e_n 20 e_w, e_e, e_s, e_n
21 21
22 D % Total discrete operator 22 D % Total discrete operator
23
24 % String, type of interface coupling
25 % Default: 'upwind'
26 % Other: 'centered'
27 coupling_type
28
29 % String, type of interpolation operators
30 % Default: 'AWW' (Almquist Wang Werpers)
31 % Other: 'MC' (Mattsson Carpenter)
32 interpolation_type
33
34
35 % Cell array, damping on upwstream and downstream sides.
36 interpolation_damping
37
38 end 23 end
39 24
40 25
41 methods 26 methods
42 function obj = Utux2D(g ,order, opSet, a, coupling_type, interpolation_type, interpolation_damping) 27 function obj = Utux2D(g ,order, opSet, a)
43 28
44 default_arg('interpolation_damping',{0,0});
45 default_arg('interpolation_type','AWW');
46 default_arg('coupling_type','upwind');
47 default_arg('a',1/sqrt(2)*[1, 1]); 29 default_arg('a',1/sqrt(2)*[1, 1]);
48 default_arg('opSet',@sbp.D2Standard); 30 default_arg('opSet',@sbp.D2Standard);
49 31
50 assert(isa(g, 'grid.Cartesian')) 32 assert(isa(g, 'grid.Cartesian'))
51 if iscell(a) 33 if iscell(a)
100 82
101 obj.m = m; 83 obj.m = m;
102 obj.h = [ops_x.h ops_y.h]; 84 obj.h = [ops_x.h ops_y.h];
103 obj.order = order; 85 obj.order = order;
104 obj.a = a; 86 obj.a = a;
105 obj.coupling_type = coupling_type;
106 obj.interpolation_type = interpolation_type;
107 obj.interpolation_damping = interpolation_damping;
108 obj.D = -(a{1}*obj.Dx + a{2}*obj.Dy); 87 obj.D = -(a{1}*obj.Dx + a{2}*obj.Dy);
109 88
110 end 89 end
111 % Closure functions return the opertors applied to the own domain to close the boundary 90 % Closure functions return the opertors applied to the own domain to close the boundary
112 % Penalty functions return the opertors to force the solution. In the case of an interface it returns the operator applied to the other doamin. 91 % Penalty functions return the opertors to force the solution. In the case of an interface it returns the operator applied to the other doamin.
128 tau = sigma*obj.a{2}*obj.e_s*obj.H_x; 107 tau = sigma*obj.a{2}*obj.e_s*obj.H_x;
129 closure = obj.Hi*tau*obj.e_s'; 108 closure = obj.Hi*tau*obj.e_s';
130 end 109 end
131 penalty = -obj.Hi*tau; 110 penalty = -obj.Hi*tau;
132 111
133 end 112 end
134 113
135 function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary,opts) 114 % opts Struct that specifies the interface coupling.
136 115 % Fields:
137 % Get neighbour boundary operator 116 % -- couplingType String, type of interface coupling
138 switch neighbour_boundary 117 % % Default: 'upwind'. Other: 'centered'
139 case {'e','E','east','East'} 118 % -- interpolation: struct of interpolation operators (empty for conforming grids)
140 e_neighbour = neighbour_scheme.e_e; 119 % -- interpolationDamping: damping on upstream and downstream sides.
141 m_neighbour = neighbour_scheme.m(2); 120 function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary,opts)
142 case {'w','W','west','West'} 121 if isempty(opts)
143 e_neighbour = neighbour_scheme.e_w; 122 [closure, penalty] = interfaceStandard(obj,boundary,neighbour_scheme,neighbour_boundary);
144 m_neighbour = neighbour_scheme.m(2); 123 else
145 case {'n','N','north','North'} 124 assertType(opts, 'struct');
146 e_neighbour = neighbour_scheme.e_n; 125 if isfield(opts, 'I_local2neighbor')
147 m_neighbour = neighbour_scheme.m(1); 126 [closure, penalty] = interfaceNonConforming(obj,boundary,neighbour_scheme,neighbour_boundary,opts);
148 case {'s','S','south','South'} 127 else
149 e_neighbour = neighbour_scheme.e_s; 128 [closure, penalty] = interfaceStandard(obj,boundary,neighbour_scheme,neighbour_boundary,opts);
150 m_neighbour = neighbour_scheme.m(1); 129 end
151 end 130 end
152 131 end
153 switch obj.coupling_type 132
154 133 function [closure, penalty] = interfaceStandard(obj,boundary,neighbour_scheme,neighbour_boundary,opts)
155 % Upwind coupling (energy dissipation) 134 default_arg('opts', struct);
156 case 'upwind' 135 default_field(opts, 'couplingType', 'upwind');
136 couplingType = opts.couplingType;
137
138 % Get neighbour boundary operator
139 switch neighbour_boundary
140 case {'e','E','east','East'}
141 e_neighbour = neighbour_scheme.e_e;
142 case {'w','W','west','West'}
143 e_neighbour = neighbour_scheme.e_w;
144 case {'n','N','north','North'}
145 e_neighbour = neighbour_scheme.e_n;
146 case {'s','S','south','South'}
147 e_neighbour = neighbour_scheme.e_s;
148 end
149
150 switch couplingType
151
152 % Upwind coupling (energy dissipation)
153 case 'upwind'
157 sigma_ds = -1; %"Downstream" penalty 154 sigma_ds = -1; %"Downstream" penalty
158 sigma_us = 0; %"Upstream" penalty 155 sigma_us = 0; %"Upstream" penalty
159 156
160 % Energy-preserving coupling (no energy dissipation) 157 % Energy-preserving coupling (no energy dissipation)
161 case 'centered' 158 case 'centered'
162 sigma_ds = -1/2; %"Downstream" penalty 159 sigma_ds = -1/2; %"Downstream" penalty
163 sigma_us = 1/2; %"Upstream" penalty 160 sigma_us = 1/2; %"Upstream" penalty
164 161
165 otherwise 162 otherwise
166 error(['Interface coupling type ' coupling_type ' is not available.']) 163 error(['Interface coupling type ' couplingType ' is not available.'])
164 end
165
166 switch boundary
167 case {'w','W','west','West'}
168 tau = sigma_ds*obj.a{1}*obj.e_w*obj.H_y;
169 closure = obj.Hi*tau*obj.e_w';
170 penalty = -obj.Hi*tau*e_neighbour';
171 case {'e','E','east','East'}
172 tau = sigma_us*obj.a{1}*obj.e_e*obj.H_y;
173 closure = obj.Hi*tau*obj.e_e';
174 penalty = -obj.Hi*tau*e_neighbour';
175 case {'s','S','south','South'}
176 tau = sigma_ds*obj.a{2}*obj.e_s*obj.H_x;
177 closure = obj.Hi*tau*obj.e_s';
178 penalty = -obj.Hi*tau*e_neighbour';
179 case {'n','N','north','North'}
180 tau = sigma_us*obj.a{2}*obj.e_n*obj.H_x;
181 closure = obj.Hi*tau*obj.e_n';
182 penalty = -obj.Hi*tau*e_neighbour';
167 end 183 end
168 184
169 % Check grid ratio for interpolation 185 end
170 switch boundary 186
171 case {'w','W','west','West','e','E','east','East'} 187 function [closure, penalty] = interfaceNonConforming(obj,boundary,neighbour_scheme,neighbour_boundary,opts)
172 m = obj.m(2); 188 default_arg('opts', struct);
173 case {'s','S','south','South','n','N','north','North'} 189 default_field(opts, 'couplingType', 'upwind');
174 m = obj.m(1); 190 default_field(opts, 'interpolationDamping', {0,0});
175 end 191 couplingType = opts.couplingType;
176 grid_ratio = m/m_neighbour; 192 interpolationDamping = opts.interpolationDamping;
177 if grid_ratio ~= 1 193
178 194 % Get neighbour boundary operator
179 [ms, index] = sort([m, m_neighbour]); 195 switch neighbour_boundary
180 orders = [obj.order, neighbour_scheme.order]; 196 case {'e','E','east','East'}
181 orders = orders(index); 197 e_neighbour = neighbour_scheme.e_e;
182 198 case {'w','W','west','West'}
183 switch obj.interpolation_type 199 e_neighbour = neighbour_scheme.e_w;
184 case 'MC' 200 case {'n','N','north','North'}
185 interpOpSet = sbp.InterpMC(ms(1),ms(2),orders(1),orders(2)); 201 e_neighbour = neighbour_scheme.e_n;
186 if grid_ratio < 1 202 case {'s','S','south','South'}
187 I_neighbour2local_us = interpOpSet.IF2C; 203 e_neighbour = neighbour_scheme.e_s;
188 I_neighbour2local_ds = interpOpSet.IF2C; 204 end
189 I_local2neighbour_us = interpOpSet.IC2F; 205
190 I_local2neighbour_ds = interpOpSet.IC2F; 206 switch couplingType
191 elseif grid_ratio > 1 207
192 I_neighbour2local_us = interpOpSet.IC2F; 208 % Upwind coupling (energy dissipation)
193 I_neighbour2local_ds = interpOpSet.IC2F; 209 case 'upwind'
194 I_local2neighbour_us = interpOpSet.IF2C; 210 sigma_ds = -1; %"Downstream" penalty
195 I_local2neighbour_ds = interpOpSet.IF2C; 211 sigma_us = 0; %"Upstream" penalty
196 end 212
197 case 'AWW' 213 % Energy-preserving coupling (no energy dissipation)
198 %String 'C2F' indicates that ICF2 is more accurate. 214 case 'centered'
199 interpOpSetF2C = sbp.InterpAWW(ms(1),ms(2),orders(1),orders(2),'F2C'); 215 sigma_ds = -1/2; %"Downstream" penalty
200 interpOpSetC2F = sbp.InterpAWW(ms(1),ms(2),orders(1),orders(2),'C2F'); 216 sigma_us = 1/2; %"Upstream" penalty
201 if grid_ratio < 1 217
202 % Local is coarser than neighbour 218 otherwise
203 I_neighbour2local_us = interpOpSetC2F.IF2C; 219 error(['Interface coupling type ' couplingType ' is not available.'])
204 I_neighbour2local_ds = interpOpSetF2C.IF2C; 220 end
205 I_local2neighbour_us = interpOpSetC2F.IC2F; 221
206 I_local2neighbour_ds = interpOpSetF2C.IC2F; 222 int_damp_us = interpolationDamping{1};
207 elseif grid_ratio > 1 223 int_damp_ds = interpolationDamping{2};
208 % Local is finer than neighbour 224
209 I_neighbour2local_us = interpOpSetF2C.IC2F; 225 % u denotes the solution in the own domain
210 I_neighbour2local_ds = interpOpSetC2F.IC2F; 226 % v denotes the solution in the neighbour domain
211 I_local2neighbour_us = interpOpSetF2C.IF2C; 227 I_local2neighbour_ds = opts.I_local2neighbor.bad;
212 I_local2neighbour_ds = interpOpSetC2F.IF2C; 228 I_local2neighbour_us = opts.I_local2neighbor.good;
213 end 229 I_neighbour2local_ds = opts.I_neighbor2local.good;
214 otherwise 230 I_neighbour2local_us = opts.I_neighbor2local.bad;
215 error(['Interpolation type ' obj.interpolation_type ... 231
216 ' is not available.' ]); 232 I_back_forth_us = I_neighbour2local_us*I_local2neighbour_us;
217 end 233 I_back_forth_ds = I_neighbour2local_ds*I_local2neighbour_ds;
218
219 else
220 % No interpolation required
221 I_neighbour2local_us = speye(m,m);
222 I_neighbour2local_ds = speye(m,m);
223 end
224
225 int_damp_us = obj.interpolation_damping{1};
226 int_damp_ds = obj.interpolation_damping{2};
227
228 I = speye(m,m);
229 I_back_forth_us = I_neighbour2local_us*I_local2neighbour_us;
230 I_back_forth_ds = I_neighbour2local_ds*I_local2neighbour_ds;
231 234
232 235
233 switch boundary 236 switch boundary
234 case {'w','W','west','West'} 237 case {'w','W','west','West'}
235 tau = sigma_ds*obj.a{1}*obj.e_w*obj.H_y; 238 tau = sigma_ds*obj.a{1}*obj.e_w*obj.H_y;
236 closure = obj.Hi*tau*obj.e_w'; 239 closure = obj.Hi*tau*obj.e_w';
237 penalty = -obj.Hi*tau*I_neighbour2local_ds*e_neighbour'; 240 penalty = -obj.Hi*tau*I_neighbour2local_ds*e_neighbour';
238 241
239 beta = int_damp_ds*obj.a{1}... 242 beta = int_damp_ds*obj.a{1}...
240 *obj.e_w*obj.H_y; 243 *obj.e_w*obj.H_y;
241 closure = closure + obj.Hi*beta*(I_back_forth_ds - I)*obj.e_w'; 244 closure = closure + obj.Hi*beta*I_back_forth_ds*obj.e_w' - obj.Hi*beta*obj.e_w';
242 case {'e','E','east','East'} 245 case {'e','E','east','East'}
243 tau = sigma_us*obj.a{1}*obj.e_e*obj.H_y; 246 tau = sigma_us*obj.a{1}*obj.e_e*obj.H_y;
244 closure = obj.Hi*tau*obj.e_e'; 247 closure = obj.Hi*tau*obj.e_e';
245 penalty = -obj.Hi*tau*I_neighbour2local_us*e_neighbour'; 248 penalty = -obj.Hi*tau*I_neighbour2local_us*e_neighbour';
246 249
247 beta = int_damp_us*obj.a{1}... 250 beta = int_damp_us*obj.a{1}...
248 *obj.e_e*obj.H_y; 251 *obj.e_e*obj.H_y;
249 closure = closure + obj.Hi*beta*(I_back_forth_us - I)*obj.e_e'; 252 closure = closure + obj.Hi*beta*I_back_forth_us*obj.e_e' - obj.Hi*beta*obj.e_e';
250 case {'s','S','south','South'} 253 case {'s','S','south','South'}
251 tau = sigma_ds*obj.a{2}*obj.e_s*obj.H_x; 254 tau = sigma_ds*obj.a{2}*obj.e_s*obj.H_x;
252 closure = obj.Hi*tau*obj.e_s'; 255 closure = obj.Hi*tau*obj.e_s';
253 penalty = -obj.Hi*tau*I_neighbour2local_ds*e_neighbour'; 256 penalty = -obj.Hi*tau*I_neighbour2local_ds*e_neighbour';
254 257
255 beta = int_damp_ds*obj.a{2}... 258 beta = int_damp_ds*obj.a{2}...
256 *obj.e_s*obj.H_x; 259 *obj.e_s*obj.H_x;
257 closure = closure + obj.Hi*beta*(I_back_forth_ds - I)*obj.e_s'; 260 closure = closure + obj.Hi*beta*I_back_forth_ds*obj.e_s' - obj.Hi*beta*obj.e_s';
258 case {'n','N','north','North'} 261 case {'n','N','north','North'}
259 tau = sigma_us*obj.a{2}*obj.e_n*obj.H_x; 262 tau = sigma_us*obj.a{2}*obj.e_n*obj.H_x;
260 closure = obj.Hi*tau*obj.e_n'; 263 closure = obj.Hi*tau*obj.e_n';
261 penalty = -obj.Hi*tau*I_neighbour2local_us*e_neighbour'; 264 penalty = -obj.Hi*tau*I_neighbour2local_us*e_neighbour';
262 265
263 beta = int_damp_us*obj.a{2}... 266 beta = int_damp_us*obj.a{2}...
264 *obj.e_n*obj.H_x; 267 *obj.e_n*obj.H_x;
265 closure = closure + obj.Hi*beta*(I_back_forth_us - I)*obj.e_n'; 268 closure = closure + obj.Hi*beta*I_back_forth_us*obj.e_n' - obj.Hi*beta*obj.e_n';
266 end 269 end
267 270
268 271
269 end 272 end
270 273