Mercurial > repos > public > sbplib

annotate +time/Rungekutta4SecondOrder.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 c6fcee3fcf1b
children
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
0
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
1 classdef Rungekutta4SecondOrder < time.Timestepper
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
2 properties
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
3 F
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
4 k
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
5 t
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
6 w
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
7 m
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
8 D
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
9 E
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
10 S
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
11 M
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
12 C
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
13 n
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
14 end
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
15
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
16
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
17 methods
549
ae905a11e32c Change away from matrix formulation in Rk4 second order form
Jonatan Werpers <jonatan@werpers.com>
parents: 345
diff changeset
18 % Solves u_tt = Du + Eu_t + S by
ae905a11e32c Change away from matrix formulation in Rk4 second order form
Jonatan Werpers <jonatan@werpers.com>
parents: 345
diff changeset
19 % Rewriting on first order form:
ae905a11e32c Change away from matrix formulation in Rk4 second order form
Jonatan Werpers <jonatan@werpers.com>
parents: 345
diff changeset
20 % w_t = M*w + C(t)
ae905a11e32c Change away from matrix formulation in Rk4 second order form
Jonatan Werpers <jonatan@werpers.com>
parents: 345
diff changeset
21 % where
ae905a11e32c Change away from matrix formulation in Rk4 second order form
Jonatan Werpers <jonatan@werpers.com>
parents: 345
diff changeset
22 % M = [
ae905a11e32c Change away from matrix formulation in Rk4 second order form
Jonatan Werpers <jonatan@werpers.com>
parents: 345
diff changeset
23 % 0, I;
ae905a11e32c Change away from matrix formulation in Rk4 second order form
Jonatan Werpers <jonatan@werpers.com>
parents: 345
diff changeset
24 % D, E;
ae905a11e32c Change away from matrix formulation in Rk4 second order form
Jonatan Werpers <jonatan@werpers.com>
parents: 345
diff changeset
25 % ]
ae905a11e32c Change away from matrix formulation in Rk4 second order form
Jonatan Werpers <jonatan@werpers.com>
parents: 345
diff changeset
26 % and
ae905a11e32c Change away from matrix formulation in Rk4 second order form
Jonatan Werpers <jonatan@werpers.com>
parents: 345
diff changeset
27 % C(t) = [
ae905a11e32c Change away from matrix formulation in Rk4 second order form
Jonatan Werpers <jonatan@werpers.com>
parents: 345
diff changeset
28 % 0;
ae905a11e32c Change away from matrix formulation in Rk4 second order form
Jonatan Werpers <jonatan@werpers.com>
parents: 345
diff changeset
29 % S(t)
ae905a11e32c Change away from matrix formulation in Rk4 second order form
Jonatan Werpers <jonatan@werpers.com>
parents: 345
diff changeset
30 % ]
222
e7e73173d44d time.Rungekutta4SecondOrder: Started adding code to allow timedependent D,E,S. Might be buggy.
Jonatan Werpers <jonatan@werpers.com>
parents: 13
diff changeset
31 % D, E, S can either all be constants or all be function handles,
e7e73173d44d time.Rungekutta4SecondOrder: Started adding code to allow timedependent D,E,S. Might be buggy.
Jonatan Werpers <jonatan@werpers.com>
parents: 13
diff changeset
32 % They can also be omitted by setting them equal to the empty matrix.
0
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
33 function obj = Rungekutta4SecondOrder(D, E, S, k, t0, v0, v0t)
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
34 obj.D = D;
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
35 obj.E = E;
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
36 obj.S = S;
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
37 obj.m = length(v0);
13
b18d3d201a71 Fixed initialization of step counter in timesteppers.
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
38 obj.n = 0;
0
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
39
222
e7e73173d44d time.Rungekutta4SecondOrder: Started adding code to allow timedependent D,E,S. Might be buggy.
Jonatan Werpers <jonatan@werpers.com>
parents: 13
diff changeset
40
e7e73173d44d time.Rungekutta4SecondOrder: Started adding code to allow timedependent D,E,S. Might be buggy.
Jonatan Werpers <jonatan@werpers.com>
parents: 13
diff changeset
41 if isa(D, 'function_handle') || isa(E, 'function_handle') || isa(S, 'function_handle')
e7e73173d44d time.Rungekutta4SecondOrder: Started adding code to allow timedependent D,E,S. Might be buggy.
Jonatan Werpers <jonatan@werpers.com>
parents: 13
diff changeset
42 default_arg('D', @(t)sparse(obj.m, obj.m));
e7e73173d44d time.Rungekutta4SecondOrder: Started adding code to allow timedependent D,E,S. Might be buggy.
Jonatan Werpers <jonatan@werpers.com>
parents: 13
diff changeset
43 default_arg('E', @(t)sparse(obj.m, obj.m));
e7e73173d44d time.Rungekutta4SecondOrder: Started adding code to allow timedependent D,E,S. Might be buggy.
Jonatan Werpers <jonatan@werpers.com>
parents: 13
diff changeset
44 default_arg('S', @(t)sparse(obj.m, 1) );
e7e73173d44d time.Rungekutta4SecondOrder: Started adding code to allow timedependent D,E,S. Might be buggy.
Jonatan Werpers <jonatan@werpers.com>
parents: 13
diff changeset
45
e7e73173d44d time.Rungekutta4SecondOrder: Started adding code to allow timedependent D,E,S. Might be buggy.
Jonatan Werpers <jonatan@werpers.com>
parents: 13
diff changeset
46 if ~isa(D, 'function_handle')
e7e73173d44d time.Rungekutta4SecondOrder: Started adding code to allow timedependent D,E,S. Might be buggy.
Jonatan Werpers <jonatan@werpers.com>
parents: 13
diff changeset
47 D = @(t)D;
e7e73173d44d time.Rungekutta4SecondOrder: Started adding code to allow timedependent D,E,S. Might be buggy.
Jonatan Werpers <jonatan@werpers.com>
parents: 13
diff changeset
48 end
345
7b5ef8b89268 Rungekutta bug fixes.
Jonatan Werpers <jonatan@werpers.com>
parents: 222
diff changeset
49 if ~isa(E, 'function_handle')
222
e7e73173d44d time.Rungekutta4SecondOrder: Started adding code to allow timedependent D,E,S. Might be buggy.
Jonatan Werpers <jonatan@werpers.com>
parents: 13
diff changeset
50 E = @(t)E;
e7e73173d44d time.Rungekutta4SecondOrder: Started adding code to allow timedependent D,E,S. Might be buggy.
Jonatan Werpers <jonatan@werpers.com>
parents: 13
diff changeset
51 end
345
7b5ef8b89268 Rungekutta bug fixes.
Jonatan Werpers <jonatan@werpers.com>
parents: 222
diff changeset
52 if ~isa(S, 'function_handle')
222
e7e73173d44d time.Rungekutta4SecondOrder: Started adding code to allow timedependent D,E,S. Might be buggy.
Jonatan Werpers <jonatan@werpers.com>
parents: 13
diff changeset
53 S = @(t)S;
e7e73173d44d time.Rungekutta4SecondOrder: Started adding code to allow timedependent D,E,S. Might be buggy.
Jonatan Werpers <jonatan@werpers.com>
parents: 13
diff changeset
54 end
0
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
55
222
e7e73173d44d time.Rungekutta4SecondOrder: Started adding code to allow timedependent D,E,S. Might be buggy.
Jonatan Werpers <jonatan@werpers.com>
parents: 13
diff changeset
56 obj.k = k;
e7e73173d44d time.Rungekutta4SecondOrder: Started adding code to allow timedependent D,E,S. Might be buggy.
Jonatan Werpers <jonatan@werpers.com>
parents: 13
diff changeset
57 obj.t = t0;
e7e73173d44d time.Rungekutta4SecondOrder: Started adding code to allow timedependent D,E,S. Might be buggy.
Jonatan Werpers <jonatan@werpers.com>
parents: 13
diff changeset
58 obj.w = [v0; v0t];
e7e73173d44d time.Rungekutta4SecondOrder: Started adding code to allow timedependent D,E,S. Might be buggy.
Jonatan Werpers <jonatan@werpers.com>
parents: 13
diff changeset
59
549
ae905a11e32c Change away from matrix formulation in Rk4 second order form
Jonatan Werpers <jonatan@werpers.com>
parents: 345
diff changeset
60 % Avoid matrix formulation because it is VERY slow
ae905a11e32c Change away from matrix formulation in Rk4 second order form
Jonatan Werpers <jonatan@werpers.com>
parents: 345
diff changeset
61 obj.F = @(w,t)[
ae905a11e32c Change away from matrix formulation in Rk4 second order form
Jonatan Werpers <jonatan@werpers.com>
parents: 345
diff changeset
62 w(obj.m+1:end);
ae905a11e32c Change away from matrix formulation in Rk4 second order form
Jonatan Werpers <jonatan@werpers.com>
parents: 345
diff changeset
63 D(t)*w(1:obj.m) + E(t)*w(obj.m+1:end) + S(t);
ae905a11e32c Change away from matrix formulation in Rk4 second order form
Jonatan Werpers <jonatan@werpers.com>
parents: 345
diff changeset
64 ];
0
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
65 else
222
e7e73173d44d time.Rungekutta4SecondOrder: Started adding code to allow timedependent D,E,S. Might be buggy.
Jonatan Werpers <jonatan@werpers.com>
parents: 13
diff changeset
66
e7e73173d44d time.Rungekutta4SecondOrder: Started adding code to allow timedependent D,E,S. Might be buggy.
Jonatan Werpers <jonatan@werpers.com>
parents: 13
diff changeset
67 default_arg('D', sparse(obj.m, obj.m));
e7e73173d44d time.Rungekutta4SecondOrder: Started adding code to allow timedependent D,E,S. Might be buggy.
Jonatan Werpers <jonatan@werpers.com>
parents: 13
diff changeset
68 default_arg('E', sparse(obj.m, obj.m));
e7e73173d44d time.Rungekutta4SecondOrder: Started adding code to allow timedependent D,E,S. Might be buggy.
Jonatan Werpers <jonatan@werpers.com>
parents: 13
diff changeset
69 default_arg('S', sparse(obj.m, 1) );
e7e73173d44d time.Rungekutta4SecondOrder: Started adding code to allow timedependent D,E,S. Might be buggy.
Jonatan Werpers <jonatan@werpers.com>
parents: 13
diff changeset
70
e7e73173d44d time.Rungekutta4SecondOrder: Started adding code to allow timedependent D,E,S. Might be buggy.
Jonatan Werpers <jonatan@werpers.com>
parents: 13
diff changeset
71 I = speye(obj.m);
e7e73173d44d time.Rungekutta4SecondOrder: Started adding code to allow timedependent D,E,S. Might be buggy.
Jonatan Werpers <jonatan@werpers.com>
parents: 13
diff changeset
72 O = sparse(obj.m,obj.m);
0
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
73
222
e7e73173d44d time.Rungekutta4SecondOrder: Started adding code to allow timedependent D,E,S. Might be buggy.
Jonatan Werpers <jonatan@werpers.com>
parents: 13
diff changeset
74 obj.M = [
e7e73173d44d time.Rungekutta4SecondOrder: Started adding code to allow timedependent D,E,S. Might be buggy.
Jonatan Werpers <jonatan@werpers.com>
parents: 13
diff changeset
75 O, I;
e7e73173d44d time.Rungekutta4SecondOrder: Started adding code to allow timedependent D,E,S. Might be buggy.
Jonatan Werpers <jonatan@werpers.com>
parents: 13
diff changeset
76 D, E;
e7e73173d44d time.Rungekutta4SecondOrder: Started adding code to allow timedependent D,E,S. Might be buggy.
Jonatan Werpers <jonatan@werpers.com>
parents: 13
diff changeset
77 ];
e7e73173d44d time.Rungekutta4SecondOrder: Started adding code to allow timedependent D,E,S. Might be buggy.
Jonatan Werpers <jonatan@werpers.com>
parents: 13
diff changeset
78 obj.C = [
e7e73173d44d time.Rungekutta4SecondOrder: Started adding code to allow timedependent D,E,S. Might be buggy.
Jonatan Werpers <jonatan@werpers.com>
parents: 13
diff changeset
79 zeros(obj.m,1);
e7e73173d44d time.Rungekutta4SecondOrder: Started adding code to allow timedependent D,E,S. Might be buggy.
Jonatan Werpers <jonatan@werpers.com>
parents: 13
diff changeset
80 S;
e7e73173d44d time.Rungekutta4SecondOrder: Started adding code to allow timedependent D,E,S. Might be buggy.
Jonatan Werpers <jonatan@werpers.com>
parents: 13
diff changeset
81 ];
0
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
82
222
e7e73173d44d time.Rungekutta4SecondOrder: Started adding code to allow timedependent D,E,S. Might be buggy.
Jonatan Werpers <jonatan@werpers.com>
parents: 13
diff changeset
83 obj.k = k;
e7e73173d44d time.Rungekutta4SecondOrder: Started adding code to allow timedependent D,E,S. Might be buggy.
Jonatan Werpers <jonatan@werpers.com>
parents: 13
diff changeset
84 obj.t = t0;
e7e73173d44d time.Rungekutta4SecondOrder: Started adding code to allow timedependent D,E,S. Might be buggy.
Jonatan Werpers <jonatan@werpers.com>
parents: 13
diff changeset
85 obj.w = [v0; v0t];
e7e73173d44d time.Rungekutta4SecondOrder: Started adding code to allow timedependent D,E,S. Might be buggy.
Jonatan Werpers <jonatan@werpers.com>
parents: 13
diff changeset
86
e7e73173d44d time.Rungekutta4SecondOrder: Started adding code to allow timedependent D,E,S. Might be buggy.
Jonatan Werpers <jonatan@werpers.com>
parents: 13
diff changeset
87 obj.F = @(w,t)(obj.M*w + obj.C);
e7e73173d44d time.Rungekutta4SecondOrder: Started adding code to allow timedependent D,E,S. Might be buggy.
Jonatan Werpers <jonatan@werpers.com>
parents: 13
diff changeset
88 end
0
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
89 end
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
90
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
91 function [v,t] = getV(obj)
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
92 v = obj.w(1:end/2);
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
93 t = obj.t;
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
94 end
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
95
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
96 function [vt,t] = getVt(obj)
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
97 vt = obj.w(end/2+1:end);
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
98 t = obj.t;
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
99 end
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
100
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
101 function obj = step(obj)
846
c6fcee3fcf1b Add generalized RungeKutta and RungeKuttaRV class which extracts its coefficients from a butcher tableau, specified on the scheme.
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents: 549
diff changeset
102 obj.w = time.rk.rungekutta_4(obj.w, obj.t, obj.k, obj.F);
0
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
103 obj.t = obj.t + obj.k;
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
104 obj.n = obj.n + 1;
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
105 end
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
106 end
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
107
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
108
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
109 methods (Static)
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
110 function k = getTimeStep(lambda)
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
111 k = rk4.get_rk4_time_step(lambda);
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
112 end
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
113 end
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
114
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
115 end