Mercurial > repos > public > sbplib

changeset 1017:2d7c1333bd6c feature/advectionRV

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Add support for using the ODE to approximate the time derivative in the residual
author Vidar Stiernström <vidar.stiernstrom@it.uu.se>
date Tue, 11 Dec 2018 16:29:21 +0100
parents 4b42999874c0
children eca4d9df9365
files +rv/+time/RungekuttaExteriorRV.m +rv/+time/RungekuttaInteriorRV.m +rv/+time/rungekuttaRV.m +rv/ResidualViscosity.m
diffstat 4 files changed, 83 insertions(+), 68 deletions(-) [+]
line wrap: on
line diff
diff -r 4b42999874c0 -r 2d7c1333bd6c +rv/+time/RungekuttaExteriorRV.m
--- a/+rv/+time/RungekuttaExteriorRV.m	Fri Dec 07 13:11:53 2018 +0100
+++ b/+rv/+time/RungekuttaExteriorRV.m	Tue Dec 11 16:29:21 2018 +0100
@@ -6,17 +6,26 @@
         v       % Solution vector
         n       % Time level
         coeffs  % The coefficents used for the RK time integration
-        RV      % Residual Viscosity
-        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
-        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.
+        
+        % 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.
+                        % dictates which accuracy the boot-strapping should start from.
+        upperBdfOrder   % Orders of the approximation of the time deriative, used for the RV evaluation.
+                        % Dictates the order of accuracy used once the boot-strapping is complete.
+        
+        % Convenience properties. Only for plotting
+        residual
+        dvdt
+        Df
     end
     methods
 
-        function obj = RungekuttaExteriorRV(F, k, t0, v0, RV, rkOrder, bdfOrders)
+        % TODO: Decide on how to compute dvdt.
+        function obj = RungekuttaExteriorRV(F, k, t0, v0, RV, DvDt, rkOrder, bdfOrders)
             obj.F = F;
             obj.k = k;
             obj.t = t0;
@@ -34,11 +43,21 @@
             %       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.lowerBdfOrder = bdfOrders.lowerBdfOrder;
+            % obj.upperBdfOrder = bdfOrders.upperBdfOrder;
+            % assert((obj.lowerBdfOrder >= 1) && (obj.upperBdfOrder <= 6));
+            % obj.v_prev = [];
+            % obj.DvDt = rv.time.BdfDerivative();
+            % obj.viscosity = zeros(size(v0));
+            % obj.residual = zeros(size(v0));
+            % obj.dvdt = zeros(size(v0));
+            % obj.Df = zeros(size(v0));
+
+            % Using the ODE:
+            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;
         end
 
         function [v, t] = getV(obj)
@@ -47,30 +66,36 @@
         end
 
         function state = getState(obj)
-            [residual, u_t, grad_f] = obj.RV.getResidual();
-            state = struct('v', obj.v, 'residual', residual, 'u_t', u_t, 'grad_f', grad_f, 'viscosity', obj.RV.getViscosity(), 't', obj.t);
+            state = struct('v', obj.v, 'residual', obj.residual, 'dvdt', obj.dvdt, 'Df', obj.Df, 'viscosity', obj.viscosity, 't', obj.t);
         end
 
         function obj = step(obj)
-            % Store current time level
-            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;
-            end    
+            % Store current time level and update v_prev
+            % 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;
+            % end
+            
+            obj.dvdt = obj.DvDt(obj.v);
+            [obj.viscosity, obj.Df] = obj.RV.evaluate(obj.v,obj.dvdt);
+            obj.residual = obj.dvdt + obj.Df;
+
             % Fix the viscosity of the RHS function F
-            F_visc = @(v,t) obj.F(v,t,obj.RV.getViscosity()); %TBD: Remove state in RV?
+            F_visc = @(v,t) obj.F(v,t,obj.viscosity);
             obj.v = time.rk.rungekutta(obj.v, obj.t, obj.k, F_visc, obj.coeffs);
             obj.t = obj.t + obj.k;
             obj.n = obj.n + 1;
-            % Calculate dvdt and update RV for the new time level
-            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
+
+            % %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;
+            % end
         end
     end
 end
\ No newline at end of file
diff -r 4b42999874c0 -r 2d7c1333bd6c +rv/+time/RungekuttaInteriorRV.m
--- a/+rv/+time/RungekuttaInteriorRV.m	Fri Dec 07 13:11:53 2018 +0100
+++ b/+rv/+time/RungekuttaInteriorRV.m	Tue Dec 11 16:29:21 2018 +0100
@@ -7,20 +7,30 @@
         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
+
+        residual
+        dvdt
+        Df
     end
 
     methods
-        function obj = RungekuttaInteriorRV(F, k, t0, v0, RV, order)
+        function obj = RungekuttaInteriorRV(F, k, t0, v0, RV, DvDt, order)
             obj.F = F;
             obj.k = k;
             obj.t = t0;
             obj.v = v0;
             obj.n = 0;
