Mercurial > repos > public > sbplib

comparison +rv/ResidualViscosity.m @ 1151:03ecf18d035f feature/rv

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Set RV artificial viscosity to zero around borders of the domain
author Vidar Stiernström <vidar.stiernstrom@it.uu.se>
date Mon, 18 Feb 2019 09:00:00 +0100
parents 922996695952
children 010bb2677230
comparison
equal deleted inserted replaced
1150:922996695952 1151:03ecf18d035f
5 Cmax % Constant controlling relative amount of upwind dissipation 5 Cmax % Constant controlling relative amount of upwind dissipation
6 Cres % Constant controling relative amount of upwind dissipation 6 Cres % Constant controling relative amount of upwind dissipation
7 h % Length scale used for scaling the viscosity. Typically grid spacing. 7 h % Length scale used for scaling the viscosity. Typically grid spacing.
8 normalization % Function used to normalize the residual such that it is amplified in the 8 normalization % Function used to normalize the residual such that it is amplified in the
9 % shocks. 9 % shocks.
10 boundaryIndices % indices of the boundary of the domain
10 end 11 end
11 12
12 methods 13 methods
13 % TBD: Decide on how to treat waveSpeed. It would be nice to just pass a constant value without 14 % TBD: Decide on how to treat waveSpeed. It would be nice to just pass a constant value without
14 % wrapping it in a function. 15 % wrapping it in a function.
15 function obj = ResidualViscosity(Df, waveSpeed, Cmax, Cres, h, normalization) 16 function obj = ResidualViscosity(g, Df, waveSpeed, Cmax, Cres, h, normalization)
16 default_arg('normalization',@(v)abs(obj.minmaxDiffNeighborhood1d(v)-norm(v-mean(v),inf))); 17 default_arg('normalization',@(v)abs(obj.minmaxDiffNeighborhood1d(v)-norm(v-mean(v),inf)));
17 %default_arg('normalization',@(v)norm(v/2,inf));
18 obj.Df = Df; 18 obj.Df = Df;
19 obj.waveSpeed = waveSpeed; 19 obj.waveSpeed = waveSpeed;
20 obj.h = h; 20 obj.h = h;
21 obj.Cmax = Cmax; 21 obj.Cmax = Cmax;
22 obj.Cres = Cres; 22 obj.Cres = Cres;
23 obj.normalization = normalization; 23 obj.normalization = normalization;
24 obj.boundaryIndices = obj.getBoundaryIndices(g);
24 end 25 end
25 26
26 function [viscosity, Df, firstOrderViscosity, residualViscosity] = evaluate(obj, v, dvdt) 27 function [viscosity, Df, firstOrderViscosity, residualViscosity] = evaluate(obj, v, dvdt)
27 Df = obj.Df(v); 28 Df = obj.Df(v);
28 firstOrderViscosity = obj.Cmax*obj.h*abs(obj.waveSpeed(v)); 29 firstOrderViscosity = obj.Cmax*obj.h*abs(obj.waveSpeed(v));
29 residualViscosity = obj.Cres*obj.h^2*abs(dvdt + Df)./obj.normalization(v); 30 residualViscosity = obj.Cres*obj.h^2*abs(dvdt + Df)./obj.normalization(v);
30 viscosity = min(firstOrderViscosity, residualViscosity); 31 viscosity = min(firstOrderViscosity, residualViscosity);
32 % TODO: If the boundary conditions are implemented correctly, this might not be needed.
33 viscosity(obj.boundaryIndices) = 0;
31 end 34 end
32 end 35 end
33 methods (Static) 36 methods (Static)
37
38 function boundaryind = getBoundaryIndices(g)
39 ind = grid.funcToMatrix(g, 1:g.N());
40 switch g.D()
41 case 1
42 boundaryind = [ind(1) ind(end)];
43 case 2
44 boundaryind = [ind(1,:) ind(end,:), ind(:,1)' ind(:,end)'];
45 case 3
46 error;
47 end
48 end
34 49
35 % function f = minmaxDiffNeighborhood(g) 50 % function f = minmaxDiffNeighborhood(g)
36 % switch g.D() 51 % switch g.D()
37 % case 1 52 % case 1
38 % f = @(u)minmaxDiffNeighborhood1d(u); 53 % f = @(u)minmaxDiffNeighborhood1d(u);