Mercurial > repos > public > sbplib

changeset 1016:4b42999874c0 feature/advectionRV

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Add lower level for boot-strapping to RungeKuttaExteriorRV - Add a lower level to RungeKuttaExteriorRV for which bootstrapping starts, e.g start bootstrapping from time level 3 using a 3rd order BDF - Clean up ResidualViscosity
author Vidar Stiernström <vidar.stiernstrom@it.uu.se>
date Fri, 07 Dec 2018 13:11:53 +0100
parents 9b7fcd5e4480
children 2d7c1333bd6c
files +rv/+time/BDFDerivative.m +rv/+time/BdfDerivative.m +rv/+time/RungekuttaExteriorRV.m +rv/ResidualViscosity.m
diffstat 4 files changed, 59 insertions(+), 79 deletions(-) [+]
line wrap: on
line diff
--- a/+rv/+time/BDFDerivative.m	Thu Dec 06 17:03:22 2018 +0100
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,31 +0,0 @@
-classdef BDFDerivative < handle
-    properties
-        coefficients
-    end
-    methods %TBD: Decide if the BDF order should be passed in construction
-            %     and only a row of coefficients stored based on the order.
-            %     There would still be an implicit dependancy on the number
-            %     of vectors in v_prev and elements in coefficients.
-            %     In addition, dt could be stored, but this would be
-            %     inflexible if different step sizes are employed.
-        function obj = BDFDerivative()
-            obj.coefficients = obj.getBDFCoefficients();
-        end
-        % Add asserts here?
-        function DvDt = evaluate(obj, v, v_prev, dt)
-            order = size(v_prev,2);
-            DvDt = (obj.coefficients(order,1)*v - sum(obj.coefficients(order,2:order+1).*v_prev,2))/dt;
-        end
-    end
-    methods(Static)
-        function c = getBDFCoefficients()
-            c = zeros(6,7);
-            c(1,1) = 1; c(1,2) = 1;
-            c(2,1) = 3/2; c(2,2) = 4/2; c(2,3) = -1/2;
-            c(3,1) = 11/6; c(3,2) = 18/6; c(3,3) = -9/6; c(3,4) = 2/6;
-            c(4,1) = 25/12; c(4,2) = 48/12; c(4,3) = -36/12; c(4,4) = 16/12; c(4,5) = -3/12;
-            c(5,1) = 137/60; c(5,2) = 300/60; c(5,3) = -300/60; c(5,4) = 200/60; c(5,5) = -75/60; c(5,6) = 12/60;
-            c(6,1) = 147/60; c(6,2) = 360/60; c(6,3) = -450/60; c(6,4) = 400/60; c(6,5) = -225/60; c(6,6) = 72/60; c(6,7) = -10/60;
-        end
-    end
-end
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/+rv/+time/BdfDerivative.m	Fri Dec 07 13:11:53 2018 +0100
@@ -0,0 +1,31 @@
+classdef BdfDerivative < handle
+    properties
+        coefficients
+    end
+    methods %TBD: Decide if the BDF order should be passed in construction
+            %     and only a row of coefficients stored based on the order.
+            %     There would still be an implicit dependancy on the number
+            %     of vectors in v_prev and elements in coefficients.
+            %     In addition, dt could be stored, but this would be
+            %     inflexible if different step sizes are employed.
+        function obj = BdfDerivative()
+            obj.coefficients = obj.getBdfCoefficients();
+        end
+        % Add asserts here?
+        function DvDt = evaluate(obj, v, v_prev, dt)
+            order = size(v_prev,2);
+            DvDt = (obj.coefficients(order,1)*v - sum(obj.coefficients(order,2:order+1).*v_prev,2))/dt;
+        end
+    end
+    methods(Static)
+        function c = getBdfCoefficients()
+            c = zeros(6,7);
+            c(1,1) = 1; c(1,2) = 1;
+            c(2,1) = 3/2; c(2,2) = 4/2; c(2,3) = -1/2;
+            c(3,1) = 11/6; c(3,2) = 18/6; c(3,3) = -9/6; c(3,4) = 2/6;
+            c(4,1) = 25/12; c(4,2) = 48/12; c(4,3) = -36/12; c(4,4) = 16/12; c(4,5) = -3/12;
+            c(5,1) = 137/60; c(5,2) = 300/60; c(5,3) = -300/60; c(5,4) = 200/60; c(5,5) = -75/60; c(5,6) = 12/60;
+            c(6,1) = 147/60; c(6,2) = 360/60; c(6,3) = -450/60; c(6,4) = 400/60; c(6,5) = -225/60; c(6,6) = 72/60; c(6,7) = -10/60;
+        end
+    end
+end
--- a/+rv/+time/RungekuttaExteriorRV.m	Thu Dec 06 17:03:22 2018 +0100
+++ b/+rv/+time/RungekuttaExteriorRV.m	Fri Dec 07 13:11:53 2018 +0100
@@ -7,16 +7,16 @@
         n       % Time level
         coeffs  % The coefficents used for the RK time integration
         RV      % Residual Viscosity
