Mercurial > repos > public > sbplib

annotate +time/ExplicitRungeKuttaSecondOrderDiscreteData.m @ 958:72cd29107a9a feature/poroelastic

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Temporary changes in multiblock.DiffOp. Change traction operators in Elastic2dvariable to be true boundary operators. But adjoint FD conv test fails for dirichlet BC so need to debug!
author Martin Almquist <malmquist@stanford.edu>
date Wed, 05 Dec 2018 18:58:10 -0800
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