Mercurial > repos > public > sbplib

changeset 1019:f029b97dbc72 feature/advectionRV

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Support upwind opSet in Utux2d
author Vidar Stiernström <vidar.stiernstrom@it.uu.se>
date Fri, 14 Dec 2018 18:30:27 +0100
parents eca4d9df9365
children 5359a61cb4d9
files +scheme/Utux2d.m
diffstat 1 files changed, 18 insertions(+), 5 deletions(-) [+]
line wrap: on
line diff
--- a/+scheme/Utux2d.m	Tue Dec 11 17:31:28 2018 +0100
+++ b/+scheme/Utux2d.m	Fri Dec 14 18:30:27 2018 +0100
@@ -16,6 +16,9 @@
         % Derivatives
         Dx, Dy
 
+        % Dissipation Operators
+        DissOpx, DissOpy
+
         % Boundary operators
         e_w, e_e, e_s, e_n
 
@@ -69,10 +72,21 @@
             obj.Hyi = kron(Ix,Hyi);
 
             % Derivatives
-            Dx = ops_x.D1;
-            Dy = ops_y.D1;
-            obj.Dx = kron(Dx,Iy);
-            obj.Dy = kron(Ix,Dy);
+            if (isequal(opSet,@sbp.D1Upwind))
+                Dx = (ops_x.Dp + ops_x.Dm)/2;
+                Dy = (ops_y.Dp + ops_y.Dm)/2;
+                DissOpx = (ops_x.Dm - ops_x.Dp)/2;
+                DissOpy = (ops_y.Dm - ops_y.Dp)/2;
+                obj.Dx = kron(Dx,Iy);
+                obj.Dy = kron(Ix,Dy);
+                obj.DissOpx = kron(DissOpx,Iy);
+                obj.DissOpy = kron(Ix,DissOpy);
+            else
+                Dx = ops_x.D1;
+                Dy = ops_y.D1;
+                obj.Dx = kron(Dx,Iy);
+                obj.Dy = kron(Ix,Dy);
+            end
 
             % Boundary operators
             obj.e_w = kr(ops_x.e_l, Iy);
@@ -85,7 +99,6 @@
             obj.order = order;
             obj.a = a;
             obj.D = -(a{1}*obj.Dx + a{2}*obj.Dy);
-
         end
         % Closure functions return the opertors applied to the own domain to close the boundary
         % Penalty functions return the opertors to force the solution. In the case of an interface it returns the operator applied to the other doamin.