Mercurial > repos > public > sbplib

comparison +rv/constructDiffOps.m @ 1026:44c3ea38097e feature/advectionRV

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Fix minor issues in rv/constructDiffOps - Fix issue where the order of the RV flux diff op could be set to less than 2 - Remove faulty assert of BCs
author Vidar Stiernström <vidar.stiernstrom@it.uu.se>
date Mon, 07 Jan 2019 16:32:40 +0100
parents defc9d0cc1f2
children 1a5c8723c9be 2d7ba44340d0
comparison
equal deleted inserted replaced
1025:ac80bedc8df7 1026:44c3ea38097e
5 [D, solutionPenalties] = constructTotalFluxDiffOp(scheme, grid, order, opSet, waveSpeed, BCs, fluxSplitting); 5 [D, solutionPenalties] = constructTotalFluxDiffOp(scheme, grid, order, opSet, waveSpeed, BCs, fluxSplitting);
6 D2 = constructSymmetricD2Operator(grid, order, opSet); 6 D2 = constructSymmetricD2Operator(grid, 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, order-2, opSet, waveSpeed, BCs, fluxSplitting); 10 [D_flux, residualPenalties] = constructTotalFluxDiffOp(scheme, grid, 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 end 17 end
18 18
19 function [D, penalties] = constructTotalFluxDiffOp(scheme, grid, order, opSet, waveSpeed, BCs, fluxSplitting) 19 function [D, penalties] = constructTotalFluxDiffOp(scheme, grid, order, opSet, waveSpeed, BCs, fluxSplitting)
20 if isequal(opSet, @sbp.D1Upwind) 20 if isequal(opSet, @sbp.D1Upwind)
21 assert(size(fluxSplitting,1) == grid.D());
22 diffOp = scheme(grid, order, opSet, waveSpeed, fluxSplitting); 21 diffOp = scheme(grid, order, opSet, waveSpeed, fluxSplitting);
23 else 22 else
24 diffOp = scheme(grid, order, opSet, waveSpeed); 23 diffOp = scheme(grid, order, opSet, waveSpeed);
25 end 24 end
26 [D, penalties] = addClosuresToDiffOp(diffOp, BCs); 25 [D, penalties] = addClosuresToDiffOp(diffOp, BCs);