Merge pull request #106 from LukasFuchs/main

update tracer averaging and added timer output

Author: LukasFuchs

Project.toml

@@ -16,3 +16,7 @@ Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
 SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
+TimerOutputs = "a759f4b9-e2f1-59dc-863e-4aeb61b1ea8f"
+
+[compat]
+TimerOutputs = "0.5.29"

docs/src/assets/RTI_Growth_Rate_Res_Test_const_NM_arith.png

@@ -1,13 +1,75 @@
 # Examples and Benchmarks
 
-This section provides various one- and two-dimensional examples and benchmark problems for each of the governing equations. The examples demonstrate how to implement different numerical solvers, apply scaling, and evaluate the advantages and limitations of various finite difference schemes.
+The `GeoModBox.jl` provides various one- and two-dimensional examples and benchmark problems for each of the governing equations. The examples demonstrate how to implement different numerical solvers, apply scaling, and evaluate the advantages and limitations of various finite difference schemes.
 
-By clicking on the title of each page, you will be directed to the corresponding Julia script in the [examples directory](https://github.com/GeoSci-FFM/GeoModBox.jl/blob/main/examples).
+By clicking on the title of each document page, you will be directed to the corresponding Julia script in the [examples directory](https://github.com/GeoSci-FFM/GeoModBox.jl/blob/main/examples).
+
+**Advection** 
+- [2-D advection with constant velocity field](./examples/Advection2D.md)
+- [Resolution test of 2-D advection](./examples/AdvectionRestest2D.md)
+
+**Heat Diffusion** 
+- [1-D continental geotherm](./examples/ContinentalGeotherm.md)
+- [Comparison of FD schemes on a Gaussian anomaly](./examples/GaussianDiffusion1D.md)
+- [1-D oceanic geotherm](./examples/OceanicGeotherm.md)
+- [2-D Gaussian anomaly using iterative backward Euler solver](./examples/BackwardEuler_DC.md)
+- [2-D Gaussian anomaly using iterative forward Euler solver](./examples/ForwardEuler_DC.md)
+- [2-D resolution test with Gaussian anomaly](./examples/GaussianDiffusion2D.md)
+- [2-D Poisson equation resolution test](./examples/PoissonRestest.md)
+- [](./examples/PoissonVariablek.md)
+
+**Thermal Convection Models** 
+- [Bottom heated, isoviscous thermal convection](./examples/BottomHeatedConvection.md)
+- [Internally heated, isoviscous thermal convection](./examples/InternallyHeatedConvection.md)
+- [Mixed heated, isoviscous thermal convection](./examples/MixedHeatedConvection.md)
+
+**Stokes Equation**
+- [1D channel flow with constant and depth-dependent viscosity](./examples/ChannelFlow1D.md)
+- [2D falling block benchmark](./examples/FallingBlockBenchmark.md)
+- [2D falling block with constant or variable viscosity (defect correction)](./examples/FallingBlockDC.md)
+- [2D Rayleigh–Taylor instability](./examples/RTI.md)
+- [2D Rayleigh–Taylor instability benchmark](./examples/RTI_growth_rate.md)
+- [2D viscous inclusion problem](./examples/ViscousInclusion.md)
+
+In the following, the runtime for each of the provided examples is listed as a reference. 
+
+| Example                           | Total Runtime                                         |
+| --------------------------------- | ----------------------------------------------------- |
+| **Advection ===**                                                                         |
+| 2D_Advection.jl                   | 1) Upwind: 5.88 s <br> 2) SLF: 6.11 s 3) Semi-lag: 15.7 s <br> 4) Tracers: 172 s |
+| 2D_Advection_ResolutionTest.jl    | 335 s                                                 |
+| **Heat Diffusion ===**                                                                    |
+| *1D* ---                                                                                  |
+| ContinentalGeotherm_1D.jl         | 530 ms                                                | 
+| Heat_1D_discretization.jl         | 2.28 s                                                |
+| OceanicGeotherm_1D.jl             | 431 ms                                                |
+| *2D* ---                                                                                  |
+| BackwardEuler.jl                  | 29.9 s                                                | 
+| ForwardEuler.jl                   | 5.89 s                                                |
+| Gaussian_Diffusion.jl             | 678 s                                                 | 
+| Poisson_RestTest.jl               | 34.7 s                                                |
+| Poisson_variable_k.jl             | 5.24 s                                                |
+| **Thermal Convection Models ===**                                                         |
+| BottomHeated.jl                   | 6.37 h                                                |
+| InternallyHeated.jl               | 5.94 h                                                |
+| MixedHeated.jl                    |                                                       |
+| **Stokes Equation ===**                                                                   |
+| *1D* ---                                                                                  |
+| ChannelFlow_1D.jl                 | 3.61 s                                                |
+| *2D* ---                                                                                  |
+| FallingBlockBenchmark.jl          | 1) Steady State: 8 s <br> 2) Time-Dependent <br> a) Upwind: 74.4 s <br> b) SLF: 170 s <br> c) Semi-lag: 175 s <br> d) Tracers: 233 s                           |
+| FallingBlockConstEta_DC.jl        | 332 ms                                                | 
+| FallingBlockVarEta_DC.jl          | 25.8 s                                                |
+| RTI.jl                            | 154 s                                                 |
+| RTI_GrowthRate.jl                 | 49.4 s                                                |
+| RTI_Growth_Rate_Res_Test.jl       | 250 s                                                 |
+| RTI_Growth_Rate_Res_Test_2.jl     | 261 s                                                 |
+| RTI_Growth_Rate_Res_Test_3.jl     | 446 s                                                 |
 
 > **Note:** In `GeoModBox.jl`, thermal and kinematic boundary conditions are explicitly implemented within the solvers. The absolute values at *ghost nodes* are computed based on the values provided in the boundary condition tuple `BC`. Each tuple specifies the `type` (either Dirichlet or Neumann) and the corresponding `val`ue at each boundary.  
