Mercurial > repos > public > sbplib

comparison +rv/+time/RungekuttaRvInstage.m @ 1185:abb1b3ab8c23 feature/rv

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Fix incorrect of RungekuttaRvInstage
author Vidar Stiernström <vidar.stiernstrom@it.uu.se>
date Fri, 05 Jul 2019 18:12:10 +0200
parents +rv/+time/RungeKuttaRvInstage.m@9ac86ccfd6a1
children
comparison
equal deleted inserted replaced
1184:ecc605453733 1185:abb1b3ab8c23
1 classdef RungekuttaRvInstage < 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 coeffs % The coefficents used for the RK time integration
9 RV % Residual Viscosity
10 DvDt % Function for computing the time deriative used for the RV evaluation
11 end
12
13 methods
14 function obj = RungekuttaRvInstage(F, k, t0, v0, RV, DvDt, order)
15 obj.F = F;
16 obj.k = k;
17 obj.t = t0;
18 obj.v = v0;
19 obj.n = 0;
20 % Extract the coefficients for the specified order
21 % used for the RK updates from the Butcher tableua.
22 [s,a,b,c] = time.rk.butcherTableau(order);
23 obj.coeffs = struct('s',s,'a',a,'b',b,'c',c);
24 obj.RV = RV;
25 obj.DvDt = DvDt;
26 end
27
28 function [v, t] = getV(obj)
29 v = obj.v;
30 t = obj.t;
31 end
32
33 function state = getState(obj)
34 dvdt = obj.DvDt(obj.v);
35 [viscosity, Df, firstOrderViscosity, residualViscosity] = obj.RV.evaluate(obj.v, dvdt);
36 state = struct('v', obj.v, 'dvdt', dvdt, 'Df', Df, 'viscosity', viscosity, 'residualViscosity', residualViscosity, 'firstOrderViscosity', firstOrderViscosity, 't', obj.t);
37 end
38
39 % Advances the solution vector one time step using the Runge-Kutta method given by
40 % obj.coeffs, updating the Residual Viscosity in each Runge-Kutta stage
41 function obj = step(obj)
42 obj.v = rv.time.rungekuttaRV(obj.v, obj.t, obj.k, obj.F, obj.RV, obj.DvDt, obj.coeffs);
43 obj.t = obj.t + obj.k;
44 obj.n = obj.n + 1;
45 end
46 end
47
48 % Takes one time step of size dt using the rungekutta method
49 % starting from v and where the function F(v,t,RV) gives the
50 % time derivatives. coeffs is a struct holding the RK coefficients
51 % for the specific method. RV is the residual viscosity which is updated
52 % in between the stages and after the updated solution is computed.
53 methods (Static)
54 function v = rungekuttaRV(v, t , dt, F, RV, DvDt, coeffs)
55 % Move one stage outside to avoid branching for updating the
56 % residual inside the loop.
57 k = zeros(length(v), coeffs.s);
58 k(:,1) = F(v,t,RV.evaluateViscosity(v,DvDt(v)));
59
60 % Compute the intermediate stages k
61 for i = 2:coeffs.s
62 u = v;
63 for j = 1:i-1
64 u = u + dt*coeffs.a(i,j)*k(:,j);
65 end
66 k(:,i) = F(u,t+coeffs.c(i)*dt, RV.evaluateViscosity(u,DvDt(u)));
67 end
68
69 % Compute the updated solution as a linear combination
70 % of the intermediate stages.
71 u = v;
72 for i = 1:coeffs.s
73 u = u + dt*coeffs.b(i)*k(:,i);
74 end
75 v = u;
76 end
77
78
79 end