Merge pull request #5714 from orionjohnston/add_particle_integration_…

…benchmark

Add particle integration benchmark

benchmarks/particle_integration_scheme/circle.prm

@@ -0,0 +1,80 @@
+set Dimension                              = 2
+set Start time                             = 0
+set End time                               = 1.0
+set Use years in output instead of seconds = false
+set CFL number                             = 1.0
+set Output directory                       = circle_euler_1.0
+set Timing output frequency                = 100
+
+subsection Geometry model
+  set Model name = box
+  subsection Box
+    set X extent = 4
+    set Y extent = 4
+    set Box origin X coordinate = -2
+    set Box origin Y coordinate = -2
+  end
+end
+
+
+subsection Prescribed Stokes solution
+  set Model name = function
+  subsection Velocity function
+    set Variable names = x,y,t
+    set Function expression = 2 * pi * (1+0.0*sin(t))*-y; 2 * pi * (1+0.0*sin(t))*x
+  end
+end
+
+set Nonlinear solver scheme                = Advection only
+
+subsection Initial temperature model
+  set Model name = function
+end
+
+
+subsection Gravity model
+  set Model name = vertical
+
+  subsection Vertical
+    set Magnitude = 0
+  end
+end
+
+
+subsection Material model
+  set Model name = simple
+end
+
+
+subsection Mesh refinement
+  set Initial global refinement                = 4
+  set Initial adaptive refinement              = 0
+  set Time steps between mesh refinement       = 0
+end
+
+############### Parameters describing what to do with the solution
+
+subsection Postprocess
+  set List of postprocessors = visualization, particles
+
+  subsection Visualization
+    set Output format                 = vtu
+    set Time between graphical output = 1
+  end
+
+  subsection Particles
+    set Number of particles = 50
+    set Time between data output = 0.1
+    set Data output format = ascii
+
+    set Particle generator name = ascii file
+    set Integration scheme = rk2
+
+    subsection Generator
+      subsection Ascii file
+        set Data directory = $ASPECT_SOURCE_DIR/benchmarks/particle_integration_scheme/
+        set Data file name = particle.dat
+      end
+    end
+  end
+end

benchmarks/particle_integration_scheme/create.sh

@@ -0,0 +1,31 @@
+#!/bin/bash
+
+for int in euler rk2 rk4; do #euler rk2 rk4
+  rm max_deviation_$int
+
+  for cfl in 1.0 0.5 0.25 0.125 0.0625 0.03125; do
+  # Prepare and run model
+  OUTPUT_DIR=circle_${int}_${cfl}
+  echo set Output directory = $OUTPUT_DIR > output_dir.prm
+  echo set CFL number = $cfl > cfl.prm
+
+  cat circle.prm output_dir.prm cfl.prm ${int}.prm | aspect --
+
+  # format and print results
+    tail -n +8 $OUTPUT_DIR/particles/particles-00000.0000.gnuplot | sort -g -k 3 > sort_temp_0
+    tail -n +8 $OUTPUT_DIR/particles/particles-00001.0000.gnuplot | sort -g -k 3 > sort_temp_1
+    paste sort_temp_0 sort_temp_1 > sort_temp
+    gawk -v cfl=$cfl 'BEGIN{max_deviation=0.0} {deviation=sqrt(($1-$4)*($1-$4)+($2-$5)*($2-$5));if(deviation>max_deviation) {max_deviation=deviation}} END{print cfl, max_deviation}' sort_temp >> max_deviation_$int
+
+
+  done
+done
+
+# make visualization
+echo set Output directory = circle_vis > output_dir.prm
+cat circle.prm vis.prm output_dir.prm | aspect --
+
+# make plot
+#gnuplot plot_convergence.plt
+
+rm output_dir.prm cfl.prm sort_temp_0 sort_temp_1 sort_temp

benchmarks/particle_integration_scheme/doc/particle_integration_scheme.md

@@ -0,0 +1,27 @@
+# Particle integration scheme benchmark
+
+*This section was contributed by Gabriel Johnston and Rene Gassmöller*
+
+This benchmark is designed to test the convergence of various particle integration schemes (Forward-Euler, RK2, RK4) as detailed in {cite}`gassmoller:etal:2018`. The first benchmark involves 100 particles moving in circular motion around a ring in a 2D box. The error is measured as the distance between the final and initial positions, which should ideally be zero. The second benchmark involves particles advected in a flow field defined by {math:numref}`eq:flow`.
+
+```{math}
+:label: eq:flow
+\mathbf{u}_x(t) = 1 + e^t \quad \text{and} \quad \mathbf{u}_y(t) = 1 + e^t.
+```
+The final analytic position is represented by {math:numref}`eq:fposition` and the distance between the initial and final position is defined by {math:numref}`eq:fidistance`.
+
+```{math}
+:label: eq:fposition
+\mathbf{x}(1) = \mathbf{x}_0 + \int_{0}^{1} \mathbf{u}(t) \, dt
+```
+
+```{math}
+:label: eq:fidistance
+d = \| \mathbf{x}_1 - \mathbf{x}_0 \|_2 = \sqrt{2} e
+```
+
+```{figure-md} fig:particle-integration-scheme
+<img src="convergence.*" width="100%" />
+
+Figure shows the convergence order of particle advection schemes in spatially variable flow (a) and temporally changing flow (b). Additionally, the figure (right) depicts the RK4 scheme using the analytical velocity solution for the halfway time instead of interpolating the discrete time steps. Figure adapted from {cite}`gassmoller:etal:2018`
+```

