diff --git a/benchmarks/annulus/doc/annulus.md b/benchmarks/annulus/doc/annulus.md index 136739890c3..e230e261c97 100644 --- a/benchmarks/annulus/doc/annulus.md +++ b/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. diff --git a/benchmarks/annulus/transient/plot_transient_convergence.py b/benchmarks/annulus/transient/plot_transient_convergence.py index e478be642b6..850e0088e50 100644 --- a/benchmarks/annulus/transient/plot_transient_convergence.py +++ b/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,7 +182,7 @@ 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 @@ -189,7 +190,6 @@ def plot_error_over_resolution(statistics, timestep, output_file): 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') @@ -197,10 +197,17 @@ def plot_error_over_resolution(statistics, timestep, output_file): 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 diff --git a/benchmarks/annulus/transient/transient_annulus.prm b/benchmarks/annulus/transient/transient_annulus.prm index f9cca81da0e..23aadcff9ed 100644 --- a/benchmarks/annulus/transient/transient_annulus.prm +++ b/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 diff --git a/benchmarks/rigid_shear/doc/rigid_shear.md b/benchmarks/rigid_shear/doc/rigid_shear.md index 4ce6ad82cfc..a28ee2a23ef 100644 --- a/benchmarks/rigid_shear/doc/rigid_shear.md +++ b/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. diff --git a/benchmarks/rigid_shear/transient/plot_transient_convergence.py b/benchmarks/rigid_shear/transient/plot_transient_convergence.py index 1461a02bd4b..31a4ea29e89 100644 --- a/benchmarks/rigid_shear/transient/plot_transient_convergence.py +++ b/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,7 +88,7 @@ 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 @@ -93,18 +96,23 @@ def plot_error_over_resolution(statistics, timestep, output_file): 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 diff --git a/doc/sphinx/references.bib b/doc/sphinx/references.bib index 563f7073f7c..ee8c9054a73 100644 --- a/doc/sphinx/references.bib +++ b/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} +} diff --git a/tests/annulus_transient.prm b/tests/annulus_transient.prm index 82f01358973..af155b4e38d 100644 --- a/tests/annulus_transient.prm +++ b/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 diff --git a/tests/annulus_transient/screen-output b/tests/annulus_transient/screen-output index ef70d62865a..18afd612f6c 100644 --- a/tests/annulus_transient/screen-output +++ b/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