Mercurial > repos > public > sbplib

comparison +scheme/Burgers1D.m @ 852:fbb8be3177c8 feature/burgers1d

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Fix bug in SAT terms for upwind operators - Set the appropriate boundary derivatives in respect to choice of upwind discretization of D2 - Distinguish between upwind discretizations e.g upwind+ vs upwind-
author Vidar Stiernström <vidar.stiernstrom@it.uu.se>
date Thu, 27 Sep 2018 09:30:21 +0200
parents ab2e5a24ddde
children cda996e64925
comparison
equal deleted inserted replaced
851:ab2e5a24ddde 852:fbb8be3177c8
32 if (strcmp(dissipation,'on')) 32 if (strcmp(dissipation,'on'))
33 DissipationOp = -1*sbp.dissipationOperator(m, order, ops.HI); 33 DissipationOp = -1*sbp.dissipationOperator(m, order, ops.HI);
34 end 34 end
35 d_l = ops.d1_l'; 35 d_l = ops.d1_l';
36 d_r = ops.d1_r'; 36 d_r = ops.d1_r';
37 case 'upwind' 37 case 'upwind-'
38 ops = sbp.D1Upwind(m, lim, order); 38 ops = sbp.D1Upwind(m, lim, order);
39 D1 = (ops.Dp + ops.Dm)/2; 39 D1 = (ops.Dp + ops.Dm)/2;
40 D2 = @(eps) ops.Dp*spdiag(eps)*ops.Dm; 40 D2 = @(eps) ops.Dp*spdiag(eps)*ops.Dm;
41 %D2 = @(eps) ops.Dm*diag(eps)*ops.Dp; 41 if (strcmp(dissipation,'on'))
42 %D2 = @(eps) (ops.Dp*diag(eps)*ops.Dm + ops.Dm*diag(eps)*ops.Dp)/2; 42 DissipationOp = (ops.Dp-ops.Dm)/2;
43 end
44 d_l = ops.e_l'*ops.Dm;
45 d_r = ops.e_r'*ops.Dm;
46 case 'upwind+'
47 ops = sbp.D1Upwind(m, lim, order);
48 D1 = (ops.Dp + ops.Dm)/2;
49 D2 = @(eps) ops.Dm*spdiag(eps)*ops.Dp;
50 if (strcmp(dissipation,'on'))
51 DissipationOp = (ops.Dp-ops.Dm)/2;
52 end
53 d_l = ops.e_l'*ops.Dp;
54 d_r = ops.e_r'*ops.Dp;
55 case 'upwind+-'
56 ops = sbp.D1Upwind(m, lim, order);
57 D1 = (ops.Dp + ops.Dm)/2;
58 D2 = @(eps) (ops.Dp*spdiag(eps)*ops.Dm + ops.Dm*spdiag(eps)*ops.Dp)/2;
43 if (strcmp(dissipation,'on')) 59 if (strcmp(dissipation,'on'))
44 DissipationOp = (ops.Dp-ops.Dm)/2; 60 DissipationOp = (ops.Dp-ops.Dm)/2;
45 end 61 end
46 d_l = ops.e_l'*D1; 62 d_l = ops.e_l'*D1;
47 d_r = ops.e_r'*D1; 63 d_r = ops.e_r'*D1;