benchmarks/particle_integration_scheme/euler.prm

@@ -0,0 +1,5 @@
+subsection Postprocess
+  subsection Particles
+    set Integration scheme = euler
+  end
+end

benchmarks/particle_integration_scheme/particle.dat

@@ -0,0 +1,100 @@
+1.0000 0.0000
+0.9980 0.0634
+0.9920 0.1266
+0.9819 0.1893
+0.9679 0.2511
+0.9501 0.3120
+0.9284 0.3717
+0.9029 0.4298
+0.8738 0.4862
+0.8413 0.5406
+0.8053 0.5929
+0.7660 0.6428
+0.7237 0.6901
+0.6785 0.7346
+0.6306 0.7761
+0.5801 0.8146
+0.5272 0.8497
+0.4723 0.8815
+0.4154 0.9096
+0.3569 0.9341
+0.2969 0.9549
+0.2358 0.9718
+0.1736 0.9848
+0.1108 0.9938
+0.0476 0.9989
+-0.0159 0.9999
+-0.0792 0.9969
+-0.1423 0.9898
+-0.2048 0.9788
+-0.2665 0.9638
+-0.3271 0.9450
+-0.3863 0.9224
+-0.4441 0.8960
+-0.5000 0.8660
+-0.5539 0.8326
+-0.6056 0.7958
+-0.6549 0.7557
+-0.7015 0.7127
+-0.7453 0.6668
+-0.7861 0.6182
+-0.8237 0.5671
+-0.8580 0.5137
+-0.8888 0.4582
+-0.9161 0.4009
+-0.9397 0.3420
+-0.9595 0.2817
+-0.9754 0.2203
+-0.9874 0.1580
+-0.9955 0.0951
+-0.9995 0.0317
+-0.9995 -0.0317
+-0.9955 -0.0951
+-0.9874 -0.1580
+-0.9754 -0.2203
+-0.9595 -0.2817
+-0.9397 -0.3420
+-0.9161 -0.4009
+-0.8888 -0.4582
+-0.8580 -0.5137
+-0.8237 -0.5671
+-0.7861 -0.6182
+-0.7453 -0.6668
+-0.7015 -0.7127
+-0.6549 -0.7557
+-0.6056 -0.7958
+-0.5539 -0.8326
+-0.5000 -0.8660
+-0.4441 -0.8960
+-0.3863 -0.9224
+-0.3271 -0.9450
+-0.2665 -0.9638
+-0.2048 -0.9788
+-0.1423 -0.9898
+-0.0792 -0.9969
+-0.0159 -0.9999
+0.0476 -0.9989
+0.1108 -0.9938
+0.1736 -0.9848
+0.2358 -0.9718
+0.2969 -0.9549
+0.3569 -0.9341
+0.4154 -0.9096
+0.4723 -0.8815
+0.5272 -0.8497
+0.5801 -0.8146
+0.6306 -0.7761
+0.6785 -0.7346
+0.7237 -0.6901
+0.7660 -0.6428
+0.8053 -0.5929
+0.8413 -0.5406
+0.8738 -0.4862
+0.9029 -0.4298
+0.9284 -0.3717
+0.9501 -0.3120
+0.9679 -0.2511
+0.9819 -0.1893
+0.9920 -0.1266
+0.9980 -0.0634
+1.0000 -0.0000

benchmarks/particle_integration_scheme/plot_convergence.plt

