Mercurial > repos > public > sbplib
changeset 1013:eb441fbdf379 feature/advectionRV
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 | 1e437c9e5132 |
| children | e547794a9407 |
| files | +rv/+time/BDFDerivative.m +rv/+time/RungekuttaExteriorRV.m |
| diffstat | 2 files changed, 90 insertions(+), 0 deletions(-) [+] |
line wrap: on
line diff
diff -r 1e437c9e5132 -r eb441fbdf379 +rv/+time/BDFDerivative.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+rv/+time/BDFDerivative.m Wed Dec 05 15:04:44 2018 +0100 @@ -0,0 +1,33 @@ +class BDFDerivative < handle + properties + dt + coefficents + end + methods + function obj = BDFDerivative(order, dt) + obj.coefficents = obj.getBDFCoefficients(order); + obj.dt = dt; + end + + function DvDt = evaluate(v, v_prev) + DvDt = (obj.coefficients(1)*v - sum(obj.coefficients(2:end).*v_prev),2)/obj.dt; + end + + function c = getBDFCoefficients(obj, order) + switch order + case 1 + c(1,1) = 1; c(1,2) = 1; + case 2 + c(1,1) = 3/2; c(1,2) = 4/2; c(1,3) = -1/2; + case 3 + c(1,1) = 11/6; c(1,2) = 18/6; c(1,3) = -9/6; c(1,4) = 2/6; + case 4 + c(1,1) = 25/12; c(1,2) = 48/12; c(1,3) = -36/12; c(1,4) = 16/12; c(1,5) = -3/12; + case 5 + c(1,1) = 137/60; c(1,2) = 300/60; c(1,3) = -300/60; c(1,4) = 200/60; c(1,5) = -75/60; c(1,6) = 12/60; + case 6 + c(1,1) = 147/60; c(1,2) = 360/60; c(1,3) = -450/60; c(1,4) = 400/60; c(1,5) = -225/60; c(1,6) = 72/60; c(1,7) = -10/60; + assert(length(c) == (order+1)); + end + end +end
diff -r 1e437c9e5132 -r eb441fbdf379 +rv/+time/RungekuttaExteriorRV.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+rv/+time/RungekuttaExteriorRV.m Wed Dec 05 15:04:44 2018 +0100 @@ -0,0 +1,57 @@ +classdef RungekuttaExteriorRV < time.Timestepper + properties + F % RHS of the ODE + k % Time step + t % Time point + v % Solution vector + n % Time level + coeffs % The coefficents used for the RK time integration + RV % Residual Viscosity + v_prev % Solution vector at previous time levels, used for the RV update + DvDt % Function for computing the time deriative used for the RV update + end + methods + + function obj = RungekuttaExteriorRV(F, k, t0, v0, RV, rkOrder, bdfOrder) + obj.F = F; + obj.k = k; + obj.t = t0; + obj.v = v0; + obj.n = 0; + % Extract the coefficients for the specified rkOrder + % used for the RK updates from the Butcher tableua. + [s,a,b,c] = time.rk.butcherTableau(rkOrder); + obj.coeffs = struct('s',s,'a',a,'b',b,'c',c); + + obj.RV = RV; + % TODO: Initialize v_prev. + obj.v_prev = zeros(length(obj.v),bdfOrder); + % TBD: For cases where h~k we could probably use rkOrder-2 here. + obj.DvDt = rv.time.BDFDerivative(bdfOrder,k); + end + + function [v, t] = getV(obj) + v = obj.v; + t = obj.t; + end + + function state = getState(obj) + [residual, u_t, grad_f] = obj.RV.getResidual(); + state = struct('v', obj.v, 'residual', residual, 'u_t', u_t, 'grad_f', grad_f, 'viscosity', obj.RV.getViscosity(), 't', obj.t); + end + + function obj = step(obj) + % Store current time level + obj.v_prev(:,2:end) = obj.v_prev(:,1:end-1); + obj.v_prev(:,1) = obj.v; + % Fix the viscosity of the RHS function F + F_visc = @(v,t)F(v,t,obj.RV.getViscosity()); %Remove state in RV? + obj.v = time.rk.rungekutta(obj.v, obj.t, obj.k, F_visc, obj.RV, obj.coeffs); + % Calculate dvdt and update RV + dvdt = obj.DvDt.evaluate(v, v_prev, lenght(v_prev)); + RV.update(v,dvdt); + obj.t = obj.t + obj.k; + obj.n = obj.n + 1; + end + end +end \ No newline at end of file