Skip to content

Commit 6725a88

Browse files
halterwhalterw
committed
speedup
1 parent b8b94ba commit 6725a88

17 files changed

+42
-25
lines changed

‎A_simple_computer_program_Matlab_v2_corrected.m

+1-1
Original file line numberDiff line numberDiff line change
@@ -35,7 +35,7 @@
3535
% Parameters for pseudo-transient iterations
3636
tol = 5e-6; err_absolute = 1; err_relative = 1;
3737
%%% Numerical parameters for iterative solver %%%
38-
CFLV = 1/5; CFLP = CFLV/nx/20; % Numerical parameter for iterative solver
38+
CFLV = 1/2.05; CFLP = CFLV/nx*2.5; % Numerical parameter for iterative solver
3939
dpt_P = zeros(nx,ny);
4040
dpt_Vx = CFLV./(max(ETA(1:end-1,2:end-1),ETA(2:end,2:end-1))/dx^2 + max(ETA_XY(:,1:end-1),ETA_XY(:,2:end))/dx/dy);
4141
dpt_Vy = CFLV./(max(ETA(2:end-1,1:end-1),ETA(2:end-1,2:end))/dy^2 + max(ETA_XY(1:end-1,:),ETA_XY(2:end,:))/dx/dy);

‎A_simple_computer_program_Octave_v2_corrected.m

+1-1
Original file line numberDiff line numberDiff line change
@@ -47,7 +47,7 @@
4747
% Parameters for pseudo-transient iterations
4848
tol = 5e-6; err_absolute = 1; err_relative = 1;
4949
%%% Numerical parameters for iterative solver %%%
50-
CFLV = 1/5; CFLP = CFLV/nx/20; % Numerical parameter for iterative solver
50+
CFLV = 1/2.05; CFLP = CFLV/nx*2.5; % Numerical parameter for iterative solver
5151
dpt_P = zeros(nx,ny);
5252
dpt_Vx = CFLV./(max(ETA(1:end-1,2:end-1),ETA(2:end,2:end-1))/dx^2 + max(ETA_XY(:,1:end-1),ETA_XY(:,2:end))/dx/dy);
5353
dpt_Vy = CFLV./(max(ETA(2:end-1,1:end-1),ETA(2:end-1,2:end))/dy^2 + max(ETA_XY(1:end-1,:),ETA_XY(2:end,:))/dx/dy);

‎Examples_Matlab/FIG06a_Rectangle_PS.m

+1-1
Original file line numberDiff line numberDiff line change
@@ -38,7 +38,7 @@
3838
% Parameters for pseudo-transient iterations
3939
tol = 1e-6; err_absolute = 1; err_relative = 1;
4040
%%% Numerical parameters for iterative solver %%%
41-
CFLV = 1/5; CFLP = CFLV/nx/20; % Numerical parameter for iterative solver
41+
CFLV = 1/2.05; CFLP = CFLV/nx/3; % Numerical parameter for iterative solver
4242
dpt_P = zeros(nx,ny);
4343
dpt_Vx = CFLV./(max(ETA(1:end-1,2:end-1),ETA(2:end,2:end-1))/dx^2 + max(ETA_XY(:,1:end-1),ETA_XY(:,2:end))/dx/dy);
4444
dpt_Vy = CFLV./(max(ETA(2:end-1,1:end-1),ETA(2:end-1,2:end))/dy^2 + max(ETA_XY(1:end-1,:),ETA_XY(2:end,:))/dx/dy);

‎Examples_Matlab/FIG06b_Rectangle_SS.m

