Mercurial > repos > public > sbplib

comparison +time/ExplicitRungeKuttaDiscreteData.m @ 1331:60c875c18de3 feature/D2_boundary_opt

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Merge with feature/poroelastic for Elastic schemes
author Vidar Stiernström <vidar.stiernstrom@it.uu.se>
date Thu, 10 Mar 2022 16:54:26 +0100
parents 7963a9cab637
children 0aefcb30cab4
comparison
equal deleted inserted replaced
1330:855871e0b852 1331:60c875c18de3
1 classdef ExplicitRungeKuttaDiscreteData < time.Timestepper
2 properties
3 D
4 S % Function handle for time-dependent data
5 data % Matrix of data vectors, one column per stage
6 k
7 t
8 v
9 m
10 n
11 order
12 a, b, c, s % Butcher tableau
13 K % Stage rates
14 U % Stage approximations
15 T % Stage times
16 end
17
18
19 methods
20 function obj = ExplicitRungeKuttaDiscreteData(D, S, data, k, t0, v0, order)
21 default_arg('order', 4);
22 default_arg('S', []);
23 default_arg('data', []);
24
25 obj.D = D;
26 obj.S = S;
27 obj.k = k;
28 obj.t = t0;
29 obj.v = v0;
30 obj.m = length(v0);
31 obj.n = 0;
32 obj.order = order;
33 obj.data = data;
34
35 switch order
36 case 4
37 [obj.a, obj.b, obj.c, obj.s] = time.rkparameters.rk4();
38 otherwise
39 error('That RK method is not available');
40 end
41
42 obj.K = zeros(obj.m, obj.s);
43 obj.U = zeros(obj.m, obj.s);
44
45 end
46
47 function [v,t,U,T,K] = getV(obj)
48 v = obj.v;
49 t = obj.t;
50 U = obj.U; % Stage approximations in previous time step.
51 T = obj.T; % Stage times in previous time step.
52 K = obj.K; % Stage rates in previous time step.
53 end
54
55 function [a,b,c,s] = getTableau(obj)
56 a = obj.a;
57 b = obj.b;
58 c = obj.c;
59 s = obj.s;
60 end
61
62 % Returns quadrature weights for stages in one time step
63 function quadWeights = getTimeStepQuadrature(obj)
64 [~, b] = obj.getTableau();
65 quadWeights = obj.k*b;
66 end
67
68 function obj = step(obj)
69 v = obj.v;
70 a = obj.a;
71 b = obj.b;
72 c = obj.c;
73 s = obj.s;
74 S = obj.S;
75 dt = obj.k;
76 K = obj.K;
77 U = obj.U;
78 D = obj.D;
79 data = obj.data;
80
81 for i = 1:s
82 U(:,i) = v;
83 for j = 1:i-1
84 U(:,i) = U(:,i) + dt*a(i,j)*K(:,j);
85 end
86
87 K(:,i) = D*U(:,i);
88 obj.T(i) = obj.t + c(i)*dt;
89
90 % Data from continuous function and discrete time-points.
91 if ~isempty(S)
92 K(:,i) = K(:,i) + S(obj.T(i));
93 end
94 if ~isempty(data)
95 K(:,i) = K(:,i) + data(:,obj.n*s + i);
96 end
97
98 end
99
100 obj.v = v + dt*K*b;
101 obj.t = obj.t + dt;
102 obj.n = obj.n + 1;
103 obj.U = U;
104 obj.K = K;
105 end
106 end
107
108
109 methods (Static)
110 function k = getTimeStep(lambda, order)
111 default_arg('order', 4);
112 switch order
113 case 4
114 k = time.rk4.get_rk4_time_step(lambda);
115 otherwise
116 error('Time-step function not available for this order');
117 end
118 end
119 end
120
121 end