Mercurial > repos > public > sbplib
annotate +time/CdiffNonlin.m @ 1213:43f1cd11e8e8 feature/poroelastic
Add physical normals to AnisotropicCurvilinear
| author | Martin Almquist <malmquist@stanford.edu> |
|---|---|
| date | Mon, 14 Oct 2019 13:54:50 -0700 |
| parents | d1f9dd55a2b0 |
| children | b5e5b195da1e |
| rev | line source |
|---|---|
|
1
5ae4f23d9130
Added CdiffNonlin timestepper. Probably fixed a bug with Cdiff. Added default arguments to Rk4SecondOrderNonlin.
Jonatan Werpers <jonatan@werpers.com>
parents:
0
diff
changeset
|
1 classdef CdiffNonlin < time.Timestepper |
| 0 | 2 properties |
| 3 D | |
| 4 E | |
| 5 S | |
| 6 k | |
| 7 t | |
| 8 v | |
| 9 v_prev | |
| 10 n | |
| 11 end | |
| 12 | |
| 13 | |
| 14 methods | |
|
13
b18d3d201a71
Fixed initialization of step counter in timesteppers.
Jonatan Werpers <jonatan@werpers.com>
parents:
4
diff
changeset
|
15 function obj = CdiffNonlin(D, E, S, k, t0,n0, v, v_prev) |
|
4
8e14b5a577a6
Attempt att optimization of CdiffNonlin.
Jonatan Werpers <jonatan@werpers.com>
parents:
1
diff
changeset
|
16 m = size(D(v),1); |
|
30
cd2e28c5ecd2
Corrected parameter treatment in nonlin timesteppers.
Jonatan Werpers <jonatan@werpers.com>
parents:
13
diff
changeset
|
17 default_arg('E',0); |
|
cd2e28c5ecd2
Corrected parameter treatment in nonlin timesteppers.
Jonatan Werpers <jonatan@werpers.com>
parents:
13
diff
changeset
|
18 default_arg('S',0); |
|
cd2e28c5ecd2
Corrected parameter treatment in nonlin timesteppers.
Jonatan Werpers <jonatan@werpers.com>
parents:
13
diff
changeset
|
19 |
|
cd2e28c5ecd2
Corrected parameter treatment in nonlin timesteppers.
Jonatan Werpers <jonatan@werpers.com>
parents:
13
diff
changeset
|
20 if isnumeric(S) |
|
cd2e28c5ecd2
Corrected parameter treatment in nonlin timesteppers.
Jonatan Werpers <jonatan@werpers.com>
parents:
13
diff
changeset
|
21 S = @(v,t)S; |
|
cd2e28c5ecd2
Corrected parameter treatment in nonlin timesteppers.
Jonatan Werpers <jonatan@werpers.com>
parents:
13
diff
changeset
|
22 end |
|
cd2e28c5ecd2
Corrected parameter treatment in nonlin timesteppers.
Jonatan Werpers <jonatan@werpers.com>
parents:
13
diff
changeset
|
23 |
|
cd2e28c5ecd2
Corrected parameter treatment in nonlin timesteppers.
Jonatan Werpers <jonatan@werpers.com>
parents:
13
diff
changeset
|
24 if isnumeric(E) |
|
cd2e28c5ecd2
Corrected parameter treatment in nonlin timesteppers.
Jonatan Werpers <jonatan@werpers.com>
parents:
13
diff
changeset
|
25 E = @(v)E; |
|
cd2e28c5ecd2
Corrected parameter treatment in nonlin timesteppers.
Jonatan Werpers <jonatan@werpers.com>
parents:
13
diff
changeset
|
26 end |
| 0 | 27 |
|
1
5ae4f23d9130
Added CdiffNonlin timestepper. Probably fixed a bug with Cdiff. Added default arguments to Rk4SecondOrderNonlin.
Jonatan Werpers <jonatan@werpers.com>
parents:
0
diff
changeset
|
28 |
|
5ae4f23d9130
Added CdiffNonlin timestepper. Probably fixed a bug with Cdiff. Added default arguments to Rk4SecondOrderNonlin.
Jonatan Werpers <jonatan@werpers.com>
parents:
0
diff
changeset
|
29 % m = size(D,1); |
|
5ae4f23d9130
Added CdiffNonlin timestepper. Probably fixed a bug with Cdiff. Added default arguments to Rk4SecondOrderNonlin.
Jonatan Werpers <jonatan@werpers.com>
parents:
0
diff
changeset
|
30 % default_arg('E',sparse(m,m)); |
|
5ae4f23d9130
Added CdiffNonlin timestepper. Probably fixed a bug with Cdiff. Added default arguments to Rk4SecondOrderNonlin.
Jonatan Werpers <jonatan@werpers.com>
parents:
0
diff
changeset
|
31 % default_arg('S',sparse(m,1)); |
|
5ae4f23d9130
Added CdiffNonlin timestepper. Probably fixed a bug with Cdiff. Added default arguments to Rk4SecondOrderNonlin.
Jonatan Werpers <jonatan@werpers.com>
parents:
0
diff
changeset
|
32 |
| 0 | 33 obj.D = D; |
| 34 obj.E = E; | |
| 35 obj.S = S; | |
| 36 obj.k = k; | |
|
1
5ae4f23d9130
Added CdiffNonlin timestepper. Probably fixed a bug with Cdiff. Added default arguments to Rk4SecondOrderNonlin.
Jonatan Werpers <jonatan@werpers.com>
parents:
0
diff
changeset
|
37 obj.t = t0; |
|
13
b18d3d201a71
Fixed initialization of step counter in timesteppers.
Jonatan Werpers <jonatan@werpers.com>
parents:
4
diff
changeset
|
38 obj.n = n0; |
| 0 | 39 obj.v = v; |
| 40 obj.v_prev = v_prev; | |
| 41 end | |
| 42 | |
| 43 function [v,t] = getV(obj) | |
| 44 v = obj.v; | |
| 45 t = obj.t; | |
| 46 end | |
| 47 | |
| 48 function [vt,t] = getVt(obj) | |
| 49 vt = (obj.v-obj.v_prev)/obj.k; % Could be improved using u_tt = f(u)) | |
| 50 t = obj.t; | |
| 51 end | |
| 52 | |
| 53 function obj = step(obj) | |
|
4
8e14b5a577a6
Attempt att optimization of CdiffNonlin.
Jonatan Werpers <jonatan@werpers.com>
parents:
1
diff
changeset
|
54 D = obj.D(obj.v); |
|
8e14b5a577a6
Attempt att optimization of CdiffNonlin.
Jonatan Werpers <jonatan@werpers.com>
parents:
1
diff
changeset
|
55 E = obj.E(obj.v); |
|
30
cd2e28c5ecd2
Corrected parameter treatment in nonlin timesteppers.
Jonatan Werpers <jonatan@werpers.com>
parents:
13
diff
changeset
|
56 S = obj.S(obj.v,obj.t); |
|
4
8e14b5a577a6
Attempt att optimization of CdiffNonlin.
Jonatan Werpers <jonatan@werpers.com>
parents:
1
diff
changeset
|
57 |
|
8e14b5a577a6
Attempt att optimization of CdiffNonlin.
Jonatan Werpers <jonatan@werpers.com>
parents:
1
diff
changeset
|
58 m = size(D,1); |
|
8e14b5a577a6
Attempt att optimization of CdiffNonlin.
Jonatan Werpers <jonatan@werpers.com>
parents:
1
diff
changeset
|
59 I = speye(m); |
|
8e14b5a577a6
Attempt att optimization of CdiffNonlin.
Jonatan Werpers <jonatan@werpers.com>
parents:
1
diff
changeset
|
60 |
|
8e14b5a577a6
Attempt att optimization of CdiffNonlin.
Jonatan Werpers <jonatan@werpers.com>
parents:
1
diff
changeset
|
61 %% Calculate for which indices we need to solve system of equations |
|
8e14b5a577a6
Attempt att optimization of CdiffNonlin.
Jonatan Werpers <jonatan@werpers.com>
parents:
1
diff
changeset
|
62 [rows,cols] = find(E); |
|
8e14b5a577a6
Attempt att optimization of CdiffNonlin.
Jonatan Werpers <jonatan@werpers.com>
parents:
1
diff
changeset
|
63 j = union(rows,cols); |
|
8e14b5a577a6
Attempt att optimization of CdiffNonlin.
Jonatan Werpers <jonatan@werpers.com>
parents:
1
diff
changeset
|
64 i = setdiff(1:m,j); |
|
8e14b5a577a6
Attempt att optimization of CdiffNonlin.
Jonatan Werpers <jonatan@werpers.com>
parents:
1
diff
changeset
|
65 |
|
8e14b5a577a6
Attempt att optimization of CdiffNonlin.
Jonatan Werpers <jonatan@werpers.com>
parents:
1
diff
changeset
|
66 |
|
8e14b5a577a6
Attempt att optimization of CdiffNonlin.
Jonatan Werpers <jonatan@werpers.com>
parents:
1
diff
changeset
|
67 %% Calculate matrices need for the timestep |
|
8e14b5a577a6
Attempt att optimization of CdiffNonlin.
Jonatan Werpers <jonatan@werpers.com>
parents:
1
diff
changeset
|
68 % Before optimization: A = 1/k^2 * I - 1/(2*k)*E; |
|
8e14b5a577a6
Attempt att optimization of CdiffNonlin.
Jonatan Werpers <jonatan@werpers.com>
parents:
1
diff
changeset
|
69 k = obj.k; |
|
31
d1f9dd55a2b0
Fixed bug in step() for CdiffNonlin.
Jonatan Werpers <jonatan@werpers.com>
parents:
30
diff
changeset
|
70 |
|
4
8e14b5a577a6
Attempt att optimization of CdiffNonlin.
Jonatan Werpers <jonatan@werpers.com>
parents:
1
diff
changeset
|
71 Aj = 1/k^2 * I(j,j) - 1/(2*k)*E(j,j); |
|
8e14b5a577a6
Attempt att optimization of CdiffNonlin.
Jonatan Werpers <jonatan@werpers.com>
parents:
1
diff
changeset
|
72 B = 2/k^2 * I + D; |
|
8e14b5a577a6
Attempt att optimization of CdiffNonlin.
Jonatan Werpers <jonatan@werpers.com>
parents:
1
diff
changeset
|
73 C = -1/k^2 * I - 1/(2*k)*E; |
|
8e14b5a577a6
Attempt att optimization of CdiffNonlin.
Jonatan Werpers <jonatan@werpers.com>
parents:
1
diff
changeset
|
74 |
|
8e14b5a577a6
Attempt att optimization of CdiffNonlin.
Jonatan Werpers <jonatan@werpers.com>
parents:
1
diff
changeset
|
75 %% Take the timestep |
|
8e14b5a577a6
Attempt att optimization of CdiffNonlin.
Jonatan Werpers <jonatan@werpers.com>
parents:
1
diff
changeset
|
76 v = obj.v; |
|
31
d1f9dd55a2b0
Fixed bug in step() for CdiffNonlin.
Jonatan Werpers <jonatan@werpers.com>
parents:
30
diff
changeset
|
77 v_prev = obj.v_prev; |
|
4
8e14b5a577a6
Attempt att optimization of CdiffNonlin.
Jonatan Werpers <jonatan@werpers.com>
parents:
1
diff
changeset
|
78 |
|
31
d1f9dd55a2b0
Fixed bug in step() for CdiffNonlin.
Jonatan Werpers <jonatan@werpers.com>
parents:
30
diff
changeset
|
79 % Want to solve the seq A*v_next = b where |
|
d1f9dd55a2b0
Fixed bug in step() for CdiffNonlin.
Jonatan Werpers <jonatan@werpers.com>
parents:
30
diff
changeset
|
80 b = (B*v + C*v_prev + S); |
|
d1f9dd55a2b0
Fixed bug in step() for CdiffNonlin.
Jonatan Werpers <jonatan@werpers.com>
parents:
30
diff
changeset
|
81 |
|
d1f9dd55a2b0
Fixed bug in step() for CdiffNonlin.
Jonatan Werpers <jonatan@werpers.com>
parents:
30
diff
changeset
|
82 % Before optimization: obj.v = A\b; |
|
d1f9dd55a2b0
Fixed bug in step() for CdiffNonlin.
Jonatan Werpers <jonatan@werpers.com>
parents:
30
diff
changeset
|
83 |
|
d1f9dd55a2b0
Fixed bug in step() for CdiffNonlin.
Jonatan Werpers <jonatan@werpers.com>
parents:
30
diff
changeset
|
84 obj.v(i) = k^2*b(i); |
|
d1f9dd55a2b0
Fixed bug in step() for CdiffNonlin.
Jonatan Werpers <jonatan@werpers.com>
parents:
30
diff
changeset
|
85 obj.v(j) = Aj\b(j); |
|
d1f9dd55a2b0
Fixed bug in step() for CdiffNonlin.
Jonatan Werpers <jonatan@werpers.com>
parents:
30
diff
changeset
|
86 |
|
4
8e14b5a577a6
Attempt att optimization of CdiffNonlin.
Jonatan Werpers <jonatan@werpers.com>
parents:
1
diff
changeset
|
87 obj.v_prev = v; |
|
8e14b5a577a6
Attempt att optimization of CdiffNonlin.
Jonatan Werpers <jonatan@werpers.com>
parents:
1
diff
changeset
|
88 |
|
8e14b5a577a6
Attempt att optimization of CdiffNonlin.
Jonatan Werpers <jonatan@werpers.com>
parents:
1
diff
changeset
|
89 %% Update state of the timestepper |
| 0 | 90 obj.t = obj.t + obj.k; |
| 91 obj.n = obj.n + 1; | |
| 92 end | |
| 93 end | |
| 94 end |