Mercurial > repos > public > sbplib
changeset 1190:79618b58b0a0 feature/rv
Refactor constructDiffops and remove closures from some of the diffops used to compute the residual.
| author | Vidar Stiernström <vidar.stiernstrom@it.uu.se> |
|---|---|
| date | Mon, 08 Jul 2019 15:12:19 +0200 |
| parents | 6cb03209f0a7 |
| children | 881afc40a3d2 |
| files | +rv/+diffops/addClosuresToDiffOp.m +rv/+diffops/constructDiffOps.m +rv/+diffops/constructFluxDiffOp.m +rv/+diffops/constructFluxDiffOpWithClosures.m +rv/+diffops/constructFluxDiffOps.m +rv/+diffops/constructSchemeDiffOp.m +rv/+diffops/constructSchemeDiffOps.m |
| diffstat | 7 files changed, 92 insertions(+), 53 deletions(-) [+] |
line wrap: on
line diff
diff -r 6cb03209f0a7 -r 79618b58b0a0 +rv/+diffops/addClosuresToDiffOp.m --- a/+rv/+diffops/addClosuresToDiffOp.m Mon Jul 08 14:50:46 2019 +0200 +++ b/+rv/+diffops/addClosuresToDiffOp.m Mon Jul 08 15:12:19 2019 +0200 @@ -1,8 +1,7 @@ -function [D, penalties] = addClosuresToDiffOp(diffOp, D, BCs) - penalties = cell(size(BCs)); - 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); +function D = addClosuresToDiffOp(D, closures) + for i = 1:size(closures,1) + for j = 1:size(closures,2) + closure = closures{i,j}; if ~isa(closure, 'function_handle') closure = @(v) closure*v; end
diff -r 6cb03209f0a7 -r 79618b58b0a0 +rv/+diffops/constructDiffOps.m --- a/+rv/+diffops/constructDiffOps.m Mon Jul 08 14:50:46 2019 +0200 +++ b/+rv/+diffops/constructDiffOps.m Mon Jul 08 15:12:19 2019 +0200 @@ -15,84 +15,120 @@ function diffOpsStruct = diffOpsStandard(scheme, g, schemeOrder, residualOrder, schemeParams, opSet, BCs) % DiffOps for stabilized solution vector - [D_scheme, penalties_scheme, D_t] = rv.diffops.constructSchemeDiffOp(scheme, g, schemeOrder, schemeParams, opSet, BCs); + [D_scheme, penalties_scheme, D_t] = rv.diffops.constructSchemeDiffOps(scheme, g, schemeOrder, schemeParams, opSet, BCs); %% DiffOps for residual viscosity - D_flux = rv.diffops.constructFluxDiffOpWithClosures(scheme, g, residualOrder, schemeParams, opSet, BCs); + D_flux = rv.diffops.constructFluxDiffOps(scheme, g, residualOrder, schemeParams, opSet, BCs); D_flux = @(v) -D_flux(v); diffOpsStruct = struct('D_scheme', D_scheme,... 'D_flux', D_flux,... 'D_t', D_t,... - 'penalties_scheme', penalties_scheme); + 'penalties_scheme', {penalties_scheme}); end function diffOpsStruct = diffOpsBdf(scheme, g, schemeOrder, residualOrder, schemeParams, opSet, BCs) %% DiffOps for solution vector - [D_scheme, penalties_scheme, ~] = rv.diffops.constructSchemeDiffOp(scheme, g, schemeOrder, schemeParams, opSet, BCs); + [D_scheme, penalties_scheme] = rv.diffops.constructSchemeDiffOps(scheme, g, schemeOrder, schemeParams, opSet, BCs); %% DiffOps for residual viscosity - [D_flux, ~] = rv.diffops.constructFluxDiffOp(scheme, g, residualOrder, schemeParams, opSet, BCs); + D_flux = rv.diffops.constructFluxDiffOps(scheme, g, residualOrder, schemeParams, opSet, BCs); D_flux = @(v) -D_flux(v); diffOpsStruct = struct('D_scheme', D_scheme,... 'D_flux', D_flux,... - 'penalties_scheme', penalties_scheme); + 'penalties_scheme', {penalties_scheme}); end +% TODO: Remove? function diffOpsStruct = diffOpsMultiStage(scheme, g, schemeOrder, residualOrder, schemeParams, opSet, BCs) % DiffOps for stabilized solution vector - [D_scheme, penalties_scheme] = rv.diffops.constructSchemeDiffOp(scheme, g, schemeOrder, schemeParams, opSet, BCs); + [D_scheme, penalties_scheme] = rv.diffops.constructSchemeDiffOps(scheme, g, schemeOrder, schemeParams, opSet, BCs); % DiffOp for unstabilized solution vector - D_unstable = rv.diffops.constructFluxDiffOpWithClosures(scheme, g, schemeOrder, schemeParams, opSet, BCs); + [D_unstable, closures] = rv.diffops.constructFluxDiffOps(scheme, g, schemeOrder, schemeParams, opSet, BCs); + D_unstable = rv.diffops.addClosuresToDiffOp(D_unstable, closures); %% DiffOps for residual viscosity - [D_t, penalties_res] = rv.diffops.constructFluxDiffOpWithClosures(scheme, g, residualOrder, schemeParams, opSet, BCs); + D_t = rv.diffops.constructFluxDiffOps(scheme, g, residualOrder, schemeParams, opSet, BCs); D_flux = @(v) -D_t(v); + % TODO: Use residual penalties as well here? diffOpsStruct = struct('D_scheme', D_scheme,... 'D_unstable', D_unstable,... 'D_flux', D_flux,... 'D_t', D_t,... - 'penalties_scheme', penalties_scheme,... - 'penalties_res', penalties_res); + 'penalties_scheme', {penalties_scheme}); end function diffOpsStruct = diffOpsMultiGrid1(scheme, g, schemeOrder, residualOrder, schemeParams, opSet, BCs) % DiffOps for stabilized solution vector - [D_scheme, penalties_scheme, D] = rv.diffops.constructSchemeDiffOp(scheme, g, schemeOrder, schemeParams, opSet, BCs); + [D_scheme, penalties_scheme, D_f] = rv.diffops.constructSchemeDiffOps(scheme, g, schemeOrder, schemeParams, opSet, BCs); % DiffOp for unstabilized solution vector - D_coarse = coarserGridDiffOp(scheme, g, residualOrder, schemeParams, opSet, BCs); + D_coarse = coarserGridDiffOpWithClosures(scheme, g, residualOrder, schemeParams, opSet, BCs); %% DiffOps for residual viscosity - [D_flux, penalties_res] = rv.diffops.constructFluxDiffOpWithClosures(scheme, g, residualOrder, schemeParams, opSet, BCs); + D_flux = rv.diffops.constructFluxDiffOps(scheme, g, residualOrder, schemeParams, opSet, BCs); D_flux = @(v) -D_flux(v); - D_t = @(v) D(v); + D_t = @(v) D_f(v); + % TODO: Use residual penalties as well here? diffOpsStruct = struct('D_scheme', D_scheme,... 'D_coarse', D_coarse,... 'D_flux', D_flux,... 'D_t', D_t,... - 'penalties_scheme', penalties_scheme,... - 'penalties_res', penalties_res); + 'penalties_scheme', {penalties_scheme}); end function diffOpsStruct = diffOpsMultiGrid2(scheme, g, schemeOrder, residualOrder, schemeParams, opSet, BCs) % DiffOps for stabilized solution vector - [D_scheme, penalties_scheme, D] = rv.diffops.constructSchemeDiffOp(scheme, g, schemeOrder, schemeParams, opSet, BCs); + [D_scheme, penalties_scheme] = rv.diffops.constructSchemeDiffOps(scheme, g, schemeOrder, schemeParams, opSet, BCs); % TODO: What orders to use here? D_coarse = coarserGridDiffOp(scheme, g, residualOrder, schemeParams, opSet, BCs); %% DiffOps for residual viscosity - D_flux = rv.diffops.constructFluxDiffOpWithClosures(scheme, g, schemeOrder, schemeParams, opSet, BCs); + D_flux = rv.diffops.constructFluxDiffOps(scheme, g, schemeOrder , schemeParams, opSet, BCs); D_flux = @(v) -D_flux(v); D_t = @(v) D_coarse(v); diffOpsStruct = struct('D_scheme', D_scheme,... 'D_flux', D_flux,... 'D_t', D_t,... - 'penalties_scheme', penalties_scheme); + 'penalties_scheme', {penalties_scheme}); +end +%% Multigrid functions: +% TODO: Implement properly. +function v_c = coarse_v(g, g_c, v) + V = grid.funcToMatrix(g,v); + v_c = reshape(V(1:2:end,1:2:end)',g_c.N(),1); end -% TODO: Only works for equidistant grids -function D_coarse = coarserGridDiffOp(scheme, g, schemeOrder, schemeParams, opSet, BCs) +% Note: Only works for equidistant grids. +% TODO: Change from using matlabs interp to using proper interpolation operators.. +function D_c = coarserGridDiffOp(scheme, g, schemeOrder, schemeParams, opSet, BCs) m = g.m(); lim = g.lim(); - m_coarse = (m-1)/2 + 1; - g_coarse = grid.equidistant(m_coarse, lim{1}); - D_coarse = rv.diffops.constructFluxDiffOpWithClosures(scheme, g_coarse, schemeOrder, schemeParams, opSet, BCs); - x = g.x{1}; - x_coarse = x(1:2:end); - % TODO: Fix interpolation - D_coarse = @(v) interp1(x_coarse,D_coarse(v(1:2:end)),x,'spline'); -end \ No newline at end of file + m_c = (m-1)/2 + 1; + switch g.D() + case 1 + g_c = grid.equidistant(m_c, lim{1}); + D_c = rv.diffops.constructFluxDiffOps(scheme, g_c, schemeOrder, schemeParams, opSet, BCs); + x = g.x{1}; + x_c = x(1:2:end); + D_c = @(v) interp1(x_c,D_c(v(1:2:end)),x,'spline'); + case 2 + g_c = grid.equidistant(m_c, lim{1}, lim{2}); + D_c = rv.diffops.constructFluxDiffOps(scheme, g_c, schemeOrder, schemeParams, opSet, BCs); + D_c = @(v) reshape(interp2(grid.funcToMatrix(g_c,D_c(coarse_v(g,g_c,v))),'spline')',g.N(),1); + end +end + +function D_c = coarserGridDiffOpWithClosures(scheme, g, schemeOrder, schemeParams, opSet, BCs) + m = g.m(); + lim = g.lim(); + m_c = (m-1)/2 + 1; + switch g.D() + case 1 + g_c = grid.equidistant(m_c, lim{1}); + [D_c, closures] = rv.diffops.constructFluxDiffOps(scheme, g_c, schemeOrder, schemeParams, opSet, BCs); + D_c = rv.diffops.addClosuresToDiffOp(D_c, closures); + x = g.x{1}; + x_c = x(1:2:end); + D_c = @(v) interp1(x_c,D_c(v(1:2:end)),x,'spline'); + case 2 + g_c = grid.equidistant(m_c, lim{1}, lim{2}); + [D_c, closures] = rv.diffops.constructFluxDiffOps(scheme, g_c, schemeOrder, schemeParams, opSet, BCs); + D_c = rv.diffops.addClosuresToDiffOp(D_c, closures); + D_c = @(v) reshape(interp2(grid.funcToMatrix(g_c,D_c(coarse_v(g,g_c,v))),'spline')',g.N(),1); + end +end +
diff -r 6cb03209f0a7 -r 79618b58b0a0 +rv/+diffops/constructFluxDiffOp.m --- a/+rv/+diffops/constructFluxDiffOp.m Mon Jul 08 14:50:46 2019 +0200 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,8 +0,0 @@ -function [D, diffOp] = constructFluxDiffOp(scheme, g, order, schemeParams, opSet, BCs) - diffOp = scheme(g, order, schemeParams{:}, opSet); - if ~isa(diffOp.D, 'function_handle') - D = @(v) diffOp.D*v; - else - D = diffOp.D; - end -end \ No newline at end of file
diff -r 6cb03209f0a7 -r 79618b58b0a0 +rv/+diffops/constructFluxDiffOpWithClosures.m --- a/+rv/+diffops/constructFluxDiffOpWithClosures.m Mon Jul 08 14:50:46 2019 +0200 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,4 +0,0 @@ -function [D, penalties] = constructFluxDiffOpWithClosures(scheme, g, order, schemeParams, opSet, BCs) - [D, diffOp] = rv.diffops.constructFluxDiffOp(scheme, g, order, schemeParams, opSet, BCs); - [D, penalties] = rv.diffops.addClosuresToDiffOp(diffOp, D, BCs); -end \ No newline at end of file
diff -r 6cb03209f0a7 -r 79618b58b0a0 +rv/+diffops/constructFluxDiffOps.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+rv/+diffops/constructFluxDiffOps.m Mon Jul 08 15:12:19 2019 +0200 @@ -0,0 +1,15 @@ +function [D_f, closures, penalties] = constructFluxDiffOps(scheme, g, order, schemeParams, opSet, BCs) + diffOp = scheme(g, order, schemeParams{:}, opSet); + if ~isa(diffOp.D, 'function_handle') + D_f = @(v) diffOp.D*v; + else + D_f = diffOp.D; + end + penalties = cell(size(BCs)); + closures = cell(size(BCs)); + for i = 1:size(BCs,1) + for j = 1:size(BCs,2) + [closures{i,j}, penalties{i,j}] = diffOp.boundary_condition(BCs{i,j}.boundary, BCs{i,j}.type); + end + end +end \ No newline at end of file
diff -r 6cb03209f0a7 -r 79618b58b0a0 +rv/+diffops/constructSchemeDiffOp.m --- a/+rv/+diffops/constructSchemeDiffOp.m Mon Jul 08 14:50:46 2019 +0200 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,6 +0,0 @@ -function [D_scheme, penalties_scheme, D] = constructSchemeDiffOp(scheme, g, schemeOrder, schemeParams, opSet, BCs) - %% DiffOps for solution vector - [D, penalties_scheme] = rv.diffops.constructFluxDiffOpWithClosures(scheme, g, schemeOrder, schemeParams, opSet, BCs); - D2 = rv.diffops.constructSymmetricD2(g, schemeOrder, opSet); - D_scheme = @(v,Viscosity)(D(v) + D2(Viscosity)*v); -end \ No newline at end of file
diff -r 6cb03209f0a7 -r 79618b58b0a0 +rv/+diffops/constructSchemeDiffOps.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+rv/+diffops/constructSchemeDiffOps.m Mon Jul 08 15:12:19 2019 +0200 @@ -0,0 +1,7 @@ +function [D_scheme, penalties, D_f] = constructSchemeDiffOps(scheme, g, order, schemeParams, opSet, BCs) + %% DiffOps for solution vector + [D_f, closures, penalties] = rv.diffops.constructFluxDiffOps(scheme, g, order, schemeParams, opSet, BCs); + D_scheme = rv.diffops.addClosuresToDiffOp(D_f, closures); + D2 = rv.diffops.constructSymmetricD2(g, order, opSet); + D_scheme = @(v,viscosity)(D_scheme(v) + D2(viscosity)*v); +end \ No newline at end of file