Mercurial > repos > public > sbplib
diff +rv/+diffops/constructDiffOps.m @ 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 | ecc605453733 |
| children | 881afc40a3d2 |
line wrap: on
line diff
--- 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 +