Mercurial > repos > public > sbplib
view +rv/+diffops/constructSymmetricD2.m @ 1225:68ee061639a1 feature/rv
Make sure matrices are sparse.
author | Vidar Stiernström <vidar.stiernstrom@it.uu.se> |
---|---|
date | Wed, 06 Nov 2019 14:52:07 +0100 |
parents | 433c89bf19e0 |
children |
line wrap: on
line source
function D2 = constructSymmetricD2(g, order, opSet) m = g.size(); ops = cell(g.D(),1); I = cell(g.D(),1); for i = 1:g.D() lim = {g.x{i}(1), g.x{i}(end)}; ops{i} = opSet(m(i), lim, order); I{i} = speye(m(i)); end switch g.D() case 1 e_r = ops{1}.e_r; e_l = ops{1}.e_l; Hi = ops{1}.HI; if isequal(opSet,@sbp.D1Upwind) B = e_r*e_r' - e_l*e_l'; Dm = sparse(ops{1}.Dm); Dp = sparse(ops{1}.Dp); M = sparse(Dm - Hi*B); D2 = @(Viscosity) M*Viscosity*Dp; else % TODO: Fix Viscosity not being vector d1_r = ops{1}.d1_r'; d1_l = ops{1}.d1_l'; D2 = @(Viscosity)ops{1}.D2(diag(Viscosity)) + Hi*(Viscosity(1,1)*e_l*d1_l - e_r*Viscosity(end,end)*d1_r); end case 2 % TODO: % Currently only implemented for upwind operators. % Remove this part once the time-dependent D2 operator is implemented for other opSets % or if it is decided that it should only be supported for upwind operators. assert(isequal(opSet,@sbp.D1Upwind)) e_e = kron(ops{1}.e_r,I{2}); e_w = kron(ops{1}.e_l,I{2}); Dm_x = kron(ops{1}.Dm,I{2}); Dp_x = kron(ops{1}.Dp,I{2}); H_x = kron(ops{1}.HI,I{2}); B_x = e_e*e_e' - e_w*e_w'; M_x = Dm_x-H_x*B_x; D2_x = @(Viscosity) M_x*Viscosity*Dp_x; e_n = kron(I{1},ops{2}.e_r); e_s = kron(I{1},ops{2}.e_l); Dm_y = kron(I{1},ops{2}.Dm); Dp_y = kron(I{1},ops{2}.Dp); H_y = kron(I{1},ops{2}.HI); B_y = e_n*e_n' - e_s*e_s'; M_y = Dm_y-H_y*B_y; D2_y = @(Viscosity) M_y*Viscosity*Dp_y; D2 = @(Viscosity) D2_x(Viscosity) + D2_y(Viscosity); otherwise error('3D not yet implemented') end end