Mercurial > repos > public > sbplib

comparison +rv/ResidualViscosity.m @ 1158:2ba63553ccfc feature/rv

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Add option to compute maximum neighbordhood residual
author Vidar Stiernström <vidar.stiernstrom@it.uu.se>
date Tue, 25 Jun 2019 14:14:25 +0200
parents 010bb2677230
children 1ad7da049b50
comparison
equal deleted inserted replaced
1157:82315fa6adb1 1158:2ba63553ccfc
1 classdef ResidualViscosity < handle 1 classdef ResidualViscosity < handle
2 properties 2 properties
3 grid % grid
3 Df % Diff op approximating the gradient of the flux f(u) 4 Df % Diff op approximating the gradient of the flux f(u)
4 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?
5 Cmax % Constant controlling relative amount of upwind dissipation 6 Cmax % Constant controlling relative amount of upwind dissipation
6 Cres % Constant controling relative amount of upwind dissipation 7 Cres % Constant controling relative amount of upwind dissipation
7 h % Length scale used for scaling the viscosity. Typically grid spacing. 8 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 9 normalization % Function used to normalize the residual such that it is amplified in the
9 % shocks. 10 % shocks and suppressed elsewhere.
10 boundaryIndices % indices of the boundary of the domain 11 Mres % Coefficients for the residual viscosity
12 Mfirst % Coefficients for the first order viscosity
13 fRes % Function handle for computing the residual.
11 end 14 end
12 15
13 methods 16 methods
14 % 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
15 % wrapping it in a function. 18 % wrapping it in a function.
16 function obj = ResidualViscosity(g, Df, waveSpeed, Cmax, Cres, h, normalization) 19 function obj = ResidualViscosity(grid, Df, waveSpeed, Cmax, Cres, h, normalization, postProcess)
17 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');
18 obj.Df = Df; 22 obj.Df = Df;
19 obj.waveSpeed = waveSpeed; 23 obj.waveSpeed = waveSpeed;
20 obj.h = h; 24 obj.h = h;
21 obj.Cmax = Cmax; 25 obj.Cmax = Cmax;
26
22 obj.Cres = Cres; 27 obj.Cres = Cres;
23 obj.normalization = normalization; 28 obj.normalization = normalization;
24 obj.boundaryIndices = obj.getBoundaryIndices(g); 29 obj.grid = grid;
30 obj.Mres = obj.Cres*obj.h^2;
31 obj.Mfirst = obj.Cmax*obj.h;
32 switch postProcess
33 case 'none'
34 obj.fRes = @(v,dvdt) obj.Mres*abs(dvdt + obj.Df(v))./obj.normalization(v);
35 case 'maximum neighbors'
36 obj.fRes = @obj.maxResidualNeighbors
37 end
25 end 38 end
26 39
27 function viscosity = evaluateViscosity(obj, v, dvdt) 40 function viscosity = evaluateViscosity(obj, v, dvdt)
28 viscosity = min(obj.Cmax*obj.h*abs(obj.waveSpeed(v)), obj.Cres*obj.h^2*abs(dvdt + obj.Df(v))./obj.normalization(v)); 41 viscosity = min(obj.Mfirst*abs(obj.waveSpeed(v)), obj.fRes(v,dvdt));
29 % TODO: If the boundary conditions are implemented correctly, this might not be needed.
30 viscosity(obj.boundaryIndices) = 0;
31 end 42 end
32 43
33 function [viscosity, Df, firstOrderViscosity, residualViscosity] = evaluate(obj, v, dvdt) 44 function [viscosity, Df, firstOrderViscosity, residualViscosity] = evaluate(obj, v, dvdt)
34 Df = obj.Df(v); 45 Df = obj.Df(v);
35 firstOrderViscosity = obj.Cmax*obj.h*abs(obj.waveSpeed(v)); 46 firstOrderViscosity = obj.Mfirst*abs(obj.waveSpeed(v));
36 residualViscosity = obj.Cres*obj.h^2*abs(dvdt + Df)./obj.normalization(v); 47 residualViscosity = obj.fRes(v,dvdt);
37 viscosity = min(firstOrderViscosity, residualViscosity); 48 viscosity = min(firstOrderViscosity, obj.fRes(v,dvdt));
38 % TODO: If the boundary conditions are implemented correctly, this might not be needed. 49 end
39 viscosity(obj.boundaryIndices) = 0; 50 function res = maxResidualNeighbors(obj,v,dvdt)
51 res = obj.Mres*abs(dvdt + obj.Df(v))./obj.normalization(v);
52 resMat = grid.funcToMatrix(obj.grid,res);
53 res = reshape(max(movmax(resMat,3,1),movmax(resMat,3,2))',obj.grid.N(),1);
40 end 54 end
41 end 55 end
42 methods (Static) 56 methods (Static)
43
44 function boundaryind = getBoundaryIndices(g)
45 ind = grid.funcToMatrix(g, 1:g.N());
46 switch g.D()
47 case 1
48 boundaryind = [ind(1) ind(end)];
49 case 2
50 boundaryind = [ind(1,:) ind(end,:), ind(:,1)' ind(:,end)'];
51 case 3
52 error;
53 end
54 end
55
56 function minmaxDiff = minmaxDiffNeighborhood1d(u) 57 function minmaxDiff = minmaxDiffNeighborhood1d(u)
57 umax = movmax(u,3); 58 umax = movmax(u,3);
58 umin = movmin(u,3); 59 umin = movmin(u,3);
59 minmaxDiff = umax - umin; 60 minmaxDiff = umax - umin;
60 end 61 end
61 function minmaxDiff = minmaxDiffNeighborhood2d(g, u) 62
62 m = g.m(); 63 function minmaxDiff = minmaxDiffNeighborhood2d(grid, u)
63 uMatrix = grid.funcToMatrix(g,u); 64 uMatrix = grid.funcToMatrix(grid,u);
64 umax = max(movmax(uMatrix,3,1),movmax(uMatrix,3,2)); 65 umax = max(movmax(uMatrix,3,1),movmax(uMatrix,3,2));
65 umin = min(movmin(uMatrix,3,1),movmin(uMatrix,3,2)); 66 umin = min(movmin(uMatrix,3,1),movmin(uMatrix,3,2));
66 minmaxDiff = umax - umin; 67 minmaxDiff = umax - umin;
67 minmaxDiff = reshape(minmaxDiff',m(1)*m(2),1); 68 minmaxDiff = reshape(minmaxDiff',g.N(),1);
68 end 69 end
69 end 70 end
70 end 71 end