+1-1
Original file line numberDiff line numberDiff line change
@@ -38,7 +38,7 @@
3838
% Parameters for pseudo-transient iterations
3939
tol = 1e-6; err_absolute = 1; err_relative = 1;
4040
%%% Numerical parameters for iterative solver %%%
41-
CFLV = 1/5; CFLP = CFLV/nx/20; % Numerical parameter for iterative solver
41+
CFLV = 1/2.05; CFLP = CFLV/nx/3; % Numerical parameter for iterative solver
4242
dpt_P = zeros(nx,ny);
4343
dpt_Vx = CFLV./(max(ETA(1:end-1,2:end-1),ETA(2:end,2:end-1))/dx^2 + max(ETA_XY(:,1:end-1),ETA_XY(:,2:end))/dx/dy);
4444
dpt_Vy = CFLV./(max(ETA(2:end-1,1:end-1),ETA(2:end-1,2:end))/dy^2 + max(ETA_XY(1:end-1,:),ETA_XY(2:end,:))/dx/dy);

‎Examples_Matlab/FIG08a_Ellipse_linVisc.m

+1-1
Original file line numberDiff line numberDiff line change
@@ -35,7 +35,7 @@
3535
% Parameters for pseudo-transient iterations
3636
tol = 1e-6; err_absolute = 1; err_relative = 1;
3737
%%% Numerical parameters for iterative solver %%%
38-
CFLV = 1/5; CFLP = CFLV/nx/20; % Numerical parameter for iterative solver
38+
CFLV = 1/2.05; CFLP = CFLV/nx/3; % Numerical parameter for iterative solver
3939
dpt_P = zeros(nx,ny);
4040
dpt_Vx = CFLV./(max(ETA(1:end-1,2:end-1),ETA(2:end,2:end-1))/dx^2 + max(ETA_XY(:,1:end-1),ETA_XY(:,2:end))/dx/dy);
4141
dpt_Vy = CFLV./(max(ETA(2:end-1,1:end-1),ETA(2:end-1,2:end))/dy^2 + max(ETA_XY(1:end-1,:),ETA_XY(2:end,:))/dx/dy);

‎Examples_Matlab/FIG08b_Ellipse_pwlaw5.m

+1-1
Original file line numberDiff line numberDiff line change
@@ -35,7 +35,7 @@
3535
% Parameters for pseudo-transient iterations
3636
tol = 1e-6; err_absolute = 1; err_relative = 1;
3737
%%% Numerical parameters for iterative solver %%%
38-
CFLV = 1/5; CFLP = CFLV/nx/20; % Numerical parameter for iterative solver
38+
CFLV = 1/2.05; CFLP = CFLV/nx/3; % Numerical parameter for iterative solver
3939
dpt_P = zeros(nx,ny);
4040
dpt_Vx = CFLV./(max(ETA(1:end-1,2:end-1),ETA(2:end,2:end-1))/dx^2 + max(ETA_XY(:,1:end-1),ETA_XY(:,2:end))/dx/dy);
4141
dpt_Vy = CFLV./(max(ETA(2:end-1,1:end-1),ETA(2:end-1,2:end))/dy^2 + max(ETA_XY(1:end-1,:),ETA_XY(2:end,:))/dx/dy);

‎Examples_Matlab/FIG09_MulitInclusions.m

+1-2
Original file line numberDiff line numberDiff line change
@@ -80,8 +80,7 @@
8080
% Parameters for pseudo-transient iterations
8181
tol = 1e-6; err_absolute = 1; err_relative = 1;
8282
%%% Numerical parameters for iterative solver %%%
83-
%CFLV = 1/5; CFLP = CFLV/nx/20; % Numerical parameter for iterative solver
84-
CFLV = 1/5; CFLP = CFLV/nx/50; % Numerical parameter for iterative solver
83+
CFLV = 1/2.10; CFLP = CFLV/nx/25; % Numerical parameter for iterative solver
8584
dpt_P = zeros(nx,ny);
8685
dpt_Vx = CFLV./(max(ETA(1:end-1,2:end-1),ETA(2:end,2:end-1))/dx^2 + max(ETA_XY(:,1:end-1),ETA_XY(:,2:end))/dx/dy);
8786
dpt_Vy = CFLV./(max(ETA(2:end-1,1:end-1),ETA(2:end-1,2:end))/dy^2 + max(ETA_XY(1:end-1,:),ETA_XY(2:end,:))/dx/dy);

