Merge pull request #39 from GeoSci-FFM/development

Finalize Exercise 9 and 10

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
@@ -231,18 +225,18 @@ for i=2:nt
 
     # Advection ===
     if FD.Method.Adv==:upwind
-        upwindc2D!(D,NC,T,Δ)
+        upwindc2D!(D.T,D.T_ex,D.vxc,D.vyc,NC,T.Δ[1],Δ.x,Δ.y)
     elseif FD.Method.Adv==:slf
-        slfc2D!(D,NC,T,Δ)   
+        slfc2D!(D.T,D.T_ex,D.T_exo,D.vxc,D.vyc,NC,T.Δ[1],Δ.x,Δ.y)
     elseif FD.Method.Adv==:semilag
-        semilagc2D!(D,[],[],x,y,T)
+        semilagc2D!(D.T,D.T_ex,D.vxc,D.vyc,[],[],x,y,T.Δ[1])
     elseif FD.Method.Adv==:tracers
         # Advect tracers ---
         AdvectTracer2D(Ma,nmark,D,x,y,T.Δ[1],Δ,NC,rkw,rkv,1)
         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
@@ -222,18 +222,18 @@ for m = 1:ns # Loop over advection schemes
             #@printf("Time step: #%04d\n ",i) 
 
             if FD.Method.Adv == "upwind"
-                upwindc2D!(D,NC,T,Δ)
+                upwindc2D!(D.T,D.T_ex,D.vxc,D.vyc,NC,T.Δ[1],Δ.x,Δ.y)
             elseif FD.Method.Adv == "slf"
-                slfc2D!(D,NC,T,Δ)   
+                slfc2D!(D.T,D.T_ex,D.T_exo,D.vxc,D.vyc,NC,T.Δ[1],Δ.x,Δ.y)
             elseif FD.Method.Adv == "semilag"
-                semilagc2D!(D,[],[],x,y,T)
+                semilagc2D!(D.T,D.T_ex,D.vxc,D.vyc,[],[],x,y,T.Δ[1])
             elseif FD.Method.Adv == "tracers"
                 # Advect tracers ---
                 AdvectTracer2D(Ma,nmark,D,x,y,T.Δ[1],Δ,NC,rkw,rkv,1)
                 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))

examples/AdvectionEquation/Results/2D_advection_circle_RigidBody_semilag.gif

@@ -111,7 +111,7 @@ where
 
 $$
 \begin{equation}
-\left. c_{W} = \frac{\partial{T}}{\partial{x}} \right\vert_{W},\ and \ \left. c_{E} = \frac{\partial{T}}{\partial{x}} \right\vert_{E}, 
+\left. c_{W} = \frac{\partial{T}}{\partial{x}} \right\vert_{W},\ \textrm{and}\ \left. c_{E} = \frac{\partial{T}}{\partial{x}} \right\vert_{E}, 
 \end{equation}
 $$
 
@@ -186,7 +186,7 @@ $$
 \end{equation}
 $$
 
-### Defection Correction Method
+### Defect Correction Method
 
   The defection correction method is an iterative solution, in which the residual of the conductive part of the *temperature equation* for an initial temperature condition is reduced by a correction term. In case, the system is linear, one iteration is sufficient enough to optain the exact solution. 
 
@@ -256,7 +256,7 @@ where
 
 $$
 \begin{equation}
-a = \frac{\kappa}{\Delta{x}^2},\ and \ b = \frac{1}{\Delta{t}}. 
+a = \frac{\kappa}{\Delta{x}^2},\ \textrm{and} \ b = \frac{1}{\Delta{t}}. 
 \end{equation}
 $$
 
@@ -377,7 +377,7 @@ where $\rho, c_{p}, T, t, k, H,$ and $y$ are the density [kg/m<sup>3</sup>], the
 &emsp;A proper conservative finite difference scheme means that the 1-D vertical heat flux and the thermal conductivity *k* are defined on the vertices and *q* is defined as:
 
 $$\begin{equation}
-\left. q_{y,m} = -k_m \frac{\partial T}{\partial y}\right\vert_{m},\ for\ m = 1:nv, 
+\left. q_{y,m} = -k_m \frac{\partial T}{\partial y}\right\vert_{m},\ \textrm{for}\ m = 1:nv, 
 \end{equation}
 $$
 
@@ -388,7 +388,7 @@ where $nv$ is the number of vertices.
 &emsp;Following the discretization as described above, one needs to solve the following equation for each centroid (in an [explicit finite difference formulation](../../../src/HeatEquation/1Dsolvers.jl)):
 
 $$\begin{equation}
-\rho_j c_{p,j} \frac{T_{j}^{n+1} - T_{j}^{n}}{\Delta{t}} = -\frac{q_{y,j+1}^{n} - q_{y,j}^{n} }{\Delta{y}} + \rho_j H_j,\ for\ j = 1:nc, 
+\rho_j c_{p,j} \frac{T_{j}^{n+1} - T_{j}^{n}}{\Delta{t}} = -\frac{q_{y,j+1}^{n} - q_{y,j}^{n} }{\Delta{y}} + \rho_j H_j,\ \textrm{for}\ j = 1:nc, 
 \end{equation}
 $$
 

examples/AdvectionEquation/Results/2D_advection_circle_RigidBody_slf.gif

@@ -1,5 +1,7 @@
 using GeoModBox.HeatEquation.TwoD, ExtendableSparse, Plots
 
+function Poisson_variable_k()
+
 # Physikalischer Parameter ---------------------------------------------- #
 P       =   (
     L           =   4e3,      #   [m]
@@ -106,6 +108,9 @@ display(q)
 savefig(p,"./examples/DiffusionEquation/2D/Results/Poisson_variable_k_01.png")
 savefig(q,"./examples/DiffusionEquation/2D/Results/Poisson_variable_k_02.png")
 # ----------------------------------------------------------------------- #
+end
+
+Poisson_variable_k()