Merge pull request #5502 from gassmoeller/update_particle_benchmark

Update particle benchmark documentation

benchmarks/annulus/doc/annulus.md

@@ -62,6 +62,6 @@ Additionally, the subdirectory
 [benchmarks/annulus/transient](https://github.com/geodynamics/aspect/tree/main/benchmarks/annulus/transient)
 contains an extension of the benchmark
 to time-dependent flow. The benchmark and its results are described in
-Gassmoeller et al. (2023), "Benchmarking the accuracy of higher order
-particle methods in geodynamic models of transient flow", see there
+{cite:t}`gassmoeller:etal:2023`, "Benchmarking the accuracy of higher
+order particle methods in geodynamic models of transient flow" see there
 for a detailed description.

benchmarks/annulus/transient/plot_transient_convergence.py

@@ -14,11 +14,12 @@
 
 refinements = ["2","3","4","5","6","7"]
 models = ["analytical_density", "compositional_field","continuous_compositional_field", "higher_order_true","higher_order_false"]
-labels = ["Analytical density", "DGQ2 field","Q2 field", "Particles RK2","Particles RK2 FOT"]
+labels = ["Density: Benchmark", "Density: FE field ($DGQ_2$)","Density: FE field ($Q_2$)", "Density: Particles (RK2)","Density: Particles (RK2 FOT)"]
 errors = ["u_L2","p_L2","rho_L2"]
 ylabels = [r"$\|\boldsymbol{u} - \boldsymbol{u}_h\|_{L_2}$",r"$\|p - p_h\|_{L_2}$",r"$\|\rho - \rho_h\|_{L_2}$"]
-markers=['o','X','P','v','s','D','<','>','^','+','x']
+markers=['v','<','>','^','o','D','<','>','^','+','x']
 colors = plt.rcParams['axes.prop_cycle'].by_key()['color']
+plt.rcParams['lines.markersize'] = 8.5
 
 h = []
 for refinement in refinements:
@@ -65,10 +66,9 @@ def plot_error_over_time(statistics, output_file):
         for model,label,marker,color in zip(models,labels,markers,colors):
             ax.loglog(statistics[model]["Time (seconds)"],statistics[model][errors[i_error]], label=label, color=color)
 
-        if (i_error == 2):
-            ax.set_ybound(lower=1e-6,upper=None)
+        ax.set_ybound(5e-8,0.9)
 
-    plt.xlabel("Time")
+    plt.xlabel("Time $t$")
     ax.legend(bbox_to_anchor=(1.05, 1), loc='upper left', borderaxespad=0.)
     plt.savefig(output_file, bbox_inches='tight',dpi=200)
     return None
@@ -80,6 +80,7 @@ def plot_error_over_time_steps(statistics, output_file):
         ax = plt.subplot(3,1,i_error+1)
         ax.set_ylabel(ylabels[i_error])
         plt.xlim(1,1e4)
+        ax.set_ybound(5e-8,0.9)
 
         for model,label,marker,color in zip(models,labels,markers,colors):
             ax.loglog(statistics[model]["Time step number"],statistics[model][errors[i_error]], label=label, color=color)
@@ -123,12 +124,12 @@ def plot_rk2_error_over_resolution(statistics, timestep, output_file):
         # The density error is not reliable, because the analytical density error
         # is 0. Therefore we cannot compute the relative error.
         if i_error != 2:
-            line3, = ax.loglog(h,error_values[models[i]][errors[i_error]], marker=markers[6], color=colors[6], label="relative RK2")
+            line3, = ax.loglog(h,error_values[models[i]][errors[i_error]], marker=markers[5], color=colors[6], label="(RK2 - Benchmark) / Benchmark")
 
-    plt.xlabel("h")
+    plt.xlabel("Cell size $h$")
 
-    additional_legend_elements = [Line2D([0], [0], color=colors[0], label='Analytical density'),
-    Line2D([0], [0], color=colors[3], label='RK2')]
+    additional_legend_elements = [Line2D([0], [0], color=colors[0], label='Density: Benchmark'),
+    Line2D([0], [0], color=colors[3], label='Density: Particles (RK2)')]
 
     # Create the figure
     ax.legend(bbox_to_anchor=(1.05, 1), loc='upper left', borderaxespad=0., handles=[line1,line2,additional_legend_elements[0],additional_legend_elements[1],line3])
@@ -165,7 +166,7 @@ def plot_rk2_error_over_time(statistics, output_file):
         if (i_error == 2):
             ax.set_ybound(lower=3e-5,upper=5e-5)
 
-    plt.xlabel("Time (seconds)")
+    plt.xlabel("Time $t$")
     plt.savefig(output_file, bbox_inches='tight',dpi=200)
     return None
 
@@ -181,26 +182,32 @@ def plot_error_over_resolution(statistics, timestep, output_file):
             for refinement in refinements:
                 error_values[model][error].append(statistics[refinement][model].iloc[timestep][error])
 
-    scale_factors = [1.0,15.0,3.0]
+    scale_factors = [1.0,15.0,7.0]
     x = np.linspace(7e-3,0.25,100)
     y = 0.1 * x
     y2 = 0.1 * x*x
     y3 = 0.1 * x*x*x
 
     for i_error in range(3):
         ax = plt.subplot(3,1,i_error+1)
-        ax.set_ybound(1e-7,0.5)
         ax.plot(x,scale_factors[i_error] * y, label='$h$', linestyle='--', color='grey')
         ax.plot(x,scale_factors[i_error] * y2, label='$h^2$', linestyle='-.', color='grey')
         ax.plot(x,scale_factors[i_error] * y3, label='$h^3$', linestyle=':', color='grey')
 
         ax.set_ylabel(ylabels[i_error])
 
         for i in range(len(models)):
+            # dont plot density error for the analytical model (it is 0)
+            if i_error == 2 and i == 0:
+                continue
             ax.loglog(h,error_values[models[i]][errors[i_error]], marker=markers[i], label=labels[i], color=colors[i])
 
-    plt.xlabel("h")
-    ax.legend(ncol=2, bbox_to_anchor=(1.05, 1), loc='upper left', borderaxespad=0.)
+        ax.set_ybound(5e-8,0.9)
+
+        if i_error == 0:
+            ax.legend(ncol=2, bbox_to_anchor=(1.05, 1), loc='upper left', borderaxespad=0.)
+
+    plt.xlabel("Cell size $h$")
     plt.savefig(output_file, bbox_inches='tight',dpi=200)
     return None
 

benchmarks/annulus/transient/transient_annulus.prm

@@ -76,5 +76,6 @@ subsection Solver parameters
     set Number of cheap Stokes solver steps = 200
     set Stokes solver type = block GMG
     set Krylov method for cheap solver steps = IDR(s)
+    set IDR(s) parameter = 4
   end
 end

benchmarks/rigid_shear/doc/rigid_shear.md

@@ -8,11 +8,12 @@ instantaneous, steady-state and transient versions of the benchmark are in the
 corresponding subdirectories. Each directory has a "run_all_models" shell
 script that will run the benchmark setups and a "create_formatted_tables.sh"
 script that formats the output into readable text files. For a more
-detailed description of the benchmark see Gassmoeller et al. (2019), "Evaluating
+detailed description of the benchmark see {cite:t}`gassmoller:etal:2019`, "Evaluating
 the Accuracy of Hybrid Finite Element/Particle-In-Cell Methods for
-Modeling Incompressible Stokes Flow", Geophysical Journal International, 219, 3.
+Modeling Incompressible Stokes Flow".
 
 Additionally, the "transient" directory contains an extension of the benchmark
 to time-dependent flow. The benchmark and its results are described in
-Gassmoeller et al. (2023), "Benchmarking the accuracy of higher order particle
-methods in geodynamic models of transient flow".
+{cite:t}`gassmoeller:etal:2023`, "Benchmarking the accuracy of higher order
+particle methods in geodynamic models of transient flow",
+see there for a more detailed description.

benchmarks/rigid_shear/transient/plot_transient_convergence.py

@@ -1,7 +1,9 @@
 #!/usr/bin/python
 
-# This script can be used to plot errors for the operator
-# splitting 'advection reaction' benchmark. 
+# This script can be used to plot errors for the
+# 'transient rigid shear' benchmark described in Gassmoeller
+# et al. (2023), "Benchmarking the accuracy of higher order
+# particle methods in geodynamic models of transient flow"
 
 
 import numpy as np
@@ -14,11 +16,13 @@
           "higher_order_true_interpolation_bilinear_least_squares", \
           "higher_order_false_interpolation_bilinear_least_squares"]
 
