Merge pull request #40 from LukasFuchs/main

Finalizing falling block benchmark

Author: LukasFuchs

README.md

@@ -25,7 +25,7 @@ The **Geod**ynamic **Mod**elling Tool**Box** is a julia package mainly used for
 - Solution of the conductive part -->
 ------------------
 ### [Heat Diffusion Equation](./examples/DiffusionEquation/1D/README.md)
-  The **GeoModBox** provides different finite difference (**FD**) schemes (e.g., [forward](./src/HeatEquation/ForwardEuler.jl) and [backward](./src/HeatEquation/BackwardEuler.jl) Euler, [Crank-Nicholson approach](./src/HeatEquation/CNA.jl), [ADI](./src/HeatEquation/ADI.jl)) to solve the diffusive part of the time-dependent or steady-state *temperature equation* including radioactive heating, in [1-D](./examples/DiffusionEquation/1D/README.md) and [2-D](./examples/DiffusionEquation/2D/README.md). The solvers are located in [src/HeatEquation](./src/HeatEquation/). So far, only *Dirichlet* and *Neumann* thermal boundary conditions are available. Currently, most of the functions assume constant thermal parameters. 
+  The **GeoModBox** provides different finite difference (**FD**) schemes (e.g., [forward](./src/HeatEquation/ForwardEuler.jl) and [backward](./src/HeatEquation/BackwardEuler.jl) Euler, [Crank-Nicholson approach](./src/HeatEquation/CNA.jl), [ADI](./src/HeatEquation/ADI.jl)) to solve the *diffusive part* of the time-dependent or steady-state *temperature equation* including radioactive heating, in [1-D](./examples/DiffusionEquation/1D/README.md) and [2-D](./examples/DiffusionEquation/2D/README.md). The solvers are located in [src/HeatEquation](./src/HeatEquation/). So far, only *Dirichlet* and *Neumann* thermal boundary conditions are available. Currently, most of the functions assume constant thermal parameters (except for the 1-D solvers). 
 
 The examples of solving the *heat diffusion equation* include, amongst others: 
 - the determination of an [oceanic](./examples/DiffusionEquation/1D/OceanicGeotherm_1D.jl) and [continental](./examples/DiffusionEquation/1D/ContinentalGeotherm_1D.jl) 1-D geotherm profile, 
@@ -44,6 +44,25 @@ All numerical schemes methods can be used in the [thermal convection code]() and
 ------------------ -->
 
 ### [Heat Advection Equation](./examples/AdvectionEquation/README.md)
+
+  To solve the *advective part* of the *temperature equation*, the ```GeoModBox.jl``` provides the following different methods: 
+- an upwind scheme,
+- the staggered -leaped frog scheme, 
+- a semi-lagrangian advection, and
+- passive tracers/markers. 
+
+The solvers for a the tracer advection method are located in [./src/Tracers](./src/Tracers/), where the remaining advection routines are located in [./src/AdvectionEquation](./src/AdvectionEquation/). The routines are structured in such a way, that any property, as long as the property is defined on the *centroids* including *ghost nodes* on all boundaries, can be advected with the first **three** advection methods listed above. Using passive tracers, one can so far choose to either advect the temperature or to advect the phases. In case of advecting phases, one can define a certain rheology ($\eta$) or density ($\rho$) associated to each phase. The phase ID is used to interpolate the corresponding property from the tracers to the *centroids*. 
+
+  A key aspect for the advection equation is the conservation of the **amplitude** and the **shape** of an anomaly. Depending on the method, numerical diffusion or interpolation effects lead to strong deviations of the initial anomaly. For more details see [here](./examples/AdvectionEquation/README.md). The ```GeoModBox.jl``` contains several [routines](./src/InitialCondition/2Dini.jl) to setup a certain initial anomaly, either for properties defined on their correspondig grid (i.e., temperature, velocity, or phase) or for tracers advecting (so far) a certain temperature or phase. Within the examples and the exercise one can choose differen initial temperature and velocity conditions.
+
+  The examples for a two dimensional advection problem include:
+- [a 2-D advection, assuming a constant velocity field (e.g., a rigid body rotation)](./examples/AdvectionEquation/2D_Advection.jl), and
+- [a resolution test of the same advection example](./examples/AdvectionEquation/2D_Advection_ResolutionTest.jl). 
+
+  The exercises include a: 
+- [1-D advection of a gaussian or block anomaly](./exercises/06_1D_Advection.ipynb), and
+- [a 2-D advection coupled with the solution of the diffusion equation](./exercises/07_2D_Energy_Equation.ipynb).
+
 <!-- 
 - Different properties can be advected! 
 ### Discretization Schemes
@@ -54,9 +73,34 @@ All numerical schemes methods can be used in the [thermal convection code]() and
 ------------------
 ------------------
 ## [Momentum Conservation Equation](./examples/StokesEquation/README.md)
+
+<!--
+- Constant viscosity
+- Variable viscosity
+- Coefficients assembly 
+- Right-hand side updates
+- Direct solution 
+- Defect correction
+-->
+
+------------------
+------------------
+## Code Structure
+
+### Initial Conditions 
+
+<!---
+- IniTemperature, IniVelocity, IniPhase,
+- Tracer Initialization
+- ...
+-->
+
+### Scaling
+
 ------------------
 ------------------
 ## [Benchmarks](./examples/)
