Mercurial > repos > public > sbplib

comparison +rv/ResidualViscosity.m @ 1016:4b42999874c0 feature/advectionRV

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Add lower level for boot-strapping to RungeKuttaExteriorRV - Add a lower level to RungeKuttaExteriorRV for which bootstrapping starts, e.g start bootstrapping from time level 3 using a 3rd order BDF - Clean up ResidualViscosity
author Vidar Stiernström <vidar.stiernstrom@it.uu.se>
date Fri, 07 Dec 2018 13:11:53 +0100
parents 9b7fcd5e4480
children 2d7c1333bd6c
comparison
equal deleted inserted replaced
1015:9b7fcd5e4480 1016:4b42999874c0
1 % class describing the viscosity
2 classdef ResidualViscosity < handle 1 classdef ResidualViscosity < handle
3 properties 2 properties
4 D % Diff op approximating the gradient of the flux f(u) 3 D % Diff op approximating the gradient of the flux f(u)
5 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?
6 Cmax % Constant controling magnitude of upwind dissipation 5 Cmax % Constant controlling relative amount of upwind dissipation
7 Cres % Constant controling magnitude residual dissipation 6 Cres % Constant controling relative amount of upwind dissipation
8 h % Length scale used for scaling the viscosity. 7 h % Length scale used for scaling the viscosity. Typically grid spacing.
9 viscosity % Stores the computed viscosity. 8 viscosity % Stores the computed viscosity.
9 normalization % Function used to normalize the residual such that it is amplified in the
10 % shocks.
10 11
11 % Convenice (for verification and plotting) TBD: Decide on if it should be kept. 12 % Convenice (for verification and plotting) TBD: Decide on if it should be kept.
12 u_t % Stores the latest approximated time derivative of the solution. 13 u_t % Stores the latest approximated time derivative of the solution.
13 grad_f % Stores the latest approximated gradient of the flux 14 grad_f % Stores the latest approximated gradient of the flux
14 residual % Stores the computed residual 15 residual % Stores the computed residual
15 end 16 end
16 17
17 methods 18 methods
18 % TODO: - Consider passing residual normalization as a function handle. 19 % TBD: Decide on how to treat waveSpeed. It would be nice to just pass a constant value without
19 % or choosing a type of normalization on construction. 20 % wrapping it in a function.
20 % Could for example be 1, norm((v-mean(v),inf) or normInfNeighborhood(v) 21 function obj = ResidualViscosity(D, waveSpeed, Cmax, Cres, h, N, normalization)
21 % but working 22 default_arg('normalization',@(v)norm(v-mean(v),inf));
22 % - Decide on how to treat waveSpeed. It would be nice to just pass a constant value without
23 % wrapping it in a function.
24 function obj = ResidualViscosity(D, waveSpeed, Cmax, Cres, h, N)
25 obj.D = D; 23 obj.D = D;
26 obj.waveSpeed = waveSpeed; 24 obj.waveSpeed = waveSpeed;
27 obj.h = h; 25 obj.h = h;
28 obj.Cmax = Cmax; 26 obj.Cmax = Cmax;
29 obj.Cres = Cres; 27 obj.Cres = Cres;
30 obj.viscosity = zeros(N,1); 28 obj.viscosity = zeros(N,1);
31 obj.u_t = zeros(N,1); 29 obj.u_t = zeros(N,1);
32 obj.grad_f = zeros(N,1); 30 obj.grad_f = zeros(N,1);
33 obj.residual = zeros(N,1); 31 obj.residual = zeros(N,1);
32 obj.normalization = normalization;
34 end 33 end
35 34
36 function obj = update(obj, v, dvdt) 35 function obj = update(obj, v, dvdt)
37 obj.u_t = dvdt; 36 obj.u_t = dvdt;
38 obj.grad_f = obj.D(v); 37 obj.grad_f = obj.D(v);
39 obj.residual = obj.u_t + obj.grad_f; 38 obj.residual = obj.u_t + obj.grad_f;
40 %obj.viscosity = min(obj.Cmax*obj.h*abs(obj.waveSpeed(v)), obj.Cres*obj.h^2*abs(obj.residual)/norm(v-mean(v),inf)); 39 obj.viscosity = min(obj.Cmax*obj.h*abs(obj.waveSpeed(v)), obj.Cres*obj.h^2*abs(obj.residual)/obj.normalization(v));
41 obj.viscosity = obj.smoothen(obj.Cres*obj.h^2*abs(obj.residual)/norm(v-mean(v),inf));
42
43 end
44
45 function smoothendVector = smoothen(obj, vector)
46 smoothendVector = vector;
47 for i = 2:length(vector)-1
48 smoothendVector(i) = (1/6)*(vector(i-1) + 4*vector(i) + vector(i+1));
49 end
50 end 40 end
51 41
52 function [residual, u_t, grad_f] = getResidual(obj) 42 function [residual, u_t, grad_f] = getResidual(obj)
53 residual = obj.residual; 43 residual = obj.residual;
54 u_t = obj.u_t; 44 u_t = obj.u_t;
57 47
58 function viscosity = getViscosity(obj) 48 function viscosity = getViscosity(obj)
59 viscosity = obj.viscosity; 49 viscosity = obj.viscosity;
60 end 50 end
61 end 51 end
62 % Remove or fix. Should be able to handle values close to zero. Should work in 2d and 3d.
63 methods (Static)
64 function R_norm = normInfNeighborhood(v)
65 n = length(v);
66 R_norm = zeros(n,1);
67 R_norm(1,1) = norm(v(1:3), inf);
68 R_norm(n,1) = norm(v(n-3:n), inf);
69 for i = 2:n-1
70 R_norm(i,1) = norm(v(i-1:i+1), inf);
71 end
72 end
73 end
74 end 52 end