Mercurial > repos > public > sbplib

comparison +scheme/Burgers2d.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 1fe48cbd379a
children
comparison
equal deleted inserted replaced
1152:010bb2677230 1221:0c906f7ab8bf
53 switch pde_form 53 switch pde_form
54 case 'skew-symmetric' 54 case 'skew-symmetric'
55 switch length(fluxSplitting) 55 switch length(fluxSplitting)
56 case 1 56 case 1
57 DissOp = DissOpx + DissOpy; 57 DissOp = DissOpx + DissOpy;
58 obj.D = @(v) -(1/3*D1*v.*v + (1/3*spdiag(v)*D1 + fluxSplitting{1}(v)*DissOp)*v); 58 obj.D = @(v) -(1/3*D1*(v.*v) + (1/3*spdiag(v)*D1 + fluxSplitting{1}(v)*DissOp)*v);
59 case 2 59 case 2
60 obj.D = @(v) -(1/3*D1*v.*v + (1/3*spdiag(v)*D1 + fluxSplitting{1}(v)*DissOpx + fluxSplitting{2}(v)*DissOpy)*v); 60 obj.D = @(v) -(1/3*D1*(v.*v) + (1/3*spdiag(v)*D1 + fluxSplitting{1}(v)*DissOpx + fluxSplitting{2}(v)*DissOpy)*v);
61 end 61 end
62 case 'conservative' 62 case 'conservative'
63 switch length(fluxSplitting) 63 switch length(fluxSplitting)
64 case 1 64 case 1
65 DissOp = DissOpx + DissOpy; 65 DissOp = DissOpx + DissOpy;
66 obj.D = @(v) -(1/2*D1*v.*v + fluxSplitting{1}(v)*DissOp*v); 66 obj.D = @(v) -(1/2*D1*spdiag(v) + fluxSplitting{1}(v)*DissOp);
67 case 2 67 case 2
68 obj.D = @(v) -(1/2*D1*v.*v + (fluxSplitting{1}(v)*DissOpx + fluxSplitting{2}(v)*DissOpy)*v); 68 obj.D = @(v) -(1/2*D1*spdiag(v) + (fluxSplitting{1}(v)*DissOpx + fluxSplitting{2}(v)*DissOpy));
69 end 69 end
70 70
71 end 71 end
72 else 72 else
73 Dx = kron(ops_x.D1,Iy); 73 Dx = kron(ops_x.D1,Iy);
74 Dy = kron(Ix,ops_y.D1); 74 Dy = kron(Ix,ops_y.D1);
75 D1 = Dx + Dy; 75 D1 = Dx + Dy;
76 switch pde_form 76 switch pde_form
77 case 'skew-symmetric' 77 case 'skew-symmetric'
78 obj.D = @(v) -(1/3*D1*v.*v + 1/3*spdiag(v)*D1*v); 78 obj.D = @(v) -(1/3*D1*(v.*v) + 1/3*spdiag(v)*D1*v);
79 case 'conservative' 79 case 'conservative'
80 obj.D = @(v) -1/2*D1*v.*v; 80 obj.D = @(v) -1/2*D1*(v.*v);
81 end 81 end
82 end 82 end
83 83
84 % Boundary operators 84 % Boundary operators
85 obj.e_w = kr(ops_x.e_l, Iy); 85 obj.e_w = kr(ops_x.e_l, Iy);
101 % Stable dirchlet-like boundary conditions (u+-abs(u))*u/3 101 % Stable dirchlet-like boundary conditions (u+-abs(u))*u/3
102 % with +- at left/right boundaries in each coordinate direction 102 % with +- at left/right boundaries in each coordinate direction
103 case {'D', 'd', 'dirichlet', 'Dirichlet'} 103 case {'D', 'd', 'dirichlet', 'Dirichlet'}
104 104
105 magnitude = 2/3; 105 magnitude = 2/3;
106 tau = @(v) s*magnitude*obj.Hi*e*H_b*spdiag((v(index)-s*abs(v(index)))/2); 106 Tau = s*magnitude*obj.Hi*e*H_b;
107 closure = @(v) tau(v)*v(index); 107 tau = @(v) Tau*spdiag((v(index)-s*abs(v(index)))/2);
108 closure = @(v) tau(v)*e';
108 penalty = @(v) -tau(v); 109 penalty = @(v) -tau(v);
109 110
110 111
111 % tau = s*e*H_b; 112 % tau = s*e*H_b;
112 % closure = @(v) (obj.Hi*tau)*(((v(index)-s*abs(v(index)))/3).*v(index)); 113 % closure = @(v) (obj.Hi*tau)*(((v(index)-s*abs(v(index)))/3).*v(index));