+
 ### [Gaussian Temperature Diffusion](./examples/DiffusionEquation/2D/Gaussian_Diffusion.jl)
 <img src="./examples/DiffusionEquation/2D/Results/Gaussian_Diffusion_CNA_nx_100_ny_100.gif" alt="drawing" width="600"/> <br>
 **Figure 1. Gaussian Diffusion.** Time-dependent, diffusive solution of a 2-D Gaussian temperature anomaly using the [Crank-Nicholson approach](./src/HeatEquation/CNA.jl) in comparison to its analytical solution. Top Left: 2-D temperature field of the numerical solution and isotherms lines of the numerical (solid black) and analytical (dashed yellow) solution. Top Right: Total deviation to the analytical solution. Bottom Left: 1-D y-profile along x=0. Bottom Right: Root Mean Square total devation of the temperature over time. 
@@ -74,6 +118,23 @@ All numerical schemes methods can be used in the [thermal convection code]() and
 
 **Figure 3. Rigid-Body-Rotation.** Time-dependent solution of a rotating circular temperature anomaly using the **upwind (first)**, **semi-lagrangian (second)**, and **tracer (third)** method. Within a circular area of our model domain the velocity is set to the velocity of a rigid rotation and outside euqal to zero. The temperature is scaled by the maximum temperature of the anomaly. 
 
+### [Falling Block](./examples/StokesEquation/2D/FallingBlockBenchmark.jl)
+
+<img src="./examples/StokesEquation/2D/Results/Falling_block_ηr_0.0_tracers.gif" alt="drawing" width="500"/> <br>
+
+**Figure 4. Isoviscous Falling Block.** Time-dependent solution of an isoviscous falling block example. The problem is solved with a solver for variable viscosities. The tracers advect the phase ID, which is used to interpolate the density and viscosity on the centroids and vertices, respectively. 
+
+<img src="./examples/StokesEquation/2D/Results/FallingBlock_SinkingVeloc_tracers.png" alt="drawing" width="500"/> <br>
+
+**Figure 5. Falling Block Sinking Velocity.** Sinking velocity of the block with respect to the viscosity ratio $\eta_r$ at the initial condition. 
+
+<img src="./examples/StokesEquation/2D/Results/FallingBlock_FinalStage_tracers.png" alt="drawing" width="500"/> <br>
+
+**Figure 6. Falling Block Benchmark.** Final tracers distribution for specific cases with $\eta_r \ge 0$. 
+
+------------------
+------------------
+
 <!--
 - Blanckenbach
 - Channel Flow
@@ -83,6 +144,3 @@ All numerical schemes methods can be used in the [thermal convection code]() and
 - Rigid Body Rotation, check! 
 - Viscous Inclusion
 -->
-
-------------------
-------------------

examples/AdvectionEquation/2D_Advection.jl

@@ -27,12 +27,6 @@ using GeoModBox.InitialCondition
 using Base.Threads
 using Printf
 
-@doc raw"""
-    Advection_2D()
-
-...
-
-"""
 function Advection_2D()
 
 @printf("Running on %d thread(s)\n", nthreads())
@@ -175,7 +169,7 @@ if FD.Method.Adv==:tracers
         wt_th   =   [similar(D.wt) for _ = 1:nthreads()], # per thread
     )
     MPC     =   merge(MPC,MPC1)
-    Ma      =   IniTracer2D(Aparam,nmx,nmy,Δ,M,NC,noise)
+    Ma      =   IniTracer2D(Aparam,nmx,nmy,Δ,M,NC,noise,0,0)
     # RK4 weights ---
     rkw     =   1.0/6.0*[1.0 2.0 2.0 1.0]   # for averaging
     rkv     =   1.0/2.0*[1.0 1.0 2.0 2.0]   # for time stepping
@@ -242,7 +236,7 @@ for i=2:nt
         CountMPC(Ma,nmark,MPC,M,x,y,Δ,NC,i)
 
         # Interpolate temperature from tracers to grid ---
-        Markers2Cells(Ma,nmark,MPC.PG_th,D.T,MPC.wt_th,D.wt,x,y,Δ,Aparam)           
+        Markers2Cells(Ma,nmark,MPC.PG_th,D.T,MPC.wt_th,D.wt,x,y,Δ,Aparam,0)           
         D.T_ex[2:end-1,2:end-1]     .= D.T
     end
 

examples/AdvectionEquation/2D_Advection_ResolutionTest.jl

@@ -167,7 +167,7 @@ for m = 1:ns # Loop over advection schemes
                 wt_th   =   [similar(D.wt) for _ = 1:nthreads()], # per thread
             )
             MPC     =   merge(MPC,MPC1)
-            Ma      =   IniTracer2D(Aparam,nmx,nmy,Δ,M,NC,noise)
+            Ma      =   IniTracer2D(Aparam,nmx,nmy,Δ,M,NC,noise,Ini.p,phase)
             # RK4 weights ---
             rkw     =   1.0/6.0*[1.0 2.0 2.0 1.0]   # for averaging
             rkv     =   1.0/2.0*[1.0 1.0 2.0 2.0]   # for time stepping
@@ -233,7 +233,7 @@ for m = 1:ns # Loop over advection schemes
                 CountMPC(Ma,nmark,MPC,M,x,y,Δ,NC,i)
 
                 # Interpolate temperature from tracers to grid ---
-                Markers2Cells(Ma,nmark,MPC.PG_th,D.T,MPC.wt_th,D.wt,x,y,Δ,Aparam)        
+                Markers2Cells(Ma,nmark,MPC.PG_th,D.T,MPC.wt_th,D.wt,x,y,Δ,Aparam,0)        
             end
 
             display(string("ΔT = ",((maximum(filter(!isnan,D.T))-D.Tmax[1])/D.Tmax[1])*100))