Mercurial > repos > public > sbplib

comparison +rv/constructDiffOps.m @ 1020:5359a61cb4d9 feature/advectionRV

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Add utility for constructing the operators used by a discretization emplying RV-stabilization - Add constructDiffOps which creates the operatores used for a scheme with RV-stabilization - Minor clean up in RungekuttaExteriorRV
author Vidar Stiernström <vidar.stiernstrom@it.uu.se>
date Fri, 14 Dec 2018 18:33:10 +0100
parents
children cc61dde120cd
comparison
equal deleted inserted replaced
1019:f029b97dbc72 1020:5359a61cb4d9
1 function [D_rv, D_flux, DvDt, solutionPenalties, residualPenalties] = constructDiffOps(scheme, grid, order, opSet, waveSpeed, BCs)
2 assert(size(BCs,1) == grid.D());
3
4 %% DiffOps for solution vector
5 [D, solutionPenalties] = constructTotalFluxDiffOp(scheme, grid, order, opSet, waveSpeed, BCs);
6 D2 = constructSymmetricD2Operator(grid, order, opSet);
7 D_rv = @(v,viscosity)(D + D2(viscosity))*v;
8
9 %% DiffOps for residual viscosity
10 [D_flux, residualPenalties] = constructTotalFluxDiffOp(scheme, grid, order-2, opSet, waveSpeed, BCs);
11 % DiffOp for flux in residual viscosity. Due to sign conventions of the implemnted schemes, we need to
12 % change the sign.
13 D_flux = -D_flux;
14 D_flux = @(v) D_flux*v;
15 % DiffOp for time derivative in residual viscosity
16 DvDt = @(v)D*v;
17 end
18
19 function [D, penalties] = constructTotalFluxDiffOp(scheme, grid, order, opSet, waveSpeed, BCs)
20 diffOp = scheme(grid, order, opSet, waveSpeed);
21 [D, penalties] = addClosuresToDiffOp(diffOp, BCs);
22 if (isequal(opSet,@sbp.D1Upwind))
23 D = D + getUpwindDissOp(diffOp, waveSpeed);
24 end
25 end
26
27 function [D, penalties] = addClosuresToDiffOp(diffOp, BCs)
28 D = diffOp.D;
29 penalties = cell(size(BCs));
30 for i = 1:size(BCs,1)
31 for j = 1:size(BCs,2)
32 [closure, penalties{i,j}] = diffOp.boundary_condition(BCs{i,j}.boundary, BCs{i,j}.type);
33 D = D + closure;
34 end
35 end
36 end
37
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)
53 % TODO:
54 % Currently only implemented for upwind operators.
55 % Remove this part once the time-dependent D2 operator is implemented for other opSets
56 assert(isequal(opSet,@sbp.D1Upwind))
57
58 m = grid.size();
59 ops = cell(grid.D(),1);
60 I = cell(grid.D(),1);
61 for i = 1:grid.D()
62 lim = {grid.x{i}(1), grid.x{i}(end)};
63 ops{i} = opSet(m(i), lim, order);
64 I{i} = speye(m(i));
65 end
66
67 % TBD: How is this generalized to a loop over dimensions or similar?
68 switch grid.D()
69 case 1
70 e_r = ops{1}.e_r;
71 e_l = ops{1}.e_l;
72 Dm = ops{1}.Dm;
73 Dp = ops{1}.Dp;
74 Hi = ops{1}.HI;
75 B = e_r*e_r' - e_l*e_l';
76 D2 = @(viscosity) Dm*spdiag(viscosity)*Dp-Hi*(B*spdiag(viscosity)*Dp);
77 case 2
78 e_e = kron(ops{1}.e_r,I{2});
79 e_w = kron(ops{1}.e_l,I{2});
80 Dm_x = kron(ops{1}.Dm,I{2});
81 Dp_x = kron(ops{1}.Dp,I{2});
82 H_x = kron(ops{1}.HI,I{2});
83 B_x = e_e*e_e' - e_w*e_w';
84 D2_x = @(viscosity) Dm_x*spdiag(viscosity)*Dp_x-H_x*(B_x*spdiag(viscosity)*Dp_x);
85
86 e_n = kron(I{1},ops{2}.e_r);
87 e_s = kron(I{1},ops{2}.e_l);
88 Dm_y = kron(I{1},ops{2}.Dm);
89 Dp_y = kron(I{1},ops{2}.Dp);
90 H_y = kron(I{1},ops{2}.HI);
91 B_y = e_n*e_n' - e_s*e_s';
92 D2_y = @(viscosity) Dm_y*spdiag(viscosity)*Dp_y-H_y*(B_y*spdiag(viscosity)*Dp_y);
93 D2 = @(viscosity)D2_x(viscosity) + D2_y(viscosity);
94 otherwise
95 error('3D not yet implemented')
96 end
97 end