-labels = ["Analytical density", "DGQ2 field","Q2 field", "Particles RK2","Particles RK2 FOT"]
+labels = ["Density: Benchmark", "Density: FE field ($DGQ_2$)","Density: FE field ($Q_2$)", "Density: Particles (RK2)","Density: Particles (RK2 FOT)"]
 errors = ["u_L2","p_L2","rho_L2"]
 ylabels = [r"$\|\boldsymbol{u} - \boldsymbol{u}_h\|_{L_2}$",r"$\|p - p_h\|_{L_2}$",r"$\|\rho - \rho_h\|_{L_2}$"]
 
-markers=['o','X','P','v','s','D','<','>','^','+','x']
+markers=['v','<','>','^','o','D','<','>','^','+','x']
+colors = plt.rcParams['axes.prop_cycle'].by_key()['color']
+plt.rcParams['lines.markersize'] = 8.5
 
 h = []
 for refinement in refinements:
@@ -61,10 +65,9 @@ def plot_error_over_time(statistics, output_file):
         for model,label,marker in zip(models,labels,markers):
             ax.loglog(statistics[model]["Time (seconds)"],statistics[model][errors[i_error]], label=label)
 
-        if (i_error == 2):
-            ax.set_ybound(lower=1e-5,upper=None)
+        ax.set_ybound(5e-8,1.0)
 
