Mercurial > repos > public > sbplib
changeset 610:b7b3c11fab4d feature/utux2D
Add interpolation to scheme
| author | Martin Almquist <malmquist@stanford.edu> |
|---|---|
| date | Sat, 14 Oct 2017 22:38:42 -0700 |
| parents | 8cbecf22075b |
| children | 2d85f17a8aec |
| files | +scheme/Utux2D.m |
| diffstat | 1 files changed, 66 insertions(+), 3 deletions(-) [+] |
line wrap: on
line diff
diff -r 8cbecf22075b -r b7b3c11fab4d +scheme/Utux2D.m --- a/+scheme/Utux2D.m Sat Oct 14 22:36:31 2017 -0700 +++ b/+scheme/Utux2D.m Sat Oct 14 22:38:42 2017 -0700 @@ -25,13 +25,19 @@ % Other: 'centered' coupling_type + % String, type of interpolation operators + % Default: 'AWW' (Almquist Wang Werpers) + % Other: 'MC' (Mattsson Carpenter) + interpolation_type + end methods - function obj = Utux2D(g ,order, opSet, a, coupling_type) + function obj = Utux2D(g ,order, opSet, a, coupling_type, interpolation_type) + default_arg('interpolation_type','AWW'); default_arg('coupling_type','upwind'); default_arg('a',1/sqrt(2)*[1, 1]); default_arg('opSet',@sbp.D2Standard); @@ -84,6 +90,7 @@ obj.order = order; obj.a = a; obj.coupling_type = coupling_type; + obj.interpolation_type = interpolation_type; obj.D = -(a(1)*obj.Dx + a(2)*obj.Dy); end @@ -117,12 +124,16 @@ switch neighbour_boundary case {'e','E','east','East'} e_neighbour = neighbour_scheme.e_e; + m_neighbour = neighbour_scheme.m(2); case {'w','W','west','West'} e_neighbour = neighbour_scheme.e_w; + m_neighbour = neighbour_scheme.m(2); case {'n','N','north','North'} e_neighbour = neighbour_scheme.e_n; + m_neighbour = neighbour_scheme.m(1); case {'s','S','south','South'} e_neighbour = neighbour_scheme.e_s; + m_neighbour = neighbour_scheme.m(1); end switch obj.coupling_type @@ -140,22 +151,74 @@ otherwise error(['Interface coupling type ' coupling_type ' is not available.']) end + + % Check grid ratio for interpolation + switch boundary + case {'w','W','west','West','e','E','east','East'} + m = obj.m(2); + case {'s','S','south','South','n','N','north','North'} + m = obj.m(1); + end + grid_ratio = m/m_neighbour; + if grid_ratio ~= 1 + + [ms, index] = sort([m, m_neighbour]); + orders = [obj.order, neighbour_scheme.order]; + orders = orders(index); + + switch obj.interpolation_type + case 'MC' + interpOpSet = sbp.InterpMC(ms(1),ms(2),orders(1),orders(2)); + if grid_ratio < 1 + I_neighbour2local_us = interpOpSet.IF2C; + I_neighbour2local_ds = interpOpSet.IF2C; + elseif grid_ratio > 1 + I_neighbour2local_us = interpOpSet.IC2F; + I_neighbour2local_ds = interpOpSet.IC2F; + end + case 'AWW' + %String 'C2F' indicates that ICF2 is more accurate. + interpOpSetF2C = sbp.InterpAWW(ms(1),ms(2),orders(1),orders(2),'F2C'); + interpOpSetC2F = sbp.InterpAWW(ms(1),ms(2),orders(1),orders(2),'C2F'); + if grid_ratio < 1 + % Local is coarser than neighbour + I_neighbour2local_us = interpOpSetC2F.IF2C; + I_neighbour2local_ds = interpOpSetF2C.IF2C; + elseif grid_ratio > 1 + % Local is finer than neighbour + I_neighbour2local_us = interpOpSetF2C.IC2F; + I_neighbour2local_ds = interpOpSetC2F.IC2F; + end + otherwise + error(['Interpolation type ' obj.interpolation_type ... + ' is not available.' ]); + end + + else + % No interpolation required + I_neighbour2local_us = speye(m,m); + I_neighbour2local_ds = speye(m,m); + end switch boundary case {'w','W','west','West'} tau = sigma_ds*obj.a(1)*obj.e_w*obj.H_y; - closure = obj.Hi*tau*obj.e_w'; + closure = obj.Hi*tau*obj.e_w'; + penalty = -obj.Hi*tau*I_neighbour2local_ds*e_neighbour'; case {'e','E','east','East'} tau = sigma_us*obj.a(1)*obj.e_e*obj.H_y; closure = obj.Hi*tau*obj.e_e'; + penalty = -obj.Hi*tau*I_neighbour2local_us*e_neighbour'; case {'s','S','south','South'} tau = sigma_ds*obj.a(2)*obj.e_s*obj.H_x; closure = obj.Hi*tau*obj.e_s'; + penalty = -obj.Hi*tau*I_neighbour2local_ds*e_neighbour'; case {'n','N','north','North'} tau = sigma_us*obj.a(2)*obj.e_n*obj.H_x; closure = obj.Hi*tau*obj.e_n'; + penalty = -obj.Hi*tau*I_neighbour2local_us*e_neighbour'; end - penalty = -obj.Hi*tau*e_neighbour'; + end