‎Examples_Matlab/FIG10_Garnet.m

+1-1
Original file line numberDiff line numberDiff line change
@@ -41,7 +41,7 @@
4141
% Parameters for pseudo-transient iterations
4242
tol = 1e-6; err_absolute = 1; err_relative = 1;
4343
%%% Numerical parameters for iterative solver %%%
44-
CFLV = 1/5; CFLP = CFLV/nx/20; % Numerical parameter for iterative solver
44+
CFLV = 1/2.10; CFLP = CFLV/nx*2; % Numerical parameter for iterative solver
4545
dpt_P = zeros(nx,ny);
4646
dpt_Vx = CFLV./(max(ETA(1:end-1,2:end-1),ETA(2:end,2:end-1))/dx^2 + max(ETA_XY(:,1:end-1),ETA_XY(:,2:end))/dx/dy);
4747
dpt_Vy = CFLV./(max(ETA(2:end-1,1:end-1),ETA(2:end-1,2:end))/dy^2 + max(ETA_XY(1:end-1,:),ETA_XY(2:end,:))/dx/dy);

‎Examples_Matlab/FIG11_Boudinage.m

+1-1
Original file line numberDiff line numberDiff line change
@@ -34,7 +34,7 @@
3434
% Parameters for pseudo-transient iterations
3535
tol = 1e-6; err_absolute = 1; err_relative = 1;
3636
%%% Numerical parameters for iterative solver %%%
37-
CFLV = 1/5; CFLP = CFLV/nx/20; % Numerical parameter for iterative solver
37+
CFLV = 1/2.05; CFLP = CFLV/nx*3; % Numerical parameter for iterative solver
3838
dpt_P = zeros(nx,ny);
3939
dpt_Vx = CFLV./(max(ETA(1:end-1,2:end-1),ETA(2:end,2:end-1))/dx^2 + max(ETA_XY(:,1:end-1),ETA_XY(:,2:end))/dx/dy);
4040
dpt_Vy = CFLV./(max(ETA(2:end-1,1:end-1),ETA(2:end-1,2:end))/dy^2 + max(ETA_XY(1:end-1,:),ETA_XY(2:end,:))/dx/dy);

‎README.md