-    plt.xlabel("Time")
+    plt.xlabel("Time $t$")
     ax.legend(bbox_to_anchor=(1.05, 1), loc='upper left', borderaxespad=0.)
     plt.savefig(output_file, bbox_inches='tight',dpi=200)
     return None
@@ -85,26 +88,31 @@ def plot_error_over_resolution(statistics, timestep, output_file):
             for refinement in refinements:
                 error_values[model][error].append(statistics[refinement][model].iloc[timestep][error])
 
-    scale_factors = [1.0,15.0,3.0]
+    scale_factors = [1.5,12.0,11.0]
     x = np.linspace(7e-3,0.25,100)
     y = 0.1 * x
     y2 = 0.1 * x*x
     y3 = 0.1 * x*x*x
 
     for i_error in range(3):
         ax = plt.subplot(3,1,i_error+1)
-        ax.set_ybound(1e-7,0.5)
         ax.plot(x,scale_factors[i_error] * y, label='$h$', linestyle='--', color='grey')
         ax.plot(x,scale_factors[i_error] * y2, label='$h^2$', linestyle='-.', color='grey')
         ax.plot(x,scale_factors[i_error] * y3, label='$h^3$', linestyle=':', color='grey')
-
         ax.set_ylabel(ylabels[i_error])
 
         for i in range(len(models)):
-            ax.loglog(h,error_values[models[i]][errors[i_error]], marker=markers[i], label=labels[i])
+            # dont plot density error for the analytical model (it is 0)
+            if i_error == 2 and i == 0:
+                continue
+            ax.loglog(h,error_values[models[i]][errors[i_error]], marker=markers[i], label=labels[i], color=colors[i])
+
+        ax.set_ybound(5e-8,1.0)
+
+        if i_error == 0:
+            ax.legend(ncol=2, bbox_to_anchor=(1.05, 1), loc='upper left', borderaxespad=0.)
 
-    plt.xlabel("h")
-    ax.legend(ncol=2, bbox_to_anchor=(1.05, 1), loc='upper left', borderaxespad=0.)
+    plt.xlabel("Cell size $h$")
     plt.savefig(output_file, bbox_inches='tight',dpi=200)
     return None
 

doc/sphinx/references.bib

