Mercurial > repos > public > sbplib

changeset 77:0b07ff8a0a12

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Fixed bugs and added outflow bc for rho. Added missing abs() added missing 1/2.
author Jonatan Werpers <jonatan@werpers.com>
date Wed, 25 Nov 2015 15:19:55 +0100
parents 5c569cbef49e
children 80948a4084f3
files +scheme/Euler1d.m
diffstat 1 files changed, 33 insertions(+), 15 deletions(-) [+]
line wrap: on
line diff
diff -r 5c569cbef49e -r 0b07ff8a0a12 +scheme/Euler1d.m
--- a/+scheme/Euler1d.m	Mon Nov 23 17:49:34 2015 +0100
+++ b/+scheme/Euler1d.m	Wed Nov 25 15:19:55 2015 +0100
@@ -170,6 +170,8 @@
                     closure = obj.boundary_condition_inflow(boundary,varargin{:});
                 case 'outflow'
                     closure = obj.boundary_condition_outflow(boundary,varargin{:});
+                    case 'outflow_rho'
+                    closure = obj.boundary_condition_outflow_rho(boundary,varargin{:});
                 case 'char'
                     closure = obj.boundary_condition_char(boundary,varargin{:});
                 otherwise
@@ -207,29 +209,27 @@
 
                 % Convert bc from ordinary form to characteristic form.
                 %   Lq = g  => w_in = Rw_ot + g_tilde
-                L = L_fun(rho,u,e);
-                g = g_fun(t);
 
                 Lambda = obj.Lambda(q_s);
 
-                R =-inv(L*Tin)*(L*Tot);
-                g_tilde = inv(L*Tin)*g;
-
                 % Setup the penalty parameter
-                tau1 = -2*Lambda(p_in,p_in);
+                tau1 = -2*abs(Lambda(p_in,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,:));
 
-                % Calculate the numerical caracteristics and apply the penalty
-                w_s = inv(T)*q_s;
-                % w_s = T\q_s;
-                % w_s = Tinv * q_s; % Med analytisk matris
-                w_in = w_s(p_in);
-                w_ot = w_s(p_ot);
+                L = L_fun(rho,u,e);
+                g = g_fun(t);
 
-                o = obj.Hi * tau * (w_in - R*w_ot - g_tilde);
+                % printExpr('s')
+                % penalty = tauHat(pt,:)*inv(L*Tin)*(L*q_s - g);
+                % tauHatPt = tauHat(pt,:);
+                % display(tauHatPt);
+                % display(penalty);
+                % pause
+
+                o = 1/2*obj.Hi * tau * inv(L*Tin)*(L*q_s - g);
             end
             closure = @closure_fun;
         end
@@ -254,7 +254,7 @@
 
                 T = obj.T(q_s);
 
-                tau1 = -2*diag(Lambda(p_in));
+                tau1 = -2*diag(abs(Lambda(p_in)));
                 tau2 = zeros(length(p_ot),length(p_in)); % Penalty only on ingoing char.
 
                 tauHat = [tau1; tau2];
@@ -267,7 +267,7 @@
                 w_s_data = w_data(t);
                 w_in_data = w_s_data(p_in);
 
-                o = obj.Hi * tau * (w_in - w_in_data);
+                o = 1/2*obj.Hi * tau * (w_in - w_in_data);
             end
 
             closure = @closure_fun;
@@ -317,6 +317,24 @@
 
         end
 
+        function closure = boundary_condition_outflow_rho(obj, boundary, rho_data)
+            [~,~,s] = obj.get_boundary_ops(boundary);
+
+             switch s
+                case -1
+                    p_in = 2;
+                case 1
+                    p_in = 3;
+            end
+
+            L = @(~,~,~)[1 0 0];
+            g = @(t)[rho_data(t)];
+
+
+            closure = boundary_condition_L(obj, boundary, L, g, p_in);
+
+        end
+
         % Set wall boundary condition v = 0.
         function closure = boundary_condition_wall(obj,boundary)
             [e_s,e_S,s] = obj.get_boundary_ops(boundary);