Mercurial > repos > public > sbplib

annotate +time/SBPInTimeSecondOrderFormImplicit.m @ 774:66eb4a2bbb72 feature/grids

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Remove default scaling of the system. The scaling doens't seem to help actual solutions. One example that fails in the flexural code. With large timesteps the solutions seems to blow up. One particular example is profilePresentation on the tdb_presentation_figures branch with k = 0.0005
author Jonatan Werpers <jonatan@werpers.com>
date Wed, 18 Jul 2018 15:42:52 -0700
parents f2f5efc3644f
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