diff --git a/benchmarks/newton_solver_benchmark_set/nonlinear_channel_flow/README.md b/benchmarks/newton_solver_benchmark_set/nonlinear_channel_flow/README.md
index 5237a686809..7acc5876a84 100644
--- a/benchmarks/newton_solver_benchmark_set/nonlinear_channel_flow/README.md
+++ b/benchmarks/newton_solver_benchmark_set/nonlinear_channel_flow/README.md
@@ -2,8 +2,25 @@
The files in [this directory](https://github.com/geodynamics/aspect/tree/main/benchmarks/newton_solver_benchmark_set/nonlinear_channel_flow)
can be used to recreate the nonlinear channel flow figures from
-{cite:t}`fraters:etal:2019`. Adjust both `metabash.sh` and `bash.sh` to reflect
-the parameter search you are interested in and run the metabash script. If you
-have run exactly the same runs as shown in the paper, you can use the gnuplot
-files in `plotfiles/` to recreate the figures in that paper. This benchmark
-requires the normal nonlinear channel flow benchmark to be compiled.
+{cite:t}`fraters:etal:2019`.
+Start by compiling the `simple_nonlinear` plugin in this directory as described (at the example of another benchmark)
+in {ref}`sec:benchmark-run` .
+After compiling the plugin adjust `metabash.sh` and/or `bash.sh` to reflect
+the parameter search you are interested in and run the script.
+Running the `bash.sh` script without changes will recreate the model runs of Fig. 1 and Fig. 2 in {cite:t}`fraters:etal:2019`
+(by default the script will assume your processor has 4 compute cores, but you can adjust the number in the script).
+After running the model series you can use the gnuplot files in `plotfiles/` to recreate the following figures.
+
+```{figure-md} fig:benchmark-newton-nonlinear-channel-flow-tractions
+
+
+Convergence history for several methods for a rheology with n = 3 where in- and outflow are described by prescribing the traction. Top row: Computations in which we switch abruptly from Picard iterations to Newton iterations. Bottom row: Same as top, but using the residual scaling method in which we gradually introduce the Newton method. Left-hand column: Unmodified Newton iterations. Right-hand column: Results where we applied the SPD stabilizations to the Newton matrix. Horizontal axes: number of the non-linear (outer) iteration. Vertical axes: non-linear residual.
+```
+
+```{figure-md} fig:benchmark-newton-nonlinear-channel-flow-velocities
+
+
+Convergence history for several methods for a rheology with n = 3 where in- and outflow are described by prescribing the velocity. Top row: Computations in which we switch abruptly from Picard iterations to Newton iterations. Bottom row: Same as top, but using the residual scaling method in which we gradually introduce the Newton method. Left-hand column: Unmodified Newton iterations. Right-hand column: Results where we applied the SPD stabilizations to the Newton matrix. Horizontal axes: number of the non-linear (outer) iteration. Vertical axes: non-linear residual.
+```
+
+The nonlinear convergence behavior is shown in {numref}`fig:benchmark-newton-nonlinear-channel-flow-tractions` and {numref}`fig:benchmark-newton-nonlinear-channel-flow-velocities`.
diff --git a/benchmarks/newton_solver_benchmark_set/nonlinear_channel_flow/bash.sh b/benchmarks/newton_solver_benchmark_set/nonlinear_channel_flow/bash.sh
index 05e53f2082c..ce03fd92dc0 100755
--- a/benchmarks/newton_solver_benchmark_set/nonlinear_channel_flow/bash.sh
+++ b/benchmarks/newton_solver_benchmark_set/nonlinear_channel_flow/bash.sh
@@ -2,30 +2,30 @@
version="b1"
build_dir="../../../build/"
-declare -a boundary_type=("v" "t")
-declare -a grid=(5)
-declare -a UFS=("false" "true") # use failsafe
-declare -a NSP=("SPD" "none") # "Symmetric" "PD" "none") #Use Newton stabilization preconditioner
-declare -a NSA=("SPD" "symmetric" "none") # "Symmetric" "PD" "none") #Use Newton stabilization A block
-declare -a agrid=(0)
+declare -a boundary_type=("v" "t") # velocity or traction boundary conditions
+declare -a grid=(4) # global refinement levels
+declare -a agrid=(0) # adaptive refinement levels
+declare -a UFS=("false") # use Newton failsafe
+declare -a NSP=("SPD") # ("Symmetric" "PD" "none") # Use Newton stabilization preconditioner
+declare -a NSA=("SPD" "none") # ("Symmetric" "PD" "none") # Use Newton stabilization A block
declare -a n=(3) #1 2 5 10 25 50 100)
-I=150
-declare -a P=(0 5 10 50 15 150 20 25 30) #0 5 10 15 20 150) # 1 2 3 4 5 10 15 25 50 100)
-declare -a LS=(0 5 10)
-declare -a ST=(1e-20)
-declare -a LT=("1e-5") # "1e-8") # "1e-8") #"1e-5") # "1e-8")
-declare -a NLT=("1e-14") # "1e-9" "1e-10" "1e-14" "1e-20")
-declare -a ABT=("1e-2")
-declare -a RSM=("false" "true") #"true" #"false"
-declare -a SF=("9e-1")
-declare -a UDS=("true")
+I=150 # maximum nonlinear iterations
+declare -a P=(5 10 15 150) # (0 5 10 15 20 25 30 50 150) # Picard iterations
+declare -a LS=(0 100) # (0 5 10 100) # Line search iterations
+declare -a ST=(1e-20) # Nonlinear newton solver switch tolerance
+declare -a LT=("1e-5") # "1e-8" # Linear solver tolerance
+declare -a NLT=("1e-14") # "1e-9" "1e-10" "1e-14" "1e-20" # Nonlinear solver tolerance
+declare -a ABT=("1e-2") # A block solver tolerance
+declare -a RSM=("false" "true") # "true" #"false" # Residual scaling method
+declare -a SF=("9e-1") # SPD safety factor
+declare -a UDS=("false") # ("false" "true") # Use deviator of strain-rate
declare -a AEW=("false") # always use Eisenstat Walker, even for picard
-AV=-1
-COMP="0" #"4e-12"
-declare -a OS=("9e-1" "1e-1" "1e-2" "1e-8") # "1e-3" "1e-4" "1e-5" "1e-6" "1e-6" "1e-7" "1e-8")
-materialmodelnameShort="SNL" #"VP2"
+AV=-1 # Viscosity averaging
+COMP="0" # "4e-12" # compressibility
+declare -a MLT=("9e-1" "1e-8") # ("9e-1" "1e-1" "1e-2" "1e-8") # maximum linear Stokes solver tolerance
+materialmodelnameShort="SNL"
materialmodelname="SNL"
-phi=0 # 0
+phi=0
declare -a vel=(0) #(25 50 125) # 50 125) #25 50 125)
declare -a BV=(1e19) #("1e23" "1e24" "5e24") #"1e23" "1e24" "5e24")
processes=4
@@ -86,7 +86,7 @@ do
do
for i_AEW in "${AEW[@]}"
do
- for i_OS in "${OS[@]}"
+ for i_MLT in "${MLT[@]}"
do
for i_RSM in "${RSM[@]}"
do
@@ -103,7 +103,7 @@ do
U="3.96109558e-10"
fi
- dirname_clean="$version""$materialmodelnameShort""_BT""$i_boundary_type""_""$SOLVER_SHORT""_ST""$i_ST""_UFS""$i_UFS""_NSP-""$i_NSP""_NSA-""$i_NSA""_C""$COMP""_g""$i_grid""_ag""$i_agrid""_AEW""$i_AEW""_UDS""$i_UDS""_SF""$i_SF""_NLT""$i_NLT""_ABT""$i_ABT""_LT""$i_LT""_mLT""$i_OS""_I""$I""_P""$i_P""_EW1""_theta1""_LS""$i_LS""_RSM""$i_RSM""_AV""$AV""_phi""$phi""_vel""$i_vel""_BV""$i_BV""_n""$i_n"
+ dirname_clean="$version""$materialmodelnameShort""_BT""$i_boundary_type""_""$SOLVER_SHORT""_ST""$i_ST""_UFS""$i_UFS""_NSP-""$i_NSP""_NSA-""$i_NSA""_C""$COMP""_g""$i_grid""_ag""$i_agrid""_AEW""$i_AEW""_UDS""$i_UDS""_SF""$i_SF""_NLT""$i_NLT""_ABT""$i_ABT""_LT""$i_LT""_mLT""$i_MLT""_I""$I""_P""$i_P""_EW1""_theta1""_LS""$i_LS""_RSM""$i_RSM""_AV""$AV""_phi""$phi""_vel""$i_vel""_BV""$i_BV""_n""$i_n"
dirname="results/""$dirname_clean"
infilename="$dirname""/input.prm"
outfilename="$dirname""/output.log"
@@ -141,13 +141,13 @@ do
-e "s/set SPD safety factor = .*/set SPD safety factor = $i_SF/g" \
-e "s/set Use deviator of strain-rate = .*/set Use deviator of strain-rate = $i_UDS/g" \
-e "s/set Max Newton line search iterations = .*/ set Max Newton line search iterations = $i_LS/g" \
- -e "s/set Maximum linear Stokes solver tolerance =.*/set Maximum linear Stokes solver tolerance = $i_OS/g" \
+ -e "s/set Maximum linear Stokes solver tolerance =.*/set Maximum linear Stokes solver tolerance = $i_MLT/g" \
-e "s/set Use Newton residual scaling method .*/ set Use Newton residual scaling method = $i_RSM/g" \
$input_name > "$infilename"
nohup mpirun -np $processes $build_dir./aspect $infilename > $outfilename 2>$errorfilename
- grep "Relative nonlinear residual" $outfilename > $outplotfilename
+ grep "Relative nonlinear residual " $outfilename > $outplotfilename
done
done
diff --git a/benchmarks/newton_solver_benchmark_set/nonlinear_channel_flow/doc/figure_t.png b/benchmarks/newton_solver_benchmark_set/nonlinear_channel_flow/doc/figure_t.png
new file mode 100644
index 00000000000..5c466200cee
Binary files /dev/null and b/benchmarks/newton_solver_benchmark_set/nonlinear_channel_flow/doc/figure_t.png differ
diff --git a/benchmarks/newton_solver_benchmark_set/nonlinear_channel_flow/doc/figure_v.png b/benchmarks/newton_solver_benchmark_set/nonlinear_channel_flow/doc/figure_v.png
new file mode 100644
index 00000000000..9ab3c7f883e
Binary files /dev/null and b/benchmarks/newton_solver_benchmark_set/nonlinear_channel_flow/doc/figure_v.png differ
diff --git a/benchmarks/newton_solver_benchmark_set/nonlinear_channel_flow/plotfiles/plot5_b1_g4_UFSfalse_NSPSPD_mlt1e-8_n3_UDSfalse_AEWfalse_BTt.gnuplot b/benchmarks/newton_solver_benchmark_set/nonlinear_channel_flow/plotfiles/plot5_b1_g4_UFSfalse_NSPSPD_mlt1e-8_n3_UDSfalse_AEWfalse_BTt.gnuplot
index 570c824fbef..af685db09be 100644
--- a/benchmarks/newton_solver_benchmark_set/nonlinear_channel_flow/plotfiles/plot5_b1_g4_UFSfalse_NSPSPD_mlt1e-8_n3_UDSfalse_AEWfalse_BTt.gnuplot
+++ b/benchmarks/newton_solver_benchmark_set/nonlinear_channel_flow/plotfiles/plot5_b1_g4_UFSfalse_NSPSPD_mlt1e-8_n3_UDSfalse_AEWfalse_BTt.gnuplot
@@ -1,4 +1,5 @@
-set terminal qt size 1500,1000 enhanced font 'Verdana,14'
+set terminal png size 1500,1000 enhanced font 'Verdana,14'
+set output 'figure_t.png'
set yrange [1e-14:1]
set xrange [0:60]
set format y '%g'
@@ -9,71 +10,70 @@ set origin 0,0.48
unset key
set title "Without stabilization of the velocity block and without the RSM"
plot \
-'results/b1SNL_BTt_NS_ST1e-20_UFSfalse_NSP-SPD_NSA-none_C0_g4_ag0_AEWfalse_UDSfalse_SF9e-1_NLT1e-14_ABT1e-2_LT1e-5_mLT1e-8_I150_P150_EW1_theta1_LS0_RSMfalse_AV-1_phi0_vel0_BV1e19_n3/plot.dat' u 10:11 w l lw 2 lc rgb 'green' t 'Defect correction Picard', \
-'results/b1SNL_BTt_NS_ST1e-20_UFSfalse_NSP-SPD_NSA-none_C0_g4_ag0_AEWfalse_UDSfalse_SF9e-1_NLT1e-14_ABT1e-2_LT1e-5_mLT1e-8_I150_P5_EW1_theta1_LS0_RSMfalse_AV-1_phi0_vel0_BV1e19_n3/plot.dat' u 10:11 w l dt 1 lw 2 lc rgb 'blue' t 'Newton solver with 5 Picard and without line search', \
-'results/b1SNL_BTt_NS_ST1e-20_UFSfalse_NSP-SPD_NSA-none_C0_g4_ag0_AEWfalse_UDSfalse_SF9e-1_NLT1e-14_ABT1e-2_LT1e-5_mLT1e-8_I150_P5_EW1_theta1_LS100_RSMfalse_AV-1_phi0_vel0_BV1e19_n3/plot.dat' u 10:11 w p dt 3 lw 2 lc rgb 'blue' pt 7 t 'Newton solver with 5 Picard and line search', \
-'results/b1SNL_BTt_NS_ST1e-20_UFSfalse_NSP-SPD_NSA-none_C0_g4_ag0_AEWfalse_UDSfalse_SF9e-1_NLT1e-14_ABT1e-2_LT1e-5_mLT1e-8_I150_P10_EW1_theta1_LS0_RSMfalse_AV-1_phi0_vel0_BV1e19_n3/plot.dat' u 10:11 w l dt 1 lw 2 lc rgb 'orange' t 'Newton solver with 10 Picard and without line search', \
-'results/b1SNL_BTt_NS_ST1e-20_UFSfalse_NSP-SPD_NSA-none_C0_g4_ag0_AEWfalse_UDSfalse_SF9e-1_NLT1e-14_ABT1e-2_LT1e-5_mLT1e-8_I150_P10_EW1_theta1_LS100_RSMfalse_AV-1_phi0_vel0_BV1e19_n3/plot.dat' u 10:11 w p dt 3 lw 2 lc rgb 'orange' pt 7 ps 1 t 'Newton solver with 10 Picard and with line search', \
-'results/b1SNL_BTt_NS_ST1e-20_UFSfalse_NSP-SPD_NSA-none_C0_g4_ag0_AEWfalse_UDSfalse_SF9e-1_NLT1e-14_ABT1e-2_LT1e-5_mLT1e-8_I150_P15_EW1_theta1_LS0_RSMfalse_AV-1_phi0_vel0_BV1e19_n3/plot.dat' u 10:11 w l dt 1 lw 2 lc rgb 'red' t 'Newton solver with 15 Picard and without line search', \
-'results/b1SNL_BTt_NS_ST1e-20_UFSfalse_NSP-SPD_NSA-none_C0_g4_ag0_AEWfalse_UDSfalse_SF9e-1_NLT1e-14_ABT1e-2_LT1e-5_mLT1e-8_I150_P15_EW1_theta1_LS100_RSMfalse_AV-1_phi0_vel0_BV1e19_n3/plot.dat' u 10:11 w p dt 3 lw 2 lc rgb 'red' pt 7 ps 1 t 'Newton solver with 15 Picard and with line search', \
-'