Mercurial > repos > public > sbplib

comparison +scheme/Utux2D.m @ 942:35701c85e356 feature/utux2D

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Make Utux2D work with new interface type etc.
author Martin Almquist <malmquist@stanford.edu>
date Tue, 04 Dec 2018 14:54:28 -0800
parents 50cafc4b9e40
children 3dd7f87c9a1b
comparison
equal deleted inserted replaced
941:cd74e177d256 942:35701c85e356
109 end 109 end
110 penalty = -obj.Hi*tau; 110 penalty = -obj.Hi*tau;
111 111
112 end 112 end
113 113
114 % opts Struct that specifies the interface coupling. 114 % type Struct that specifies the interface coupling.
115 % Fields: 115 % Fields:
116 % -- couplingType String, type of interface coupling 116 % -- couplingType String, type of interface coupling
117 % % Default: 'upwind'. Other: 'centered' 117 % % Default: 'upwind'. Other: 'centered'
118 % -- interpolation: struct of interpolation operators (empty for conforming grids) 118 % -- interpolation: type of interpolation, default 'none'
119 % -- interpolationDamping: damping on upstream and downstream sides. 119 % -- interpolationDamping: damping on upstream and downstream sides, when using interpolation.
120 function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary,opts) 120 % Default {0,0} gives zero damping.
121 if isempty(opts) 121 function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary,type)
122 [closure, penalty] = interfaceStandard(obj,boundary,neighbour_scheme,neighbour_boundary); 122
123 else 123 defaultType.couplingType = 'upwind';
124 assertType(opts, 'struct'); 124 defaultType.interpolation = 'none';
125 if isfield(opts, 'I_local2neighbor') 125 defaultType.interpolationDamping = {0,0};
126 [closure, penalty] = interfaceNonConforming(obj,boundary,neighbour_scheme,neighbour_boundary,opts); 126 default_struct('type', defaultType);
127 else 127
128 [closure, penalty] = interfaceStandard(obj,boundary,neighbour_scheme,neighbour_boundary,opts); 128 switch type.interpolation
129 end 129 case {'none', ''}
130 end 130 [closure, penalty] = interfaceStandard(obj,boundary,neighbour_scheme,neighbour_boundary,type);
131 end 131 case {'op','OP'}
132 132 [closure, penalty] = interfaceNonConforming(obj,boundary,neighbour_scheme,neighbour_boundary,type);
133 function [closure, penalty] = interfaceStandard(obj,boundary,neighbour_scheme,neighbour_boundary,opts) 133 otherwise
134 default_arg('opts', struct); 134 error('Unknown type of interpolation: %s ', type.interpolation);
135 default_field(opts, 'couplingType', 'upwind'); 135 end
136 couplingType = opts.couplingType; 136 end
137
138 function [closure, penalty] = interfaceStandard(obj,boundary,neighbour_scheme,neighbour_boundary,type)
139 couplingType = type.couplingType;
137 140
138 % Get neighbour boundary operator 141 % Get neighbour boundary operator
139 switch neighbour_boundary 142 switch neighbour_boundary
140 case {'e','E','east','East'} 143 case {'e','E','east','East'}
141 e_neighbour = neighbour_scheme.e_e; 144 e_neighbour = neighbour_scheme.e_e;
158 case 'centered' 161 case 'centered'
159 sigma_ds = -1/2; %"Downstream" penalty 162 sigma_ds = -1/2; %"Downstream" penalty
160 sigma_us = 1/2; %"Upstream" penalty 163 sigma_us = 1/2; %"Upstream" penalty
161 164
162 otherwise 165 otherwise
163 error(['Interface coupling type ' couplingType ' is not available.']) 166 error(['Interface coupling type ' couplingType ' is not available.'])
164 end 167 end
165 168
166 switch boundary 169 switch boundary
167 case {'w','W','west','West'} 170 case {'w','W','west','West'}
168 tau = sigma_ds*obj.a{1}*obj.e_w*obj.H_y; 171 tau = sigma_ds*obj.a{1}*obj.e_w*obj.H_y;
182 penalty = -obj.Hi*tau*e_neighbour'; 185 penalty = -obj.Hi*tau*e_neighbour';
183 end 186 end
184 187
185 end 188 end
186 189
187 function [closure, penalty] = interfaceNonConforming(obj,boundary,neighbour_scheme,neighbour_boundary,opts) 190 function [closure, penalty] = interfaceNonConforming(obj,boundary,neighbour_scheme,neighbour_boundary,type)
188 default_arg('opts', struct); 191
189 default_field(opts, 'couplingType', 'upwind'); 192 % User can request special interpolation operators by specifying type.interpOpSet
190 default_field(opts, 'interpolationDamping', {0,0}); 193 default_field(type, 'interpOpSet', @sbp.InterpOpsOP);
191 couplingType = opts.couplingType; 194
192 interpolationDamping = opts.interpolationDamping; 195 interpOpSet = type.interpOpSet;
196 couplingType = type.couplingType;
197 interpolationDamping = type.interpolationDamping;
193 198
194 % Get neighbour boundary operator 199 % Get neighbour boundary operator
195 switch neighbour_boundary 200 switch neighbour_boundary
196 case {'e','E','east','East'} 201 case {'e','E','east','East'}
197 e_neighbour = neighbour_scheme.e_e; 202 e_neighbour = neighbour_scheme.e_e;
222 int_damp_us = interpolationDamping{1}; 227 int_damp_us = interpolationDamping{1};
223 int_damp_ds = interpolationDamping{2}; 228 int_damp_ds = interpolationDamping{2};
224 229
225 % u denotes the solution in the own domain 230 % u denotes the solution in the own domain
226 % v denotes the solution in the neighbour domain 231 % v denotes the solution in the neighbour domain
227 I_local2neighbour_ds = opts.I_local2neighbor.bad; 232 % Find the number of grid points along the interface
228 I_local2neighbour_us = opts.I_local2neighbor.good; 233 switch boundary
229 I_neighbour2local_ds = opts.I_neighbor2local.good; 234 case {'w','e'}
230 I_neighbour2local_us = opts.I_neighbor2local.bad; 235 m_u = obj.m(2);
236 case {'s','n'}
237 m_u = obj.m(1);
238 end
239 m_v = size(e_neighbour, 2);
240
241 % Build interpolation operators
242 intOps = interpOpSet(m_u, m_v, obj.order, neighbour_scheme.order);
243 Iu2v = intOps.Iu2v;
244 Iv2u = intOps.Iv2u;
245
246 I_local2neighbour_ds = intOps.Iu2v.bad;
247 I_local2neighbour_us = intOps.Iu2v.good;
248 I_neighbour2local_ds = intOps.Iv2u.good;
249 I_neighbour2local_us = intOps.Iv2u.bad;
231 250
232 I_back_forth_us = I_neighbour2local_us*I_local2neighbour_us; 251 I_back_forth_us = I_neighbour2local_us*I_local2neighbour_us;
233 I_back_forth_ds = I_neighbour2local_ds*I_local2neighbour_ds; 252 I_back_forth_ds = I_neighbour2local_ds*I_local2neighbour_ds;
234 253
235 254
236 switch boundary 255 switch boundary
237 case {'w','W','west','West'} 256 case {'w','W','west','West'}
238 tau = sigma_ds*obj.a{1}*obj.e_w*obj.H_y; 257 tau = sigma_ds*obj.a{1}*obj.e_w*obj.H_y;
239 closure = obj.Hi*tau*obj.e_w'; 258 closure = obj.Hi*tau*obj.e_w';
240 penalty = -obj.Hi*tau*I_neighbour2local_ds*e_neighbour'; 259 penalty = -obj.Hi*tau*I_neighbour2local_ds*e_neighbour';
241 260
242 beta = int_damp_ds*obj.a{1}... 261 beta = int_damp_ds*obj.a{1}...
243 *obj.e_w*obj.H_y; 262 *obj.e_w*obj.H_y;
244 closure = closure + obj.Hi*beta*I_back_forth_ds*obj.e_w' - obj.Hi*beta*obj.e_w'; 263 closure = closure + obj.Hi*beta*I_back_forth_ds*obj.e_w' - obj.Hi*beta*obj.e_w';
245 case {'e','E','east','East'} 264 case {'e','E','east','East'}
246 tau = sigma_us*obj.a{1}*obj.e_e*obj.H_y; 265 tau = sigma_us*obj.a{1}*obj.e_e*obj.H_y;
247 closure = obj.Hi*tau*obj.e_e'; 266 closure = obj.Hi*tau*obj.e_e';
248 penalty = -obj.Hi*tau*I_neighbour2local_us*e_neighbour'; 267 penalty = -obj.Hi*tau*I_neighbour2local_us*e_neighbour';
249 268
250 beta = int_damp_us*obj.a{1}... 269 beta = int_damp_us*obj.a{1}...
251 *obj.e_e*obj.H_y; 270 *obj.e_e*obj.H_y;
252 closure = closure + obj.Hi*beta*I_back_forth_us*obj.e_e' - obj.Hi*beta*obj.e_e'; 271 closure = closure + obj.Hi*beta*I_back_forth_us*obj.e_e' - obj.Hi*beta*obj.e_e';
253 case {'s','S','south','South'} 272 case {'s','S','south','South'}
254 tau = sigma_ds*obj.a{2}*obj.e_s*obj.H_x; 273 tau = sigma_ds*obj.a{2}*obj.e_s*obj.H_x;
255 closure = obj.Hi*tau*obj.e_s'; 274 closure = obj.Hi*tau*obj.e_s';
256 penalty = -obj.Hi*tau*I_neighbour2local_ds*e_neighbour'; 275 penalty = -obj.Hi*tau*I_neighbour2local_ds*e_neighbour';
257 276
258 beta = int_damp_ds*obj.a{2}... 277 beta = int_damp_ds*obj.a{2}...
259 *obj.e_s*obj.H_x; 278 *obj.e_s*obj.H_x;
260 closure = closure + obj.Hi*beta*I_back_forth_ds*obj.e_s' - obj.Hi*beta*obj.e_s'; 279 closure = closure + obj.Hi*beta*I_back_forth_ds*obj.e_s' - obj.Hi*beta*obj.e_s';
261 case {'n','N','north','North'} 280 case {'n','N','north','North'}
262 tau = sigma_us*obj.a{2}*obj.e_n*obj.H_x; 281 tau = sigma_us*obj.a{2}*obj.e_n*obj.H_x;
263 closure = obj.Hi*tau*obj.e_n'; 282 closure = obj.Hi*tau*obj.e_n';
264 penalty = -obj.Hi*tau*I_neighbour2local_us*e_neighbour'; 283 penalty = -obj.Hi*tau*I_neighbour2local_us*e_neighbour';
265 284
266 beta = int_damp_us*obj.a{2}... 285 beta = int_damp_us*obj.a{2}...
267 *obj.e_n*obj.H_x; 286 *obj.e_n*obj.H_x;
268 closure = closure + obj.Hi*beta*I_back_forth_us*obj.e_n' - obj.Hi*beta*obj.e_n'; 287 closure = closure + obj.Hi*beta*I_back_forth_us*obj.e_n' - obj.Hi*beta*obj.e_n';
269 end 288 end
270 289
271 290
272 end 291 end
273 292