Mercurial > repos > public > sbplib

comparison +rv/+time/RungekuttaExteriorRV.m @ 1013:eb441fbdf379 feature/advectionRV

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Draft implementation of RungeKutta with exterior RV updates - Draft class RungeKuttaExteriorRV which updates the residual using a BDF approximation of the derivtive outside of the RK stages. - Draft class for BDF derivative
author Vidar Stiernström <vidar.stiernstrom@it.uu.se>
date Wed, 05 Dec 2018 15:04:44 +0100
parents
children e547794a9407
comparison
equal deleted inserted replaced
1012:1e437c9e5132 1013:eb441fbdf379
1 classdef RungekuttaExteriorRV < 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 v_prev % Solution vector at previous time levels, used for the RV update
11 DvDt % Function for computing the time deriative used for the RV update
12 end
13 methods
14
15 function obj = RungekuttaExteriorRV(F, k, t0, v0, RV, rkOrder, bdfOrder)
16 obj.F = F;
17 obj.k = k;
18 obj.t = t0;
19 obj.v = v0;
20 obj.n = 0;
21 % Extract the coefficients for the specified rkOrder
22 % used for the RK updates from the Butcher tableua.
23 [s,a,b,c] = time.rk.butcherTableau(rkOrder);
24 obj.coeffs = struct('s',s,'a',a,'b',b,'c',c);
25
26 obj.RV = RV;
27 % TODO: Initialize v_prev.
28 obj.v_prev = zeros(length(obj.v),bdfOrder);
29 % TBD: For cases where h~k we could probably use rkOrder-2 here.
30 obj.DvDt = rv.time.BDFDerivative(bdfOrder,k);
31 end
32
33 function [v, t] = getV(obj)
34 v = obj.v;
35 t = obj.t;
36 end
37
38 function state = getState(obj)
39 [residual, u_t, grad_f] = obj.RV.getResidual();
40 state = struct('v', obj.v, 'residual', residual, 'u_t', u_t, 'grad_f', grad_f, 'viscosity', obj.RV.getViscosity(), 't', obj.t);
41 end
42
43 function obj = step(obj)
44 % Store current time level
45 obj.v_prev(:,2:end) = obj.v_prev(:,1:end-1);
46 obj.v_prev(:,1) = obj.v;
47 % Fix the viscosity of the RHS function F
48 F_visc = @(v,t)F(v,t,obj.RV.getViscosity()); %Remove state in RV?
49 obj.v = time.rk.rungekutta(obj.v, obj.t, obj.k, F_visc, obj.RV, obj.coeffs);
50 % Calculate dvdt and update RV
51 dvdt = obj.DvDt.evaluate(v, v_prev, lenght(v_prev));
52 RV.update(v,dvdt);
53 obj.t = obj.t + obj.k;
54 obj.n = obj.n + 1;
55 end
56 end
57 end