Mercurial > repos > public > sbplib

comparison +rv/constructDiffOps.m @ 1156:8c0e2b50f018 feature/rv

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Attempt to fix oscilations when not using bdfs by adding lower-order 2nd derivative correction term to the Residual
author Vidar Stiernström <vidar.stiernstrom@it.uu.se>
date Tue, 25 Jun 2019 13:06:42 +0200
parents 3108963cc42c
children
comparison
equal deleted inserted replaced
1155:336ee37a0617 1156:8c0e2b50f018
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, schemeOrder, residualOrder, schemeParams, opSet, BCs)
2 %% DiffOps for solution vector 2 %% DiffOps for solution vector
3 [D, solutionPenalties] = constructFluxDiffOp(scheme, g, order, schemeParams, opSet, BCs); 3 [D, solutionPenalties] = constructFluxDiffOp(scheme, g, schemeOrder, schemeParams, opSet, BCs);
4 D2 = constructSymmetricD2Operator(g, order, opSet); 4 D2 = constructSymmetricD2Operator(g, schemeOrder, opSet);
5 D_rv = @(v,Viscosity)(D(v) + D2(Viscosity)*v); 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] = constructFluxDiffOp(scheme, g, max(order-2,2), schemeParams, opSet, BCs); 8 [D_flux, residualPenalties] = constructFluxDiffOp(scheme, g, residualOrder, schemeParams, opSet, BCs);
9 % DiffOp for flux in residual viscosity. Due to sign conventions of the implemnted schemes, we need to 9 % TODO: Construct D_flux without closures when using bdfs.
10 % diffOp = scheme(g, residualOrder, schemeParams{:}, opSet);
11 % if ~isa(diffOp.D, 'function_handle')
12 % D_flux = @(v) diffOp.D*v;
13 % else
14 % D_flux = diffOp.D;
15 % end
16
17 % DiffOp for flux in residual viscosity. Due to sign conventions of the implemented schemes, we need to
10 % change the sign. 18 % change the sign.
11 D_flux = @(v) -D_flux(v); 19 D_flux = @(v) -D_flux(v);
20 lim = {g.x{1}(1), g.x{1}(end)};
21 o = 2;
22 h = g.scaling();
23 gam = 1;
24 if isequal(opSet,@sbp.D1Upwind)
25 ops = sbp.D1Upwind(g.m(), lim, o);
26 Hi = ops.HI;
27 e_r = ops.e_r;
28 e_l = ops.e_l;
29 B = e_r*e_r' - e_l*e_l';
30 M = ops.Dm - Hi*B;
31 D2 = M*ops.Dp;
32 else
33 ops = sbp.D2Standard(g.m(), lim, o);
34 D2 = ops.D2;
35 end
36 D2_low = gam*h^(residualOrder)*D2;
37 D2_low = @(v) D2_low*v;
12 % DiffOp for time derivative in residual viscosity 38 % DiffOp for time derivative in residual viscosity
13 DvDt = D; 39 DvDt = @(v) D(v)+D2_low(v);
40 %DvDt = @(v) D(v);
14 end 41 end
15 42
16 function [D, penalties] = constructFluxDiffOp(scheme, g, order, schemeParams, opSet, BCs) 43 function [D, penalties] = constructFluxDiffOp(scheme, g, order, schemeParams, opSet, BCs)
17 diffOp = scheme(g, order, schemeParams{:}, opSet); 44 diffOp = scheme(g, order, schemeParams{:}, opSet);
18 [D, penalties] = addClosuresToDiffOp(diffOp, BCs); 45 [D, penalties] = addClosuresToDiffOp(diffOp, BCs);
35 end 62 end
36 end 63 end
37 end 64 end
38 65
39 function D2 = constructSymmetricD2Operator(g, order, opSet) 66 function D2 = constructSymmetricD2Operator(g, order, opSet)
40
41
42 m = g.size(); 67 m = g.size();
43 ops = cell(g.D(),1); 68 ops = cell(g.D(),1);
44 I = cell(g.D(),1); 69 I = cell(g.D(),1);
45 for i = 1:g.D() 70 for i = 1:g.D()
46 lim = {g.x{i}(1), g.x{i}(end)}; 71 lim = {g.x{i}(1), g.x{i}(end)};
59 Dp = ops{1}.Dp; 84 Dp = ops{1}.Dp;
60 M = Dm - Hi*B; 85 M = Dm - Hi*B;
61 D2 = @(Viscosity) M*Viscosity*Dp; 86 D2 = @(Viscosity) M*Viscosity*Dp;
62 else 87 else
63 % TODO: Fix closure for D2Variable 88 % TODO: Fix closure for D2Variable
64 D2 = @(Viscosity)ops{1}.D2(Viscosity); 89 % TODO: Fix Viscosity not being vector
90 D2 = @(Viscosity)ops{1}.D2(diag(Viscosity));
65 end 91 end
66 case 2 92 case 2
67 % TODO: 93 % TODO:
68 % Currently only implemented for upwind operators. 94 % Currently only implemented for upwind operators.
69 % Remove this part once the time-dependent D2 operator is implemented for other opSets 95 % Remove this part once the time-dependent D2 operator is implemented for other opSets