Mercurial > repos > public > sbplib
comparison +rv/+diffops/constructSymmetricD2.m @ 1163:65a577db5ca0 feature/rv
Move all functions in constructDiffOps to subpackage and refactor
| author | Vidar Stiernström <vidar.stiernstrom@it.uu.se> |
|---|---|
| date | Thu, 27 Jun 2019 11:04:34 +0200 |
| parents | +rv/constructDiffOps.m@856bd6291d17 |
| children | 433c89bf19e0 |
comparison
equal
deleted
inserted
replaced
| 1162:0ec06ca3fc36 | 1163:65a577db5ca0 |
|---|---|
| 1 function D2 = constructSymmetricD2(g, order, opSet) | |
| 2 m = g.size(); | |
| 3 ops = cell(g.D(),1); | |
| 4 I = cell(g.D(),1); | |
| 5 for i = 1:g.D() | |
| 6 lim = {g.x{i}(1), g.x{i}(end)}; | |
| 7 ops{i} = opSet(m(i), lim, order); | |
| 8 I{i} = speye(m(i)); | |
| 9 end | |
| 10 | |
| 11 switch g.D() | |
| 12 case 1 | |
| 13 e_r = ops{1}.e_r; | |
| 14 e_l = ops{1}.e_l; | |
| 15 Hi = ops{1}.HI; | |
| 16 B = e_r*e_r' - e_l*e_l'; | |
| 17 if isequal(opSet,@sbp.D1Upwind) | |
| 18 Dm = ops{1}.Dm; | |
| 19 Dp = ops{1}.Dp; | |
| 20 M = Dm - Hi*B; | |
| 21 D2 = @(Viscosity) M*Viscosity*Dp; | |
| 22 else | |
| 23 % TODO: Fix closure for D2Variable | |
| 24 % TODO: Fix Viscosity not being vector | |
| 25 D2 = @(Viscosity)ops{1}.D2(diag(Viscosity)); | |
| 26 end | |
| 27 case 2 | |
| 28 % TODO: | |
| 29 % Currently only implemented for upwind operators. | |
| 30 % Remove this part once the time-dependent D2 operator is implemented for other opSets | |
| 31 % or if it is decided that it should only be supported for upwind operators. | |
| 32 assert(isequal(opSet,@sbp.D1Upwind)) | |
| 33 e_e = kron(ops{1}.e_r,I{2}); | |
| 34 e_w = kron(ops{1}.e_l,I{2}); | |
| 35 Dm_x = kron(ops{1}.Dm,I{2}); | |
| 36 Dp_x = kron(ops{1}.Dp,I{2}); | |
| 37 H_x = kron(ops{1}.HI,I{2}); | |
| 38 B_x = e_e*e_e' - e_w*e_w'; | |
| 39 M_x = Dm_x-H_x*B_x; | |
| 40 D2_x = @(Viscosity) M_x*Viscosity*Dp_x; | |
| 41 | |
| 42 e_n = kron(I{1},ops{2}.e_r); | |
| 43 e_s = kron(I{1},ops{2}.e_l); | |
| 44 Dm_y = kron(I{1},ops{2}.Dm); | |
| 45 Dp_y = kron(I{1},ops{2}.Dp); | |
| 46 H_y = kron(I{1},ops{2}.HI); | |
| 47 B_y = e_n*e_n' - e_s*e_s'; | |
| 48 M_y = Dm_y-H_y*B_y; | |
| 49 D2_y = @(Viscosity) M_y*Viscosity*Dp_y; | |
| 50 D2 = @(Viscosity) D2_x(Viscosity) + D2_y(Viscosity); | |
| 51 otherwise | |
| 52 error('3D not yet implemented') | |
| 53 end | |
| 54 end |