Mercurial > repos > public > sbplib

comparison +rv/constructDiffOps.m @ 1154:3108963cc42c feature/rv

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Improve efficiency of diffOps in Burgers2d, the artificial diffusion operator in rv.constructDiffOps and the RungekuttaExteriorRv time-steppers
author Vidar Stiernström <vidar.stiernstrom@it.uu.se>
date Wed, 06 Mar 2019 09:45:52 +0100
parents 52f59d27b40f
children 8c0e2b50f018 82315fa6adb1
comparison
equal deleted inserted replaced
1153:635386c073b9 1154:3108963cc42c
1 function [D_rv, D_flux, DvDt, solutionPenalties, residualPenalties] = constructDiffOps(scheme, g, order, schemeParams, opSet, BCs) 1 function [D_rv, D_flux, DvDt, solutionPenalties, residualPenalties] = constructDiffOps(scheme, g, order, schemeParams, opSet, BCs)
2 %% DiffOps for solution vector 2 %% DiffOps for solution vector
3 [D, solutionPenalties] = constructTotalFluxDiffOp(scheme, g, order, schemeParams, opSet, BCs); 3 [D, solutionPenalties] = constructFluxDiffOp(scheme, g, order, schemeParams, opSet, BCs);
4 D2 = constructSymmetricD2Operator(g, order, opSet); 4 D2 = constructSymmetricD2Operator(g, order, opSet);
5 D_rv = @(v,viscosity)(D(v) + D2(v, viscosity)); 5 D_rv = @(v,Viscosity)(D(v) + D2(Viscosity)*v);
6 6
7 %% DiffOps for residual viscosity 7 %% DiffOps for residual viscosity
8 [D_flux, residualPenalties] = constructTotalFluxDiffOp(scheme, g, max(order-2,2), schemeParams, opSet, BCs); 8 [D_flux, residualPenalties] = constructFluxDiffOp(scheme, g, max(order-2,2), schemeParams, opSet, BCs);
9 % 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
10 % change the sign. 10 % change the sign.
11 D_flux = @(v) -D_flux(v); 11 D_flux = @(v) -D_flux(v);
12 % DiffOp for time derivative in residual viscosity 12 % DiffOp for time derivative in residual viscosity
13 DvDt = D; 13 DvDt = D;
14 end 14 end
15 15
16 function [D, penalties] = constructTotalFluxDiffOp(scheme, g, order, schemeParams, opSet, BCs) 16 function [D, penalties] = constructFluxDiffOp(scheme, g, order, schemeParams, opSet, BCs)
17 diffOp = scheme(g, order, schemeParams{:}, opSet); 17 diffOp = scheme(g, order, schemeParams{:}, opSet);
18 [D, penalties] = addClosuresToDiffOp(diffOp, BCs); 18 [D, penalties] = addClosuresToDiffOp(diffOp, BCs);
19 end 19 end
20 20
21 function [D, penalties] = addClosuresToDiffOp(diffOp, BCs) 21 function [D, penalties] = addClosuresToDiffOp(diffOp, BCs)
46 lim = {g.x{i}(1), g.x{i}(end)}; 46 lim = {g.x{i}(1), g.x{i}(end)};
47 ops{i} = opSet(m(i), lim, order); 47 ops{i} = opSet(m(i), lim, order);
48 I{i} = speye(m(i)); 48 I{i} = speye(m(i));
49 end 49 end
50 50
51 % TBD: How is this generalized to a loop over dimensions or similar?
52 switch g.D() 51 switch g.D()
53 case 1 52 case 1
54
55 e_r = ops{1}.e_r; 53 e_r = ops{1}.e_r;
56 e_l = ops{1}.e_l; 54 e_l = ops{1}.e_l;
57 Hi = ops{1}.HI; 55 Hi = ops{1}.HI;
58 B = e_r*e_r' - e_l*e_l'; 56 B = e_r*e_r' - e_l*e_l';
59 if isequal(opSet,@sbp.D1Upwind) 57 if isequal(opSet,@sbp.D1Upwind)
60 Dm = ops{1}.Dm; 58 Dm = ops{1}.Dm;
61 Dp = ops{1}.Dp; 59 Dp = ops{1}.Dp;
62 D2 = @(viscosity) Dm*spdiag(viscosity)*Dp-Hi*(B*spdiag(viscosity)*Dp); 60 M = Dm - Hi*B;
61 D2 = @(Viscosity) M*Viscosity*Dp;
63 else 62 else
64 D2 = @(viscosity)ops{1}.D2(viscosity); 63 % TODO: Fix closure for D2Variable
64 D2 = @(Viscosity)ops{1}.D2(Viscosity);
65 end 65 end
66 case 2 66 case 2
67 % TODO: 67 % TODO:
68 % Currently only implemented for upwind operators. 68 % Currently only implemented for upwind operators.
69 % Remove this part once the time-dependent D2 operator is implemented for other opSets 69 % Remove this part once the time-dependent D2 operator is implemented for other opSets
73 e_w = kron(ops{1}.e_l,I{2}); 73 e_w = kron(ops{1}.e_l,I{2});
74 Dm_x = kron(ops{1}.Dm,I{2}); 74 Dm_x = kron(ops{1}.Dm,I{2});
75 Dp_x = kron(ops{1}.Dp,I{2}); 75 Dp_x = kron(ops{1}.Dp,I{2});
76 H_x = kron(ops{1}.HI,I{2}); 76 H_x = kron(ops{1}.HI,I{2});
77 B_x = e_e*e_e' - e_w*e_w'; 77 B_x = e_e*e_e' - e_w*e_w';
78 D2_x = @(viscosity) Dm_x*spdiag(viscosity)*Dp_x-H_x*(B_x*spdiag(viscosity)*Dp_x); 78 M_x = Dm_x-H_x*B_x;
79 D2_x = @(Viscosity) M_x*Viscosity*Dp_x;
79 80
80 e_n = kron(I{1},ops{2}.e_r); 81 e_n = kron(I{1},ops{2}.e_r);
81 e_s = kron(I{1},ops{2}.e_l); 82 e_s = kron(I{1},ops{2}.e_l);
82 Dm_y = kron(I{1},ops{2}.Dm); 83 Dm_y = kron(I{1},ops{2}.Dm);
83 Dp_y = kron(I{1},ops{2}.Dp); 84 Dp_y = kron(I{1},ops{2}.Dp);
84 H_y = kron(I{1},ops{2}.HI); 85 H_y = kron(I{1},ops{2}.HI);
85 B_y = e_n*e_n' - e_s*e_s'; 86 B_y = e_n*e_n' - e_s*e_s';
86 D2_y = @(viscosity) Dm_y*spdiag(viscosity)*Dp_y-H_y*(B_y*spdiag(viscosity)*Dp_y); 87 M_y = Dm_y-H_y*B_y;
87 D2 = @(viscosity)D2_x(viscosity) + D2_y(viscosity); 88 D2_y = @(Viscosity) M_y*Viscosity*Dp_y;
89 D2 = @(Viscosity) D2_x(Viscosity) + D2_y(Viscosity);
88 otherwise 90 otherwise
89 error('3D not yet implemented') 91 error('3D not yet implemented')
90 end 92 end
91 D2 = @(v, viscosity) D2(viscosity)*v;
92 end 93 end