Mercurial > repos > public > sbplib

comparison +rv/constructDiffOps.m @ 1221:0c906f7ab8bf rv_diffOp_test

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Attempt to reduce unnessecary operations when using the RV diffOps. Does not seem to improve efficiency at this stage.
author Vidar Stiernström <vidar.stiernstrom@it.uu.se>
date Tue, 05 Mar 2019 10:56:29 +0100
parents 52f59d27b40f
children
comparison
equal deleted inserted replaced
1152:010bb2677230 1221:0c906f7ab8bf
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] = constructTotalFluxDiffOp(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
6 if isa(D, 'function_handle')
7 D_rv = @(v,viscosity)(D(v) + D2(viscosity))*v;
8 [D_flux, residualPenalties] = constructTotalFluxDiffOp(scheme, g, max(order-2,2), schemeParams, opSet, BCs);
9 D_flux = @(v) -D_flux(v)*v;
10 DvDt = @(v)D(v)*v;
11 else
12 D_rv = @(v,viscosity)(D + D2(viscosity))*v;
13 [D_flux, residualPenalties] = constructTotalFluxDiffOp(scheme, g, max(order-2,2), schemeParams, opSet, BCs);
14 D_flux = @(v) -D_flux*v;
15 DvDt = @(v)D*v;
16 end
6 17
7 %% DiffOps for residual viscosity 18 %% DiffOps for residual viscosity
8 [D_flux, residualPenalties] = constructTotalFluxDiffOp(scheme, g, max(order-2,2), schemeParams, opSet, BCs); 19
9 % DiffOp for flux in residual viscosity. Due to sign conventions of the implemnted schemes, we need to 20 % DiffOp for flux in residual viscosity. Due to sign conventions of the implemnted schemes, we need to
10 % change the sign. 21 % change the sign.
11 D_flux = @(v) -D_flux(v); 22
12 % DiffOp for time derivative in residual viscosity 23 % DiffOp for time derivative in residual viscosity
13 DvDt = D; 24
14 end 25 end
15 26
16 function [D, penalties] = constructTotalFluxDiffOp(scheme, g, order, schemeParams, opSet, BCs) 27 function [D, penalties] = constructTotalFluxDiffOp(scheme, g, order, schemeParams, opSet, BCs)
17 diffOp = scheme(g, order, schemeParams{:}, opSet); 28 diffOp = scheme(g, order, schemeParams{:}, opSet);
18 [D, penalties] = addClosuresToDiffOp(diffOp, BCs); 29 if isa(diffOp.D, 'function_handle')
30 [D, penalties] = addClosuresToDiffOpFunction(diffOp, BCs);
31 else
32 [D, penalties] = addClosuresToDiffOp(diffOp, BCs);
33 end
19 end 34 end
20 35
21 function [D, penalties] = addClosuresToDiffOp(diffOp, BCs) 36 function [D, penalties] = addClosuresToDiffOp(diffOp, BCs)
22 if ~isa(diffOp.D, 'function_handle')
23 D = @(v) diffOp.D*v;
24 else
25 D = diffOp.D;
26 end
27 penalties = cell(size(BCs)); 37 penalties = cell(size(BCs));
38 D = diffOp.D;
28 for i = 1:size(BCs,1) 39 for i = 1:size(BCs,1)
29 for j = 1:size(BCs,2) 40 for j = 1:size(BCs,2)
30 [closure, penalties{i,j}] = diffOp.boundary_condition(BCs{i,j}.boundary, BCs{i,j}.type); 41 [closure, penalties{i,j}] = diffOp.boundary_condition(BCs{i,j}.boundary, BCs{i,j}.type);
31 if ~isa(closure, 'function_handle') 42 D = D + closure;
32 closure = @(v) closure*v; 43 end
33 end 44 end
45
46 end
47
48 function [D, penalties] = addClosuresToDiffOpFunction(diffOp, BCs)
49 penalties = cell(size(BCs));
50 D = diffOp.D;
51 for i = 1:size(BCs,1)
52 for j = 1:size(BCs,2)
53 [closure, penalties{i,j}] = diffOp.boundary_condition(BCs{i,j}.boundary, BCs{i,j}.type);
34 D = @(v) D(v) + closure(v); 54 D = @(v) D(v) + closure(v);
35 end 55 end
36 end 56 end
57
37 end 58 end
38 59
39 function D2 = constructSymmetricD2Operator(g, order, opSet) 60 function D2 = constructSymmetricD2Operator(g, order, opSet)
40
41
42 m = g.size(); 61 m = g.size();
43 ops = cell(g.D(),1); 62 ops = cell(g.D(),1);
44 I = cell(g.D(),1); 63 I = cell(g.D(),1);
45 for i = 1:g.D() 64 for i = 1:g.D()
46 lim = {g.x{i}(1), g.x{i}(end)}; 65 lim = {g.x{i}(1), g.x{i}(end)};
47 ops{i} = opSet(m(i), lim, order); 66 ops{i} = opSet(m(i), lim, order);
48 I{i} = speye(m(i)); 67 I{i} = speye(m(i));
49 end 68 end
50 69
51 % TBD: How is this generalized to a loop over dimensions or similar?
52 switch g.D() 70 switch g.D()
53 case 1 71 case 1
54
55 e_r = ops{1}.e_r; 72 e_r = ops{1}.e_r;
56 e_l = ops{1}.e_l; 73 e_l = ops{1}.e_l;
57 Hi = ops{1}.HI; 74 Hi = ops{1}.HI;
58 B = e_r*e_r' - e_l*e_l'; 75 B = e_r*e_r' - e_l*e_l';
59 if isequal(opSet,@sbp.D1Upwind) 76 if isequal(opSet,@sbp.D1Upwind)
60 Dm = ops{1}.Dm; 77 Dm = ops{1}.Dm;
61 Dp = ops{1}.Dp; 78 Dp = ops{1}.Dp;
62 D2 = @(viscosity) Dm*spdiag(viscosity)*Dp-Hi*(B*spdiag(viscosity)*Dp); 79 HiB = Hi*B;
80 D2 = @(viscosity) Dm*spdiag(viscosity)*Dp-HiB*spdiag(viscosity)*Dp;
63 else 81 else
64 D2 = @(viscosity)ops{1}.D2(viscosity); 82 D2 = @(viscosity)ops{1}.D2(viscosity);
65 end 83 end
66 case 2 84 case 2
67 % TODO: 85 % TODO:
73 e_w = kron(ops{1}.e_l,I{2}); 91 e_w = kron(ops{1}.e_l,I{2});
74 Dm_x = kron(ops{1}.Dm,I{2}); 92 Dm_x = kron(ops{1}.Dm,I{2});
75 Dp_x = kron(ops{1}.Dp,I{2}); 93 Dp_x = kron(ops{1}.Dp,I{2});
76 H_x = kron(ops{1}.HI,I{2}); 94 H_x = kron(ops{1}.HI,I{2});
77 B_x = e_e*e_e' - e_w*e_w'; 95 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); 96 H_xB = H_x*B_x;
97 D2_x = @(viscosity) (Dm_x*spdiag(viscosity)-H_xB*spdiag(viscosity))*Dp_x;
79 98
80 e_n = kron(I{1},ops{2}.e_r); 99 e_n = kron(I{1},ops{2}.e_r);
81 e_s = kron(I{1},ops{2}.e_l); 100 e_s = kron(I{1},ops{2}.e_l);
82 Dm_y = kron(I{1},ops{2}.Dm); 101 Dm_y = kron(I{1},ops{2}.Dm);
83 Dp_y = kron(I{1},ops{2}.Dp); 102 Dp_y = kron(I{1},ops{2}.Dp);
84 H_y = kron(I{1},ops{2}.HI); 103 H_y = kron(I{1},ops{2}.HI);
85 B_y = e_n*e_n' - e_s*e_s'; 104 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); 105 H_yB = H_y*B_y;
106 D2_y = @(viscosity) (Dm_y*spdiag(viscosity)-H_yB*spdiag(viscosity))*Dp_y;
87 D2 = @(viscosity)D2_x(viscosity) + D2_y(viscosity); 107 D2 = @(viscosity)D2_x(viscosity) + D2_y(viscosity);
88 otherwise 108 otherwise
89 error('3D not yet implemented') 109 error('3D not yet implemented')
90 end 110 end
91 D2 = @(v, viscosity) D2(viscosity)*v; 111 D2 = @(viscosity) D2(viscosity);
92 end 112 end