+33-15
Original file line numberDiff line numberDiff line change
@@ -23,9 +23,15 @@ ps = 1; % ps = 1 models pure shear; ps = 0 models simple shear
2323
```
2424

2525
where in Fig. 6a we activate pure shear, whereas in Fig. 6b we activate simple shear.
26+
<figure>
27+
<img src="figures/Fig_06a_output.png" width="800">
28+
<figcaption>Fig. 6a setup: hard rectangular inclusion, pure shear</figcaption>
29+
</figure>
30+
<figure>
31+
<img src="figures/Fig_06b_output.png" width="800">
32+
<figcaption>Fig. 6b: hard rectangular inclusion, simple shear</figcaption>
33+
</figure>
2634

27-
<img src="figures/Fig_06a_output.png" width="800">
28-
<img src="figures/Fig_06b_output.png" width="800">
2935

3036

3137

@@ -43,43 +49,55 @@ n_exp = 5;
4349

4450
where in Fig. 8b we chose a power-law exponent > 1.
4551

46-
<img src="figures/Fig_08a_output.png" width="800">
47-
<img src="figures/Fig_08b_output.png" width="800">
52+
<figure>
53+
<img src="figures/Fig_08a_output.png" width="800">
54+
<figcaption>Fig. 8a: soft elliptical inclusion, linear viscous</figcaption>
55+
</figure>
56+
<figure>
57+
<img src="figures/Fig_08b_output.png" width="800">
58+
<figcaption>Fig. 8b: soft elliptical inclusion, power-law viscous</figcaption>
59+
</figure>
60+
4861

4962
### Multiple inclusions
5063

51-
Fig. 9 illustrates the interaction between multiple inclusions. This model configuration is calculated on a higher resolution (901 x 601) and can take multiple hours to fully converge.
64+
Fig. 9 illustrates the interaction between multiple inclusions. This model configuration is calculated on a higher resolution (901 x 601) and takes multiple hours to fully converge. Results will be added later.
5265

53-
<img src="figures/Fig_9_output.png" width="800">
5466

5567
### Garnet
5668

5769
In Fig. 10 we want to investigate whether we can infer the pressure field inside and around a rigid garnet porphyroblast. To define the shape of the garnet we created a polygon using Matlab's [ginput](https://ch.mathworks.com/help/matlab/ref/ginput.html) function on a photo of a real rock.
5870

59-
<img src="figures/Fig_10_output.png" width="800">
71+
<figure>
72+
<img src="figures/Fig_10_output.png" width="800">
73+
<figcaption>Fig. 10: garnet inclusion</figcaption>
74+
</figure>
75+
6076

6177
Note the error convergence behaviour. In this example, the (absolute) error reaches a plateau before fully converging to the bottom. This behaviour is typical for the pseudo-transient method.
6278

6379
### Boudinage
6480

6581
In Fig. 11 we show the pressure and stress distribution inside and between 2 separating boudin blocks. Feel free to modify the code visualization to confirm that the horizontal total stress (Sxx) is indeed continuous across the weak gap, as discussed in Halter et al. 2022.
6682

67-
<img src="figures/Fig_11_output.png" width="800">
6883

84+
<figure>
85+
<img src="figures/Fig_11_output.png" width="800">
86+
<figcaption>Fig. 11: boudinage</figcaption>
87+
</figure>
6988

7089
### Computation time
7190

7291
An overview of calculation times obtained on a commercial laptop (Lenovo ThinkPad P1 gen 3).
7392

7493
| Model configuration | Resolution (nx,ny) | Tolerance | Computation time (seconds) |
7594
|-------------------------------------|--------------------|-----------|----------------------------|
76-
| Fig. 6a - Rectangle pure shear | (201,201) | 1e-6 | 3397 |
77-
| Fig. 6b - Rectangle simple shear | (201,201) | 1e-6 | 3714 |
78-
| Fig. 8a - Ellipse linear viscous | (201,201) | 1e-6 | 4220 |
79-
| Fig. 8b - Ellipes power-law viscous | (201,201) | 1e-6 | 4184 |
80-
| Fig. 9 - Multiple inclusions | (901,601) | 1e-6 | TBD |
81-
| Fig. 10 - Garnet | (301,186) | 1e-6 | 7373 |
82-
| Fig. 11 - Boudinage | (201,201) | 1e-6 | 4259 |
95+
| Fig. 6a - Rectangle pure shear | (201,201) | 1e-6 | 450 |
96+
| Fig. 6b - Rectangle simple shear | (201,201) | 1e-6 | 450 |
97+
| Fig. 8a - Ellipse linear viscous | (201,201) | 1e-6 | 207 |
98+
| Fig. 8b - Ellipes power-law viscous | (201,201) | 1e-6 | 162 |
99+
| Fig. 10 - Garnet | (301,186) | 1e-6 | 803 |
100+
| Fig. 11 - Boudinage | (201,201) | 1e-6 | 28 |
83101

84102
Note that there exist other codes using different numerical approaches (e.g., [MDOODZ7.0](https://github.com/tduretz/MDOODZ7.0)) which are capable of solving the same problem much faster. However, as emphasized in Halter et al. 2022, the beauty in our code lies in it's simplicity and readability, making coding accessible and transparent.
85103

‎figures/Fig_06a_output.png

545 Bytes
Loading

‎figures/Fig_06b_output.png

790 Bytes
Loading

‎figures/Fig_08a_output.png

203 Bytes
Loading

‎figures/Fig_08b_output.png

182 Bytes
Loading

‎figures/Fig_10_output.png

1.98 KB
Loading

‎figures/Fig_11_output.png

3.13 KB
Loading

‎figures/Fig_9_output.png

-222 KB
Binary file not shown.

0 commit comments

Comments
 (0)