@@ -0,0 +1,63 @@
+set terminal svg size 1500,900 font "Arial,20" enhanced background rgb 'white'
+set output 'convergence.svg'
+
+set multiplot layout 1,2
+
+# Plot A
+set size ratio 1.0
+set origin 0.0,0.0
+set lmargin at screen 0.1
+set rmargin at screen 0.45
+set tmargin at screen 0.9
+set bmargin at screen 0.25
+
+set xrange[2.5e-2:1.2]
+set yrange[1e-14:1.0]
+set xlabel "CFL Number" font "Arial,18"
+set ylabel "L2 errors" font "Arial,18"
+set logscale xy
+
+f(x) = 0.1 * x
+g(x) = 0.005 * x*x
+h(x) = 1e-6 * x*x*x*x
+set key font "Arial,16" outside center bottom samplen 1 maxrows 3 box
+
+set label 1 "A)" at screen 0.25, screen 0.95 font "Arial,22"
+
+plot "spatial/max_deviation_euler" using 1:2 with linespoints linetype 1 lw 3 pointtype 7 ps 1.0 linecolor rgb "#0A5F02" title "Forward Euler", \
+     "spatial/max_deviation_rk2" using 1:2 with linespoints linetype 1 lw 3 pointtype 9 ps 1.0 linecolor rgb "#1c567a" title "RK2", \
+     "spatial/max_deviation_rk4" using 1:2 with linespoints linetype 1 lw 3 pointtype 11 ps 1.0 linecolor rgb "#814292" title "RK4", \
+     f(x) title 'linear' with lines dashtype 2 lw 3 linecolor rgb "gray", \
+     g(x) title 'quadratic' with lines dashtype 3 lw 3 linecolor rgb "gray", \
+     h(x) title 'quartic' with lines dashtype 4 lw 3 linecolor rgb "gray"
+
+# Plot B
+set size ratio 1.0
+set origin 0.5,0.0
+set lmargin at screen 0.55
+set rmargin at screen 0.9
+set tmargin at screen 0.9
+set bmargin at screen 0.25
+
+set xrange[2.5e-2:1.2]
+set yrange[1e-14:1.0]
+set xlabel "CFL Number" font "Arial,18"
+set ylabel "L2 errors" font "Arial,18"
+set logscale xy
+
+f(x) = 0.05 * x
+g(x) = 2e-3 * x*x
+h(x) = 1e-7 * x*x*x*x
+set key font "Arial,16" outside center bottom samplen 1 maxrows 4 box
+
+set label 2 "B)" at screen 0.75, screen 0.95 font "Arial,22"
+
+plot "time_exponential/max_deviation_euler" using 1:2 with linespoints linetype 1 lw 3 pointtype 7 ps 1.0 linecolor rgb "#0A5F02" notitle , \
+     "time_exponential/max_deviation_rk2" using 1:2 with linespoints linetype 1 lw 3 pointtype 8 ps 1.0 linecolor rgb "#1c567a" notitle , \
+     "time_exponential/max_deviation_rk4" using 1:2 with linespoints linetype 1 lw 3 pointtype 11 ps 1.0 linecolor rgb "#814292" notitle, \
+     "time_exponential/max_deviation_rk4_analytic" using 1:2 with linespoints linetype 2 lw 3 pointtype 13 ps 1.0 linecolor rgb "#814292" title "RK4 with analytic halfstep", \
+     f(x) notitle with lines dashtype 2 lw 3 linecolor rgb "gray", \
+     g(x) notitle with lines dashtype 3 lw 3 linecolor rgb "gray", \
+     h(x) notitle with lines dashtype 4 lw 3 linecolor rgb "gray"
+
+unset multiplot

benchmarks/particle_integration_scheme/rk2.prm

@@ -0,0 +1,5 @@
+subsection Postprocess
+  subsection Particles
+    set Integration scheme = rk2
+  end
+end

benchmarks/particle_integration_scheme/rk4.prm

@@ -0,0 +1,5 @@
+subsection Postprocess
+  subsection Particles
+    set Integration scheme = rk4
+  end
+end

doc/modules/changes/20240601_orionjohnston

@@ -0,0 +1,3 @@
+New: Benchmark from Gassmoller et al. 2018 for particle integration schemes.
+<br>
+(Gabriel Johnston, 2024/06/01)

doc/sphinx/user/benchmarks/index.md

@@ -50,6 +50,7 @@ benchmarks/inclusion/doc/inclusion.md
 benchmarks/burstedde/doc/burstedde.md
 benchmarks/slab_detachment/doc/slab_detachment.md
 benchmarks/hollow_sphere/doc/hollow_sphere.md
+benchmarks/particle_integration_scheme/doc/particle_integration_scheme.md
 benchmarks/annulus/doc/annulus.md
 benchmarks/rigid_shear/doc/rigid_shear.md
 ../cookbooks/cookbooks/stokes/doc/stokes.md

tests/benchmark_particle_integration_scheme.prm

@@ -0,0 +1,7 @@
+# A benchmark of particle integration schemes in Gassmoller 2018
+#  See the benchmark folder for more information.
+set Dimension = 2
+
+include $ASPECT_SOURCE_DIR/benchmarks/particle_integration_scheme/circle.prm
+
+set End time = 0.1