Mercurial > repos > public > sbplib

changeset 61:183d9349b4c1

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Refactored euler1d to have real methods.
author Jonatan Werpers <jonatan@werpers.com>
date Sat, 14 Nov 2015 19:26:30 -0800
parents e707cd9e6d0a
children 00344139cd7d
files +scheme/Euler1d.m
diffstat 1 files changed, 51 insertions(+), 54 deletions(-) [+]
line wrap: on
line diff
diff -r e707cd9e6d0a -r 183d9349b4c1 +scheme/Euler1d.m
--- a/+scheme/Euler1d.m	Sat Nov 14 15:43:01 2015 -0800
+++ b/+scheme/Euler1d.m	Sat Nov 14 19:26:30 2015 -0800
@@ -11,12 +11,6 @@
         M % Derivative norm
         gamma
 
-        T
-        p
-        c
-        Lambda
-
-
         H % Discrete norm
         Hi
         e_l, e_r, e_L, e_R;
@@ -76,64 +70,20 @@
             % F=[rho*u, rho*u^2+p, (e+p)*u] ^T
             % p=(gamma-1)*(e-rho*u^2/2);
 
+
             %Solving on form q_t + F_x = 0
-            function o = F(Q)
-                % Flux: f = [q2; q2.^2/q1 + p(q); (q3+p(q))*q2/q1];
-                o = [Q(2,:); Q(2,:).^2./Q(1,:) + p(Q); (Q(3,:)+p(Q)).*Q(2,:)./Q(1,:)];
-            end
-
-            % Equation of state
-            function o = p(Q)
-                % Pressure p = (gamma-1)*(q3-q2.^2/q1/2)
-                o = (gamma-1)*( Q(3,:)-1/2*Q(2,:).^2./Q(1,:) );
-            end
-            obj.p = @p;
-
-            function o = c(Q)
-                % Speed of light c = sqrt(obj.gamma*p/rho);
-                o = sqrt(gamma*p(Q)./Q(1,:));
-            end
-            obj.c = @c;
-
-            function o = Lambda(q)
-                u = q(2)/q(1);
-                c = obj.c(q);
-                L = [u, u+c, u-c];
-                o = diag(L);
-            end
-            obj.Lambda = @Lambda;
-
-            function o = T(q)
-                % T is the transformation such that A = T*Lambda*inv(T)
-                % where Lambda=diag(u, u+c, u-c)
-                rho = q(1);
-                u = q(2)/q(1);
-                e = q(3);
-
-                c = sqrt(gamma*p(q)/rho);
-
-                sqrt2gamm = sqrt(2*(gamma-1));
-
-                o = [
-                     sqrt2gamm*rho      , rho                               , rho                               ;
-                     sqrt2gamm*rho*u    , rho*(u+c)                         , rho*(u-c)                         ;
-                     sqrt2gamm*rho*u^2/2, e+(gamma-1)*(e-rho*u^2/2)+rho*u*c , e+(gamma-1)*(e-rho*u^2/2)-rho*u*c ;
-                ];
-                % Devide columns by stuff to get rid of extra comp?
-            end
-            obj.T = @T;
 
             function o = Fx(q)
                 Q = reshape(q,3,m);
-                o = reshape(F(Q),3*m,1);
+                o = reshape(obj.F(Q),3*m,1);
                 o = D1*o;
             end
 
             function o = Fx_disp(q)
                 Q = reshape(q,3,m);
-                f = reshape(F(Q),3*m,1);
+                f = reshape(obj.F(Q),3*m,1);
 
-                c = c(Q);
+                c = obj.c(Q);
                 lambda_max = c+abs(Q(2,:)./Q(1,:));
                 % lambda_max = max(lambda_max);
 
@@ -157,6 +107,53 @@
             obj.gamma = gamma;
         end
 
+        % Flux function
+        function o = F(obj, Q)
+            % Flux: f = [q2; q2.^2/q1 + p(q); (q3+p(q))*q2/q1];
+            p = obj.p(Q);
+            o = [Q(2,:); Q(2,:).^2./Q(1,:) + p; (Q(3,:)+p).*Q(2,:)./Q(1,:)];
+        end
+
+        % Equation of state
+        function o = p(obj, Q)
+            % Pressure p = (gamma-1)*(q3-q2.^2/q1/2)
+            o = (obj.gamma-1)*( Q(3,:)-1/2*Q(2,:).^2./Q(1,:) );
+        end
+
+        % Speed of sound
+        function o = c(obj, Q)
+            % Speed of light c = sqrt(obj.gamma*p/rho);
+            o = sqrt(obj.gamma*obj.p(Q)./Q(1,:));
+        end
+
+        % Eigen value matrix
+        function o = Lambda(q)
+            u = q(2)/q(1);
+            c = obj.c(q);
+            L = [u, u+c, u-c];
+            o = diag(L);
+        end
+
+        % Diagonalization transformation
+        function o = T(obj, q)
+            % T is the transformation such that A = T*Lambda*inv(T)
+            % where Lambda=diag(u, u+c, u-c)
+            rho = q(1);
+            u = q(2)/q(1);
+            e = q(3);
+            gamma = obj.gamma;
+
+            c = sqrt(gamma*obj.p(q)/rho);
+
+            sqrt2gamm = sqrt(2*(gamma-1));
+
+            o = [
+                 sqrt2gamm*rho      , rho                               , rho                               ;
+                 sqrt2gamm*rho*u    , rho*(u+c)                         , rho*(u-c)                         ;
+                 sqrt2gamm*rho*u^2/2, e+(gamma-1)*(e-rho*u^2/2)+rho*u*c , e+(gamma-1)*(e-rho*u^2/2)-rho*u*c ;
+            ];
+            % Devide columns by stuff to get rid of extra comp?
+        end
 
         % Enforces the boundary conditions
         %  w+ = R*w- + g(t)