Mercurial > repos > public > sbplib

comparison +rv/ResidualViscosity.m @ 1150:922996695952 feature/rv

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Implement local residual normalization functions
author Vidar Stiernström <vidar.stiernstrom@it.uu.se>
date Thu, 24 Jan 2019 10:23:52 +0100
parents a8ee5eca0e6c
children 03ecf18d035f
comparison
equal deleted inserted replaced
1149:1fe48cbd379a 1150:922996695952
11 11
12 methods 12 methods
13 % TBD: Decide on how to treat waveSpeed. It would be nice to just pass a constant value without 13 % TBD: Decide on how to treat waveSpeed. It would be nice to just pass a constant value without
14 % wrapping it in a function. 14 % wrapping it in a function.
15 function obj = ResidualViscosity(Df, waveSpeed, Cmax, Cres, h, normalization) 15 function obj = ResidualViscosity(Df, waveSpeed, Cmax, Cres, h, normalization)
16 default_arg('normalization',@(v)norm(v/2,inf)); 16 default_arg('normalization',@(v)abs(obj.minmaxDiffNeighborhood1d(v)-norm(v-mean(v),inf)));
17 %default_arg('normalization',@(v)norm(v/2,inf));
17 obj.Df = Df; 18 obj.Df = Df;
18 obj.waveSpeed = waveSpeed; 19 obj.waveSpeed = waveSpeed;
19 obj.h = h; 20 obj.h = h;
20 obj.Cmax = Cmax; 21 obj.Cmax = Cmax;
21 obj.Cres = Cres; 22 obj.Cres = Cres;
27 firstOrderViscosity = obj.Cmax*obj.h*abs(obj.waveSpeed(v)); 28 firstOrderViscosity = obj.Cmax*obj.h*abs(obj.waveSpeed(v));
28 residualViscosity = obj.Cres*obj.h^2*abs(dvdt + Df)./obj.normalization(v); 29 residualViscosity = obj.Cres*obj.h^2*abs(dvdt + Df)./obj.normalization(v);
29 viscosity = min(firstOrderViscosity, residualViscosity); 30 viscosity = min(firstOrderViscosity, residualViscosity);
30 end 31 end
31 end 32 end
33 methods (Static)
34
35 % function f = minmaxDiffNeighborhood(g)
36 % switch g.D()
37 % case 1
38 % f = @(u)minmaxDiffNeighborhood1d(u);
39 % case 2
40 % f = @(u)minmaxDiffNeighborhood2d(g,u);
41 % end
42 % end
43
44 function minmaxDiff = minmaxDiffNeighborhood1d(u)
45 umax = movmax(u,3);
46 umin = movmin(u,3);
47 minmaxDiff = umax - umin;
48 end
49 function minmaxDiff = minmaxDiffNeighborhood2d(g, u)
50 m = g.m();
51 uMatrix = grid.funcToMatrix(g,u);
52 umax = max(movmax(uMatrix,3,1),movmax(uMatrix,3,2));
53 umin = min(movmin(uMatrix,3,1),movmin(uMatrix,3,2));
54 minmaxDiff = umax - umin;
55 minmaxDiff = reshape(minmaxDiff',m(1)*m(2),1);
56 end
57 end
32 end 58 end