-> For constant velocity boundary conditions, additional values must be defined in `val` for the boundary nodes (e.g., `BC.val.vxW`, `BC.val.vxE`). These additional values are required to directly solve the momentum equation and update the right-hand side of the linear system. Furthermore, these values must be assigned to the initial boundary nodes of the respective velocity fields.  
+> For constant velocity boundary conditions, additional values must be defined in `val` for the boundary nodes (e.g., `BC.val.vxW`, `BC.val.vxE`). These additional values are required to **directly** solve the momentum equation using a the backslash operator and to update the right-hand side of the linear system. Furthermore, if non-zero, these values must be assigned to the initial boundary nodes of the respective velocity fields.  
 > For more details on the implementation of constant velocity boundaries, refer to the documentation of the [viscous inclusion](./examples/ViscousInclusion.md) example or the [initial velocity setup](Ini.md).
 
 > **Note:** By default, the results of time-dependent examples in `GeoModBox.jl` are stored as GIF animations. To visualize solutions at specific time steps without generating a GIF, set the parameter `save_fig = 0`. In this case, individual plots are not saved, so caution is advised when running problems that require multiple time step iterations.
 
-> **Note:** Some examples use *named tuples* to define constants and parameters. Alternatively, *mutable structures* can be used—particularly useful when parameters need to be modified after initialization (e.g., for scaling purposes). A full transition from *named tuples* to *mutable structures* is planned for future versions of `GeoModBox.jl`.
+> **Note:** Some examples use *named tuples* to define constants and parameters. Alternatively, *mutable structures* can be used—particularly useful when parameters need to be modified after initialization (e.g., for scaling purposes). A full transition from *named tuples* to *mutable structures* is planned for future versions of `GeoModBox.jl`.

docs/src/man/Examples.md

@@ -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,3 +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_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 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/RTI_growth_rate.md

@@ -26,17 +26,19 @@ using GeoModBox.AdvectionEquation.TwoD, GeoModBox.Tracers.TwoD
 using GeoModBox.InitialCondition
 using Base.Threads
 using Printf
+using TimerOutputs
 
 function Advection_2D()
-
+to      =   TimerOutput()
 @printf("Running on %d thread(s)\n", nthreads())
 
 save_fig    =   1
 
+@timeit to "Ini" begin
 # Define Numerical Scheme ============================================ #
 # Advection ---
 #   1) upwind, 2) slf, 3) semilag, 4) tracers
-FD          =   (Method     = (Adv=:tracers,),)
+FD          =   (Method     = (Adv=:upwind,),)
 # -------------------------------------------------------------------- #
 # Define Initial Condition =========================================== #
 # Temperature --- 
@@ -126,7 +128,9 @@ D       =   (
     Tmean   =   [0.0],
 )
 # -------------------------------------------------------------------- #
+end
 # Initial Conditions ================================================= #
+@timeit to "IniCondition" begin
 # Temperature ---
 IniTemperature!(Ini.T,M,NC,D,x,y)
 if FD.Method.Adv==:slf
@@ -145,6 +149,7 @@ IniVelocity!(Ini.V,D,BC,NV,Δ,M,x,y)            # [ m/s ]
     end
 end
 @. D.vc        = sqrt(D.vxc^2 + D.vyc^2)
