Mercurial > repos > public > sbplib

comparison +rv/ResidualViscosity.m @ 1012:1e437c9e5132 feature/advectionRV

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Create residual viscosity package +rv and generalize the ResidualViscosity class - Generalize residual viscosity, by passing user-defined flux and calculating the time derivative outside of the update. - Create separate RungekuttaRV specifically using interior RV updates - Separate the artifical dissipation operator from the scheme AdvectionRV1D so that the same scheme can be reused for creating the diff op used by the ResidualViscosity class
author Vidar Stiernström <vidar.stiernstrom@it.uu.se>
date Wed, 05 Dec 2018 13:44:10 +0100
parents
children e547794a9407
comparison
equal deleted inserted replaced
1011:e0560bc4fb7d 1012:1e437c9e5132
1 % class describing the viscosity
2 classdef ResidualViscosity < handle
3 properties
4 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?
6 Cmax % Constant controling magnitude of upwind dissipation
7 Cres % Constant controling magnitude residual dissipation
8 h % Length scale used for scaling the viscosity.
9 viscosity % Stores the computed viscosity.
10
11 % 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 grad_f % Stores the latest approximated gradient of the flux
14 residual % Stores the computed residual
15 end
16
17 methods
18 % TODO: - Consider passing residual normalization as a function handle.
19 % or choosing a type of normalization on construction.
20 % Could for example be 1, norm((v-mean(v),inf) or normInfNeighborhood(v)
21 % but working
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;
26 obj.waveSpeed = waveSpeed;
27 obj.h = h;
28 obj.Cmax = Cmax;
29 obj.Cres = Cres;
30 obj.viscosity = zeros(N,1);
31 obj.u_t = zeros(N,1);
32 obj.grad_f = zeros(N,1);
33 obj.residual = zeros(N,1);
34 end
35
36 function obj = update(obj, v, dvdt)
37 obj.u_t = dvdt;
38 obj.grad_f = obj.D(v);
39 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));
41 end
42
43 function [residual, u_t, grad_f] = getResidual(obj)
44 residual = obj.residual;
45 u_t = obj.u_t;
46 grad_f = obj.grad_f;
47 end
48
49 function viscosity = getViscosity(obj)
50 viscosity = obj.viscosity;
51 end
52 end
53 % Remove or fix. Should be able to handle values close to zero. Should work in 2d and 3d.
54 methods (Static)
55 function R_norm = normInfNeighborhood(v)
56 n = length(v);
57 R_norm = zeros(n,1);
58 R_norm(1,1) = norm(v(1:3), inf);
59 R_norm(n,1) = norm(v(n-3:n), inf);
60 for i = 2:n-1
61 R_norm(i,1) = norm(v(i-1:i+1), inf);
62 end
63 end
64 end
65 end