Mercurial > repos > public > sbplib
changeset 1188:7173b6fd4063 feature/rv
Rename class property grid to g, to avoid confusion with grid package
| author | Vidar Stiernström <vidar.stiernstrom@it.uu.se> |
|---|---|
| date | Mon, 08 Jul 2019 14:50:19 +0200 |
| parents | 5aa3049a4212 |
| children | 6cb03209f0a7 |
| files | +rv/ResidualViscosity.m |
| diffstat | 1 files changed, 17 insertions(+), 15 deletions(-) [+] |
line wrap: on
line diff
--- a/+rv/ResidualViscosity.m Mon Jul 08 14:48:56 2019 +0200 +++ b/+rv/ResidualViscosity.m Mon Jul 08 14:50:19 2019 +0200 @@ -1,6 +1,6 @@ classdef ResidualViscosity < handle properties - grid % grid + g % grid Df % Diff op approximating the gradient of the flux f(u) waveSpeed % Wave speed at each grid point, e.g f'(u). %TBD: Better naming? Cmax % Constant controlling relative amount of upwind dissipation @@ -16,7 +16,7 @@ methods % TBD: Decide on how to treat waveSpeed. It would be nice to just pass a constant value without % wrapping it in a function. - function obj = ResidualViscosity(grid, Df, waveSpeed, Cmax, Cres, h, normalization, postProcess) + function obj = ResidualViscosity(g, Df, waveSpeed, Cmax, Cres, h, normalization, postProcess) default_arg('normalization',@(v)abs(obj.minmaxDiffNeighborhood1d(v)-norm(v-mean(v),inf))); default_arg('postProcess','maximum neighbors'); obj.Df = Df; @@ -26,17 +26,19 @@ obj.Cres = Cres; obj.normalization = normalization; - obj.grid = grid; + obj.g = g; obj.Mres = obj.Cres*obj.h^2; obj.Mfirst = obj.Cmax*obj.h; switch postProcess case 'none' obj.fRes = @(v,dvdt) obj.Mres*abs(dvdt + obj.Df(v))./obj.normalization(v); case 'filter' - F = obj.shapiroFilter(obj.grid, order); - obj.fRes = @(v,dvdt) F*obj.Mres*abs(dvdt + obj.Df(v))./obj.normalization(v); + order = 4; + F = obj.shapiroFilter(obj.g, order); + obj.Mres = F*obj.Mres; + obj.fRes = @(v,dvdt) obj.Mres*abs(dvdt + obj.Df(v))./obj.normalization(v); case 'maximum neighbors' - switch grid.D() + switch g.D() case 1 obj.fRes = @(v,dvdt) movmax(obj.Mres*abs(dvdt + obj.Df(v))./obj.normalization(v),3); case 2 @@ -58,8 +60,8 @@ function res = maxResidualNeighbors2d(obj,v,dvdt) res = obj.Mres*abs(dvdt + obj.Df(v))./obj.normalization(v); - resMat = grid.funcToMatrix(obj.grid,res); - res = reshape(max(movmax(resMat,3,1),movmax(resMat,3,2))',obj.grid.N(),1); + resMat = grid.funcToMatrix(obj.g,res); + res = reshape(max(movmax(resMat,3,1),movmax(resMat,3,2))',obj.g.N(),1); end end methods (Static) @@ -69,23 +71,23 @@ minmaxDiff = umax - umin; end - function minmaxDiff = minmaxDiffNeighborhood2d(grid, u) - uMatrix = grid.funcToMatrix(grid,u); + function minmaxDiff = minmaxDiffNeighborhood2d(g, u) + uMatrix = grid.funcToMatrix(g,u); umax = max(movmax(uMatrix,3,1),movmax(uMatrix,3,2)); umin = min(movmin(uMatrix,3,1),movmin(uMatrix,3,2)); minmaxDiff = umax - umin; minmaxDiff = reshape(minmaxDiff',g.N(),1); end - function F = shapiroFilter(grid, order) + function F = shapiroFilter(g, order) switch order case 2 - F = spdiags(repmat([1 2 1],grid.m(1),1), -1:1, grid.m(1), grid.m(1)); + F = spdiags(repmat([1 2 1],g.m(1),1), -1:1, g.m(1), g.m(1)); case 4 - F = spdiags(repmat([-1 4 10 4 -1],grid.m(1),1), -2:2, grid.m(1), grid.m(1)); + F = spdiags(repmat([-1 4 10 4 -1],g.m(1),1), -2:2, g.m(1), g.m(1)); case 6 - F = spdiags(repmat([1 -6 15 44 15 -6 1],grid.m(1),1), -3:3, grid.m(1), grid.m(1)); + F = spdiags(repmat([1 -6 15 44 15 -6 1],g.m(1),1), -3:3, g.m(1), g.m(1)); case 8 - F = spdiags(repmat([-1 8 -28 56 186 56 -28 8 -1],grid.m(1),1), -4:4, grid.m(1), grid.m(1)); + F = spdiags(repmat([-1 8 -28 56 186 56 -28 8 -1],g.m(1),1), -4:4, g.m(1), g.m(1)); end F = 1/2^order * F; % Fx = spdiags(repmat(1/(mid+2)*[1 mid 1],g.m(1),1), -1:1, g.m(1), g.m(1));