Mercurial > repos > public > sbplib

changeset 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 78c75c95b7dd
children 1a5c8723c9be 037f203b9bf5 1f0f2276e678
files +rv/+time/RungekuttaExteriorRv.m +rv/+time/RungekuttaExteriorRvBdf.m +rv/+time/RungekuttaInteriorRv.m +rv/ResidualViscosity.m
diffstat 4 files changed, 32 insertions(+), 29 deletions(-) [+]
line wrap: on
line diff
--- a/+rv/+time/RungekuttaExteriorRv.m	Fri Jan 11 15:52:48 2019 +0100
+++ b/+rv/+time/RungekuttaExteriorRv.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
 
@@ -38,8 +39,7 @@
             obj.RV = RV;
             obj.DvDt = DvDt;
             obj.dvdt = obj.DvDt(obj.v);
-            [obj.viscosity, obj.Df] = RV.evaluate(obj.v,obj.dvdt);
-            obj.residual = obj.dvdt + obj.Df;
+            [obj.viscosity,  obj.Df, obj.firstOrderViscosity, obj.residualViscosity] = RV.evaluate(obj.v,obj.dvdt);
         end
 
         function [v, t] = getV(obj)
@@ -48,13 +48,12 @@
         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)            
             obj.dvdt = obj.DvDt(obj.v);
-            [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);
 
             % Fix the viscosity of the RHS function F
             F_visc = @(v,t) obj.F(v,t,obj.viscosity);
--- 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
--- a/+rv/+time/RungekuttaInteriorRv.m	Fri Jan 11 15:52:48 2019 +0100
+++ b/+rv/+time/RungekuttaInteriorRv.m	Thu Jan 17 10:25:06 2019 +0100
@@ -7,12 +7,14 @@
         n       % Time level
         coeffs  % The coefficents used for the RK time integration
         RV      % Residual Viscosity
-        viscosity       % Viscosity vector
-        DvDt            % Function for computing the time deriative used for the RV evaluation
+        DvDt    % Function for computing the time deriative used for the RV evaluation
 
-        residual
-        dvdt
-        Df
+        % Convenience properties. Only for plotting
+        viscosity % Total viscosity
+        residualViscosity % Residual viscosity
+        firstOrderViscosity % first order viscosity
+        dvdt % Evaluated time derivative in residual
+        Df % Evaluated flux in residual
     end
 
     methods
@@ -29,8 +31,7 @@
             obj.RV = RV;
             obj.DvDt = DvDt;
             obj.dvdt = obj.DvDt(obj.v);
-            [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
 
         function [v, t] = getV(obj)
@@ -39,7 +40,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)
@@ -47,8 +48,7 @@
             obj.t = obj.t + obj.k;
             obj.n = obj.n + 1;
             obj.dvdt = obj.DvDt(obj.v);
-            [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
\ No newline at end of file
--- a/+rv/ResidualViscosity.m	Fri Jan 11 15:52:48 2019 +0100
+++ b/+rv/ResidualViscosity.m	Thu Jan 17 10:25:06 2019 +0100
@@ -22,9 +22,12 @@
             obj.normalization = normalization;
         end
 
-        function [viscosity, Df] = evaluate(obj, v, dvdt)
+        function [viscosity, Df, firstOrderViscosity, residualViscosity] = evaluate(obj, v, dvdt)
             Df = obj.Df(v);
-            viscosity = min(obj.Cmax*obj.h*abs(obj.waveSpeed(v)), obj.Cres*obj.h^2*abs(dvdt + Df)/obj.normalization(v));
+            firstOrderViscosity = obj.Cmax*obj.h*abs(obj.waveSpeed(v));
+            residualViscosity = obj.Cres*obj.h^2*abs(dvdt + Df)/obj.normalization(v);
+            viscosity = min(firstOrderViscosity, residualViscosity);
+            %viscosity = obj.Cmax*obj.h*abs(obj.waveSpeed(v));
         end
     end
 end
\ No newline at end of file