Mercurial > repos > public > sbplib

comparison +rv/constructDiffOps.m @ 1021:cc61dde120cd feature/advectionRV

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Add upwind dissipation to the operator inside Utux2d
author Vidar Stiernström <vidar.stiernstrom@it.uu.se>
date Wed, 19 Dec 2018 20:00:27 +0100
parents 5359a61cb4d9
children defc9d0cc1f2
comparison
equal deleted inserted replaced
1020:5359a61cb4d9 1021:cc61dde120cd
1 function [D_rv, D_flux, DvDt, solutionPenalties, residualPenalties] = constructDiffOps(scheme, grid, order, opSet, waveSpeed, BCs) 1 function [D_rv, D_flux, DvDt, solutionPenalties, residualPenalties] = constructDiffOps(scheme, grid, order, opSet, waveSpeed, BCs, fluxSplitting)
2 default_arg('fluxSplitting',[]);
2 assert(size(BCs,1) == grid.D()); 3 assert(size(BCs,1) == grid.D());
3 4
4 %% DiffOps for solution vector 5 %% DiffOps for solution vector
5 [D, solutionPenalties] = constructTotalFluxDiffOp(scheme, grid, order, opSet, waveSpeed, BCs); 6 [D, solutionPenalties] = constructTotalFluxDiffOp(scheme, grid, order, opSet, waveSpeed, BCs, fluxSplitting);
6 D2 = constructSymmetricD2Operator(grid, order, opSet); 7 D2 = constructSymmetricD2Operator(grid, order, opSet);
7 D_rv = @(v,viscosity)(D + D2(viscosity))*v; 8 D_rv = @(v,viscosity)(D + D2(viscosity))*v;
8 9
9 %% DiffOps for residual viscosity 10 %% DiffOps for residual viscosity
10 [D_flux, residualPenalties] = constructTotalFluxDiffOp(scheme, grid, order-2, opSet, waveSpeed, BCs); 11 [D_flux, residualPenalties] = constructTotalFluxDiffOp(scheme, grid, order-2, opSet, waveSpeed, BCs, fluxSplitting);
11 % DiffOp for flux in residual viscosity. Due to sign conventions of the implemnted schemes, we need to 12 % DiffOp for flux in residual viscosity. Due to sign conventions of the implemnted schemes, we need to
12 % change the sign. 13 % change the sign.
13 D_flux = -D_flux; 14 D_flux = -D_flux;
14 D_flux = @(v) D_flux*v; 15 D_flux = @(v) D_flux*v;
15 % DiffOp for time derivative in residual viscosity 16 % DiffOp for time derivative in residual viscosity
16 DvDt = @(v)D*v; 17 DvDt = @(v)D*v;
17 end 18 end
18 19
19 function [D, penalties] = constructTotalFluxDiffOp(scheme, grid, order, opSet, waveSpeed, BCs) 20 function [D, penalties] = constructTotalFluxDiffOp(scheme, grid, order, opSet, waveSpeed, BCs, fluxSplitting)
20 diffOp = scheme(grid, order, opSet, waveSpeed); 21 if isequal(opSet, @sbp.D1Upwind)
22 assert(size(fluxSplitting,1) == grid.D());
23 diffOp = scheme(grid, order, opSet, waveSpeed, fluxSplitting);
24 else
25 diffOp = scheme(grid, order, opSet, waveSpeed);
26 end
21 [D, penalties] = addClosuresToDiffOp(diffOp, BCs); 27 [D, penalties] = addClosuresToDiffOp(diffOp, BCs);
22 if (isequal(opSet,@sbp.D1Upwind))
23 D = D + getUpwindDissOp(diffOp, waveSpeed);
24 end
25 end 28 end
26 29
27 function [D, penalties] = addClosuresToDiffOp(diffOp, BCs) 30 function [D, penalties] = addClosuresToDiffOp(diffOp, BCs)
28 D = diffOp.D; 31 D = diffOp.D;
29 penalties = cell(size(BCs)); 32 penalties = cell(size(BCs));
33 D = D + closure; 36 D = D + closure;
34 end 37 end
35 end 38 end
36 end 39 end
37 40
38 % TODO: Pass flux splitting components here instead of hard-coding the splitting.
39 % TBD: Decide if this should be done in the scheme instead of by the user.
40 function DissOp = getUpwindDissOp(diffOp, waveSpeed)
41 switch diffOp.grid.D()
42 case 1
43 DissOp = -(abs(waveSpeed(1))*diffOp.Diss);
44 case 2
45 DissOp = -(abs(waveSpeed(1))*diffOp.DissOpx + abs(waveSpeed(2))*diffOp.DissOpy);
46 otherwise
47 error('3D not yet implemented')
48 end
49
50 end
51
52 function D2 = constructSymmetricD2Operator(grid, order, opSet) 41 function D2 = constructSymmetricD2Operator(grid, order, opSet)
53 % TODO: 42 % TODO:
54 % Currently only implemented for upwind operators. 43 % Currently only implemented for upwind operators.
55 % Remove this part once the time-dependent D2 operator is implemented for other opSets 44 % Remove this part once the time-dependent D2 operator is implemented for other opSets
45 % or if it is decided that it should only be supported for upwind operators.
56 assert(isequal(opSet,@sbp.D1Upwind)) 46 assert(isequal(opSet,@sbp.D1Upwind))
57 47
58 m = grid.size(); 48 m = grid.size();
59 ops = cell(grid.D(),1); 49 ops = cell(grid.D(),1);
60 I = cell(grid.D(),1); 50 I = cell(grid.D(),1);