+end
 # -------------------------------------------------------------------- #
 # Time =============================================================== #
 T   =   ( 
@@ -159,6 +164,7 @@ nt          =   ceil(Int,T.tmax[1]/T.Δ[1])
 # -------------------------------------------------------------------- #
 # Tracer Advection =================================================== #
 if FD.Method.Adv==:tracers 
+    @timeit to "TracerIni" begin
     # Tracer Initialization ---
     nmx,nmy     =   3,3
     noise       =   1
@@ -186,6 +192,7 @@ if FD.Method.Adv==:tracers
     end
     # Count marker per cell ---
     CountMPC(Ma,nmark,MPC,M,x,y,Δ,NC,NV,1)
+    end
 end
 # -------------------------------------------------------------------- #
 # Visualize initial condition ======================================== #
@@ -228,6 +235,7 @@ for i=2:nt
     @printf("Time step: #%04d\n ",i)
 
     # Advection ===
+    @timeit to "Advection" begin
     if FD.Method.Adv==:upwind
         upwindc2D!(D.T,D.T_ex,D.vxc,D.vyc,NC,T.Δ[1],Δ.x,Δ.y)
     elseif FD.Method.Adv==:slf
@@ -243,6 +251,7 @@ for i=2:nt
         Markers2Cells(Ma,nmark,MAVG.PC_th,D.T_ex,MAVG.wte_th,D.wte,x,y,Δ,Aparam,0)           
         D.T     .=  D.T_ex[2:end-1,2:end-1]
     end
+    end
 
     display(string("ΔT = ",((maximum(filter(!isnan,D.T))-D.Tmax[1])/D.Tmax[1])*100))
 
@@ -292,6 +301,7 @@ elseif save_fig == 0
     display(plot(p))
 end
 # -------------------------------------------------------------------- #
+display(to)
 end # Function end
 
 Advection_2D()

examples/AdvectionEquation/2D_Advection.jl

@@ -3,6 +3,7 @@ using GeoModBox.AdvectionEquation.TwoD
 using GeoModBox.InitialCondition, GeoModBox.Tracers.TwoD
 using Base.Threads
 using Printf
+using TimerOutputs
 
 @doc raw"""
     Advection_2D_ResTest()
@@ -11,7 +12,7 @@ using Printf
 
 """
 function Advection_2D_ResTest()
-
+to      =   TimerOutput()
 @printf("Running on %d thread(s)\n", nthreads())
 
 nrnxny      =   5
@@ -35,6 +36,7 @@ St      = (
 #   1) RigidBody, 2) ShearCell
 Ini         =   (T=:circle,V=:RigidBody,) 
 # -------------------------------------------------------------------- #
+@timeit to "AdvectionScLoop" begin
 for m = 1:ns # Loop over advection schemes
     # Define Numerical Scheme ======================================== #
     FD          =   (Method     = (Adv=Scheme[m],),)    
@@ -57,6 +59,7 @@ for m = 1:ns # Loop over advection schemes
     )
     BC  =   ()  # dummy
     # ---------------------------------------------------------------- #
+    @timeit to "ResolutionLoop" begin
     for l = 1:nrnxny # Loop over differnet resolutions
         # Numerical Constants ======================================== #
         NC  =   (
@@ -124,6 +127,7 @@ for m = 1:ns # Loop over advection schemes
             Tmean   =   [0.0],
         )
         # Initial Condition ========================================== #
+        @timeit to "IniCondition" begin
         # Temperature ---
         IniTemperature!(Ini.T,M,NC,D,x,y)
         if FD.Method.Adv == "slf"
@@ -142,6 +146,7 @@ for m = 1:ns # Loop over advection schemes
             end
         end
         @. D.vc        = sqrt(D.vxc^2 + D.vyc^2)
+        end
         # ------------------------------------------------------------ #
         # Time ======================================================= #
         T   =   ( 
@@ -156,6 +161,7 @@ for m = 1:ns # Loop over advection schemes
         # ------------------------------------------------------------ #
         # Tracer Advection =========================================== #
         if FD.Method.Adv == "tracers"
+            @timeit to "TracerIni" begin
             # Tracer Initialization ---
             nmx,nmy     =   3,3
             noise       =   1
@@ -184,6 +190,7 @@ for m = 1:ns # Loop over advection schemes
             end
             # Count marker per cell ---
             CountMPC(Ma,nmark,MPC,M,x,y,Δ,NC,NV,1)
+            end
         end
         # ------------------------------------------------------------ #
         # Visualize initial condition -------------------------------- #
@@ -223,6 +230,7 @@ for m = 1:ns # Loop over advection schemes
         for i=2:nt
             #@printf("Time step: #%04d\n ",i) 
 
+            @timeit to "Advection" begin
             if FD.Method.Adv == "upwind"
                 upwindc2D!(D.T,D.T_ex,D.vxc,D.vyc,NC,T.Δ[1],Δ.x,Δ.y)
             elseif FD.Method.Adv == "slf"
@@ -239,7 +247,7 @@ for m = 1:ns # Loop over advection schemes
                 Markers2Cells(Ma,nmark,MAVG.PC_th,D.T_ex,MAVG.wte_th,D.wte,x,y,Δ,Aparam,0)           
                 D.T     .=  D.T_ex[2:end-1,2:end-1]
             end
-            
+            end
             display(string("ΔT = ",((maximum(filter(!isnan,D.T))-D.Tmax[1])/D.Tmax[1])*100))
 
             # Plot Solution ---
@@ -295,9 +303,9 @@ for m = 1:ns # Loop over advection schemes
         # ------------------------------------------------------------ #
 
     end # End resolution loop
-
+    end
 end # End method loop
-
+end
 q   =   plot(0,0,layout=(1,3))
 for m=1:ns    
     plot!(q,St.nxny[m,:],St.Δ[m,:],
@@ -331,7 +339,7 @@ if save_fig == 1 || save_fig == -1
                         Ini.V,"_ResTest.png"))
 end
 # --------------------------------------------------------------------- #
-
+display(to)
 end # Function end
 
 Advection_2D_ResTest()