Mercurial > repos > public > sbplib

annotate +time/SBPInTimeSecondOrderFormImplicit.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 66eb4a2bbb72
children 8894e9c49e40
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
460
e0caae9ef6ed Add SBPinTime for linear DAE formulations (BUGS!)
Jonatan Werpers <jonatan@werpers.com>
parents: 399
diff changeset
1 classdef SBPInTimeSecondOrderFormImplicit < time.Timestepper
399
c92d2f8319c2 Add wrapper for time stepping second order problems.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
2 properties
460
e0caae9ef6ed Add SBPinTime for linear DAE formulations (BUGS!)
Jonatan Werpers <jonatan@werpers.com>
parents: 399
diff changeset
3 A, B, C, f
e0caae9ef6ed Add SBPinTime for linear DAE formulations (BUGS!)
Jonatan Werpers <jonatan@werpers.com>
parents: 399
diff changeset
4 AA, BB, ff
399
c92d2f8319c2 Add wrapper for time stepping second order problems.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
5
c92d2f8319c2 Add wrapper for time stepping second order problems.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
6 n
c92d2f8319c2 Add wrapper for time stepping second order problems.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
7 t
c92d2f8319c2 Add wrapper for time stepping second order problems.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
8 k
c92d2f8319c2 Add wrapper for time stepping second order problems.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
9
c92d2f8319c2 Add wrapper for time stepping second order problems.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
10 firstOrderTimeStepper
c92d2f8319c2 Add wrapper for time stepping second order problems.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
11 end
c92d2f8319c2 Add wrapper for time stepping second order problems.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
12
c92d2f8319c2 Add wrapper for time stepping second order problems.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
13 methods
463
45a3924140f4 Move f(t) to the RHS
Jonatan Werpers <jonatan@werpers.com>
parents: 460
diff changeset
14 % Solves A*u_tt + B*u_t + C*u = f(t)
399
c92d2f8319c2 Add wrapper for time stepping second order problems.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
15 % A, B can either both be constants or both be function handles,
c92d2f8319c2 Add wrapper for time stepping second order problems.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
16 % They can also be omitted by setting them equal to the empty matrix.
744
94bd0f3293c8 SBPInTimeSecondOrderForm: Scale the system acording to C by default
Jonatan Werpers <jonatan@werpers.com>
parents: 464
diff changeset
17 function obj = SBPInTimeSecondOrderFormImplicit(A, B, C, f, k, t0, v0, v0t, do_scaling, TYPE, order, blockSize)
460
e0caae9ef6ed Add SBPinTime for linear DAE formulations (BUGS!)
Jonatan Werpers <jonatan@werpers.com>
parents: 399
diff changeset
18 default_arg('f', []);
399
c92d2f8319c2 Add wrapper for time stepping second order problems.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
19 default_arg('TYPE', []);
c92d2f8319c2 Add wrapper for time stepping second order problems.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
20 default_arg('order', []);
c92d2f8319c2 Add wrapper for time stepping second order problems.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
21 default_arg('blockSize',[]);
774
66eb4a2bbb72 Remove default scaling of the system.
Jonatan Werpers <jonatan@werpers.com>
parents: 763
diff changeset
22 default_arg('do_scaling', false);
399
c92d2f8319c2 Add wrapper for time stepping second order problems.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
23
c92d2f8319c2 Add wrapper for time stepping second order problems.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
24 m = length(v0);
c92d2f8319c2 Add wrapper for time stepping second order problems.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
25
763
f2f5efc3644f time.SBPInTimeSecondOrderFormImplicit: Fix default value for A
Jonatan Werpers <jonatan@werpers.com>
parents: 744
diff changeset
26 default_arg('A', speye(m, m));
399
c92d2f8319c2 Add wrapper for time stepping second order problems.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
27 default_arg('B', sparse(m, m));
460
e0caae9ef6ed Add SBPinTime for linear DAE formulations (BUGS!)
Jonatan Werpers <jonatan@werpers.com>
parents: 399
diff changeset
28 default_arg('C', sparse(m, m));
399
c92d2f8319c2 Add wrapper for time stepping second order problems.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
29
c92d2f8319c2 Add wrapper for time stepping second order problems.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
30 I = speye(m);
c92d2f8319c2 Add wrapper for time stepping second order problems.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
31 O = sparse(m,m);
c92d2f8319c2 Add wrapper for time stepping second order problems.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
32
460
e0caae9ef6ed Add SBPinTime for linear DAE formulations (BUGS!)
Jonatan Werpers <jonatan@werpers.com>
parents: 399
diff changeset
33 % Rewrite to
e0caae9ef6ed Add SBPinTime for linear DAE formulations (BUGS!)
Jonatan Werpers <jonatan@werpers.com>
parents: 399
diff changeset
34 % AA*w_t = BB*w + ff(t);
e0caae9ef6ed Add SBPinTime for linear DAE formulations (BUGS!)
Jonatan Werpers <jonatan@werpers.com>
parents: 399
diff changeset
35
e0caae9ef6ed Add SBPinTime for linear DAE formulations (BUGS!)
Jonatan Werpers <jonatan@werpers.com>
parents: 399
diff changeset
36 obj.AA = [
e0caae9ef6ed Add SBPinTime for linear DAE formulations (BUGS!)
Jonatan Werpers <jonatan@werpers.com>
parents: 399
diff changeset
37 I, O;
e0caae9ef6ed Add SBPinTime for linear DAE formulations (BUGS!)
Jonatan Werpers <jonatan@werpers.com>
parents: 399
diff changeset
38 O, A;
399
c92d2f8319c2 Add wrapper for time stepping second order problems.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
39 ];
460
e0caae9ef6ed Add SBPinTime for linear DAE formulations (BUGS!)
Jonatan Werpers <jonatan@werpers.com>
parents: 399
diff changeset
40 obj.BB = [
e0caae9ef6ed Add SBPinTime for linear DAE formulations (BUGS!)
Jonatan Werpers <jonatan@werpers.com>
parents: 399
diff changeset
41 O, I;
464
d4b999585af1 Fix math in conversion to first order form system
Jonatan Werpers <jonatan@werpers.com>
parents: 463
diff changeset
42 -C, -B;
399
c92d2f8319c2 Add wrapper for time stepping second order problems.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
43 ];
c92d2f8319c2 Add wrapper for time stepping second order problems.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
44
460
e0caae9ef6ed Add SBPinTime for linear DAE formulations (BUGS!)
Jonatan Werpers <jonatan@werpers.com>
parents: 399
diff changeset
45 if ~isempty(f)
e0caae9ef6ed Add SBPinTime for linear DAE formulations (BUGS!)
Jonatan Werpers <jonatan@werpers.com>
parents: 399
diff changeset
46 obj.ff = @(t)[
e0caae9ef6ed Add SBPinTime for linear DAE formulations (BUGS!)
Jonatan Werpers <jonatan@werpers.com>
parents: 399
diff changeset
47 sparse(m,1);
463
45a3924140f4 Move f(t) to the RHS
Jonatan Werpers <jonatan@werpers.com>
parents: 460
diff changeset
48 f(t);
460
e0caae9ef6ed Add SBPinTime for linear DAE formulations (BUGS!)
Jonatan Werpers <jonatan@werpers.com>
parents: 399
diff changeset
49 ];
e0caae9ef6ed Add SBPinTime for linear DAE formulations (BUGS!)
Jonatan Werpers <jonatan@werpers.com>
parents: 399
diff changeset
50 else
e0caae9ef6ed Add SBPinTime for linear DAE formulations (BUGS!)
Jonatan Werpers <jonatan@werpers.com>
parents: 399
diff changeset
51 obj.ff = @(t) sparse(2*m,1);
e0caae9ef6ed Add SBPinTime for linear DAE formulations (BUGS!)
Jonatan Werpers <jonatan@werpers.com>
parents: 399
diff changeset
52 end
e0caae9ef6ed Add SBPinTime for linear DAE formulations (BUGS!)
Jonatan Werpers <jonatan@werpers.com>
parents: 399
diff changeset
53
399
c92d2f8319c2 Add wrapper for time stepping second order problems.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
54 w0 = [v0; v0t];
c92d2f8319c2 Add wrapper for time stepping second order problems.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
55
c92d2f8319c2 Add wrapper for time stepping second order problems.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
56 obj.k = k;
c92d2f8319c2 Add wrapper for time stepping second order problems.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
57 obj.t = t0;
c92d2f8319c2 Add wrapper for time stepping second order problems.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
58 obj.n = 0;
c92d2f8319c2 Add wrapper for time stepping second order problems.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
59
744
94bd0f3293c8 SBPInTimeSecondOrderForm: Scale the system acording to C by default
Jonatan Werpers <jonatan@werpers.com>
parents: 464
diff changeset
60 if do_scaling
94bd0f3293c8 SBPInTimeSecondOrderForm: Scale the system acording to C by default
Jonatan Werpers <jonatan@werpers.com>
parents: 464
diff changeset
61 scaling = [ones(m,1); sqrt(diag(C))];
94bd0f3293c8 SBPInTimeSecondOrderForm: Scale the system acording to C by default
Jonatan Werpers <jonatan@werpers.com>
parents: 464
diff changeset
62 obj.firstOrderTimeStepper = time.SBPInTimeScaled(obj.AA, obj.BB, obj.ff, obj.k, obj.t, w0, scaling, TYPE, order, blockSize);
94bd0f3293c8 SBPInTimeSecondOrderForm: Scale the system acording to C by default
Jonatan Werpers <jonatan@werpers.com>
parents: 464
diff changeset
63 else
94bd0f3293c8 SBPInTimeSecondOrderForm: Scale the system acording to C by default
Jonatan Werpers <jonatan@werpers.com>
parents: 464
diff changeset
64 obj.firstOrderTimeStepper = time.SBPInTimeImplicitFormulation(obj.AA, obj.BB, obj.ff, obj.k, obj.t, w0, TYPE, order, blockSize);
94bd0f3293c8 SBPInTimeSecondOrderForm: Scale the system acording to C by default
Jonatan Werpers <jonatan@werpers.com>
parents: 464
diff changeset
65 end
399
c92d2f8319c2 Add wrapper for time stepping second order problems.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
66 end
c92d2f8319c2 Add wrapper for time stepping second order problems.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
67
c92d2f8319c2 Add wrapper for time stepping second order problems.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
68 function [v,t] = getV(obj)
c92d2f8319c2 Add wrapper for time stepping second order problems.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
69 w = obj.firstOrderTimeStepper.getV();
c92d2f8319c2 Add wrapper for time stepping second order problems.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
70 v = w(1:end/2);
c92d2f8319c2 Add wrapper for time stepping second order problems.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
71 t = obj.t;
c92d2f8319c2 Add wrapper for time stepping second order problems.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
72 end
c92d2f8319c2 Add wrapper for time stepping second order problems.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
73
c92d2f8319c2 Add wrapper for time stepping second order problems.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
74 function [vt,t] = getVt(obj)
c92d2f8319c2 Add wrapper for time stepping second order problems.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
75 w = obj.firstOrderTimeStepper.getV();
c92d2f8319c2 Add wrapper for time stepping second order problems.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
76 vt = w(end/2+1:end);
c92d2f8319c2 Add wrapper for time stepping second order problems.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
77 t = obj.t;
c92d2f8319c2 Add wrapper for time stepping second order problems.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
78 end
c92d2f8319c2 Add wrapper for time stepping second order problems.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
79
c92d2f8319c2 Add wrapper for time stepping second order problems.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
80 function obj = step(obj)
c92d2f8319c2 Add wrapper for time stepping second order problems.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
81 obj.firstOrderTimeStepper.step();
c92d2f8319c2 Add wrapper for time stepping second order problems.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
82 obj.t = obj.firstOrderTimeStepper.t;
c92d2f8319c2 Add wrapper for time stepping second order problems.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
83 obj.n = obj.firstOrderTimeStepper.n;
c92d2f8319c2 Add wrapper for time stepping second order problems.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
84 end
c92d2f8319c2 Add wrapper for time stepping second order problems.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
85 end
c92d2f8319c2 Add wrapper for time stepping second order problems.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
86 end