Mercurial > repos > public > sbplib
annotate +rv/+time/RungekuttaRv.m @ 1182:f35ff0861d5a feature/rv
Add standard RungekuttaRv
author | Vidar Stiernström <vidar.stiernstrom@it.uu.se> |
---|---|
date | Fri, 05 Jul 2019 17:47:13 +0200 |
parents | |
children |
rev | line source |
---|---|
1182
f35ff0861d5a
Add standard RungekuttaRv
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff
changeset
|
1 classdef RungekuttaRv < time.Timestepper |
f35ff0861d5a
Add standard RungekuttaRv
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff
changeset
|
2 properties |
f35ff0861d5a
Add standard RungekuttaRv
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff
changeset
|
3 F % RHS of the ODE |
f35ff0861d5a
Add standard RungekuttaRv
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff
changeset
|
4 k % Time step |
f35ff0861d5a
Add standard RungekuttaRv
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff
changeset
|
5 t % Time point |
f35ff0861d5a
Add standard RungekuttaRv
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff
changeset
|
6 v % Solution vector |
f35ff0861d5a
Add standard RungekuttaRv
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff
changeset
|
7 n % Time level |
f35ff0861d5a
Add standard RungekuttaRv
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff
changeset
|
8 rkScheme % The particular RK scheme used for time integration |
f35ff0861d5a
Add standard RungekuttaRv
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff
changeset
|
9 RV % Residual Viscosity operator |
f35ff0861d5a
Add standard RungekuttaRv
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff
changeset
|
10 DvDt % Function for computing the time deriative used for the RV evaluation |
f35ff0861d5a
Add standard RungekuttaRv
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff
changeset
|
11 end |
f35ff0861d5a
Add standard RungekuttaRv
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff
changeset
|
12 methods |
f35ff0861d5a
Add standard RungekuttaRv
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff
changeset
|
13 |
f35ff0861d5a
Add standard RungekuttaRv
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff
changeset
|
14 function obj = RungekuttaRv(F, k, t0, v0, RV, DvDt, order) |
f35ff0861d5a
Add standard RungekuttaRv
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff
changeset
|
15 obj.F = F; |
f35ff0861d5a
Add standard RungekuttaRv
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff
changeset
|
16 obj.k = k; |
f35ff0861d5a
Add standard RungekuttaRv
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff
changeset
|
17 obj.t = t0; |
f35ff0861d5a
Add standard RungekuttaRv
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff
changeset
|
18 obj.v = v0; |
f35ff0861d5a
Add standard RungekuttaRv
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff
changeset
|
19 obj.n = 0; |
f35ff0861d5a
Add standard RungekuttaRv
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff
changeset
|
20 |
f35ff0861d5a
Add standard RungekuttaRv
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff
changeset
|
21 if (order == 4) % Use specialized RK4 scheme |
f35ff0861d5a
Add standard RungekuttaRv
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff
changeset
|
22 obj.rkScheme = @time.rk.rungekutta_4; |
f35ff0861d5a
Add standard RungekuttaRv
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff
changeset
|
23 else |
f35ff0861d5a
Add standard RungekuttaRv
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff
changeset
|
24 % Extract the coefficients for the specified order |
f35ff0861d5a
Add standard RungekuttaRv
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff
changeset
|
25 % used for the RK updates from the Butcher tableua. |
f35ff0861d5a
Add standard RungekuttaRv
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff
changeset
|
26 [s,a,b,c] = time.rk.butcherTableau(order); |
f35ff0861d5a
Add standard RungekuttaRv
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff
changeset
|
27 coeffs = struct('s',s,'a',a,'b',b,'c',c); |
f35ff0861d5a
Add standard RungekuttaRv
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff
changeset
|
28 obj.rkScheme = @(v,t,dt,F) time.rk.rungekutta(v, t , dt, F, coeffs); |
f35ff0861d5a
Add standard RungekuttaRv
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff
changeset
|
29 end |
f35ff0861d5a
Add standard RungekuttaRv
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff
changeset
|
30 |
f35ff0861d5a
Add standard RungekuttaRv
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff
changeset
|
31 obj.RV = RV; |
f35ff0861d5a
Add standard RungekuttaRv
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff
changeset
|
32 obj.DvDt = DvDt; |
f35ff0861d5a
Add standard RungekuttaRv
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff
changeset
|
33 end |
f35ff0861d5a
Add standard RungekuttaRv
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff
changeset
|
34 |
f35ff0861d5a
Add standard RungekuttaRv
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff
changeset
|
35 function [v, t] = getV(obj) |
f35ff0861d5a
Add standard RungekuttaRv
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff
changeset
|
36 v = obj.v; |
f35ff0861d5a
Add standard RungekuttaRv
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff
changeset
|
37 t = obj.t; |
f35ff0861d5a
Add standard RungekuttaRv
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff
changeset
|
38 end |
f35ff0861d5a
Add standard RungekuttaRv
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff
changeset
|
39 |
f35ff0861d5a
Add standard RungekuttaRv
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff
changeset
|
40 function state = getState(obj) |
f35ff0861d5a
Add standard RungekuttaRv
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff
changeset
|
41 dvdt = obj.DvDt(obj.v); |
f35ff0861d5a
Add standard RungekuttaRv
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff
changeset
|
42 [viscosity, Df, firstOrderViscosity, residualViscosity] = obj.RV.evaluate(obj.v, dvdt); |
f35ff0861d5a
Add standard RungekuttaRv
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff
changeset
|
43 state = struct('v', obj.v, 'dvdt', dvdt, 'Df', Df, 'viscosity', viscosity, 'residualViscosity', residualViscosity, 'firstOrderViscosity', firstOrderViscosity, 't', obj.t); |
f35ff0861d5a
Add standard RungekuttaRv
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff
changeset
|
44 end |
f35ff0861d5a
Add standard RungekuttaRv
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff
changeset
|
45 |
f35ff0861d5a
Add standard RungekuttaRv
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff
changeset
|
46 % Advances the solution vector one time step using the Runge-Kutta method given by |
f35ff0861d5a
Add standard RungekuttaRv
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff
changeset
|
47 % obj.coeffs, using a fixed residual viscosity for the Runge-Kutta substeps |
f35ff0861d5a
Add standard RungekuttaRv
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff
changeset
|
48 function obj = step(obj) |
f35ff0861d5a
Add standard RungekuttaRv
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff
changeset
|
49 % Fix the viscosity of the stabilized RHS |
f35ff0861d5a
Add standard RungekuttaRv
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff
changeset
|
50 m = length(obj.v); |
f35ff0861d5a
Add standard RungekuttaRv
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff
changeset
|
51 F = @(v,t) obj.F(v,t,spdiags(obj.RV.evaluateViscosity(obj.v, obj.DvDt(obj.v)),0,m,m)); |
f35ff0861d5a
Add standard RungekuttaRv
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff
changeset
|
52 obj.v = obj.rkScheme(obj.v, obj.t, obj.k, F); |
f35ff0861d5a
Add standard RungekuttaRv
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff
changeset
|
53 obj.t = obj.t + obj.k; |
f35ff0861d5a
Add standard RungekuttaRv
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff
changeset
|
54 obj.n = obj.n + 1; |
f35ff0861d5a
Add standard RungekuttaRv
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff
changeset
|
55 end |
f35ff0861d5a
Add standard RungekuttaRv
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff
changeset
|
56 end |
f35ff0861d5a
Add standard RungekuttaRv
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff
changeset
|
57 end |