Mercurial > repos > public > sbplib

comparison +rv/constructDiffOps.m @ 1157:82315fa6adb1 feature/rv

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Add residual order as an input argument to constructDiffOps
author Vidar Stiernström <vidar.stiernstrom@it.uu.se>
date Tue, 25 Jun 2019 13:24:01 +0200
parents 3108963cc42c
children 76e3bb7836cf
comparison
equal deleted inserted replaced
1155:336ee37a0617 1157:82315fa6adb1
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);
12 % DiffOp for time derivative in residual viscosity 20 % DiffOp for time derivative in residual viscosity
13 DvDt = D; 21 DvDt = D;
14 end 22 end
35 end 43 end
36 end 44 end
37 end 45 end
38 46
39 function D2 = constructSymmetricD2Operator(g, order, opSet) 47 function D2 = constructSymmetricD2Operator(g, order, opSet)
40
41
42 m = g.size(); 48 m = g.size();
43 ops = cell(g.D(),1); 49 ops = cell(g.D(),1);
44 I = cell(g.D(),1); 50 I = cell(g.D(),1);
45 for i = 1:g.D() 51 for i = 1:g.D()
46 lim = {g.x{i}(1), g.x{i}(end)}; 52 lim = {g.x{i}(1), g.x{i}(end)};
59 Dp = ops{1}.Dp; 65 Dp = ops{1}.Dp;
60 M = Dm - Hi*B; 66 M = Dm - Hi*B;
61 D2 = @(Viscosity) M*Viscosity*Dp; 67 D2 = @(Viscosity) M*Viscosity*Dp;
62 else 68 else
63 % TODO: Fix closure for D2Variable 69 % TODO: Fix closure for D2Variable
64 D2 = @(Viscosity)ops{1}.D2(Viscosity); 70 % TODO: Fix Viscosity not being vector
71 D2 = @(Viscosity)ops{1}.D2(diag(Viscosity));
65 end 72 end
66 case 2 73 case 2
67 % TODO: 74 % TODO:
68 % Currently only implemented for upwind operators. 75 % Currently only implemented for upwind operators.
69 % Remove this part once the time-dependent D2 operator is implemented for other opSets 76 % Remove this part once the time-dependent D2 operator is implemented for other opSets