@@ -11930,7 +11930,8 @@ @article{gassmoller:etal:2019
   number={3},
   pages={1915--1938},
   year={2019},
-  publisher={Oxford University Press}
+  publisher={Oxford University Press},
+  doi = {10.1093/gji/ggz405}
 }
 
 @article{Liu2019,
@@ -12179,3 +12180,14 @@ @article{priestley:etal:2018
   year={2018},
   publisher={Wiley Online Library}
 }
+
+@article{gassmoeller:etal:2023,
+  author = {Gassm\"{o}ller, R. and Dannberg, J. and Bangerth, W. and Puckett, E. G. and Thieulot, C.},
+  title = {Benchmarking the accuracy of higher order particle methods in geodynamic models of transient flow},
+  journal = {EGUsphere},
+  volume = {2023},
+  year = {2023},
+  pages = {1--29},
+  url = {https://egusphere.copernicus.org/preprints/2023/egusphere-2023-2765/},
+  doi = {10.5194/egusphere-2023-2765}
+}

tests/annulus_transient.prm

@@ -83,5 +83,6 @@ subsection Solver parameters
     set Number of cheap Stokes solver steps = 200
     set Stokes solver type = block GMG
     set Krylov method for cheap solver steps = IDR(s)
+    set IDR(s) parameter = 4
   end
 end

tests/annulus_transient/screen-output

@@ -8,49 +8,49 @@ Number of degrees of freedom: 4,560 (1,728+240+864+1,728)
 *** Timestep 0:  t=0 seconds, dt=0 seconds
    Skipping temperature solve because RHS is zero.
    Solving density_field system ... 0 iterations.
-   Solving Stokes system... 8+0 iterations.
+   Solving Stokes system... 5+0 iterations.
 
    Postprocessing:
      Writing graphical output:                                                        output-annulus_transient/solution/solution-00000
      RMS, max velocity:                                                               1.92 m/s, 3.72 m/s
      Angular momentum, Moment of inertia, Angular velocity, Surface angular velocity: -2.36e+04 kg*m^2/s, 2.36e+04 kg*m^2, -1 1/s, -1 1/s
-     Pressure at top/bottom of domain:                                                2.043e-16 Pa, 1000 Pa
+     Pressure at top/bottom of domain:                                                -1.099e-16 Pa, 1000 Pa
      Computing dynamic topography                                                     
      Errors u_L2, p_L2, rho_L2, topo_L2:                                              1.106354e-02, 3.426667e-01, 9.250441e-04, 1.037264e-03
 
 *** Timestep 1:  t=0.0196519 seconds, dt=0.0196519 seconds
    Skipping temperature solve because RHS is zero.
    Solving density_field system ... 4 iterations.
-   Solving Stokes system... 8+0 iterations.
+   Solving Stokes system... 5+0 iterations.
 
    Postprocessing:
      RMS, max velocity:                                                               1.94 m/s, 3.76 m/s
      Angular momentum, Moment of inertia, Angular velocity, Surface angular velocity: -2.4e+04 kg*m^2/s, 2.36e+04 kg*m^2, -1.02 1/s, -1.02 1/s
-     Pressure at top/bottom of domain:                                                -1.16e-16 Pa, 1000 Pa
+     Pressure at top/bottom of domain:                                                -2.761e-17 Pa, 1000 Pa
      Computing dynamic topography                                                     
      Errors u_L2, p_L2, rho_L2, topo_L2:                                              1.105923e-02, 3.426761e-01, 4.622173e-03, 1.037260e-03
 
 *** Timestep 2:  t=0.0392301 seconds, dt=0.0195783 seconds
    Skipping temperature solve because RHS is zero.
    Solving density_field system ... 3 iterations.
-   Solving Stokes system... 6+0 iterations.
+   Solving Stokes system... 4+0 iterations.
 
    Postprocessing:
      RMS, max velocity:                                                               1.97 m/s, 3.79 m/s
      Angular momentum, Moment of inertia, Angular velocity, Surface angular velocity: -2.45e+04 kg*m^2/s, 2.36e+04 kg*m^2, -1.04 1/s, -1.04 1/s
-     Pressure at top/bottom of domain:                                                1.546e-16 Pa, 1000 Pa
+     Pressure at top/bottom of domain:                                                5.301e-17 Pa, 1000 Pa
      Computing dynamic topography                                                     
      Errors u_L2, p_L2, rho_L2, topo_L2:                                              1.105671e-02, 3.426825e-01, 6.204944e-03, 1.037259e-03
 
 *** Timestep 3:  t=0.05 seconds, dt=0.0107699 seconds
    Skipping temperature solve because RHS is zero.
    Solving density_field system ... 3 iterations.
-   Solving Stokes system... 6+0 iterations.
+   Solving Stokes system... 4+0 iterations.
 
    Postprocessing:
      RMS, max velocity:                                                               1.98 m/s, 3.82 m/s
      Angular momentum, Moment of inertia, Angular velocity, Surface angular velocity: -2.48e+04 kg*m^2/s, 2.36e+04 kg*m^2, -1.05 1/s, -1.05 1/s
-     Pressure at top/bottom of domain:                                                -1.502e-16 Pa, 1000 Pa
+     Pressure at top/bottom of domain:                                                2.65e-17 Pa, 1000 Pa
      Computing dynamic topography                                                     
      Errors u_L2, p_L2, rho_L2, topo_L2:                                              1.105562e-02, 3.426849e-01, 6.461361e-03, 1.037259e-03