Mercurial > repos > public > sbplib
view +rv/ResidualViscosity.m @ 1176:ebec2b86f539 feature/rv
Update comments for RungekuttaRvMultiStage/Grid
author | Vidar Stiernström <vidar.stiernstrom@it.uu.se> |
---|---|
date | Fri, 28 Jun 2019 13:50:04 +0200 |
parents | fc2631ba4da5 |
children | beecb580c5bf |
line wrap: on
line source
classdef ResidualViscosity < handle properties grid % grid Df % Diff op approximating the gradient of the flux f(u) waveSpeed % Wave speed at each grid point, e.g f'(u). %TBD: Better naming? Cmax % Constant controlling relative amount of upwind dissipation Cres % Constant controling relative amount of upwind dissipation h % Length scale used for scaling the viscosity. Typically grid spacing. normalization % Function used to normalize the residual such that it is amplified in the % shocks and suppressed elsewhere. Mres % Coefficients for the residual viscosity Mfirst % Coefficients for the first order viscosity fRes % Function handle for computing the residual. end methods % TBD: Decide on how to treat waveSpeed. It would be nice to just pass a constant value without % wrapping it in a function. function obj = ResidualViscosity(grid, Df, waveSpeed, Cmax, Cres, h, normalization, postProcess) default_arg('normalization',@(v)abs(obj.minmaxDiffNeighborhood1d(v)-norm(v-mean(v),inf))); default_arg('postProcess','maximum neighbors'); obj.Df = Df; obj.waveSpeed = waveSpeed; obj.h = h; obj.Cmax = Cmax; obj.Cres = Cres; obj.normalization = normalization; obj.grid = grid; obj.Mres = obj.Cres*obj.h^2; obj.Mfirst = obj.Cmax*obj.h; switch postProcess case 'none' obj.fRes = @(v,dvdt) obj.Mres*abs(dvdt + obj.Df(v))./obj.normalization(v); case 'maximum neighbors' switch grid.D() case 1 obj.fRes = @(v,dvdt) movmax(obj.Mres*abs(dvdt + obj.Df(v))./obj.normalization(v),3) case 2 obj.fRes = @obj.maxResidualNeighbors2d end end end function viscosity = evaluateViscosity(obj, v, dvdt) viscosity = min(obj.Mfirst*abs(obj.waveSpeed(v)), obj.fRes(v,dvdt)); end function [viscosity, Df, firstOrderViscosity, residualViscosity] = evaluate(obj, v, dvdt) Df = obj.Df(v); firstOrderViscosity = obj.Mfirst*abs(obj.waveSpeed(v)); residualViscosity = obj.fRes(v,dvdt); viscosity = min(firstOrderViscosity, residualViscosity); end function res = maxResidualNeighbors2d(obj,v,dvdt) res = obj.Mres*abs(dvdt + obj.Df(v))./obj.normalization(v); resMat = grid.funcToMatrix(obj.grid,res); res = reshape(max(movmax(resMat,3,1),movmax(resMat,3,2))',obj.grid.N(),1); end end methods (Static) function minmaxDiff = minmaxDiffNeighborhood1d(u) umax = movmax(u,3); umin = movmin(u,3); minmaxDiff = umax - umin; end function minmaxDiff = minmaxDiffNeighborhood2d(grid, u) uMatrix = grid.funcToMatrix(grid,u); umax = max(movmax(uMatrix,3,1),movmax(uMatrix,3,2)); umin = min(movmin(uMatrix,3,1),movmin(uMatrix,3,2)); minmaxDiff = umax - umin; minmaxDiff = reshape(minmaxDiff',g.N(),1); end end end