Mercurial > repos > public > sbplib

changeset 851:ab2e5a24ddde feature/burgers1d

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
- Fix bug when constructing closure for narrow stencils - Update the residual outside of the RK time steps. At least for the inner convergence rate, updating the residual inside does not seem to be required.
author Vidar Stiernström <vidar.stiernstrom@it.uu.se>
date Fri, 21 Sep 2018 15:33:15 +0200
parents 4e5e53d6336c
children fbb8be3177c8
files +scheme/Burgers1D.m +time/RungekuttaRV.m
diffstat 2 files changed, 7 insertions(+), 4 deletions(-) [+]
line wrap: on
line diff
--- a/+scheme/Burgers1D.m	Fri Sep 21 09:14:38 2018 +0200
+++ b/+scheme/Burgers1D.m	Fri Sep 21 15:33:15 2018 +0200
@@ -43,8 +43,8 @@
                     if (strcmp(dissipation,'on'))
                         DissipationOp = (ops.Dp-ops.Dm)/2;
                     end
-                    d_l = D1;
-                    d_r = D1;
+                    d_l = ops.e_l'*D1;
+                    d_r = ops.e_r'*D1;
                 otherwise
                     error('Other operator types not yet supported', operator_type);
             end
@@ -93,7 +93,7 @@
                 % 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*(((v(i_b)-s*abs(v(i_b)))/3)*(v(i_b)) - e'*((obj.params.eps(i_b) + viscosity(i_b))*d*v));
+                    closure = @(v, viscosity) p*(((v(i_b)-s*abs(v(i_b)))/3)*(v(i_b)) - ((obj.params.eps(i_b) + viscosity(i_b))*d*v));
                     switch class(data)
                         case 'double'
                             penalty = s*p*data;
--- a/+time/RungekuttaRV.m	Fri Sep 21 09:14:38 2018 +0200
+++ b/+time/RungekuttaRV.m	Fri Sep 21 15:33:15 2018 +0200
@@ -34,9 +34,12 @@
         end
 
         function obj = step(obj)
-            obj.v = time.rk.rungekuttaRV(obj.v, obj.t, obj.k, obj.F, obj.RV, obj.coeffs);
+            F = @(v,t) obj.F(v,t,obj.RV.getViscosity());
+            v_p = obj.v;
+            obj.v = time.rk.rungekutta(obj.v, obj.t, obj.k, F, obj.coeffs);
             obj.t = obj.t + obj.k;
             obj.n = obj.n + 1;
+            obj.RV.update(obj.v,v_p,obj.k);
         end
     end