Mercurial > repos > public > sbplib

diff +rv/+time/RungekuttaExteriorRvBdf.m @ 1031:2ef20d00b386 feature/advectionRV

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
For easier comparison, return both the first order and residual viscosity when evaluating the residual. Add the first order and residual viscosity to the state of the RungekuttaRV time steppers
author Vidar Stiernström <vidar.stiernstrom@it.uu.se>
date Thu, 17 Jan 2019 10:25:06 +0100
parents dce08a74e0ad
children 010bb2677230
line wrap: on
line diff
--- a/+rv/+time/RungekuttaExteriorRvBdf.m	Fri Jan 11 15:52:48 2019 +0100
+++ b/+rv/+time/RungekuttaExteriorRvBdf.m	Thu Jan 17 10:25:06 2019 +0100
@@ -9,7 +9,6 @@
         
         % Properties related to the residual viscositys
         RV              % Residual Viscosity operator
-        viscosity       % Viscosity vector
         v_prev          % Solution vector at previous time levels, used for the RV evaluation
         DvDt            % Function for computing the time deriative used for the RV evaluation
         lowerBdfOrder   % Orders of the approximation of the time deriative, used for the RV evaluation.
@@ -18,9 +17,11 @@
                         % Dictates the order of accuracy used once the boot-strapping is complete.
         
         % Convenience properties. Only for plotting
-        residual
-        dvdt
-        Df
+        viscosity % Total viscosity
+        residualViscosity % Residual viscosity
+        firstOrderViscosity % first order viscosity
+        dvdt % Evaluated time derivative in residual
+        Df % Evaluated flux in residual
     end
     methods
         function obj = RungekuttaExteriorRvBdf(F, k, t0, v0, RV, rkOrder, bdfOrders)
@@ -47,7 +48,8 @@
             obj.v_prev = [];
             obj.DvDt = rv.time.BdfDerivative();
             obj.viscosity = zeros(size(v0));
-            obj.residual = zeros(size(v0));
+            obj.firstOrderViscosity = zeros(size(v0));
+            obj.residualViscosity = zeros(size(v0));
             obj.dvdt = zeros(size(v0));
             obj.Df = zeros(size(v0));
         end
@@ -58,7 +60,7 @@
         end
 
         function state = getState(obj)
-            state = struct('v', obj.v, 'residual', obj.residual, 'dvdt', obj.dvdt, 'Df', obj.Df, 'viscosity', obj.viscosity, 't', obj.t);
+            state = struct('v', obj.v, 'dvdt', obj.dvdt, 'Df', obj.Df, 'viscosity', obj.viscosity, 'residualViscosity', obj.residualViscosity, 'firstOrderViscosity', obj.firstOrderViscosity, 't', obj.t);
         end
 
         function obj = step(obj)
@@ -81,8 +83,7 @@
             %Calculate dvdt and evaluate RV for the new time level
             if ((numStoredStages >=  obj.lowerBdfOrder) && (numStoredStages <= obj.upperBdfOrder))
                 obj.dvdt = obj.DvDt.evaluate(obj.v, obj.v_prev, obj.k);
-                [obj.viscosity, obj.Df] = obj.RV.evaluate(obj.v,obj.dvdt);
-                obj.residual = obj.dvdt + obj.Df;
+                [obj.viscosity, obj.Df, obj.firstOrderViscosity, obj.residualViscosity] = obj.RV.evaluate(obj.v,obj.dvdt);
             end
         end
     end