Mercurial > repos > public > sbplib

comparison +rv/+time/RungekuttaInteriorRV.m @ 1017:2d7c1333bd6c feature/advectionRV

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Add support for using the ODE to approximate the time derivative in the residual
author Vidar Stiernström <vidar.stiernstrom@it.uu.se>
date Tue, 11 Dec 2018 16:29:21 +0100
parents 1e437c9e5132
children
comparison
equal deleted inserted replaced
1016:4b42999874c0 1017:2d7c1333bd6c
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 viscosity % Viscosity vector
11 DvDt % Function for computing the time deriative used for the RV evaluation
12
13 residual
14 dvdt
15 Df
10 end 16 end
11 17
12 methods 18 methods
13 function obj = RungekuttaInteriorRV(F, k, t0, v0, RV, order) 19 function obj = RungekuttaInteriorRV(F, k, t0, v0, RV, DvDt, order)
14 obj.F = F; 20 obj.F = F;
15 obj.k = k; 21 obj.k = k;
16 obj.t = t0; 22 obj.t = t0;
17 obj.v = v0; 23 obj.v = v0;
18 obj.n = 0; 24 obj.n = 0;
19 obj.RV = RV;
20 % Extract the coefficients for the specified order 25 % Extract the coefficients for the specified order
21 % used for the RK updates from the Butcher tableua. 26 % used for the RK updates from the Butcher tableua.
22 [s,a,b,c] = time.rk.butcherTableau(order); 27 [s,a,b,c] = time.rk.butcherTableau(order);
23 obj.coeffs = struct('s',s,'a',a,'b',b,'c',c); 28 obj.coeffs = struct('s',s,'a',a,'b',b,'c',c);
29 obj.RV = RV;
30 obj.DvDt = DvDt;
31 obj.dvdt = obj.DvDt(obj.v);
32 [obj.viscosity, obj.Df] = obj.RV.evaluate(obj.v,obj.dvdt);
33 obj.residual = obj.dvdt + obj.Df;
24 end 34 end
25 35
26 function [v, t] = getV(obj) 36 function [v, t] = getV(obj)
27 v = obj.v; 37 v = obj.v;
28 t = obj.t; 38 t = obj.t;
29 end 39 end
30 40
31 function state = getState(obj) 41 function state = getState(obj)
32 [residual, u_t, grad_f] = obj.RV.getResidual(); 42 state = struct('v', obj.v, 'residual', obj.residual, 'dvdt', obj.dvdt, 'Df', obj.Df, 'viscosity', obj.viscosity, 't', obj.t);
33 state = struct('v', obj.v, 'residual', residual, 'u_t', u_t, 'grad_f', grad_f, 'viscosity', obj.RV.getViscosity(), 't', obj.t);
34 end 43 end
35 44
36 function obj = step(obj) 45 function obj = step(obj)
37 obj.v = rv.time.rungekuttaRV(obj.v, obj.t, obj.k, obj.F, obj.RV, obj.coeffs); 46 obj.v = rv.time.rungekuttaRV(obj.v, obj.t, obj.k, obj.F, obj.RV, obj.DvDt, obj.coeffs);
38 obj.t = obj.t + obj.k; 47 obj.t = obj.t + obj.k;
39 obj.n = obj.n + 1; 48 obj.n = obj.n + 1;
49 obj.dvdt = obj.DvDt(obj.v);
50 [obj.viscosity, obj.Df] = obj.RV.evaluate(obj.v,obj.dvdt);
51 obj.residual = obj.dvdt + obj.Df;
40 end 52 end
41 end 53 end
42 end 54 end