Mercurial > repos > public > sbplib

comparison +rv/+time/RungekuttaExteriorRV.m @ 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
comparison
equal deleted inserted replaced
1015:9b7fcd5e4480 1016:4b42999874c0
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
11 v_prev % Solution vector at previous time levels, used for the RV update 10 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 11 DvDt % Function for computing the time deriative used for the RV update
13 12 lowerBdfOrder % Orders of the approximation of the time deriative, used for the RV update.
14 13 % dictates which accuracy the boot-strapping should start from.
15 dudt 14 upperBdfOrder % Orders of the approximation of the time deriative, used for the RV update.
15 % Dictates the order of accuracy used once the boot-strapping is complete.
16 end 16 end
17 methods 17 methods
18 18
19 function obj = RungekuttaExteriorRV(F, k, t0, v0, RV, rkOrder, bdfOrder, dudt) 19 function obj = RungekuttaExteriorRV(F, k, t0, v0, RV, rkOrder, bdfOrders)
20 obj.F = F; 20 obj.F = F;
21 obj.k = k; 21 obj.k = k;
22 obj.t = t0; 22 obj.t = t0;
23 obj.v = v0; 23 obj.v = v0;
24 obj.n = 0; 24 obj.n = 0;
26 % used for the RK updates from the Butcher tableua. 26 % used for the RK updates from the Butcher tableua.
27 [s,a,b,c] = time.rk.butcherTableau(rkOrder); 27 [s,a,b,c] = time.rk.butcherTableau(rkOrder);
28 obj.coeffs = struct('s',s,'a',a,'b',b,'c',c); 28 obj.coeffs = struct('s',s,'a',a,'b',b,'c',c);
29 29
30 obj.RV = RV; 30 obj.RV = RV;
31 % TBD: For cases where h~k we could probably use rkOrder-2 here.
32 % TBD: Decide on if the initialization of the previous stages used by 31 % TBD: Decide on if the initialization of the previous stages used by
33 % the BDF should be done here, or if it should be checked for each 32 % the BDF should be done here, or if it should be checked for each
34 % step taken. 33 % step taken.
35 assert((bdfOrder >= 1) && (bdfOrder <= 6)); 34 % If it is moved here, then multiple branching stages can be removed in step()
36 obj.bdfOrder = bdfOrder; 35 % but this will effectively result in a plotted simulation starting from n = upperBdfOrder.
36 % In addition, the properties lowerBdfOrder and upperBdfOrder can be removed.
37 obj.lowerBdfOrder = bdfOrders.lowerBdfOrder;
38 obj.upperBdfOrder = bdfOrders.upperBdfOrder;
39 assert((obj.lowerBdfOrder >= 1) && (obj.upperBdfOrder <= 6));
37 obj.v_prev = []; 40 obj.v_prev = [];
38 obj.DvDt = rv.time.BDFDerivative(); 41 obj.DvDt = rv.time.BdfDerivative();
39
40 obj.dudt = dudt;
41 end 42 end
42 43
43 function [v, t] = getV(obj) 44 function [v, t] = getV(obj)
44 v = obj.v; 45 v = obj.v;
45 t = obj.t; 46 t = obj.t;
50 state = struct('v', obj.v, 'residual', residual, 'u_t', u_t, 'grad_f', grad_f, 'viscosity', obj.RV.getViscosity(), 't', obj.t); 51 state = struct('v', obj.v, 'residual', residual, 'u_t', u_t, 'grad_f', grad_f, 'viscosity', obj.RV.getViscosity(), 't', obj.t);
51 end 52 end
52 53
53 function obj = step(obj) 54 function obj = step(obj)
54 % Store current time level 55 % Store current time level
55 if (size(obj.v_prev,2) < obj.bdfOrder) 56 numStoredStages = size(obj.v_prev,2);
57 if (numStoredStages < obj.upperBdfOrder)
56 obj.v_prev = [obj.v, obj.v_prev]; 58 obj.v_prev = [obj.v, obj.v_prev];
59 numStoredStages = numStoredStages+1;
57 else 60 else
58 obj.v_prev(:,2:end) = obj.v_prev(:,1:end-1); 61 obj.v_prev(:,2:end) = obj.v_prev(:,1:end-1);
59 obj.v_prev(:,1) = obj.v; 62 obj.v_prev(:,1) = obj.v;
60 end 63 end
61 % Fix the viscosity of the RHS function F 64 % Fix the viscosity of the RHS function F
62 F_visc = @(v,t) obj.F(v,t,obj.RV.getViscosity()); %TBD: Remove state in RV? 65 F_visc = @(v,t) obj.F(v,t,obj.RV.getViscosity()); %TBD: Remove state in RV?
63 obj.v = time.rk.rungekutta(obj.v, obj.t, obj.k, F_visc, obj.coeffs); 66 obj.v = time.rk.rungekutta(obj.v, obj.t, obj.k, F_visc, obj.coeffs);
64 obj.t = obj.t + obj.k; 67 obj.t = obj.t + obj.k;
65 obj.n = obj.n + 1; 68 obj.n = obj.n + 1;
66 % Calculate dvdt and update RV for the new time level 69 % Calculate dvdt and update RV for the new time level
67 dvdt = obj.DvDt.evaluate(obj.v, obj.v_prev, obj.k); 70 if ((numStoredStages >= obj.lowerBdfOrder) && (numStoredStages <= obj.upperBdfOrder))
68 if ((size(obj.v_prev,2) >= 4) && (size(obj.v_prev,2) <= obj.bdfOrder)) 71 dvdt = obj.DvDt.evaluate(obj.v, obj.v_prev, obj.k);
69 obj.RV.update(obj.v,dvdt); 72 obj.RV.update(obj.v,dvdt);
70 end 73 end
71 %obj.RV.update(obj.v,obj.dudt(obj.t));
72 end 74 end
73 end 75 end
74 end 76 end