Mercurial > repos > public > sbplib

comparison +rv/constructDiffOps.m @ 1037:2d7ba44340d0 feature/burgers1d

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Pass scheme specific parameters as cell array. This will enabale constructDiffOps to be more general. In addition, allow for schemes returning function handles as diffOps, which is currently how non-linear schemes such as Burgers1d are implemented.
author Vidar Stiernström <vidar.stiernstrom@it.uu.se>
date Fri, 18 Jan 2019 09:02:02 +0100
parents 44c3ea38097e
children 52f59d27b40f
comparison
equal deleted inserted replaced
1036:8a9393084b30 1037:2d7ba44340d0
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, schemeParams, opSet, BCs)
2 default_arg('fluxSplitting',[]);
3
4 %% DiffOps for solution vector 2 %% DiffOps for solution vector
5 [D, solutionPenalties] = constructTotalFluxDiffOp(scheme, grid, order, opSet, waveSpeed, BCs, fluxSplitting); 3 [D, solutionPenalties] = constructTotalFluxDiffOp(scheme, g, order, schemeParams, opSet, BCs);
6 D2 = constructSymmetricD2Operator(grid, order, opSet); 4 D2 = constructSymmetricD2Operator(g, order, opSet);
7 D_rv = @(v,viscosity)(D + D2(viscosity))*v; 5 D_rv = @(v,viscosity)(D(v) + D2(v, viscosity));
8 6
9 %% DiffOps for residual viscosity 7 %% DiffOps for residual viscosity
10 [D_flux, residualPenalties] = constructTotalFluxDiffOp(scheme, grid, max(order-2,2), opSet, waveSpeed, BCs, fluxSplitting); 8 [D_flux, residualPenalties] = constructTotalFluxDiffOp(scheme, g, max(order-2,2), schemeParams, opSet, BCs);
11 % DiffOp for flux in residual viscosity. Due to sign conventions of the implemnted schemes, we need to 9 % DiffOp for flux in residual viscosity. Due to sign conventions of the implemnted schemes, we need to
12 % change the sign. 10 % change the sign.
13 D_flux = -D_flux; 11 D_flux = @(v) -D_flux(v);
14 D_flux = @(v) D_flux*v;
15 % DiffOp for time derivative in residual viscosity 12 % DiffOp for time derivative in residual viscosity
16 DvDt = @(v)D*v; 13 DvDt = D;
17 end 14 end
18 15
19 function [D, penalties] = constructTotalFluxDiffOp(scheme, grid, order, opSet, waveSpeed, BCs, fluxSplitting) 16 function [D, penalties] = constructTotalFluxDiffOp(scheme, g, order, schemeParams, opSet, BCs)
20 if isequal(opSet, @sbp.D1Upwind) 17 diffOp = scheme(g, order, schemeParams{:}, opSet);
21 diffOp = scheme(grid, order, opSet, waveSpeed, fluxSplitting);
22 else
23 diffOp = scheme(grid, order, opSet, waveSpeed);
24 end
25 [D, penalties] = addClosuresToDiffOp(diffOp, BCs); 18 [D, penalties] = addClosuresToDiffOp(diffOp, BCs);
26 end 19 end
27 20
28 function [D, penalties] = addClosuresToDiffOp(diffOp, BCs) 21 function [D, penalties] = addClosuresToDiffOp(diffOp, BCs)
29 D = diffOp.D; 22 if ~isa(diffOp.D, 'function_handle')
23 D = @(v) diffOp.D*v
24 else
25 D = diffOp.D;
26 end
30 penalties = cell(size(BCs)); 27 penalties = cell(size(BCs));
31 for i = 1:size(BCs,1) 28 for i = 1:size(BCs,1)
32 for j = 1:size(BCs,2) 29 for j = 1:size(BCs,2)
33 [closure, penalties{i,j}] = diffOp.boundary_condition(BCs{i,j}.boundary, BCs{i,j}.type); 30 [closure, penalties{i,j}] = diffOp.boundary_condition(BCs{i,j}.boundary, BCs{i,j}.type);
34 D = D + closure; 31 if ~isa(closure, 'function_handle')
32 closure = @(v) closure*v;
33 end
34 D = @(v) D(v) + closure(v);
35 end 35 end
36 end 36 end
37 end 37 end
38 38
39 function D2 = constructSymmetricD2Operator(grid, order, opSet) 39 function D2 = constructSymmetricD2Operator(g, order, opSet)
40 % TODO:
41 % Currently only implemented for upwind operators.
42 % Remove this part once the time-dependent D2 operator is implemented for other opSets
43 % or if it is decided that it should only be supported for upwind operators.
44 assert(isequal(opSet,@sbp.D1Upwind))
45 40
46 m = grid.size(); 41
47 ops = cell(grid.D(),1); 42 m = g.size();
48 I = cell(grid.D(),1); 43 ops = cell(g.D(),1);
49 for i = 1:grid.D() 44 I = cell(g.D(),1);
50 lim = {grid.x{i}(1), grid.x{i}(end)}; 45 for i = 1:g.D()
46 lim = {g.x{i}(1), g.x{i}(end)};
51 ops{i} = opSet(m(i), lim, order); 47 ops{i} = opSet(m(i), lim, order);
52 I{i} = speye(m(i)); 48 I{i} = speye(m(i));
53 end 49 end
54 50
55 % TBD: How is this generalized to a loop over dimensions or similar? 51 % TBD: How is this generalized to a loop over dimensions or similar?
56 switch grid.D() 52 switch g.D()
57 case 1 53 case 1
54
58 e_r = ops{1}.e_r; 55 e_r = ops{1}.e_r;
59 e_l = ops{1}.e_l; 56 e_l = ops{1}.e_l;
60 Dm = ops{1}.Dm;
61 Dp = ops{1}.Dp;
62 Hi = ops{1}.HI; 57 Hi = ops{1}.HI;
63 B = e_r*e_r' - e_l*e_l'; 58 B = e_r*e_r' - e_l*e_l';
64 D2 = @(viscosity) Dm*spdiag(viscosity)*Dp-Hi*(B*spdiag(viscosity)*Dp); 59 if isequal(opSet,@sbp.D1Upwind)
60 Dm = ops{1}.Dm;
61 Dp = ops{1}.Dp;
62 D2 = @(viscosity) Dm*spdiag(viscosity)*Dp-Hi*(B*spdiag(viscosity)*Dp);
63 else
64 D2 = @(viscosity)ops{1}.D2(viscosity);
65 end
65 case 2 66 case 2
67 % TODO:
68 % Currently only implemented for upwind operators.
69 % Remove this part once the time-dependent D2 operator is implemented for other opSets
70 % or if it is decided that it should only be supported for upwind operators.
71 assert(isequal(opSet,@sbp.D1Upwind))
66 e_e = kron(ops{1}.e_r,I{2}); 72 e_e = kron(ops{1}.e_r,I{2});
67 e_w = kron(ops{1}.e_l,I{2}); 73 e_w = kron(ops{1}.e_l,I{2});
68 Dm_x = kron(ops{1}.Dm,I{2}); 74 Dm_x = kron(ops{1}.Dm,I{2});
69 Dp_x = kron(ops{1}.Dp,I{2}); 75 Dp_x = kron(ops{1}.Dp,I{2});
70 H_x = kron(ops{1}.HI,I{2}); 76 H_x = kron(ops{1}.HI,I{2});
80 D2_y = @(viscosity) Dm_y*spdiag(viscosity)*Dp_y-H_y*(B_y*spdiag(viscosity)*Dp_y); 86 D2_y = @(viscosity) Dm_y*spdiag(viscosity)*Dp_y-H_y*(B_y*spdiag(viscosity)*Dp_y);
81 D2 = @(viscosity)D2_x(viscosity) + D2_y(viscosity); 87 D2 = @(viscosity)D2_x(viscosity) + D2_y(viscosity);
82 otherwise 88 otherwise
83 error('3D not yet implemented') 89 error('3D not yet implemented')
84 end 90 end
91 D2 = @(v, viscosity) D2(viscosity)*v;
85 end 92 end