Mercurial > repos > public > sbplib

annotate +time/ExplicitRungeKuttaSecondOrderDiscreteData.m @ 1347:ac54767ae1fb feature/poroelastic tip

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Add interface, not fully compatible.
author Martin Almquist <martin.almquist@it.uu.se>
date Tue, 30 Apr 2024 14:58:35 +0200
parents 335b8b1366d4
children
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
858
335b8b1366d4 Add RK for second order and discrete data.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
1 classdef ExplicitRungeKuttaSecondOrderDiscreteData < time.Timestepper
335b8b1366d4 Add RK for second order and discrete data.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
2 properties
335b8b1366d4 Add RK for second order and discrete data.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
3 k
335b8b1366d4 Add RK for second order and discrete data.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
4 t
335b8b1366d4 Add RK for second order and discrete data.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
5 w
335b8b1366d4 Add RK for second order and discrete data.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
6 m
335b8b1366d4 Add RK for second order and discrete data.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
7 D
335b8b1366d4 Add RK for second order and discrete data.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
8 E
335b8b1366d4 Add RK for second order and discrete data.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
9 M
335b8b1366d4 Add RK for second order and discrete data.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
10 C_cont % Continuous part (function handle) of forcing on first order form.
335b8b1366d4 Add RK for second order and discrete data.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
11 C_discr% Discrete part (matrix) of forcing on first order form.
335b8b1366d4 Add RK for second order and discrete data.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
12 n
335b8b1366d4 Add RK for second order and discrete data.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
13 order
335b8b1366d4 Add RK for second order and discrete data.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
14 tsImplementation % Time stepper object, RK first order form,
335b8b1366d4 Add RK for second order and discrete data.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
15 % which this wraps around.
335b8b1366d4 Add RK for second order and discrete data.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
16 end
335b8b1366d4 Add RK for second order and discrete data.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
17
335b8b1366d4 Add RK for second order and discrete data.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
18
335b8b1366d4 Add RK for second order and discrete data.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
19 methods
335b8b1366d4 Add RK for second order and discrete data.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
20 % Solves u_tt = Du + Eu_t + S by
335b8b1366d4 Add RK for second order and discrete data.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
21 % Rewriting on first order form:
335b8b1366d4 Add RK for second order and discrete data.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
22 % w_t = M*w + C(t)
335b8b1366d4 Add RK for second order and discrete data.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
23 % where
335b8b1366d4 Add RK for second order and discrete data.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
24 % M = [
335b8b1366d4 Add RK for second order and discrete data.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
25 % 0, I;
335b8b1366d4 Add RK for second order and discrete data.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
26 % D, E;
335b8b1366d4 Add RK for second order and discrete data.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
27 % ]
335b8b1366d4 Add RK for second order and discrete data.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
28 % and
335b8b1366d4 Add RK for second order and discrete data.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
29 % C(t) = [
335b8b1366d4 Add RK for second order and discrete data.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
30 % 0;
335b8b1366d4 Add RK for second order and discrete data.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
31 % S(t)
335b8b1366d4 Add RK for second order and discrete data.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
32 % ]
335b8b1366d4 Add RK for second order and discrete data.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
33 % D, E, should be matrices (or empty for zero)
335b8b1366d4 Add RK for second order and discrete data.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
34 % They can also be omitted by setting them equal to the empty matrix.
335b8b1366d4 Add RK for second order and discrete data.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
35 % S = S_cont + S_discr, where S_cont is a function handle
335b8b1366d4 Add RK for second order and discrete data.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
36 % S_discr a matrix of data vectors, one column per stage.
335b8b1366d4 Add RK for second order and discrete data.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
37 function obj = ExplicitRungeKuttaSecondOrderDiscreteData(D, E, S_cont, S_discr, k, t0, v0, v0t, order)
335b8b1366d4 Add RK for second order and discrete data.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
38 default_arg('order', 4);
335b8b1366d4 Add RK for second order and discrete data.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
39 default_arg('S_cont', []);
335b8b1366d4 Add RK for second order and discrete data.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
40 default_arg('S_discr', []);
335b8b1366d4 Add RK for second order and discrete data.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
41 obj.D = D;
335b8b1366d4 Add RK for second order and discrete data.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
42 obj.E = E;
335b8b1366d4 Add RK for second order and discrete data.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
43 obj.m = length(v0);
335b8b1366d4 Add RK for second order and discrete data.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
44 obj.n = 0;
335b8b1366d4 Add RK for second order and discrete data.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
45
335b8b1366d4 Add RK for second order and discrete data.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
46 default_arg('D', sparse(obj.m, obj.m) );
335b8b1366d4 Add RK for second order and discrete data.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
47 default_arg('E', sparse(obj.m, obj.m) );
335b8b1366d4 Add RK for second order and discrete data.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
48
335b8b1366d4 Add RK for second order and discrete data.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
49 obj.k = k;
335b8b1366d4 Add RK for second order and discrete data.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
50 obj.t = t0;
335b8b1366d4 Add RK for second order and discrete data.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
51 obj.w = [v0; v0t];
335b8b1366d4 Add RK for second order and discrete data.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
52
335b8b1366d4 Add RK for second order and discrete data.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
53 I = speye(obj.m);
335b8b1366d4 Add RK for second order and discrete data.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
54 O = sparse(obj.m,obj.m);
335b8b1366d4 Add RK for second order and discrete data.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
55
335b8b1366d4 Add RK for second order and discrete data.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
56 obj.M = [
335b8b1366d4 Add RK for second order and discrete data.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
57 O, I;
335b8b1366d4 Add RK for second order and discrete data.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
58 D, E;
335b8b1366d4 Add RK for second order and discrete data.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
59 ];
335b8b1366d4 Add RK for second order and discrete data.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
60
335b8b1366d4 Add RK for second order and discrete data.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
61 % Build C_cont
335b8b1366d4 Add RK for second order and discrete data.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
62 if ~isempty(S_cont)
335b8b1366d4 Add RK for second order and discrete data.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
63 obj.C_cont = @(t)[
335b8b1366d4 Add RK for second order and discrete data.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
64 sparse(obj.m,1);
335b8b1366d4 Add RK for second order and discrete data.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
65 S_cont(t)
335b8b1366d4 Add RK for second order and discrete data.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
66 ];
335b8b1366d4 Add RK for second order and discrete data.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
67 else
335b8b1366d4 Add RK for second order and discrete data.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
68 obj.C_cont = [];
335b8b1366d4 Add RK for second order and discrete data.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
69 end
335b8b1366d4 Add RK for second order and discrete data.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
70
335b8b1366d4 Add RK for second order and discrete data.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
71 % Build C_discr
335b8b1366d4 Add RK for second order and discrete data.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
72 if ~isempty(S_discr)
335b8b1366d4 Add RK for second order and discrete data.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
73 [~, nt] = size(S_discr);
335b8b1366d4 Add RK for second order and discrete data.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
74 obj.C_discr = [sparse(obj.m, nt);
335b8b1366d4 Add RK for second order and discrete data.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
75 S_discr
335b8b1366d4 Add RK for second order and discrete data.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
76 ];
335b8b1366d4 Add RK for second order and discrete data.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
77 else
335b8b1366d4 Add RK for second order and discrete data.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
78 obj.C_discr = [];
335b8b1366d4 Add RK for second order and discrete data.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
79 end
335b8b1366d4 Add RK for second order and discrete data.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
80 obj.tsImplementation = time.ExplicitRungeKuttaDiscreteData(obj.M, obj.C_cont, obj.C_discr,...
335b8b1366d4 Add RK for second order and discrete data.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
81 k, obj.t, obj.w, order);
335b8b1366d4 Add RK for second order and discrete data.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
82 end
335b8b1366d4 Add RK for second order and discrete data.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
83
335b8b1366d4 Add RK for second order and discrete data.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
84 function [v,t,U,T,K] = getV(obj)
335b8b1366d4 Add RK for second order and discrete data.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
85 [w,t,U,T,K] = obj.tsImplementation.getV();
335b8b1366d4 Add RK for second order and discrete data.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
86
335b8b1366d4 Add RK for second order and discrete data.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
87 v = w(1:end/2);
335b8b1366d4 Add RK for second order and discrete data.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
88 U = U(1:end/2, :); % Stage approximations in previous time step.
335b8b1366d4 Add RK for second order and discrete data.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
89 K = K(1:end/2, :); % Stage rates in previous time step.
335b8b1366d4 Add RK for second order and discrete data.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
90 % T: Stage times in previous time step.
335b8b1366d4 Add RK for second order and discrete data.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
91 end
335b8b1366d4 Add RK for second order and discrete data.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
92
335b8b1366d4 Add RK for second order and discrete data.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
93 function [vt,t,U,T,K] = getVt(obj)
335b8b1366d4 Add RK for second order and discrete data.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
94 [w,t,U,T,K] = obj.tsImplementation.getV();
335b8b1366d4 Add RK for second order and discrete data.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
95
335b8b1366d4 Add RK for second order and discrete data.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
96 vt = w(end/2+1:end);
335b8b1366d4 Add RK for second order and discrete data.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
97 U = U(end/2+1:end, :); % Stage approximations in previous time step.
335b8b1366d4 Add RK for second order and discrete data.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
98 K = K(end/2+1:end, :); % Stage rates in previous time step.
335b8b1366d4 Add RK for second order and discrete data.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
99 % T: Stage times in previous time step.
335b8b1366d4 Add RK for second order and discrete data.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
100 end
335b8b1366d4 Add RK for second order and discrete data.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
101
335b8b1366d4 Add RK for second order and discrete data.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
102 function [a,b,c,s] = getTableau(obj)
335b8b1366d4 Add RK for second order and discrete data.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
103 [a,b,c,s] = obj.tsImplementation.getTableau();
335b8b1366d4 Add RK for second order and discrete data.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
104 end
335b8b1366d4 Add RK for second order and discrete data.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
105
335b8b1366d4 Add RK for second order and discrete data.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
106 % Returns quadrature weights for stages in one time step
335b8b1366d4 Add RK for second order and discrete data.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
107 function quadWeights = getTimeStepQuadrature(obj)
335b8b1366d4 Add RK for second order and discrete data.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
108 [~, b] = obj.getTableau();
335b8b1366d4 Add RK for second order and discrete data.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
109 quadWeights = obj.k*b;
335b8b1366d4 Add RK for second order and discrete data.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
110 end
335b8b1366d4 Add RK for second order and discrete data.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
111
335b8b1366d4 Add RK for second order and discrete data.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
112 % Use RK for first order form to step
335b8b1366d4 Add RK for second order and discrete data.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
113 function obj = step(obj)
335b8b1366d4 Add RK for second order and discrete data.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
114 obj.tsImplementation.step();
335b8b1366d4 Add RK for second order and discrete data.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
115 [v, t] = obj.tsImplementation.getV();
335b8b1366d4 Add RK for second order and discrete data.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
116 obj.w = v;
335b8b1366d4 Add RK for second order and discrete data.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
117 obj.t = t;
335b8b1366d4 Add RK for second order and discrete data.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
118 obj.n = obj.n + 1;
335b8b1366d4 Add RK for second order and discrete data.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
119 end
335b8b1366d4 Add RK for second order and discrete data.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
120 end
335b8b1366d4 Add RK for second order and discrete data.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
121
335b8b1366d4 Add RK for second order and discrete data.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
122 methods (Static)
335b8b1366d4 Add RK for second order and discrete data.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
123 function k = getTimeStep(lambda, order)
335b8b1366d4 Add RK for second order and discrete data.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
124 default_arg('order', 4);
335b8b1366d4 Add RK for second order and discrete data.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
125 k = obj.tsImplementation.getTimeStep(lambda, order);
335b8b1366d4 Add RK for second order and discrete data.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
126 end
335b8b1366d4 Add RK for second order and discrete data.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
127 end
335b8b1366d4 Add RK for second order and discrete data.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
128
335b8b1366d4 Add RK for second order and discrete data.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
129 end