Mercurial > repos > public > sbplib

comparison +rv/+time/RungekuttaExteriorRV.m @ 1014:e547794a9407 feature/advectionRV

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Add boot-strapping to RungeKuttaExteriorRV - Higher order BDF approximations are successively used as increasing number of time levels are obtained.
author Vidar Stiernström <vidar.stiernstrom@it.uu.se>
date Thu, 06 Dec 2018 11:30:47 +0100
parents eb441fbdf379
children 9b7fcd5e4480
comparison
equal deleted inserted replaced
1013:eb441fbdf379 1014:e547794a9407
5 t % Time point 5 t % Time point
6 v % Solution vector 6 v % Solution vector
7 n % Time level 7 n % Time level
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 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
11 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
12 end 13 end
13 methods 14 methods
14 15
22 % used for the RK updates from the Butcher tableua. 23 % used for the RK updates from the Butcher tableua.
23 [s,a,b,c] = time.rk.butcherTableau(rkOrder); 24 [s,a,b,c] = time.rk.butcherTableau(rkOrder);
24 obj.coeffs = struct('s',s,'a',a,'b',b,'c',c); 25 obj.coeffs = struct('s',s,'a',a,'b',b,'c',c);
25 26
26 obj.RV = RV; 27 obj.RV = RV;
27 % TODO: Initialize v_prev.
28 obj.v_prev = zeros(length(obj.v),bdfOrder);
29 % TBD: For cases where h~k we could probably use rkOrder-2 here. 28 % TBD: For cases where h~k we could probably use rkOrder-2 here.
30 obj.DvDt = rv.time.BDFDerivative(bdfOrder,k); 29 % 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
31 % step taken.
32 assert(bdfOrder >= 2);
33 obj.bdfOrder = bdfOrder;
34 obj.v_prev = [];
35 obj.DvDt = rv.time.BDFDerivative();
31 end 36 end
32 37
33 function [v, t] = getV(obj) 38 function [v, t] = getV(obj)
34 v = obj.v; 39 v = obj.v;
35 t = obj.t; 40 t = obj.t;
39 [residual, u_t, grad_f] = obj.RV.getResidual(); 44 [residual, u_t, grad_f] = obj.RV.getResidual();
40 state = struct('v', obj.v, 'residual', residual, 'u_t', u_t, 'grad_f', grad_f, 'viscosity', obj.RV.getViscosity(), 't', obj.t); 45 state = struct('v', obj.v, 'residual', residual, 'u_t', u_t, 'grad_f', grad_f, 'viscosity', obj.RV.getViscosity(), 't', obj.t);
41 end 46 end
42 47
43 function obj = step(obj) 48 function obj = step(obj)
44 % Store current time level 49 % Store current time level
45 obj.v_prev(:,2:end) = obj.v_prev(:,1:end-1); 50 if (size(obj.v_prev,2) < obj.bdfOrder)
46 obj.v_prev(:,1) = obj.v; 51 obj.v_prev = [obj.v, obj.v_prev];
52 else
53 obj.v_prev(:,2:end) = obj.v_prev(:,1:end-1);
54 obj.v_prev(:,1) = obj.v;
55 end
47 % Fix the viscosity of the RHS function F 56 % Fix the viscosity of the RHS function F
48 F_visc = @(v,t)F(v,t,obj.RV.getViscosity()); %Remove state in RV? 57 F_visc = @(v,t) obj.F(v,t,obj.RV.getViscosity()); %TBD: Remove state in RV?
49 obj.v = time.rk.rungekutta(obj.v, obj.t, obj.k, F_visc, obj.RV, obj.coeffs); 58 obj.v = time.rk.rungekutta(obj.v, obj.t, obj.k, F_visc, obj.coeffs);
50 % Calculate dvdt and update RV 59 % Calculate dvdt and update RV
51 dvdt = obj.DvDt.evaluate(v, v_prev, lenght(v_prev)); 60 dvdt = obj.DvDt.evaluate(obj.v, obj.v_prev, obj.k);
52 RV.update(v,dvdt); 61 obj.RV.update(obj.v,dvdt);
53 obj.t = obj.t + obj.k; 62 obj.t = obj.t + obj.k;
54 obj.n = obj.n + 1; 63 obj.n = obj.n + 1;
55 end 64 end
56 end 65 end
57 end 66 end