Mercurial > repos > public > sbplib

comparison +rv/ResidualViscosity.m @ 1017:2d7c1333bd6c feature/advectionRV

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Add support for using the ODE to approximate the time derivative in the residual
author Vidar Stiernström <vidar.stiernstrom@it.uu.se>
date Tue, 11 Dec 2018 16:29:21 +0100
parents 4b42999874c0
children 2ef20d00b386
comparison
equal deleted inserted replaced
1016:4b42999874c0 1017:2d7c1333bd6c
1 classdef ResidualViscosity < handle 1 classdef ResidualViscosity < handle
2 properties 2 properties
3 D % Diff op approximating the gradient of the flux f(u) 3 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? 4 waveSpeed % Wave speed at each grid point, e.g f'(u). %TBD: Better naming?
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 viscosity % Stores the computed viscosity.
9 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
10 % shocks. 9 % shocks.
11
12 % Convenice (for verification and plotting) TBD: Decide on if it should be kept.
13 u_t % Stores the latest approximated time derivative of the solution.
14 grad_f % Stores the latest approximated gradient of the flux
15 residual % Stores the computed residual
16 end 10 end
17 11
18 methods 12 methods
19 % 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
20 % wrapping it in a function. 14 % wrapping it in a function.
21 function obj = ResidualViscosity(D, waveSpeed, Cmax, Cres, h, N, normalization) 15 function obj = ResidualViscosity(Df, waveSpeed, Cmax, Cres, h, normalization)
22 default_arg('normalization',@(v)norm(v-mean(v),inf)); 16 default_arg('normalization',@(v)norm(v-mean(v),inf));
23 obj.D = D; 17 obj.Df = Df;
24 obj.waveSpeed = waveSpeed; 18 obj.waveSpeed = waveSpeed;
25 obj.h = h; 19 obj.h = h;
26 obj.Cmax = Cmax; 20 obj.Cmax = Cmax;
27 obj.Cres = Cres; 21 obj.Cres = Cres;
28 obj.viscosity = zeros(N,1);
29 obj.u_t = zeros(N,1);
30 obj.grad_f = zeros(N,1);
31 obj.residual = zeros(N,1);
32 obj.normalization = normalization; 22 obj.normalization = normalization;
33 end 23 end
34 24
35 function obj = update(obj, v, dvdt) 25 function [viscosity, Df] = evaluate(obj, v, dvdt)
36 obj.u_t = dvdt; 26 Df = obj.Df(v);
37 obj.grad_f = obj.D(v); 27 viscosity = min(obj.Cmax*obj.h*abs(obj.waveSpeed(v)), obj.Cres*obj.h^2*abs(dvdt + Df)/obj.normalization(v));
38 obj.residual = obj.u_t + obj.grad_f;
39 obj.viscosity = min(obj.Cmax*obj.h*abs(obj.waveSpeed(v)), obj.Cres*obj.h^2*abs(obj.residual)/obj.normalization(v));
40 end
41
42 function [residual, u_t, grad_f] = getResidual(obj)
43 residual = obj.residual;
44 u_t = obj.u_t;
45 grad_f = obj.grad_f;
46 end
47
48 function viscosity = getViscosity(obj)
49 viscosity = obj.viscosity;
50 end 28 end
51 end 29 end
52 end 30 end