diff +scheme/Euler1d.m @ 57:9a647dcccbdd

Added pausing option to noname.animate. Added characteristic bc to Euler1d.
author Jonatan Werpers <jonatan@werpers.com>
date Thu, 12 Nov 2015 17:31:33 -0800
parents 2194cd385419
children 24103284e09d
line wrap: on
line diff
--- a/+scheme/Euler1d.m	Thu Nov 12 12:03:45 2015 -0800
+++ b/+scheme/Euler1d.m	Thu Nov 12 17:31:33 2015 -0800
@@ -166,12 +166,54 @@
                     closure = obj.boundary_condition_inflow(boundary,varargin{:});
                 case 'outflow'
                     closure = obj.boundary_condition_outflow(boundary,varargin{:});
+                case 'char'
+                    closure = obj.boundary_condition_char(boundary,varargin{:});
                 otherwise
                     error('Unsupported bc type: %s', type);
             end
 
         end
 
+        function closure = boundary_condition_char(obj,boundary,w_data)
+            [e_s,e_S,s] = obj.get_boundary_ops(boundary);
+
+            function o = closure_fun(q,t)
+                q_s = e_S'*q;
+                rho = q_s(1);
+                u = q_s(2)/rho;
+                e = q_s(3);
+
+                c = obj.c(q_s);
+
+                Lambda = [u, u+c, u-c];
+
+                p_in = find(s*Lambda < 0);
+                p_ot = find(s*Lambda >= 0);
+                p = [p_in p_ot];
+                pt(p) = 1:length(p);
+
+                T = obj.T(q_s);
+
+                tau1 = -2*diag(Lambda(p_in));
+                tau2 = zeros(length(p_ot),length(p_in)); % Penalty only on ingoing char.
+
+                tauHat = [tau1; tau2];
+
+                tau = -s*e_S*sparse(T*tauHat(pt,:));
+
+                w_s = inv(T)*q_s;
+                w_in = w_s(p_in);
+
+                w_s_data = w_data(t);
+                w_in_data = w_s_data(p_in);
+
+                o = obj.Hi * tau * (w_in - w_in_data);
+            end
+
+            closure = @closure_fun;
+
+        end
+
         function closure = boundary_condition_inflow(obj, boundary, p_data, v_data)
             [e_s,e_S,s] = obj.get_boundary_ops(boundary);
 
@@ -222,6 +264,7 @@
 
 
                 tauHat = [tau1; tau2];
+
                 tau = -s*e_S*sparse(T*tauHat(pt,:));
 
                 w_s = inv(T)*q_s;