Merge branch 'Corrections'

Author: LukasFuchs

docs/src/assets/RTI_Growth_Rate_Res_Test.png

@@ -419,8 +419,7 @@ end
 # --------------------------------------------------------------------- #
 ```
 
-<!-- ![AdvResFinal](../../assets/2D_advection_circle_RigidBody_ResTest.png) -->
-![AdvResFinalA](../../../../examples/AdvectionEquation/Results/2D_advection_circle_RigidBody_ResTest.png)
+![AdvResFinal](../../assets/2D_advection_circle_RigidBody_ResTest.png)
 
 **Figure 2. Advection Resolution Test.**  
 Deviation, maximum, and mean temperature for each advection scheme across increasing grid resolutions. 

docs/src/assets/RTI_Growth_Rate_Res_Test_const_NM_arith.png

@@ -1,8 +1,22 @@
 # [RTI - Growth Rate Benchmark](https://github.com/GeoSci-FFM/GeoModBox.jl/blob/main/examples/StokesEquation/2D/RTI_GrowthRate.jl)
 
-This script performs a benchmark for the growth rate of a Rayleigh–Taylor instability, following Gerya (2009). The benchmark is based on the analytical solution by Ramberg (1968) and is used to assess the accuracy of the velocity field in a purely gravity-driven flow. 
+This script benchmarks the **growth rate of a Rayleigh–Taylor instability**, following *Gerya (2009)*.  
+The benchmark is based on the analytical solution by *Ramberg (1968)* and is used to evaluate the accuracy of the velocity field in a purely gravity-driven flow.  
 
-While small-amplitude perturbations can be analyzed theoretically—facilitated by the use of tracers and bilinear interpolation of density onto centroids—a relatively large perturbation amplitude is employed here to accommodate practical implementation constraints. 
+The script calculates the **diapiric growth rate** at the tip of the perturbation for different **perturbation amplitudes** ($\delta A$), **wavelengths** ($\lambda$), and **viscosity ratios** ($\eta_r$).  
+The numerical solution is plotted alongside the analytical one, which is arbitrarily scaled by certain constants for visualization purposes, following *Gerya (2009)*.  
+
+The amplitude of the cosine perturbation is defined as:
+
+$\begin{equation}
+\delta A = \cos\left( 2\pi \frac{x_m - L/2}{\lambda} \right),
+\end{equation}$
+
+where $x_m$ is the x-coordinate of the marker, and $L$ is the length of the model domain.  
+The perturbation is applied to the tracers by adding the perturbation amplitude to the y-coordinate of an initially equally distributed or randomly perturbed tracer field.  
+
+Using a **bilinear interpolation scheme** to map tracer properties onto the regular numerical grid, one can test how accurately the Stokes solver reproduces the analytical growth rate of the diapiric instability.  
+However, care must be taken when defining the initial tracer positions or when averaging tracer properties onto the grid to avoid numerical artifacts.  
 
 --- 
 
@@ -34,7 +48,7 @@ Pl  =   (
 )
 ```
 
