comparison +scheme/AdvectionRV1D.m @ 1012:1e437c9e5132 feature/advectionRV

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 f753bada1a46
children
comparison
equal deleted inserted replaced
1011:e0560bc4fb7d 1012:1e437c9e5132
6 D % Non-stabalized scheme operator 6 D % Non-stabalized scheme operator
7 H % Discrete norm 7 H % Discrete norm
8 Hi % Norm inverse 8 Hi % Norm inverse
9 e_l 9 e_l
10 e_r 10 e_r
11
12 D2_visc % Artificial viscosity operator
11 end 13 end
12 14
13 methods 15 methods
14 function obj = AdvectionRV1D(grid, operator_type, order, waveSpeed) 16 function obj = AdvectionRV1D(grid, operator_type, order, waveSpeed)
15 m = grid.size(); 17 m = grid.size();
17 switch operator_type 19 switch operator_type
18 case 'upwind+' 20 case 'upwind+'
19 ops = sbp.D1Upwind(m, lim, order); 21 ops = sbp.D1Upwind(m, lim, order);
20 D1 = (ops.Dp + ops.Dm)/2; 22 D1 = (ops.Dp + ops.Dm)/2;
21 B = ops.e_r*ops.e_r' - ops.e_l*ops.e_l'; 23 B = ops.e_r*ops.e_r' - ops.e_l*ops.e_l';
22 D2 = @(viscosity) ops.Dm*spdiag(viscosity)*ops.Dp-ops.HI*(B*spdiag(viscosity)*ops.Dp); 24 obj.D2_visc = @(viscosity) ops.Dm*spdiag(viscosity)*ops.Dp-ops.HI*(B*spdiag(viscosity)*ops.Dp);
25 % max(abs()) or just abs()?
23 DissipationOp = spdiag(abs(waveSpeed))*(ops.Dp-ops.Dm)/2; 26 DissipationOp = spdiag(abs(waveSpeed))*(ops.Dp-ops.Dm)/2;
24 otherwise 27 otherwise
25 error('Other operator types not yet supported', operator_type); 28 error('Other operator types not yet supported', operator_type);
26 end 29 end
27 % max(abs()) or just abs()? 30
28 obj.D = @(viscosity) (-D1 + D2(viscosity) + DissipationOp); 31 obj.D = -D1 + DissipationOp;
29 32
30 obj.grid = grid; 33 obj.grid = grid;
31 obj.order = order; 34 obj.order = order;
32 35
33 obj.H = ops.H; 36 obj.H = ops.H;