Mercurial > repos > public > sbplib

changeset 848:c8ea9bbdc62c feature/burgers1d

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Increase efficiency of computing closures
author Vidar Stiernström <vidar.stiernstrom@it.uu.se>
date Thu, 20 Sep 2018 19:02:01 +0200
parents 1c6f1595bb94
children 5b180c76578e
files +scheme/Burgers1D.m
diffstat 1 files changed, 5 insertions(+), 3 deletions(-) [+]
line wrap: on
line diff
--- a/+scheme/Burgers1D.m	Thu Sep 20 18:36:45 2018 +0200
+++ b/+scheme/Burgers1D.m	Thu Sep 20 19:02:01 2018 +0200
@@ -88,12 +88,12 @@
         function [closure, penalty] = boundary_condition(obj,boundary,type,data)
             default_arg('type','robin');
             default_arg('data',0);
-            [e, s, d] = obj.get_boundary_ops(boundary);
+            [e, s, d, i_b] = obj.get_boundary_ops(boundary);
             switch type
                 % Stable robin-like boundary conditions ((u+-abs(u))*u/3 - eps*u_x)) with +- at left/right boundary
                 case {'R','robin'}
                     p = s*obj.Hi*e;
-                    closure = @(v, viscosity) p*(e'*((v-s*abs(v))/3)*(e'*v) - e'*((obj.params.eps + viscosity).*d*v));
+                    closure = @(v, viscosity) p*(((v(i_b)-s*abs(v(i_b)))/3)*(v(i_b)) - e'*((obj.params.eps(i_b) + viscosity(i_b))*d*v));
                     switch class(data)
                         case 'double'
                             penalty = s*p*data;
@@ -110,16 +110,18 @@
 
         % Ruturns the boundary ops and sign for the boundary specified by the string boundary.
         % The right boundary is considered the positive boundary
-        function [e, s, d] = get_boundary_ops(obj,boundary)
+        function [e, s, d, i_b] = get_boundary_ops(obj,boundary)
             switch boundary
                 case 'l'
                     e = obj.e_l;
                     s = -1;
                     d = obj.d_l;
+                    i_b = 1;
                 case 'r'
                     e = obj.e_r;
                     s = 1;
                     d = obj.d_r;
+                    i_b = length(e);
                 otherwise
                     error('No such boundary: boundary = %s',boundary);
             end