-        bdfOrder
         v_prev  % Solution vector at previous time levels, used for the RV update
         DvDt    % Function for computing the time deriative used for the RV update
-
-
-        dudt
+        lowerBdfOrder % Orders of the approximation of the time deriative, used for the RV update.
+                      % dictates which accuracy the boot-strapping should start from.
+        upperBdfOrder % Orders of the approximation of the time deriative, used for the RV update.
+                      % Dictates the order of accuracy used once the boot-strapping is complete.
     end
     methods
 
-        function obj = RungekuttaExteriorRV(F, k, t0, v0, RV, rkOrder, bdfOrder, dudt)
+        function obj = RungekuttaExteriorRV(F, k, t0, v0, RV, rkOrder, bdfOrders)
             obj.F = F;
             obj.k = k;
             obj.t = t0;
@@ -28,16 +28,17 @@
             obj.coeffs = struct('s',s,'a',a,'b',b,'c',c);
         
             obj.RV = RV;
-            %  TBD: For cases where h~k we could probably use rkOrder-2 here.
             %  TBD: Decide on if the initialization of the previous stages used by
             %       the BDF should be done here, or if it should be checked for each
             %       step taken.
-            assert((bdfOrder >= 1) && (bdfOrder <= 6));
-            obj.bdfOrder = bdfOrder;
+            %       If it is moved here, then multiple branching stages can be removed in step()
+            %       but this will effectively result in a plotted simulation starting from n = upperBdfOrder.
+            %       In addition, the properties lowerBdfOrder and upperBdfOrder can be removed.
+            obj.lowerBdfOrder = bdfOrders.lowerBdfOrder;
+            obj.upperBdfOrder = bdfOrders.upperBdfOrder;
+            assert((obj.lowerBdfOrder >= 1) && (obj.upperBdfOrder <= 6));
             obj.v_prev = [];
-            obj.DvDt = rv.time.BDFDerivative();
-
-            obj.dudt = dudt;
+            obj.DvDt = rv.time.BdfDerivative();
         end
 
         function [v, t] = getV(obj)
@@ -52,8 +53,10 @@
 
         function obj = step(obj)
             % Store current time level
-            if (size(obj.v_prev,2) < obj.bdfOrder)
+            numStoredStages = size(obj.v_prev,2);
+            if (numStoredStages < obj.upperBdfOrder)
                 obj.v_prev = [obj.v, obj.v_prev];
+                numStoredStages = numStoredStages+1;
             else
                 obj.v_prev(:,2:end) = obj.v_prev(:,1:end-1);
                 obj.v_prev(:,1) = obj.v;
@@ -64,11 +67,10 @@
             obj.t = obj.t + obj.k;
             obj.n = obj.n + 1;
             % Calculate dvdt and update RV for the new time level
-            dvdt = obj.DvDt.evaluate(obj.v, obj.v_prev, obj.k);
-            if ((size(obj.v_prev,2) >= 4) && (size(obj.v_prev,2) <= obj.bdfOrder))
+            if ((numStoredStages >=  obj.lowerBdfOrder) && (numStoredStages <= obj.upperBdfOrder))
+                dvdt = obj.DvDt.evaluate(obj.v, obj.v_prev, obj.k);
                 obj.RV.update(obj.v,dvdt);
             end
