Mercurial > repos > public > sbplib

comparison +rv/+time/RungekuttaInteriorRv.m @ 1031:2ef20d00b386 feature/advectionRV

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
For easier comparison, return both the first order and residual viscosity when evaluating the residual. Add the first order and residual viscosity to the state of the RungekuttaRV time steppers
author Vidar Stiernström <vidar.stiernstrom@it.uu.se>
date Thu, 17 Jan 2019 10:25:06 +0100
parents 78c75c95b7dd
children 010bb2677230
comparison
equal deleted inserted replaced
1030:78c75c95b7dd 1031:2ef20d00b386
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 10 DvDt % Function for computing the time deriative used for the RV evaluation
11 DvDt % Function for computing the time deriative used for the RV evaluation
12 11
13 residual 12 % Convenience properties. Only for plotting
14 dvdt 13 viscosity % Total viscosity
15 Df 14 residualViscosity % Residual viscosity
15 firstOrderViscosity % first order viscosity
16 dvdt % Evaluated time derivative in residual
17 Df % Evaluated flux in residual
16 end 18 end
17 19
18 methods 20 methods
19 function obj = RungekuttaInteriorRv(F, k, t0, v0, RV, DvDt, order) 21 function obj = RungekuttaInteriorRv(F, k, t0, v0, RV, DvDt, order)
20 obj.F = F; 22 obj.F = F;
27 [s,a,b,c] = time.rk.butcherTableau(order); 29 [s,a,b,c] = time.rk.butcherTableau(order);
28 obj.coeffs = struct('s',s,'a',a,'b',b,'c',c); 30 obj.coeffs = struct('s',s,'a',a,'b',b,'c',c);
29 obj.RV = RV; 31 obj.RV = RV;
30 obj.DvDt = DvDt; 32 obj.DvDt = DvDt;
31 obj.dvdt = obj.DvDt(obj.v); 33 obj.dvdt = obj.DvDt(obj.v);
32 [obj.viscosity, obj.Df] = obj.RV.evaluate(obj.v,obj.dvdt); 34 [obj.viscosity, obj.Df, obj.firstOrderViscosity, obj.residualViscosity] = obj.RV.evaluate(obj.v,obj.dvdt);
33 obj.residual = obj.dvdt + obj.Df;
34 end 35 end
35 36
36 function [v, t] = getV(obj) 37 function [v, t] = getV(obj)
37 v = obj.v; 38 v = obj.v;
38 t = obj.t; 39 t = obj.t;
39 end 40 end
40 41
41 function state = getState(obj) 42 function state = getState(obj)
42 state = struct('v', obj.v, 'residual', obj.residual, 'dvdt', obj.dvdt, 'Df', obj.Df, 'viscosity', obj.viscosity, 't', obj.t); 43 state = struct('v', obj.v, 'dvdt', obj.dvdt, 'Df', obj.Df, 'viscosity', obj.viscosity, 'residualViscosity', obj.residualViscosity, 'firstOrderViscosity', obj.firstOrderViscosity, 't', obj.t);
43 end 44 end
44 45
45 function obj = step(obj) 46 function obj = step(obj)
46 obj.v = rv.time.rungekuttaRV(obj.v, obj.t, obj.k, obj.F, obj.RV, obj.DvDt, obj.coeffs); 47 obj.v = rv.time.rungekuttaRV(obj.v, obj.t, obj.k, obj.F, obj.RV, obj.DvDt, obj.coeffs);
47 obj.t = obj.t + obj.k; 48 obj.t = obj.t + obj.k;
48 obj.n = obj.n + 1; 49 obj.n = obj.n + 1;
49 obj.dvdt = obj.DvDt(obj.v); 50 obj.dvdt = obj.DvDt(obj.v);
50 [obj.viscosity, obj.Df] = obj.RV.evaluate(obj.v,obj.dvdt); 51 [obj.viscosity, obj.Df, obj.firstOrderViscosity, obj.residualViscosity] = obj.RV.evaluate(obj.v,obj.dvdt);
51 obj.residual = obj.dvdt + obj.Df;
52 end 52 end
53 end 53 end
54 end 54 end