Mercurial > repos > public > sbplib

comparison +scheme/Utux2D.m @ 610:b7b3c11fab4d feature/utux2D

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Add interpolation to scheme
author Martin Almquist <malmquist@stanford.edu>
date Sat, 14 Oct 2017 22:38:42 -0700
parents 0f9d20dbb7ce
children 2d85f17a8aec
comparison
equal deleted inserted replaced
609:8cbecf22075b 610:b7b3c11fab4d
23 % String, type of interface coupling 23 % String, type of interface coupling
24 % Default: 'upwind' 24 % Default: 'upwind'
25 % Other: 'centered' 25 % Other: 'centered'
26 coupling_type 26 coupling_type
27 27
28 % String, type of interpolation operators
29 % Default: 'AWW' (Almquist Wang Werpers)
30 % Other: 'MC' (Mattsson Carpenter)
31 interpolation_type
32
28 33
29 end 34 end
30 35
31 36
32 methods 37 methods
33 function obj = Utux2D(g ,order, opSet, a, coupling_type) 38 function obj = Utux2D(g ,order, opSet, a, coupling_type, interpolation_type)
34 39
40 default_arg('interpolation_type','AWW');
35 default_arg('coupling_type','upwind'); 41 default_arg('coupling_type','upwind');
36 default_arg('a',1/sqrt(2)*[1, 1]); 42 default_arg('a',1/sqrt(2)*[1, 1]);
37 default_arg('opSet',@sbp.D2Standard); 43 default_arg('opSet',@sbp.D2Standard);
38 assert(isa(g, 'grid.Cartesian')) 44 assert(isa(g, 'grid.Cartesian'))
39 45
82 obj.m = m; 88 obj.m = m;
83 obj.h = [ops_x.h ops_y.h]; 89 obj.h = [ops_x.h ops_y.h];
84 obj.order = order; 90 obj.order = order;
85 obj.a = a; 91 obj.a = a;
86 obj.coupling_type = coupling_type; 92 obj.coupling_type = coupling_type;
93 obj.interpolation_type = interpolation_type;
87 obj.D = -(a(1)*obj.Dx + a(2)*obj.Dy); 94 obj.D = -(a(1)*obj.Dx + a(2)*obj.Dy);
88 95
89 end 96 end
90 % Closure functions return the opertors applied to the own domain to close the boundary 97 % Closure functions return the opertors applied to the own domain to close the boundary
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. 98 % Penalty functions return the opertors to force the solution. In the case of an interface it returns the operator applied to the other doamin.
115 122
116 % Get neighbour boundary operator 123 % Get neighbour boundary operator
117 switch neighbour_boundary 124 switch neighbour_boundary
118 case {'e','E','east','East'} 125 case {'e','E','east','East'}
119 e_neighbour = neighbour_scheme.e_e; 126 e_neighbour = neighbour_scheme.e_e;
127 m_neighbour = neighbour_scheme.m(2);
120 case {'w','W','west','West'} 128 case {'w','W','west','West'}
121 e_neighbour = neighbour_scheme.e_w; 129 e_neighbour = neighbour_scheme.e_w;
130 m_neighbour = neighbour_scheme.m(2);
122 case {'n','N','north','North'} 131 case {'n','N','north','North'}
123 e_neighbour = neighbour_scheme.e_n; 132 e_neighbour = neighbour_scheme.e_n;
133 m_neighbour = neighbour_scheme.m(1);
124 case {'s','S','south','South'} 134 case {'s','S','south','South'}
125 e_neighbour = neighbour_scheme.e_s; 135 e_neighbour = neighbour_scheme.e_s;
136 m_neighbour = neighbour_scheme.m(1);
126 end 137 end
127 138
128 switch obj.coupling_type 139 switch obj.coupling_type
129 140
130 % Upwind coupling (energy dissipation) 141 % Upwind coupling (energy dissipation)
138 sigma_us = 1/2; %"Upstream" penalty 149 sigma_us = 1/2; %"Upstream" penalty
139 150
140 otherwise 151 otherwise
141 error(['Interface coupling type ' coupling_type ' is not available.']) 152 error(['Interface coupling type ' coupling_type ' is not available.'])
142 end 153 end
154
155 % Check grid ratio for interpolation
156 switch boundary
157 case {'w','W','west','West','e','E','east','East'}
158 m = obj.m(2);
159 case {'s','S','south','South','n','N','north','North'}
160 m = obj.m(1);
161 end
162 grid_ratio = m/m_neighbour;
163 if grid_ratio ~= 1
164
165 [ms, index] = sort([m, m_neighbour]);
166 orders = [obj.order, neighbour_scheme.order];
167 orders = orders(index);
168
169 switch obj.interpolation_type
170 case 'MC'
171 interpOpSet = sbp.InterpMC(ms(1),ms(2),orders(1),orders(2));
172 if grid_ratio < 1
173 I_neighbour2local_us = interpOpSet.IF2C;
174 I_neighbour2local_ds = interpOpSet.IF2C;
175 elseif grid_ratio > 1
176 I_neighbour2local_us = interpOpSet.IC2F;
177 I_neighbour2local_ds = interpOpSet.IC2F;
178 end
179 case 'AWW'
180 %String 'C2F' indicates that ICF2 is more accurate.
181 interpOpSetF2C = sbp.InterpAWW(ms(1),ms(2),orders(1),orders(2),'F2C');
182 interpOpSetC2F = sbp.InterpAWW(ms(1),ms(2),orders(1),orders(2),'C2F');
183 if grid_ratio < 1
184 % Local is coarser than neighbour
185 I_neighbour2local_us = interpOpSetC2F.IF2C;
186 I_neighbour2local_ds = interpOpSetF2C.IF2C;
187 elseif grid_ratio > 1
188 % Local is finer than neighbour
189 I_neighbour2local_us = interpOpSetF2C.IC2F;
190 I_neighbour2local_ds = interpOpSetC2F.IC2F;
191 end
192 otherwise
193 error(['Interpolation type ' obj.interpolation_type ...
194 ' is not available.' ]);
195 end
196
197 else
198 % No interpolation required
199 I_neighbour2local_us = speye(m,m);
200 I_neighbour2local_ds = speye(m,m);
201 end
143 202
144 switch boundary 203 switch boundary
145 case {'w','W','west','West'} 204 case {'w','W','west','West'}
146 tau = sigma_ds*obj.a(1)*obj.e_w*obj.H_y; 205 tau = sigma_ds*obj.a(1)*obj.e_w*obj.H_y;
147 closure = obj.Hi*tau*obj.e_w'; 206 closure = obj.Hi*tau*obj.e_w';
207 penalty = -obj.Hi*tau*I_neighbour2local_ds*e_neighbour';
148 case {'e','E','east','East'} 208 case {'e','E','east','East'}
149 tau = sigma_us*obj.a(1)*obj.e_e*obj.H_y; 209 tau = sigma_us*obj.a(1)*obj.e_e*obj.H_y;
150 closure = obj.Hi*tau*obj.e_e'; 210 closure = obj.Hi*tau*obj.e_e';
211 penalty = -obj.Hi*tau*I_neighbour2local_us*e_neighbour';
151 case {'s','S','south','South'} 212 case {'s','S','south','South'}
152 tau = sigma_ds*obj.a(2)*obj.e_s*obj.H_x; 213 tau = sigma_ds*obj.a(2)*obj.e_s*obj.H_x;
153 closure = obj.Hi*tau*obj.e_s'; 214 closure = obj.Hi*tau*obj.e_s';
215 penalty = -obj.Hi*tau*I_neighbour2local_ds*e_neighbour';
154 case {'n','N','north','North'} 216 case {'n','N','north','North'}
155 tau = sigma_us*obj.a(2)*obj.e_n*obj.H_x; 217 tau = sigma_us*obj.a(2)*obj.e_n*obj.H_x;
156 closure = obj.Hi*tau*obj.e_n'; 218 closure = obj.Hi*tau*obj.e_n';
157 end 219 penalty = -obj.Hi*tau*I_neighbour2local_us*e_neighbour';
158 penalty = -obj.Hi*tau*e_neighbour'; 220 end
221
159 222
160 end 223 end
161 224
162 function N = size(obj) 225 function N = size(obj)
163 N = obj.m; 226 N = obj.m;