-            %obj.RV.update(obj.v,obj.dudt(obj.t));
         end
     end
 end
\ No newline at end of file
--- a/+rv/ResidualViscosity.m	Thu Dec 06 17:03:22 2018 +0100
+++ b/+rv/ResidualViscosity.m	Fri Dec 07 13:11:53 2018 +0100
@@ -1,12 +1,13 @@
-% class describing the viscosity 
 classdef ResidualViscosity < handle
     properties
         D % Diff op approximating the gradient of the flux f(u)
         waveSpeed % Wave speed at each grid point, e.g f'(u). %TBD: Better naming?
-        Cmax % Constant controling magnitude of upwind dissipation
-        Cres % Constant controling magnitude residual dissipation
-        h % Length scale used for scaling the viscosity.
+        Cmax % Constant controlling relative amount of upwind dissipation
+        Cres % Constant controling relative amount of upwind dissipation
+        h % Length scale used for scaling the viscosity. Typically grid spacing.
         viscosity % Stores the computed viscosity.
+        normalization % Function used to normalize the residual such that it is amplified in the
+                      % shocks.
 
         % Convenice (for verification and plotting) TBD: Decide on if it should be kept.
         u_t % Stores the latest approximated time derivative of the solution.
@@ -15,13 +16,10 @@
     end
 
     methods
-        % TODO: - Consider passing residual normalization as a function handle.
-        %         or choosing a type of normalization on construction.
-        %         Could for example be 1, norm((v-mean(v),inf) or normInfNeighborhood(v)
-        %         but working
-        %       - Decide on how to treat waveSpeed. It would be nice to just pass a constant value without
-        %         wrapping it in a function.
-        function obj = ResidualViscosity(D, waveSpeed, Cmax, Cres, h, N)
+        %  TBD: Decide on how to treat waveSpeed. It would be nice to just pass a constant value without
+        %       wrapping it in a function.
+        function obj = ResidualViscosity(D, waveSpeed, Cmax, Cres, h, N, normalization)
+            default_arg('normalization',@(v)norm(v-mean(v),inf));
             obj.D = D;
             obj.waveSpeed = waveSpeed;
             obj.h = h;
@@ -31,22 +29,14 @@
             obj.u_t = zeros(N,1);
             obj.grad_f = zeros(N,1);
             obj.residual = zeros(N,1);
+            obj.normalization = normalization;
         end
 
         function obj = update(obj, v, dvdt)
             obj.u_t = dvdt;
             obj.grad_f = obj.D(v);
             obj.residual = obj.u_t + obj.grad_f;
-            %obj.viscosity = min(obj.Cmax*obj.h*abs(obj.waveSpeed(v)), obj.Cres*obj.h^2*abs(obj.residual)/norm(v-mean(v),inf));
-            obj.viscosity = obj.smoothen(obj.Cres*obj.h^2*abs(obj.residual)/norm(v-mean(v),inf));
-
-        end
-
-        function smoothendVector = smoothen(obj, vector)
-            smoothendVector = vector;
-            for i = 2:length(vector)-1
-                smoothendVector(i) = (1/6)*(vector(i-1) + 4*vector(i) + vector(i+1));
-            end
+            obj.viscosity = min(obj.Cmax*obj.h*abs(obj.waveSpeed(v)), obj.Cres*obj.h^2*abs(obj.residual)/obj.normalization(v));
         end
 
         function [residual, u_t, grad_f] = getResidual(obj)
@@ -59,16 +49,4 @@
             viscosity = obj.viscosity;
         end
     end
-    % Remove or fix. Should be able to handle values close to zero. Should work in 2d and 3d.
-    methods (Static)
-        function R_norm = normInfNeighborhood(v)
-            n = length(v);
-            R_norm = zeros(n,1);
-            R_norm(1,1) = norm(v(1:3), inf);
-            R_norm(n,1) = norm(v(n-3:n), inf);
-            for i = 2:n-1
-                R_norm(i,1) = norm(v(i-1:i+1), inf);
-            end
-        end        
-    end
 end
\ No newline at end of file