Mercurial > repos > public > sbplib

comparison +rv/+time/RungekuttaRvBdf.m @ 1169:d02e5b8a0b24 feature/rv

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Rename RungekuttaRV time steppers. Add RungekuttaRVMultiStage time stepper
author Vidar Stiernström <vidar.stiernstrom@it.uu.se>
date Fri, 28 Jun 2019 13:13:17 +0200
parents +rv/+time/RungekuttaExteriorRvBdf.m@3108963cc42c
children a4c00628a39d
comparison
equal deleted inserted replaced
1168:af3c4eb0cbbd 1169:d02e5b8a0b24
1 classdef RungekuttaRvBdf < time.Timestepper
2 properties
3 F % RHS of the ODE
4 k % Time step
5 t % Time point
6 v % Solution vector
7 n % Time level
8 rkScheme % The particular RK scheme used for time integration
9
10
11 % Properties related to the residual viscositys
12 RV % Residual Viscosity operator
13 v_prev % Solution vector at previous time levels, used for the RV evaluation
14 DvDt % Function for computing the time deriative used for the RV evaluation
15 lowerBdfOrder % Orders of the approximation of the time deriative, used for the RV evaluation.
16 % dictates which accuracy the boot-strapping should start from.
17 upperBdfOrder % Orders of the approximation of the time deriative, used for the RV evaluation.
18 % Dictates the order of accuracy used once the boot-strapping is complete.
19
20
21 end
22 methods
23 function obj = RungekuttaRvBdf(F, k, t0, v0, RV, rkOrder, bdfOrders)
24 obj.F = F;
25 obj.k = k;
26 obj.t = t0;
27 obj.v = v0;
28 obj.n = 0;
29 obj.RV = RV;
30 obj.lowerBdfOrder = bdfOrders.lowerBdfOrder;
31 obj.upperBdfOrder = bdfOrders.upperBdfOrder;
32 assert((obj.lowerBdfOrder >= 1) && (obj.upperBdfOrder <= 6));
33 obj.v_prev = [];
34 obj.DvDt = rv.time.BdfDerivative();
35
36 if (rkOrder == 4) % Use specialized RK4 scheme
37 obj.rkScheme = @time.rk.rungekutta_4;
38 else
39 % Extract the coefficients for the specified order
40 % used for the RK updates from the Butcher tableua.
41 [s,a,b,c] = time.rk.butcherTableau(rkOrder);
42 coeffs = struct('s',s,'a',a,'b',b,'c',c);
43 obj.rkScheme = @(v,t,dt,F) time.rk.rungekutta(v, t , dt, F, coeffs);
44 end
45
46 end
47
48 function [v, t] = getV(obj)
49 v = obj.v;
50 t = obj.t;
51 end
52
53 function state = getState(obj)
54 if (size(obj.v_prev,2) >= obj.lowerBdfOrder)
55 dvdt = obj.DvDt.evaluate(obj.v, obj.v_prev, obj.k);
56 [viscosity, Df, firstOrderViscosity, residualViscosity] = obj.RV.evaluate(obj.v, dvdt);
57 else
58 viscosity = zeros(size(obj.v));
59 dvdt = zeros(size(obj.v));
60 Df = zeros(size(obj.v));
61 firstOrderViscosity = zeros(size(obj.v));
62 residualViscosity = zeros(size(obj.v));
63 end
64 state = struct('v', obj.v, 'dvdt', dvdt, 'Df', Df, 'viscosity', viscosity, 'residualViscosity', residualViscosity, 'firstOrderViscosity', firstOrderViscosity, 't', obj.t);
65 end
66
67 function obj = step(obj)
68 nStoredStages = size(obj.v_prev,2);
69
70 %Calculate viscosity for the new time level
71 if (nStoredStages >= obj.lowerBdfOrder)
72 viscosity = obj.RV.evaluateViscosity(obj.v, obj.DvDt.evaluate(obj.v, obj.v_prev, obj.k));
73 else
74 viscosity = zeros(size(obj.v));
75 end
76
77 % Store current time level and update v_prev
78 if (nStoredStages < obj.upperBdfOrder)
79 obj.v_prev = [obj.v, obj.v_prev];
80 else
81 obj.v_prev(:,2:end) = obj.v_prev(:,1:end-1);
82 obj.v_prev(:,1) = obj.v;
83 end
84
85 % Fix the viscosity of the RHS function F
86 m = length(viscosity);
87 F_visc = @(v,t) obj.F(v, t, spdiags(viscosity,0,m,m));
88 obj.v = obj.rkScheme(obj.v, obj.t, obj.k, F_visc);
89 obj.t = obj.t + obj.k;
90 obj.n = obj.n + 1;
91 end
92 end
93 end