-The parameters for an initial cosinusoidal tracer perturbation are defined for a two-layer model. For benchmarking purposes, a range of wavelengths is specified. 
+The parameters for an initial cosinusoidal tracer perturbation are defined for a two-layer model. For benchmarking purposes, a range of wavelengths is specified. One can define, if the density is interpolated from the tracers to the centroids or the vertices. 
 
 ```Julia
 # Define Initial Condition ========================================== #
@@ -480,6 +494,6 @@ end
 
 **Figure 2. RTI Growth Rate.** Growth rate of an initial cosinusoidal perturbation in a two-layer system across various wavelengths $\lambda$. The growth rate is arbitrarily scaled using $b_1$ and $b_2$ for visualization, following the approach of Gerya (2000). The lines are the analytical solutions for different viscosity ratios $\eta_r$ and the markers show the corresponding numerical results for models with decreasing amplitudes (black - scaled by 15, red - scaled by 150, yellow - scaled by 1500). The rising velocity is numerically calculated following the approach shown in Figure 1. 
 
-<!-- ![RTI_Growth_Rate_Res_test](../../assets/RTI_Growth_Rate_Res_Test.png)
+![RTI_Growth_Rate_Res_test](../../assets/RTI_Growth_Rate_Res_Test_const_NM_arith.png)
 
-**Figure 3. RTI Resolution Test.** Resolution test for the RTI growth rate using a fixed layer thickness (1500 km), a fixed wavelength $\lambda = 4000 \text{ km}$ and fixed numerical resolution of 50x50 centroids for different numbers of markers per cell $\left(nm_x,nm_y\right)$, different perturbation amplitudes $\left(\delta{A}\right)$ (colored markers), different viscosity ratios $\eta_r$, and without (top row) and with (bottom row) additional noise ontop of the initial marker position (before assigning the layer phases). -->
+**Figure 3. RTI Resolution Test.** Resolution test for the RTI growth rate using a fixed layer thickness (1500 km), a fixed wavelength $\lambda = 4000 \text{ km}$ and fixed number of markers per cell (5x5) for different horizontal and vertical grid resolutions $\left(nc_x,nc_y\right)$, different perturbation amplitudes $\left(\delta{A}\right)$ (colored markers), different viscosity ratios $\eta_r$, and without (top row) and with (bottom row) additional noise ontop of the initial marker position (before assigning the layer phases). The viscosity is interpolate from the tracers to the centroids and the vertices using an arithmetic mean. 

docs/src/man/examples/AdvectionRestest2D.md

@@ -13,6 +13,7 @@ function RTI_GrowthRate()
     @timeit to "Ini" begin
     plot_fields     =:no
     save_fig        = 1
+    avgm            =:arith     # Averaging Method for η - default arith
     Pl  =   (
         qinc    =   5, 
         qsc     =   100*(60*60*24*365.25)*5e1,
@@ -259,7 +260,7 @@ function RTI_GrowthRate()
                         end
                     end     
                     # Interpolate Viscosity ---
-                    Markers2Cells(Ma,nmark,MAVG.PC_th,D.ηce,MAVG.wte_th,D.wte,x,y,Δ,Aparam,η)
+                    Markers2Cells(Ma,nmark,MAVG.PC_th,D.ηce,MAVG.wte_th,D.wte,x,y,Δ,Aparam,η;avgm)
                     D.ηc    .=   D.ηce[2:end-1,2:end-1]
                     Markers2Vertices(Ma,nmark,MAVG.PV_th,D.ηv,MAVG.wtv_th,D.wtv,x,y,Δ,Aparam,η)
                     end
@@ -331,8 +332,6 @@ function RTI_GrowthRate()
                     end
 
                     if l == 1 && k == 1 && n == 1 && o == 1
-                        # @show l,k,n,o
-                        # @show (n-1)*size(ηᵣ,2)+o
                         scatter!(q,(1/(maximum(nm)+2)/(maximum(nm)+2),
                                     PP.Kₐ[o]),
                                     ms=ms[1],markershape=:circle,
@@ -342,8 +341,8 @@ function RTI_GrowthRate()
                                     title=string("ηᵣ = ",ηᵣ[o]),
                                     layout=(size(addnoise,2),size(ηᵣ,2)),
                                     subplot=((n-1)*size(ηᵣ,2)+o),
-                                    xlims=(1/(maximum(nm)+3)/(maximum(nm)+3), .5),
-                                    ylims=(PP.Kₐ[o]-0.5*PP.Kₐ[o], PP.Kₐ[o]+0.5*PP.Kₐ[o])
+                                    # xlims=(1/(maximum(nm)+3)/(maximum(nm)+3), .5),
+                                    # ylims=(PP.Kₐ[o]-0.5*PP.Kₐ[o], PP.Kₐ[o]+0.5*PP.Kₐ[o])
                                     )
                     elseif l == 1 && k == 1 
                         scatter!(q,(1/(maximum(nm)+2)/(maximum(nm)+2),
@@ -354,8 +353,8 @@ function RTI_GrowthRate()
                                     title=string("ηᵣ = ",ηᵣ[o]),
                                     layout=(size(addnoise,2),size(ηᵣ,2)),
                                     subplot=((n-1)*size(ηᵣ,2)+o),
-                                    xlims=(1/(maximum(nm)+3)/(maximum(nm)+3), .5),
-                                    ylims=(PP.Kₐ[o]-0.5*PP.Kₐ[o], PP.Kₐ[o]+0.5*PP.Kₐ[o])
+                                    # xlims=(1/(maximum(nm)+3)/(maximum(nm)+3), .5),
+                                    # ylims=(PP.Kₐ[o]-0.5*PP.Kₐ[o], PP.Kₐ[o]+0.5*PP.Kₐ[o])
                                     )
                     end
                     if plot_fields==:yes
@@ -425,7 +424,7 @@ function RTI_GrowthRate()
     end # Loop addnoise - n
     end
     if save_fig == 1
-        savefig(q,string("./examples/StokesEquation/2D/Results/RTI_Growth_Rate_Res_Test.png"))
+        savefig(q,string("./examples/StokesEquation/2D/Results/RTI_Growth_Rate_Res_Test_",avgm,".png"))
     else
         display(q)
     end

docs/src/man/examples/RTI_growth_rate.md

@@ -13,6 +13,7 @@ function RTI_GrowthRate()
     @timeit to "Ini" begin
     plot_fields     =:no
     save_fig        = 1
+    avgm            =:geom     # Averaging Method for η - default arith
     Pl  =   (
         qinc    =   5, 
         qsc     =   100*(60*60*24*365.25)*5e1,
@@ -261,7 +262,7 @@ function RTI_GrowthRate()
                         end
                     end     
                     # Interpolate Viscosity ---
-                    Markers2Cells(Ma,nmark,MAVG.PC_th,D.ηce,MAVG.wte_th,D.wte,x,y,Δ,Aparam,η)
+                    Markers2Cells(Ma,nmark,MAVG.PC_th,D.ηce,MAVG.wte_th,D.wte,x,y,Δ,Aparam,η;avgm)
                     D.ηc    .=   D.ηce[2:end-1,2:end-1]
                     Markers2Vertices(Ma,nmark,MAVG.PV_th,D.ηv,MAVG.wtv_th,D.wtv,x,y,Δ,Aparam,η)
                     end
@@ -383,8 +384,8 @@ function RTI_GrowthRate()
                                     xlabel="1/nmx/nmy",ylabel="ε [ % ]",
                                     xscale=:log10, yscale=:log10,
                                     title=string("ηᵣ = ",ηᵣ[o]),
-                                    xlims=(1/(maximum(nm)+3)/(maximum(nm)+3), .5),
-                                    ylims=(1e-1, 1e2),
+                                    # xlims=(1/(maximum(nm)+3)/(maximum(nm)+3), .5),
+                                    # ylims=(1e-1, 1e2),
                                     layout=(size(addnoise,2),size(ηᵣ,2)),
                                     subplot=((n-1)*size(ηᵣ,2)+o))
                     else
@@ -394,8 +395,8 @@ function RTI_GrowthRate()
                                     xlabel="1/nmx/nmy",ylabel="ε [ % ]",
                                     xscale=:log10,yscale=:log10,
                                     title=string("ηᵣ = ",ηᵣ[o]),
-                                    xlims=(1/(maximum(nm)+3)/(maximum(nm)+3), .5),
-                                    ylims=(1e-1, 1e2),
+                                    # xlims=(1/(maximum(nm)+3)/(maximum(nm)+3), .5),
+                                    # ylims=(1e-1, 1e2),
                                     layout=(size(addnoise,2),size(ηᵣ,2)),
                                     subplot=((n-1)*size(ηᵣ,2)+o))
                     end
@@ -408,7 +409,7 @@ function RTI_GrowthRate()
     end # Loop addnoise - n
     end
     if save_fig == 1
-        savefig(q,string("./examples/StokesEquation/2D/Results/RTI_Growth_Rate_Res_Test_const_NC.png"))
+        savefig(q,string("./examples/StokesEquation/2D/Results/RTI_Growth_Rate_Res_Test_const_NC_",avgm,".png"))
     else
         display(q)
     end

examples/StokesEquation/2D/RTI_Growth_Rate_Res_Test.jl

@@ -13,6 +13,7 @@ function RTI_GrowthRate()
     @timeit to "Ini" begin
     plot_fields     =:no
     save_fig        = 1
+    avgm            =:geom     # Averaging Method for η - default arith
     Pl  =   (
         qinc    =   5, 
         qsc     =   100*(60*60*24*365.25)*5e1,
@@ -261,7 +262,7 @@ function RTI_GrowthRate()
                         end
                     end     
                     # Interpolate Viscosity ---
-                    Markers2Cells(Ma,nmark,MAVG.PC_th,D.ηce,MAVG.wte_th,D.wte,x,y,Δ,Aparam,η)
+                    Markers2Cells(Ma,nmark,MAVG.PC_th,D.ηce,MAVG.wte_th,D.wte,x,y,Δ,Aparam,η;avgm)
                     D.ηc    .=   D.ηce[2:end-1,2:end-1]
                     Markers2Vertices(Ma,nmark,MAVG.PV_th,D.ηv,MAVG.wtv_th,D.wtv,x,y,Δ,Aparam,η)
                     end
@@ -383,8 +384,8 @@ function RTI_GrowthRate()
                                     xlabel="1/ncx/ncy",ylabel="ε [ % ]",
                                     xscale=:log10, yscale=:log10,
                                     title=string("ηᵣ = ",ηᵣ[o]),
-                                    xlims=(1/(maximum(nc)+20)/(maximum(nc)+20), .1),
-                                    ylims=(1e-2, 1e2),
+                                    # xlims=(1/(maximum(nc)+20)/(maximum(nc)+20), .1),
+                                    # ylims=(1e-2, 1e2),
                                     layout=(size(addnoise,2),size(ηᵣ,2)),
                                     subplot=((n-1)*size(ηᵣ,2)+o))
                     else
@@ -394,8 +395,8 @@ function RTI_GrowthRate()
                                     xlabel="1/ncx/ncy",ylabel="ε [ % ]",
                                     xscale=:log10,yscale=:log10,
                                     title=string("ηᵣ = ",ηᵣ[o]),
-                                    xlims=(1/(maximum(nc)+20)/(maximum(nc)+20), .1),
-                                    ylims=(1e-2, 1e2),
+                                    # xlims=(1/(maximum(nc)+20)/(maximum(nc)+20), .1),
+                                    # ylims=(1e-2, 1e2),
                                     layout=(size(addnoise,2),size(ηᵣ,2)),
                                     subplot=((n-1)*size(ηᵣ,2)+o))
                     end
@@ -408,7 +409,7 @@ function RTI_GrowthRate()
     end # Loop addnoise - n
     end
     if save_fig == 1
-        savefig(q,string("./examples/StokesEquation/2D/Results/RTI_Growth_Rate_Res_Test_const_NM.png"))
+        savefig(q,string("./examples/StokesEquation/2D/Results/RTI_Growth_Rate_Res_Test_const_NM_",avgm,".png"))
     else
         display(q)
     end