Mercurial > repos > public > sbplib
changeset 1221:0c906f7ab8bf rv_diffOp_test
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 | 010bb2677230 |
| children | 55463e3c1e4a |
| files | +rv/constructDiffOps.m +scheme/Burgers2d.m |
| diffstat | 2 files changed, 50 insertions(+), 29 deletions(-) [+] |
line wrap: on
line diff
--- a/+rv/constructDiffOps.m Tue Mar 05 10:53:34 2019 +0100 +++ b/+rv/constructDiffOps.m Tue Mar 05 10:56:29 2019 +0100 @@ -2,43 +2,62 @@ %% DiffOps for solution vector [D, solutionPenalties] = constructTotalFluxDiffOp(scheme, g, order, schemeParams, opSet, BCs); D2 = constructSymmetricD2Operator(g, order, opSet); - D_rv = @(v,viscosity)(D(v) + D2(v, viscosity)); + + if isa(D, 'function_handle') + D_rv = @(v,viscosity)(D(v) + D2(viscosity))*v; + [D_flux, residualPenalties] = constructTotalFluxDiffOp(scheme, g, max(order-2,2), schemeParams, opSet, BCs); + D_flux = @(v) -D_flux(v)*v; + DvDt = @(v)D(v)*v; + else + D_rv = @(v,viscosity)(D + D2(viscosity))*v; + [D_flux, residualPenalties] = constructTotalFluxDiffOp(scheme, g, max(order-2,2), schemeParams, opSet, BCs); + D_flux = @(v) -D_flux*v; + DvDt = @(v)D*v; + end %% DiffOps for residual viscosity - [D_flux, residualPenalties] = constructTotalFluxDiffOp(scheme, g, max(order-2,2), schemeParams, opSet, BCs); + % DiffOp for flux in residual viscosity. Due to sign conventions of the implemnted schemes, we need to % change the sign. - D_flux = @(v) -D_flux(v); + % DiffOp for time derivative in residual viscosity - DvDt = D; + end function [D, penalties] = constructTotalFluxDiffOp(scheme, g, order, schemeParams, opSet, BCs) diffOp = scheme(g, order, schemeParams{:}, opSet); - [D, penalties] = addClosuresToDiffOp(diffOp, BCs); + if isa(diffOp.D, 'function_handle') + [D, penalties] = addClosuresToDiffOpFunction(diffOp, BCs); + else + [D, penalties] = addClosuresToDiffOp(diffOp, BCs); + end end function [D, penalties] = addClosuresToDiffOp(diffOp, BCs) - if ~isa(diffOp.D, 'function_handle') - D = @(v) diffOp.D*v; - else - D = diffOp.D; - end penalties = cell(size(BCs)); + D = diffOp.D; for i = 1:size(BCs,1) for j = 1:size(BCs,2) [closure, penalties{i,j}] = diffOp.boundary_condition(BCs{i,j}.boundary, BCs{i,j}.type); - if ~isa(closure, 'function_handle') - closure = @(v) closure*v; - end + D = D + closure; + end + end + +end + +function [D, penalties] = addClosuresToDiffOpFunction(diffOp, BCs) + penalties = cell(size(BCs)); + D = diffOp.D; + for i = 1:size(BCs,1) + for j = 1:size(BCs,2) + [closure, penalties{i,j}] = diffOp.boundary_condition(BCs{i,j}.boundary, BCs{i,j}.type); D = @(v) D(v) + closure(v); end end + end function D2 = constructSymmetricD2Operator(g, order, opSet) - - m = g.size(); ops = cell(g.D(),1); I = cell(g.D(),1); @@ -48,10 +67,8 @@ I{i} = speye(m(i)); end - % TBD: How is this generalized to a loop over dimensions or similar? switch g.D() case 1 - e_r = ops{1}.e_r; e_l = ops{1}.e_l; Hi = ops{1}.HI; @@ -59,7 +76,8 @@ if isequal(opSet,@sbp.D1Upwind) Dm = ops{1}.Dm; Dp = ops{1}.Dp; - D2 = @(viscosity) Dm*spdiag(viscosity)*Dp-Hi*(B*spdiag(viscosity)*Dp); + HiB = Hi*B; + D2 = @(viscosity) Dm*spdiag(viscosity)*Dp-HiB*spdiag(viscosity)*Dp; else D2 = @(viscosity)ops{1}.D2(viscosity); end @@ -75,7 +93,8 @@ 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'; - D2_x = @(viscosity) Dm_x*spdiag(viscosity)*Dp_x-H_x*(B_x*spdiag(viscosity)*Dp_x); + H_xB = H_x*B_x; + D2_x = @(viscosity) (Dm_x*spdiag(viscosity)-H_xB*spdiag(viscosity))*Dp_x; e_n = kron(I{1},ops{2}.e_r); e_s = kron(I{1},ops{2}.e_l); @@ -83,10 +102,11 @@ 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'; - D2_y = @(viscosity) Dm_y*spdiag(viscosity)*Dp_y-H_y*(B_y*spdiag(viscosity)*Dp_y); + H_yB = H_y*B_y; + D2_y = @(viscosity) (Dm_y*spdiag(viscosity)-H_yB*spdiag(viscosity))*Dp_y; D2 = @(viscosity)D2_x(viscosity) + D2_y(viscosity); otherwise error('3D not yet implemented') end - D2 = @(v, viscosity) D2(viscosity)*v; + D2 = @(viscosity) D2(viscosity); end \ No newline at end of file
--- a/+scheme/Burgers2d.m Tue Mar 05 10:53:34 2019 +0100 +++ b/+scheme/Burgers2d.m Tue Mar 05 10:56:29 2019 +0100 @@ -55,17 +55,17 @@ switch length(fluxSplitting) case 1 DissOp = DissOpx + DissOpy; - obj.D = @(v) -(1/3*D1*v.*v + (1/3*spdiag(v)*D1 + fluxSplitting{1}(v)*DissOp)*v); + obj.D = @(v) -(1/3*D1*(v.*v) + (1/3*spdiag(v)*D1 + fluxSplitting{1}(v)*DissOp)*v); case 2 - obj.D = @(v) -(1/3*D1*v.*v + (1/3*spdiag(v)*D1 + fluxSplitting{1}(v)*DissOpx + fluxSplitting{2}(v)*DissOpy)*v); + obj.D = @(v) -(1/3*D1*(v.*v) + (1/3*spdiag(v)*D1 + fluxSplitting{1}(v)*DissOpx + fluxSplitting{2}(v)*DissOpy)*v); end case 'conservative' switch length(fluxSplitting) case 1 DissOp = DissOpx + DissOpy; - obj.D = @(v) -(1/2*D1*v.*v + fluxSplitting{1}(v)*DissOp*v); + obj.D = @(v) -(1/2*D1*spdiag(v) + fluxSplitting{1}(v)*DissOp); case 2 - obj.D = @(v) -(1/2*D1*v.*v + (fluxSplitting{1}(v)*DissOpx + fluxSplitting{2}(v)*DissOpy)*v); + obj.D = @(v) -(1/2*D1*spdiag(v) + (fluxSplitting{1}(v)*DissOpx + fluxSplitting{2}(v)*DissOpy)); end end @@ -75,9 +75,9 @@ D1 = Dx + Dy; switch pde_form case 'skew-symmetric' - obj.D = @(v) -(1/3*D1*v.*v + 1/3*spdiag(v)*D1*v); + obj.D = @(v) -(1/3*D1*(v.*v) + 1/3*spdiag(v)*D1*v); case 'conservative' - obj.D = @(v) -1/2*D1*v.*v; + obj.D = @(v) -1/2*D1*(v.*v); end end @@ -103,8 +103,9 @@ case {'D', 'd', 'dirichlet', 'Dirichlet'} magnitude = 2/3; - tau = @(v) s*magnitude*obj.Hi*e*H_b*spdiag((v(index)-s*abs(v(index)))/2); - closure = @(v) tau(v)*v(index); + Tau = s*magnitude*obj.Hi*e*H_b; + tau = @(v) Tau*spdiag((v(index)-s*abs(v(index)))/2); + closure = @(v) tau(v)*e'; penalty = @(v) -tau(v);