Mercurial > repos > public > sbplib
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; |