Mercurial > repos > public > sbplib

comparison +rv/ResidualViscosity.m @ 1159:1ad7da049b50 feature/rv

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Fix bug using 2d maximum neighborhood in 1d
author Vidar Stiernström <vidar.stiernstrom@it.uu.se>
date Tue, 25 Jun 2019 15:04:53 +0200
parents 2ba63553ccfc
children 856bd6291d17
comparison
equal deleted inserted replaced
1158:2ba63553ccfc 1159:1ad7da049b50
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(grid, 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','none');
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
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 'maximum neighbors' 35 case 'maximum neighbors'
36 obj.fRes = @obj.maxResidualNeighbors 36 switch grid.D()
37 case 1
38 obj.fRes = @(v,dvdt) movmax(obj.Mres*abs(dvdt + obj.Df(v))./obj.normalization(v),3)
39 case 2
40 obj.fRes = @obj.maxResidualNeighbors2d
41 end
37 end 42 end
38 end 43 end
39 44
40 function viscosity = evaluateViscosity(obj, v, dvdt) 45 function viscosity = evaluateViscosity(obj, v, dvdt)
41 viscosity = min(obj.Mfirst*abs(obj.waveSpeed(v)), obj.fRes(v,dvdt)); 46 firstOrderViscosity = obj.Mfirst*abs(obj.waveSpeed(v));
47 viscosity = min(firstOrderViscosity, obj.fRes(v,dvdt));
42 end 48 end
43 49
44 function [viscosity, Df, firstOrderViscosity, residualViscosity] = evaluate(obj, v, dvdt) 50 function [viscosity, Df, firstOrderViscosity, residualViscosity] = evaluate(obj, v, dvdt)
45 Df = obj.Df(v); 51 Df = obj.Df(v);
46 firstOrderViscosity = obj.Mfirst*abs(obj.waveSpeed(v)); 52 firstOrderViscosity = obj.Mfirst*abs(obj.waveSpeed(v));
47 residualViscosity = obj.fRes(v,dvdt); 53 residualViscosity = obj.fRes(v,dvdt);
48 viscosity = min(firstOrderViscosity, obj.fRes(v,dvdt)); 54 viscosity = min(firstOrderViscosity, residualViscosity);
49 end 55 end
50 function res = maxResidualNeighbors(obj,v,dvdt) 56
57 function res = maxResidualNeighbors2d(obj,v,dvdt)
51 res = obj.Mres*abs(dvdt + obj.Df(v))./obj.normalization(v); 58 res = obj.Mres*abs(dvdt + obj.Df(v))./obj.normalization(v);
52 resMat = grid.funcToMatrix(obj.grid,res); 59 resMat = grid.funcToMatrix(obj.grid,res);
53 res = reshape(max(movmax(resMat,3,1),movmax(resMat,3,2))',obj.grid.N(),1); 60 res = reshape(max(movmax(resMat,3,1),movmax(resMat,3,2))',obj.grid.N(),1);
54 end 61 end
55 end 62 end