Mercurial > repos > public > sbplib

comparison +rv/ResidualViscosity.m @ 1193:921595039ab8 feature/rv

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Add options for postprocessing RV
author Vidar Stiernström <vidar.stiernstrom@it.uu.se>
date Mon, 29 Jul 2019 16:44:21 +0200
parents 7173b6fd4063
children bd5383809917
comparison
equal deleted inserted replaced
1192:b3c47a716d57 1193:921595039ab8
12 Mfirst % Coefficients for the first order viscosity 12 Mfirst % Coefficients for the first order viscosity
13 fRes % Function handle for computing the residual. 13 fRes % Function handle for computing the residual.
14 end 14 end
15 15
16 methods 16 methods
17 % TODO: pass opt struct with waveSpeed, normalization etc.
17 % TBD: Decide on how to treat waveSpeed. It would be nice to just pass a constant value without 18 % TBD: Decide on how to treat waveSpeed. It would be nice to just pass a constant value without
18 % wrapping it in a function. 19 % wrapping it in a function.
19 function obj = ResidualViscosity(g, Df, waveSpeed, Cmax, Cres, h, normalization, postProcess) 20 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))); 21 %default_arg('normalization',@(v)abs(obj.minmaxDiffNeighborhood1d(v)-norm(v-mean(v),inf)));
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
28 obj.normalization = normalization; 28 obj.normalization = normalization;
29 obj.g = g; 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 {'filt', 'filter'}
36 order = 4; 36 % TBD: Keep?
37 order =
37 F = obj.shapiroFilter(obj.g, order); 38 F = obj.shapiroFilter(obj.g, order);
38 obj.Mres = F*obj.Mres; 39 obj.Mres = F*obj.Mres;
39 obj.fRes = @(v,dvdt) obj.Mres*abs(dvdt + obj.Df(v))./obj.normalization(v); 40 obj.fRes = @(v,dvdt) obj.Mres*abs(dvdt + obj.Df(v))./obj.normalization(v);
40 case 'maximum neighbors' 41 case {'max', 'maximum neighbors'}
41 switch g.D() 42 switch g.D()
42 case 1 43 case 1
43 obj.fRes = @(v,dvdt) movmax(obj.Mres*abs(dvdt + obj.Df(v))./obj.normalization(v),3); 44 obj.fRes = @(v,dvdt) movmax(obj.Mres*abs(dvdt + obj.Df(v))./obj.normalization(v),3);
44 case 2 45 case 2
45 obj.fRes = @obj.maxResidualNeighbors2d; 46 obj.fRes = @obj.maxResidualNeighbors2d;