Mercurial > repos > public > sbplib

changeset 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
files +scheme/Burgers1D.m
diffstat 1 files changed, 19 insertions(+), 3 deletions(-) [+]
line wrap: on
line diff
--- a/+scheme/Burgers1D.m	Fri Sep 21 15:33:15 2018 +0200
+++ b/+scheme/Burgers1D.m	Thu Sep 27 09:30:21 2018 +0200
@@ -34,12 +34,28 @@
                     end
                     d_l = ops.d1_l';
                     d_r = ops.d1_r';
-                case 'upwind'
+                case 'upwind-'
                     ops = sbp.D1Upwind(m, lim, order);
                     D1 = (ops.Dp + ops.Dm)/2;
                     D2 = @(eps) ops.Dp*spdiag(eps)*ops.Dm;
-                    %D2 = @(eps) ops.Dm*diag(eps)*ops.Dp;
-                    %D2 = @(eps) (ops.Dp*diag(eps)*ops.Dm + ops.Dm*diag(eps)*ops.Dp)/2;
+                    if (strcmp(dissipation,'on'))
+                        DissipationOp = (ops.Dp-ops.Dm)/2;
+                    end
+                    d_l = ops.e_l'*ops.Dm;
+                    d_r = ops.e_r'*ops.Dm;
+                case 'upwind+'
+                    ops = sbp.D1Upwind(m, lim, order);
+                    D1 = (ops.Dp + ops.Dm)/2;
+                    D2 = @(eps) ops.Dm*spdiag(eps)*ops.Dp;
+                    if (strcmp(dissipation,'on'))
+                        DissipationOp = (ops.Dp-ops.Dm)/2;
+                    end
+                    d_l = ops.e_l'*ops.Dp;
+                    d_r = ops.e_r'*ops.Dp;
+                case 'upwind+-'
+                    ops = sbp.D1Upwind(m, lim, order);
+                    D1 = (ops.Dp + ops.Dm)/2;
+                    D2 = @(eps) (ops.Dp*spdiag(eps)*ops.Dm + ops.Dm*spdiag(eps)*ops.Dp)/2;
                     if (strcmp(dissipation,'on'))
                         DissipationOp = (ops.Dp-ops.Dm)/2;
                     end