Mercurial > repos > public > sbplib

comparison +rv/+time/RungekuttaExteriorRV.m @ 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
comparison
equal deleted inserted replaced
1014:e547794a9407 1015:9b7fcd5e4480
8 coeffs % The coefficents used for the RK time integration 8 coeffs % The coefficents used for the RK time integration
9 RV % Residual Viscosity 9 RV % Residual Viscosity
10 bdfOrder 10 bdfOrder
11 v_prev % Solution vector at previous time levels, used for the RV update 11 v_prev % Solution vector at previous time levels, used for the RV update
12 DvDt % Function for computing the time deriative used for the RV update 12 DvDt % Function for computing the time deriative used for the RV update
13
14
15 dudt
13 end 16 end
14 methods 17 methods
15 18
16 function obj = RungekuttaExteriorRV(F, k, t0, v0, RV, rkOrder, bdfOrder) 19 function obj = RungekuttaExteriorRV(F, k, t0, v0, RV, rkOrder, bdfOrder, dudt)
17 obj.F = F; 20 obj.F = F;
18 obj.k = k; 21 obj.k = k;
19 obj.t = t0; 22 obj.t = t0;
20 obj.v = v0; 23 obj.v = v0;
21 obj.n = 0; 24 obj.n = 0;
27 obj.RV = RV; 30 obj.RV = RV;
28 % TBD: For cases where h~k we could probably use rkOrder-2 here. 31 % TBD: For cases where h~k we could probably use rkOrder-2 here.
29 % TBD: Decide on if the initialization of the previous stages used by 32 % TBD: Decide on if the initialization of the previous stages used by
30 % the BDF should be done here, or if it should be checked for each 33 % the BDF should be done here, or if it should be checked for each
31 % step taken. 34 % step taken.
32 assert(bdfOrder >= 2); 35 assert((bdfOrder >= 1) && (bdfOrder <= 6));
33 obj.bdfOrder = bdfOrder; 36 obj.bdfOrder = bdfOrder;
34 obj.v_prev = []; 37 obj.v_prev = [];
35 obj.DvDt = rv.time.BDFDerivative(); 38 obj.DvDt = rv.time.BDFDerivative();
39
40 obj.dudt = dudt;
36 end 41 end
37 42
38 function [v, t] = getV(obj) 43 function [v, t] = getV(obj)
39 v = obj.v; 44 v = obj.v;
40 t = obj.t; 45 t = obj.t;
54 obj.v_prev(:,1) = obj.v; 59 obj.v_prev(:,1) = obj.v;
55 end 60 end
56 % Fix the viscosity of the RHS function F 61 % Fix the viscosity of the RHS function F
57 F_visc = @(v,t) obj.F(v,t,obj.RV.getViscosity()); %TBD: Remove state in RV? 62 F_visc = @(v,t) obj.F(v,t,obj.RV.getViscosity()); %TBD: Remove state in RV?
58 obj.v = time.rk.rungekutta(obj.v, obj.t, obj.k, F_visc, obj.coeffs); 63 obj.v = time.rk.rungekutta(obj.v, obj.t, obj.k, F_visc, obj.coeffs);
59 % Calculate dvdt and update RV
60 dvdt = obj.DvDt.evaluate(obj.v, obj.v_prev, obj.k);
61 obj.RV.update(obj.v,dvdt);
62 obj.t = obj.t + obj.k; 64 obj.t = obj.t + obj.k;
63 obj.n = obj.n + 1; 65 obj.n = obj.n + 1;
66 % Calculate dvdt and update RV for the new time level
67 dvdt = obj.DvDt.evaluate(obj.v, obj.v_prev, obj.k);
68 if ((size(obj.v_prev,2) >= 4) && (size(obj.v_prev,2) <= obj.bdfOrder))
69 obj.RV.update(obj.v,dvdt);
70 end
71 %obj.RV.update(obj.v,obj.dudt(obj.t));
64 end 72 end
65 end 73 end
66 end 74 end