-            obj.RV = RV;
             % Extract the coefficients for the specified order
             % used for the RK updates from the Butcher tableua.
             [s,a,b,c] = time.rk.butcherTableau(order);
             obj.coeffs = struct('s',s,'a',a,'b',b,'c',c);
+            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;
         end
 
         function [v, t] = getV(obj)
@@ -29,14 +39,16 @@
         end
 
         function state = getState(obj)
-            [residual, u_t, grad_f] = obj.RV.getResidual();
-            state = struct('v', obj.v, 'residual', residual, 'u_t', u_t, 'grad_f', grad_f, 'viscosity', obj.RV.getViscosity(), 't', obj.t);
+            state = struct('v', obj.v, 'residual', obj.residual, 'dvdt', obj.dvdt, 'Df', obj.Df, 'viscosity', obj.viscosity, 't', obj.t);
         end
 
         function obj = step(obj)
-            obj.v = rv.time.rungekuttaRV(obj.v, obj.t, obj.k, obj.F, obj.RV, obj.coeffs);
+            obj.v = rv.time.rungekuttaRV(obj.v, obj.t, obj.k, obj.F, obj.RV, obj.DvDt, obj.coeffs);
             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;
         end
     end
 end
\ No newline at end of file
diff -r 4b42999874c0 -r 2d7c1333bd6c +rv/+time/rungekuttaRV.m
--- a/+rv/+time/rungekuttaRV.m	Fri Dec 07 13:11:53 2018 +0100
+++ b/+rv/+time/rungekuttaRV.m	Tue Dec 11 16:29:21 2018 +0100
@@ -3,11 +3,11 @@
 % time derivatives. coeffs is a struct holding the RK coefficients
 % for the specific method. RV is the residual viscosity which is updated
 % in between the stages and after the updated solution is computed.
-function v = rungekuttaRV(v, t , dt, F, RV, coeffs)
+function v = rungekuttaRV(v, t , dt, F, RV, DvDt, coeffs)
     % Move one stage outside to avoid branching for updating the
     % residual inside the loop.
     k = zeros(length(v), coeffs.s);
-    k(:,1) = F(v,t,RV.getViscosity());
+    k(:,1) = F(v,t,RV.evaluate(v,DvDt(v)));
 
     % Compute the intermediate stages k
     for i = 2:coeffs.s
@@ -15,8 +15,8 @@
         for j = 1:i-1
             u = u + dt*coeffs.a(i,j)*k(:,j);
         end
-        RV.update(0.5*(u+v),(u-v)/(coeffs.c(i)*dt)); % Crank-Nicholson for time discretization
-        k(:,i) = F(u,t+coeffs.c(i)*dt, RV.getViscosity());
+        %RV.update(0.5*(u+v),(u-v)/(coeffs.c(i)*dt)); % Crank-Nicholson for time discretization
+        k(:,i) = F(u,t+coeffs.c(i)*dt, RV.evaluate(u,DvDt(u)));
     end
 
     % Compute the updated solution as a linear combination
@@ -25,6 +25,6 @@
     for i = 1:coeffs.s
         u = u + dt*coeffs.b(i)*k(:,i);
     end
-    RV.update(0.5*(u+v),(u-v)/dt); % Crank-Nicholson for time discretization
+    %RV.update(0.5*(u+v),(u-v)/dt); % Crank-Nicholson for time discretization
     v = u;
 end
diff -r 4b42999874c0 -r 2d7c1333bd6c +rv/ResidualViscosity.m
--- a/+rv/ResidualViscosity.m	Fri Dec 07 13:11:53 2018 +0100
+++ b/+rv/ResidualViscosity.m	Tue Dec 11 16:29:21 2018 +0100
@@ -1,52 +1,30 @@
 classdef ResidualViscosity < handle
     properties
-        D % Diff op approximating the gradient of the flux f(u)
+        Df % 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 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.
-        grad_f % Stores the latest approximated gradient of the flux
-        residual % Stores the computed residual
     end
 
     methods
         %  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)
+        function obj = ResidualViscosity(Df, waveSpeed, Cmax, Cres, h, normalization)
             default_arg('normalization',@(v)norm(v-mean(v),inf));
-            obj.D = D;
+            obj.Df = Df;
             obj.waveSpeed = waveSpeed;
             obj.h = h;
             obj.Cmax = Cmax;
             obj.Cres = Cres;
-            obj.viscosity = zeros(N,1);
-            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)/obj.normalization(v));
-        end
-
-        function [residual, u_t, grad_f] = getResidual(obj)
-            residual = obj.residual;
-            u_t = obj.u_t;
-            grad_f = obj.grad_f;
-        end
-
-        function viscosity = getViscosity(obj)
-            viscosity = obj.viscosity;
+        function [viscosity, Df] = 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));
         end
     end
 end
\ No newline at end of file