Mercurial > repos > public > sbplib

comparison +rv/constructDiffOps.m @ 1032:1a5c8723c9be rv-interpolation

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Created a new branch for attempting to improve rv, via interpolation
author Vidar Stiernström <vidar.stiernstrom@it.uu.se>
date Thu, 17 Jan 2019 10:42:14 +0100
parents 44c3ea38097e
children
comparison
equal deleted inserted replaced
1031:2ef20d00b386 1032:1a5c8723c9be
1 function [D_rv, D_flux, DvDt, solutionPenalties, residualPenalties] = constructDiffOps(scheme, grid, order, opSet, waveSpeed, BCs, fluxSplitting) 1 function [D_rv, D_flux, DvDt, solutionPenalties, residualPenalties] = constructDiffOps(scheme, g, order, opSet, waveSpeed, BCs, fluxSplitting)
2 default_arg('fluxSplitting',[]); 2 default_arg('fluxSplitting',[]);
3 3
4 %% DiffOps for solution vector 4 %% DiffOps for solution vector
5 [D, solutionPenalties] = constructTotalFluxDiffOp(scheme, grid, order, opSet, waveSpeed, BCs, fluxSplitting); 5 [D, solutionPenalties] = constructTotalFluxDiffOp(scheme, g, order, opSet, waveSpeed, BCs, fluxSplitting);
6 D2 = constructSymmetricD2Operator(grid, order, opSet); 6 D2 = constructSymmetricD2Operator(g, order, opSet);
7 D_rv = @(v,viscosity)(D + D2(viscosity))*v; 7 D_rv = @(v,viscosity)(D + D2(viscosity))*v;
8 8
9 %% DiffOps for residual viscosity 9 %% DiffOps for residual viscosity
10 [D_flux, residualPenalties] = constructTotalFluxDiffOp(scheme, grid, max(order-2,2), opSet, waveSpeed, BCs, fluxSplitting); 10 [D_flux, residualPenalties] = constructTotalFluxDiffOp(scheme, g, max(order-2,2), opSet, waveSpeed, BCs, fluxSplitting);
11 % DiffOp for flux in residual viscosity. Due to sign conventions of the implemnted schemes, we need to 11 % DiffOp for flux in residual viscosity. Due to sign conventions of the implemnted schemes, we need to
12 % change the sign. 12 % change the sign.
13 D_flux = -D_flux; 13 D_flux = -D_flux;
14 D_flux = @(v) D_flux*v; 14 D_flux = @(v) D_flux*v;
15 % DiffOp for time derivative in residual viscosity 15 % DiffOp for time derivative in residual viscosity
16 DvDt = @(v)D*v; 16 DvDt = @(v)D*v;
17
18 %% DiffOps for residual viscosity
19 % xl = g.getBoundary('l');
20 % xr = g.getBoundary('r');
21 % xlim = {xl, xr};
22 % m_f = 2*(g.m()-1)+1;
23 % g_f = grid.equidistant(m_f, xlim);
24
25 % [DvDt, residualPenalties] = constructTotalFluxDiffOp(scheme, g_f, 9, opSet, waveSpeed, BCs, fluxSplitting);
26 % % DiffOp for flux in residual viscosity. Due to sign conventions of the implemnted schemes, we need to
27 % % change the sign.
28 % D_flux = -D;
29 % D_flux = @(v) D_flux*v;
30 % % DiffOp for time derivative in residual viscosity
31 % DvDt = @(v)applyInterpDiffOp(g.points(), v, g_f.points(), DvDt);
17 end 32 end
33
34 % function u = applyInterpDiffOp(x_c, u_c, x_f, D_f)
35 % u_f = D_f*interp1(x_c,u_c,x_f,'linear');
36 % u = u_f(ismember(x_f,x_c));
37 % end
18 38
19 function [D, penalties] = constructTotalFluxDiffOp(scheme, grid, order, opSet, waveSpeed, BCs, fluxSplitting) 39 function [D, penalties] = constructTotalFluxDiffOp(scheme, grid, order, opSet, waveSpeed, BCs, fluxSplitting)
20 if isequal(opSet, @sbp.D1Upwind) 40 if isequal(opSet, @sbp.D1Upwind)
21 diffOp = scheme(grid, order, opSet, waveSpeed, fluxSplitting); 41 diffOp = scheme(grid, order, opSet, waveSpeed, fluxSplitting);
22 else 42 else