Mercurial > repos > public > sbplib

changeset 1015:9b7fcd5e4480 feature/advectionRV

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Debug ResidualViscosity - Pass exact time derivative to RungeKuttaExteriorRV and use that for evaluating the residual - Start bootstrapping from later time level with higher order bdf
author Vidar Stiernström <vidar.stiernstrom@it.uu.se>
date Thu, 06 Dec 2018 17:03:22 +0100
parents e547794a9407
children 4b42999874c0
files +rv/+time/RungekuttaExteriorRV.m +rv/ResidualViscosity.m
diffstat 2 files changed, 16 insertions(+), 6 deletions(-) [+]
line wrap: on
line diff
--- a/+rv/+time/RungekuttaExteriorRV.m	Thu Dec 06 11:30:47 2018 +0100
+++ b/+rv/+time/RungekuttaExteriorRV.m	Thu Dec 06 17:03:22 2018 +0100
@@ -10,10 +10,13 @@
         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
     end
     methods
 
-        function obj = RungekuttaExteriorRV(F, k, t0, v0, RV, rkOrder, bdfOrder)
+        function obj = RungekuttaExteriorRV(F, k, t0, v0, RV, rkOrder, bdfOrder, dudt)
             obj.F = F;
             obj.k = k;
             obj.t = t0;
@@ -29,10 +32,12 @@
             %  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 >= 2);
+            assert((bdfOrder >= 1) && (bdfOrder <= 6));
             obj.bdfOrder = bdfOrder;
             obj.v_prev = [];
             obj.DvDt = rv.time.BDFDerivative();
+
+            obj.dudt = dudt;
         end
 
         function [v, t] = getV(obj)
@@ -56,11 +61,14 @@
             % Fix the viscosity of the RHS function F
             F_visc = @(v,t) obj.F(v,t,obj.RV.getViscosity()); %TBD: Remove state in RV?
             obj.v = time.rk.rungekutta(obj.v, obj.t, obj.k, F_visc, obj.coeffs);
-            % Calculate dvdt and update RV
-            dvdt = obj.DvDt.evaluate(obj.v, obj.v_prev, obj.k);
-            obj.RV.update(obj.v,dvdt);
             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))
+                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 11:30:47 2018 +0100
+++ b/+rv/ResidualViscosity.m	Thu Dec 06 17:03:22 2018 +0100
@@ -37,7 +37,9 @@
             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 = 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)