Mercurial > repos > public > sbplib

comparison +rv/+time/RungekuttaExteriorRV.m @ 1029:dce08a74e0ad feature/advectionRV

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Create a separate class of RungekuttaExteriorRV which uses BDFs for computing the time derivative. Remove BDFs from RungekuttaExteriorRV
author Vidar Stiernström <vidar.stiernstrom@it.uu.se>
date Fri, 11 Jan 2019 15:47:10 +0100
parents 5359a61cb4d9
children
comparison
equal deleted inserted replaced
1028:5df155ededcd 1029:dce08a74e0ad
22 dvdt 22 dvdt
23 Df 23 Df
24 end 24 end
25 methods 25 methods
26 26
27 % TODO: Decide on how to compute dvdt. 27 function obj = RungekuttaExteriorRV(F, k, t0, v0, RV, DvDt, rkOrder)
28 function obj = RungekuttaExteriorRV(F, k, t0, v0, RV, DvDt, rkOrder, bdfOrders)
29 obj.F = F; 28 obj.F = F;
30 obj.k = k; 29 obj.k = k;
31 obj.t = t0; 30 obj.t = t0;
32 obj.v = v0; 31 obj.v = v0;
33 obj.n = 0; 32 obj.n = 0;
35 % used for the RK updates from the Butcher tableua. 34 % used for the RK updates from the Butcher tableua.
36 [s,a,b,c] = time.rk.butcherTableau(rkOrder); 35 [s,a,b,c] = time.rk.butcherTableau(rkOrder);
37 obj.coeffs = struct('s',s,'a',a,'b',b,'c',c); 36 obj.coeffs = struct('s',s,'a',a,'b',b,'c',c);
38 37
39 obj.RV = RV; 38 obj.RV = RV;
40 % TBD: Decide on if the initialization of the previous stages used by
41 % the BDF should be done here, or if it should be checked for each
42 % step taken.
43 % If it is moved here, then multiple branching stages can be removed in step()
44 % but this will effectively result in a plotted simulation starting from n = upperBdfOrder.
45 % In addition, the properties lowerBdfOrder and upperBdfOrder can be removed.
46 % obj.lowerBdfOrder = bdfOrders.lowerBdfOrder;
47 % obj.upperBdfOrder = bdfOrders.upperBdfOrder;
48 % assert((obj.lowerBdfOrder >= 1) && (obj.upperBdfOrder <= 6));
49 % obj.v_prev = [];
50 % obj.DvDt = rv.time.BdfDerivative();
51 % obj.viscosity = zeros(size(v0));
52 % obj.residual = zeros(size(v0));
53 % obj.dvdt = zeros(size(v0));
54 % obj.Df = zeros(size(v0));
55
56 % Using the ODE:
57 obj.DvDt = DvDt; 39 obj.DvDt = DvDt;
58 obj.dvdt = obj.DvDt(obj.v); 40 obj.dvdt = obj.DvDt(obj.v);
59 [obj.viscosity, obj.Df] = RV.evaluate(obj.v,obj.dvdt); 41 [obj.viscosity, obj.Df] = RV.evaluate(obj.v,obj.dvdt);
60 obj.residual = obj.dvdt + obj.Df; 42 obj.residual = obj.dvdt + obj.Df;
61 end 43 end
67 49
68 function state = getState(obj) 50 function state = getState(obj)
69 state = struct('v', obj.v, 'residual', obj.residual, 'dvdt', obj.dvdt, 'Df', obj.Df, 'viscosity', obj.viscosity, 't', obj.t); 51 state = struct('v', obj.v, 'residual', obj.residual, 'dvdt', obj.dvdt, 'Df', obj.Df, 'viscosity', obj.viscosity, 't', obj.t);
70 end 52 end
71 53
72 function obj = step(obj) 54 function obj = step(obj)
73 % % Store current time level and update v_prev
74 % numStoredStages = size(obj.v_prev,2);
75 % if (numStoredStages < obj.upperBdfOrder)
76 % obj.v_prev = [obj.v, obj.v_prev];
77 % numStoredStages = numStoredStages+1;
78 % else
79 % obj.v_prev(:,2:end) = obj.v_prev(:,1:end-1);
80 % obj.v_prev(:,1) = obj.v;
81 % end
82
83 obj.dvdt = obj.DvDt(obj.v); 55 obj.dvdt = obj.DvDt(obj.v);
84 [obj.viscosity, obj.Df] = obj.RV.evaluate(obj.v,obj.dvdt); 56 [obj.viscosity, obj.Df] = obj.RV.evaluate(obj.v,obj.dvdt);
85 obj.residual = obj.dvdt + obj.Df; 57 obj.residual = obj.dvdt + obj.Df;
86 58
87 % Fix the viscosity of the RHS function F 59 % Fix the viscosity of the RHS function F
88 F_visc = @(v,t) obj.F(v,t,obj.viscosity); 60 F_visc = @(v,t) obj.F(v,t,obj.viscosity);
89 obj.v = time.rk.rungekutta(obj.v, obj.t, obj.k, F_visc, obj.coeffs); 61 obj.v = time.rk.rungekutta(obj.v, obj.t, obj.k, F_visc, obj.coeffs);
90 obj.t = obj.t + obj.k; 62 obj.t = obj.t + obj.k;
91 obj.n = obj.n + 1; 63 obj.n = obj.n + 1;
92
93 % %Calculate dvdt and evaluate RV for the new time level
94 % if ((numStoredStages >= obj.lowerBdfOrder) && (numStoredStages <= obj.upperBdfOrder))
95 % obj.dvdt = obj.DvDt.evaluate(obj.v, obj.v_prev, obj.k);
96 % [obj.viscosity, obj.Df] = obj.RV.evaluate(obj.v,obj.dvdt);
97 % obj.residual = obj.dvdt + obj.Df;
98 % end
99 end 64 end
100 end 65 end
101 end 66 end