Mercurial > repos > public > sbplib
diff +rv/+diffops/constructDiffOps.m @ 1184:ecc605453733 feature/rv
Move specialized constructDiffOps* into constructDiffOps.m. Add multigrid RV using the standard RungekuttaRV time stepper
| author | Vidar Stiernström <vidar.stiernstrom@it.uu.se> |
|---|---|
| date | Fri, 05 Jul 2019 17:51:11 +0200 |
| parents | e3d8f24b2c1c |
| children | 79618b58b0a0 |
line wrap: on
line diff
--- a/+rv/+diffops/constructDiffOps.m Fri Jul 05 17:48:54 2019 +0200 +++ b/+rv/+diffops/constructDiffOps.m Fri Jul 05 17:51:11 2019 +0200 @@ -1,9 +1,98 @@ function diffOpsStruct = constructDiffOps(method, varargin) - switch method - case {'bdf','backwards-difference-formula'} - diffOpsStruct = rv.diffops.constructDiffOpsBdf(varargin{:}); - case {'ms','multi-stage'} - diffOpsStruct = rv.diffops.constructDiffOpsMultiStage(varargin{:}); - case {'mg','multi-grid'} - diffOpsStruct = rv.diffops.constructDiffOpsMultiGrid(varargin{:}); + switch method + case {'standard'} + diffOpsStruct = diffOpsStandard(varargin{:}); + case {'bdf','backwards-difference-formula'} + diffOpsStruct = diffOpsBdf(varargin{:}); + case {'ms','multi-stage'} + diffOpsStruct = diffOpsMultiStage(varargin{:}); + case {'mg1','multi-grid1'} + diffOpsStruct = diffOpsMultiGrid1(varargin{:}); + case {'mg2','multi-grid2'} + diffOpsStruct = diffOpsMultiGrid2(varargin{:}); + end +end + +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); + %% DiffOps for residual viscosity + D_flux = rv.diffops.constructFluxDiffOpWithClosures(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); +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); + %% DiffOps for residual viscosity + [D_flux, ~] = rv.diffops.constructFluxDiffOp(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); +end + +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); + % DiffOp for unstabilized solution vector + D_unstable = rv.diffops.constructFluxDiffOpWithClosures(scheme, g, schemeOrder, schemeParams, opSet, BCs); + %% DiffOps for residual viscosity + [D_t, penalties_res] = rv.diffops.constructFluxDiffOpWithClosures(scheme, g, residualOrder, schemeParams, opSet, BCs); + D_flux = @(v) -D_t(v); + 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); +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); + % DiffOp for unstabilized solution vector + D_coarse = coarserGridDiffOp(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 = @(v) -D_flux(v); + D_t = @(v) D(v); + 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); +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); + % 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 = @(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); +end + +% TODO: Only works for equidistant grids +function D_coarse = 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