Mercurial > repos > public > sbplib

comparison +rv/ResidualViscosity.m @ 1188:7173b6fd4063 feature/rv

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
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 3364a51f0d9e
children 921595039ab8
comparison
equal deleted inserted replaced
1187:5aa3049a4212 1188:7173b6fd4063
1 classdef ResidualViscosity < handle 1 classdef ResidualViscosity < handle
2 properties 2 properties
3 grid % grid 3 g % grid
4 Df % Diff op approximating the gradient of the flux f(u) 4 Df % Diff op approximating the gradient of the flux f(u)
5 waveSpeed % Wave speed at each grid point, e.g f'(u). %TBD: Better naming? 5 waveSpeed % Wave speed at each grid point, e.g f'(u). %TBD: Better naming?
6 Cmax % Constant controlling relative amount of upwind dissipation 6 Cmax % Constant controlling relative amount of upwind dissipation
7 Cres % Constant controling relative amount of upwind dissipation 7 Cres % Constant controling relative amount of upwind dissipation
8 h % Length scale used for scaling the viscosity. Typically grid spacing. 8 h % Length scale used for scaling the viscosity. Typically grid spacing.
14 end 14 end
15 15
16 methods 16 methods
17 % TBD: Decide on how to treat waveSpeed. It would be nice to just pass a constant value without 17 % TBD: Decide on how to treat waveSpeed. It would be nice to just pass a constant value without
18 % wrapping it in a function. 18 % wrapping it in a function.
19 function obj = ResidualViscosity(grid, Df, waveSpeed, Cmax, Cres, h, normalization, postProcess) 19 function obj = ResidualViscosity(g, Df, waveSpeed, Cmax, Cres, h, normalization, postProcess)
20 default_arg('normalization',@(v)abs(obj.minmaxDiffNeighborhood1d(v)-norm(v-mean(v),inf))); 20 default_arg('normalization',@(v)abs(obj.minmaxDiffNeighborhood1d(v)-norm(v-mean(v),inf)));
21 default_arg('postProcess','maximum neighbors'); 21 default_arg('postProcess','maximum neighbors');
22 obj.Df = Df; 22 obj.Df = Df;
23 obj.waveSpeed = waveSpeed; 23 obj.waveSpeed = waveSpeed;
24 obj.h = h; 24 obj.h = h;
25 obj.Cmax = Cmax; 25 obj.Cmax = Cmax;
26 26
27 obj.Cres = Cres; 27 obj.Cres = Cres;
28 obj.normalization = normalization; 28 obj.normalization = normalization;
29 obj.grid = grid; 29 obj.g = g;
30 obj.Mres = obj.Cres*obj.h^2; 30 obj.Mres = obj.Cres*obj.h^2;
31 obj.Mfirst = obj.Cmax*obj.h; 31 obj.Mfirst = obj.Cmax*obj.h;
32 switch postProcess 32 switch postProcess
33 case 'none' 33 case 'none'
34 obj.fRes = @(v,dvdt) obj.Mres*abs(dvdt + obj.Df(v))./obj.normalization(v); 34 obj.fRes = @(v,dvdt) obj.Mres*abs(dvdt + obj.Df(v))./obj.normalization(v);
35 case 'filter' 35 case 'filter'
36 F = obj.shapiroFilter(obj.grid, order); 36 order = 4;
37 obj.fRes = @(v,dvdt) F*obj.Mres*abs(dvdt + obj.Df(v))./obj.normalization(v); 37 F = obj.shapiroFilter(obj.g, order);
38 obj.Mres = F*obj.Mres;
39 obj.fRes = @(v,dvdt) obj.Mres*abs(dvdt + obj.Df(v))./obj.normalization(v);
38 case 'maximum neighbors' 40 case 'maximum neighbors'
39 switch grid.D() 41 switch g.D()
40 case 1 42 case 1
41 obj.fRes = @(v,dvdt) movmax(obj.Mres*abs(dvdt + obj.Df(v))./obj.normalization(v),3); 43 obj.fRes = @(v,dvdt) movmax(obj.Mres*abs(dvdt + obj.Df(v))./obj.normalization(v),3);
42 case 2 44 case 2
43 obj.fRes = @obj.maxResidualNeighbors2d; 45 obj.fRes = @obj.maxResidualNeighbors2d;
44 end 46 end
56 viscosity = min(firstOrderViscosity, residualViscosity); 58 viscosity = min(firstOrderViscosity, residualViscosity);
57 end 59 end
58 60
59 function res = maxResidualNeighbors2d(obj,v,dvdt) 61 function res = maxResidualNeighbors2d(obj,v,dvdt)
60 res = obj.Mres*abs(dvdt + obj.Df(v))./obj.normalization(v); 62 res = obj.Mres*abs(dvdt + obj.Df(v))./obj.normalization(v);
61 resMat = grid.funcToMatrix(obj.grid,res); 63 resMat = grid.funcToMatrix(obj.g,res);
62 res = reshape(max(movmax(resMat,3,1),movmax(resMat,3,2))',obj.grid.N(),1); 64 res = reshape(max(movmax(resMat,3,1),movmax(resMat,3,2))',obj.g.N(),1);
63 end 65 end
64 end 66 end
65 methods (Static) 67 methods (Static)
66 function minmaxDiff = minmaxDiffNeighborhood1d(u) 68 function minmaxDiff = minmaxDiffNeighborhood1d(u)
67 umax = movmax(u,3); 69 umax = movmax(u,3);
68 umin = movmin(u,3); 70 umin = movmin(u,3);
69 minmaxDiff = umax - umin; 71 minmaxDiff = umax - umin;
70 end 72 end
71 73
72 function minmaxDiff = minmaxDiffNeighborhood2d(grid, u) 74 function minmaxDiff = minmaxDiffNeighborhood2d(g, u)
73 uMatrix = grid.funcToMatrix(grid,u); 75 uMatrix = grid.funcToMatrix(g,u);
74 umax = max(movmax(uMatrix,3,1),movmax(uMatrix,3,2)); 76 umax = max(movmax(uMatrix,3,1),movmax(uMatrix,3,2));
75 umin = min(movmin(uMatrix,3,1),movmin(uMatrix,3,2)); 77 umin = min(movmin(uMatrix,3,1),movmin(uMatrix,3,2));
76 minmaxDiff = umax - umin; 78 minmaxDiff = umax - umin;
77 minmaxDiff = reshape(minmaxDiff',g.N(),1); 79 minmaxDiff = reshape(minmaxDiff',g.N(),1);
78 end 80 end
79 function F = shapiroFilter(grid, order) 81 function F = shapiroFilter(g, order)
80 switch order 82 switch order
81 case 2 83 case 2
82 F = spdiags(repmat([1 2 1],grid.m(1),1), -1:1, grid.m(1), grid.m(1)); 84 F = spdiags(repmat([1 2 1],g.m(1),1), -1:1, g.m(1), g.m(1));
83 case 4 85 case 4
84 F = spdiags(repmat([-1 4 10 4 -1],grid.m(1),1), -2:2, grid.m(1), grid.m(1)); 86 F = spdiags(repmat([-1 4 10 4 -1],g.m(1),1), -2:2, g.m(1), g.m(1));
85 case 6 87 case 6
86 F = spdiags(repmat([1 -6 15 44 15 -6 1],grid.m(1),1), -3:3, grid.m(1), grid.m(1)); 88 F = spdiags(repmat([1 -6 15 44 15 -6 1],g.m(1),1), -3:3, g.m(1), g.m(1));
87 case 8 89 case 8
88 F = spdiags(repmat([-1 8 -28 56 186 56 -28 8 -1],grid.m(1),1), -4:4, grid.m(1), grid.m(1)); 90 F = spdiags(repmat([-1 8 -28 56 186 56 -28 8 -1],g.m(1),1), -4:4, g.m(1), g.m(1));
89 end 91 end
90 F = 1/2^order * F; 92 F = 1/2^order * F;
91 % Fx = spdiags(repmat(1/(mid+2)*[1 mid 1],g.m(1),1), -1:1, g.m(1), g.m(1)); 93 % Fx = spdiags(repmat(1/(mid+2)*[1 mid 1],g.m(1),1), -1:1, g.m(1), g.m(1));
92 % Fx(1,:) = 0; 94 % Fx(1,:) = 0;
93 % Fx(1,1:2) = 1/(mid+1)*[mid 1]; 95 % Fx(1,1:2) = 1/(mid+1)*[mid 1];