Mercurial > repos > public > sbplib

changeset 1148:0a5503a08a36 feature/rv

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Implement the correct Dirichlet conditions for Burgers1d
author Vidar Stiernström <vidar.stiernstrom@it.uu.se>
date Thu, 24 Jan 2019 09:01:10 +0100
parents c322a9454d75
children 1fe48cbd379a
files +scheme/Burgers1d.m
diffstat 1 files changed, 23 insertions(+), 24 deletions(-) [+]
line wrap: on
line diff
--- a/+scheme/Burgers1d.m	Sat Jan 19 17:22:09 2019 +0100
+++ b/+scheme/Burgers1d.m	Thu Jan 24 09:01:10 2019 +0100
@@ -30,30 +30,23 @@
             if (isequal(opSet, @sbp.D1Upwind))
                 obj.D1 = (ops.Dp + ops.Dm)/2;
                 DissOp = (ops.Dm - ops.Dp)/2;
+                switch pde_form
+                    case 'skew-symmetric'
+                        obj.D = @(v) -(1/3*obj.D1*v.*v + (1/3*spdiag(v)*obj.D1  + fluxSplitting(v)*DissOp)*v);
+                    case 'conservative'
+                        obj.D = @(v) -(1/2*obj.D1*v.*v + fluxSplitting(v)*DissOp*v);
+                end
             else 
                 obj.D1 = ops.D1;
+                switch pde_form
+                    case 'skew-symmetric'
+                        obj.D = @(v) -(1/3*obj.D1*v.*v + 1/3*spdiag(v)*obj.D1*v);
+                    case 'conservative'
+                        obj.D = @(v) -1/2*obj.D1*v.*v;
+                end
             end
-            
-            switch pde_form
-                case 'skew-symmetric'
-                    if (isequal(opSet, @sbp.D1Upwind))
-                        D = @(v) -(1/3*obj.D1*v.*v + (1/3*spdiag(v)*obj.D1  + fluxSplitting(v)*DissOp)*v);
-                    else
-                        D = @(v) -(1/3*obj.D1*v.*v + 1/3*spdiag(v)*obj.D1*v);
-                    end
-                case 'conservative'
-                    if (isequal(opSet, @sbp.D1Upwind))
-                        D = @(v) -(1/2*obj.D1*v.*v + fluxSplitting(v)*DissOp*v);
-                    else
-                        D = @(v) -(1/2*obj.D1*v.*v);
-                    end
-                otherwise
-                    error('Not supported', pde_form);
-            end
-
             obj.grid = g;
 
-            obj.D = D;
             obj.H =  ops.H;
             obj.Hi = ops.HI;
 
@@ -70,14 +63,20 @@
         %       boundary            is a string specifying the boundary e.g. 'l','r' or 'e','w','n','s'.
         %       type                is a string specifying the type of boundary condition if there are several.
         function [closure, penalty] = boundary_condition(obj, boundary, type)
-            default_arg('type','robin');
+            default_arg('type','dirichlet');
             [e, index, s] = obj.get_boundary_ops(boundary);
             switch type
-                % Stable dirhclet-like boundary conditions ((u+-abs(u))*u/3)) with +- at left/right boundary
+                % Stable dirchlet-like boundary conditions (u+-abs(u))*u/3
+                % with +- at left/right boundaries
                 case {'D', 'd', 'dirichlet', 'Dirichlet'}
-                    tau = s*e;
-                    closure = @(v) obj.Hi*tau*(((v(index)-s*abs(v(index)))/3)*v(index));
-                    penalty = -obj.Hi*tau;
+                    % tau = s*e;
+                    % closure = @(v) obj.Hi*tau*(((v(index)-s*abs(v(index)))/3)*v(index));
+                    % penalty = -obj.Hi*tau;
+
+                    magnitude = 2/3;
+                    tau = @(v) s*magnitude*obj.Hi*e*(v(index)-s*abs(v(index)))/2;
+                    closure = @(v) tau(v)*v(index);
+                    penalty = @(v) -tau(v);
                 otherwise
                     error('No such boundary condition: type = %s',type);
             end