changeset 93:5c41941b9e5e

Euler1D: added function to detect type of flow at a baoundary.
author Jonatan Werpers <jonatan@werpers.com>
date Tue, 01 Dec 2015 16:21:19 +0100
parents ed7c7d651428
children a38c243991d0
files +scheme/Euler1d.m
diffstat 1 files changed, 42 insertions(+), 0 deletions(-) [+]
line wrap: on
line diff
diff -r ed7c7d651428 -r 5c41941b9e5e +scheme/Euler1d.m
--- a/+scheme/Euler1d.m	Mon Nov 30 17:58:33 2015 +0100
+++ b/+scheme/Euler1d.m	Tue Dec 01 16:21:19 2015 +0100
@@ -17,6 +17,14 @@
 
     end
 
+    properties (Constant)
+        SUBSONIC_INFLOW = 1;
+        SUBSONIC_OUTFLOW = -1;
+        NO_FLOW = 0;
+        SUPERSONIC_INFLOW = 2;
+        SUPERSONIC_OUTFLOW = -2;
+    end
+
     methods
         function obj = Euler1d(m,xlim,order,gama,opsGen,do_upwind)
             default_arg('opsGen',@sbp.Ordinary);
@@ -155,6 +163,40 @@
             % Devide columns by stuff to get rid of extra comp?
         end
 
+        function fs = flowStateL(obj, q)
+            q_l = obj.e_L'*q;
+            c = obj.c(q);
+            v = q_l(2,:)/q_l(1,:);
+
+            if v > c
+                fs = scheme.Euler1d.SUPERSONIC_INFLOW;
+            elseif v > 0
+                fs = scheme.Euler1d.SUBSONIC_INFLOW;
+            elseif v > -c
+                fs = scheme.Euler1d.SUBSONIC_OUTFLOW;
+            else
+                fs = scheme.Euler1d.SUPERSONIC_INFLOW;
+            end
+        end
+
+        % returns positiv values for inlfow, negative for outflow.
+        %  +-1 for subsonic
+        function fs = flowStateR(obj, q)
+            q_r = obj.e_R'*q;
+            c = obj.c(q);
+            v = q_r(2,:)/q_r(1,:);
+
+            if v < -c
+                fs = scheme.Euler1d.SUPERSONIC_INFLOW;
+            elseif v < 0
+                fs = scheme.Euler1d.SUBSONIC_INFLOW;
+            elseif v < c
+                fs = scheme.Euler1d.SUBSONIC_OUTFLOW;
+            else
+                fs = scheme.Euler1d.SUPERSONIC_INFLOW;
+            end
+        end
+
         % Enforces the boundary conditions
         %  w+ = R*w- + g(t)
         function closure = boundary_condition(obj,boundary, type, varargin)