add examples of usage of GeoParams in basic PT codes

GeoParams_PT_examples/.vscode/settings.json

@@ -0,0 +1,2 @@
+{
+}

GeoParams_PT_examples/Project.toml

@@ -0,0 +1,11 @@
+name = "GeoParams_PT_examples"
+uuid = "499c9791-4cbf-4718-8130-fce925879af7"
+authors = ["Duretz Thibault <[email protected]>"]
+version = "0.1.0"
+
+[deps]
+GeoParams = "e018b62d-d9de-4a26-8697-af89c310ae38"
+LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
+Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
+Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"

GeoParams_PT_examples/src/Poisson2D_PT_v0.jl

@@ -0,0 +1,113 @@
+# Initialisation
+using Plots, Printf, Statistics, LinearAlgebra
+Dat = Float64  # Precision (double=Float64 or single=Float32)
+# Macros
+@views    av(A) = 0.25*(A[1:end-1,1:end-1].+A[2:end,1:end-1].+A[1:end-1,2:end].+A[2:end,2:end])
+@views av_xa(A) =  0.5*(A[1:end-1,:].+A[2:end,:])
+@views av_ya(A) =  0.5*(A[:,1:end-1].+A[:,2:end]) 
+@views   avh(A) = ( 0.25./A[1:end-1,1:end-1] .+ 0.25./A[1:end-1,2:end-0] .+ 0.25./A[2:end-0,1:end-1] .+ 0.25./A[2:end-0,2:end-0]).^(-1)
+# 2D Poisson routine
+@views function Poisson2D()
+    # Physics
+    Lx, Ly   = 6.0, 6.0  # domain size
+    T_West   = 0.0
+    T_East   = 0.0
+    T_South  = 0.0
+    T_North  = 0.0
+    slope    = -15.
+    k1       = 1
+    k2       = 1e4
+    # Numerics
+    ncx, ncy = 51, 51    # numerical grid resolution
+    ε        = 1e-6      # nonlinear tolerence
+    iterMax  = 1e4       # max number of iters
+    nout     = 500       # check frequency
+    # Iterative parameters -------------------------------------------
+    Reopt    = 0.625*π
+    cfl      = 0.32
+    ρ        = cfl*Reopt/ncx
+    nsm      = 5
+    # Reopt    = 0.625*π /3
+    # cfl      = 0.3*4
+    # ρ        = cfl*Reopt/ncx
+    # nsm      = 10
+    # Preprocessing
+    Δx, Δy  = Lx/ncx, Ly/ncy
+    # Array initialisation
+    T         =   zeros(Dat, ncx+2, ncy+2)
+    ∂T∂x      =   zeros(Dat, ncx+1, ncy  )
+    ∂T∂y      =   zeros(Dat, ncx  , ncy+1)
+    qx        =   zeros(Dat, ncx+1, ncy  )
+    qy        =   zeros(Dat, ncx  , ncy+1)
+    RT        =   zeros(Dat, ncx  , ncy  )
+    dTdτ      =   zeros(Dat, ncx  , ncy  ) 
+    Δτc       =   zeros(Dat, ncx  , ncy  )       
+    kv        = k1*ones(Dat, ncx+1, ncy+1)
+    kx        =   zeros(Dat, ncx+1, ncy+0)
+    ky        =   zeros(Dat, ncx+0, ncy+1)
+    H         =    ones(Dat, ncx  , ncy  )
+    # Initialisation
+    xc, yc    = LinRange(-Lx/2+Δx/2, Lx/2-Δx/2, ncx+0), LinRange(-Ly/2+Δy/2, Ly/2-Δy/2, ncy+0)
+    xce, yce  = LinRange(-Lx/2-Δx/2, Lx/2+Δx/2, ncx+2), LinRange(-Ly/2-Δy/2, Ly/2+Δy/2, ncy+2)
+    xv, yv    = LinRange(-Lx/2, Lx/2, ncx+1), LinRange(-Ly/2, Ly/2, ncy+1)
+    (xv2,yv2) = ([x for x=xv,y=yv], [y for x=xv,y=yv])
+    below = yv2 .< xv2.*tand.(slope)
+    kv[below] .= k2
+    kx .= ( 0.5./kv[:,1:end-1] .+ 0.5./kv[:,2:end-0]).^(-1)
+    ky .= ( 0.5./kv[1:end-1,:] .+ 0.5./kv[2:end-0,:]).^(-1)
+    kc  = ( 0.25./kv[1:end-1,1:end-1] .+ 0.25./kv[1:end-1,2:end-0] .+ 0.25./kv[2:end-0,1:end-1] .+ 0.25./kv[2:end-0,2:end-0]).^(-1)
+    for it = 1:nsm
+        kv[2:end-1,2:end-1] .= av(kc)
+        kc = av(kv)
+    end
+    # Time loop
+    it=1; iter=1; err=2*ε; err_evo1=[]; err_evo2=[];
+    while (err>ε && iter<=iterMax)
+        # BCs
+        T[:,1]   .= 2*T_South .- T[:,2]     # S
+        T[:,end] .= 2*T_North .- T[:,end-1] # N
+        T[1,:]   .= 2*T_West  .- T[2,:]     # W
+        T[end,:] .= 2*T_East  .- T[end-1,:] # E
+        # Kinematics
+        ∂T∂x .= diff(T[:,2:end-1], dims=1)./Δx 
+        ∂T∂y .= diff(T[2:end-1,:], dims=2)./Δy 
+        # Stresses
+        qx .= -kx.*∂T∂x
+        qy .= -ky.*∂T∂y
+        # Residuals
+        RT .= .-diff(qx, dims=1)./Δx .- diff(qy, dims=2)./Δy + H
+        # PT time step -----------------------------------------------
+        Δτc  .= ρ*min(Δx,Δy)^2 ./ kc ./ 4.1 * cfl 
+        # Calculate rate update --------------------------------------
+        dTdτ          .= (1-ρ) .* dTdτ .+ RT
+        # Update velocity and pressure -------------------------------
+        T[2:end-1,2:end-1]  .+= Δτc ./ ρ .* dTdτ
+        # convergence check
+        if mod(iter, nout)==0
+            norm_RT = norm(RT)/sqrt(length(RT))
+            err = maximum([norm_RT])
+            push!(err_evo1, err); push!(err_evo2, itg)
+            @printf("it = %03d, iter = %04d, err = %1.3e norm[RT=%1.3e] \n", it, itg, err, norm_RT)
+        end
+        iter+=1; global itg=iter
+    end
+    # Plotting
+    p1 = heatmap(xce, yce,         T', aspect_ratio=1, xlims=(-Lx/2, Lx/2), ylims=(-Ly/2, Ly/2), c=:turbo, title="T")
+    p2 = heatmap( xc, yv, log10.(ky)', aspect_ratio=1, xlims=(-Lx/2, Lx/2), ylims=(-Ly/2, Ly/2), c=:turbo, title="ky")
+    p3 = heatmap( xc, yc, log10.(kc)', aspect_ratio=1, xlims=(-Lx/2, Lx/2), ylims=(-Ly/2, Ly/2), c=:turbo, title="kc")
+    y_int = xv.*tand.(slope)
+    p4 = plot(xv, y_int, aspect_ratio=1, xlims=(-Lx/2, Lx/2), ylims=(-Ly/2, Ly/2), color=:black, label=:none)
+    # lines
+    for j=1:ncy+1
+        y_line = yv[j]*ones(size(xv))
+        p4 = plot!(xv, y_line, aspect_ratio=1, xlims=(-Lx/2, Lx/2), ylims=(-Ly/2, Ly/2), color=:gray, label=:none)
+    end
+    for i=1:ncx+1
+        x_line = xv[i]*ones(size(yv))
+        p4 = plot!(x_line, yv, aspect_ratio=1, xlims=(-Lx/2, Lx/2), ylims=(-Ly/2, Ly/2), color=:gray, label=:none)
+    end
+    display(plot(p1, p2, p3, p4))
+    return
+end
+
+Poisson2D()

GeoParams_PT_examples/src/Stokes2D_VE_bench_GeoParams.jl

@@ -0,0 +1,152 @@
+# Initialisation
+using Plots, Printf, Statistics, LinearAlgebra, GeoParams
+Dat = Float64  # Precision (double=Float64 or single=Float32)
+# Macros
+@views    av(A) = 0.25*(A[1:end-1,1:end-1].+A[2:end,1:end-1].+A[1:end-1,2:end].+A[2:end,2:end])
+@views av_xa(A) =  0.5*(A[1:end-1,:].+A[2:end,:])
+@views av_ya(A) =  0.5*(A[:,1:end-1].+A[:,2:end]) 
+# Rheology
+@views function UpdateStressGeoParams!( τxx, τyy, τxy, ε̇xx, ε̇yy, ε̇xy, τxx0, τyy0, τxy0, MatParam, Δt )
+    ε̇iic  = ones(size(ε̇xx))
+    ε̇iiv  = ones(size(ε̇xy))
+    τii0c = sqrt.(1//2 .*(τxx0.^2 .+ τyy0.^2) .+ av(τxy0).^2)
+    τii0v = zeros(size(ε̇xy))
+    τii0v[2:end-1,2:end-1] = sqrt.(1//2 .*( av(τxx0).^2 .+ av(τyy0).^2) .+ τxy0[2:end-1,2:end-1].^2)
+    Phase = 1
+    # Centroids
+    for i in eachindex(τxx)
+        v      = MatParam[Phase].CompositeRheology[1]
+        args   = (; τII_old = τii0c[i], dt=Δt)
+        η_eff  = computeViscosity_εII(v, ε̇iic[i], args)
+        τxx[i] = 2*η_eff*ε̇xx[i]
+        τyy[i] = 2*η_eff*ε̇yy[i]
+    end
+    # Vertices
+    for i in eachindex(τxy)
+        v      = MatParam[Phase].CompositeRheology[1]
+        args   = (; τII_old = τii0v[i], dt=Δt)
+        η_eff  = computeViscosity_εII(v, ε̇iiv[i], args)
+        τxy[i] = 2*η_eff*ε̇xy[i] 
+    end
+end
+# 2D Stokes routine
+@views function Stokes2D_VE_benchmark()
+    # Physics
+    Lx, Ly  = 1.0, 1.0  # domain size
+    ξ       = 10.0      # Maxwell relaxation time
+    η0      = 1.0       # viscous viscosity
+    G       = 1.0       # elastic shear modulus
+    εbg     = 1.0       # background strain-rate
+
+    MatParam = (SetMaterialParams(Name="Matrix", Phase=1,
+                CompositeRheology = CompositeRheology(ConstantElasticity(G=G),LinearViscous(η=η0))),)
+
+    # Numerics
+    nt       = 1        # number of time steps
+    ncx, ncy = 31, 31    # numerical grid resolution
+    ε        = 1e-6      # nonlinear tolerence
+    iterMax  = 1e5       # max number of iters
+    nout     = 200       # check frequency
+    # Iterative parameters -------------------------------------------
+    Reopt    = 5π
+    cfl      = 0.5
+    ρ        = 1#cfl*Reopt/ncx
+    # Preprocessing
+    Δx, Δy  = Lx/ncx, Ly/ncy
+    Δt      = η0/(G*ξ + 1e-15)
+    # Array initialisation
+    Pt       = zeros(Dat, ncx  ,ncy  )
+    ∇V       = zeros(Dat, ncx  ,ncy  )
+    Vx       = zeros(Dat, ncx+1,ncy+2)
+    Vy       = zeros(Dat, ncx+2,ncy+1)
+    ε̇xx      = zeros(Dat, ncx  ,ncy  )
+    ε̇yy      = zeros(Dat, ncx  ,ncy  )
+    ε̇xy      = zeros(Dat, ncx+1,ncy+1)    
+    τxx      = zeros(Dat, ncx  ,ncy  )
+    τyy      = zeros(Dat, ncx  ,ncy  )
+    τxy      = zeros(Dat, ncx+1,ncy+1)
+    τxx0     = zeros(Dat, ncx  ,ncy  )
+    τyy0     = zeros(Dat, ncx  ,ncy  )
+    τxy0     = zeros(Dat, ncx+1,ncy+1)
+    Rx       = zeros(Dat, ncx+1,ncy  )
+    Ry       = zeros(Dat, ncx  ,ncy+1)
+    Rp       = zeros(Dat, ncx  ,ncy  )
+    dVxdτ    = zeros(Dat, ncx+1,ncy  ) 
+    dVydτ    = zeros(Dat, ncx  ,ncy+1)
+    dPdτ     = zeros(Dat, ncx  ,ncy  )
+    Δτv      = zeros(Dat, ncx+1,ncy+1)       
+    Δτvx     = zeros(Dat, ncx+1,ncy  )         
+    Δτvy     = zeros(Dat, ncx  ,ncy+1)         
+    κΔτp     = zeros(Dat, ncx  ,ncy  )           
+    Rog      = zeros(Dat, ncx  ,ncy  )
+    ηc       = η0*ones(Dat, ncx, ncy)
+    ηv       = η0*ones(Dat, ncx+1, ncy+1)
+    # Initialisation
+    xc, yc    = LinRange( Δx/2, Lx-Δx/2, ncx+0), LinRange( Δy/2, Ly-Δy/2, ncy+0)
+    xce, yce  = LinRange(-Δx/2, Lx+Δx/2, ncx+2), LinRange(-Δy/2, Ly+Δy/2, ncy+2)
+    xv, yv    = LinRange(0.0, Lx, ncx+1), LinRange(0.0, Ly, ncy+1)
+    (Xvx,Yvx) = ([x for x=xv,y=yce], [y for x=xv,y=yce])
+    (Xvy,Yvy) = ([x for x=xce,y=yv], [y for x=xce,y=yv])
+    Vx     .=   εbg.*Xvx
+    Vy     .= .-εbg.*Yvy
+    # Time loop
+    t=0.0; evo_t=[]; evo_τxx=[]
+    for it = 1:nt
+        iter=1; err=2*ε; err_evo1=[]; err_evo2=[]; 
+        τxx0.=τxx; τyy0.=τyy; τxy0.=τxy
+        while (err>ε && iter<=iterMax)
+            # BCs
+            Vx[:,1]   .= Vx[:,2]     # S
+            Vx[:,end] .= Vx[:,end-1] # N
+            Vy[1,:]   .= Vy[2,:]     # W
+            Vy[end,:] .= Vy[end-1,:] # E
+            # Kinematics
+            ∇V    .= diff(Vx[:,2:end-1], dims=1)./Δx .+ diff(Vy[2:end-1,:], dims=2)./Δy
+            ε̇xx   .= diff(Vx[:,2:end-1], dims=1)./Δx .- 1.0/3.0*∇V
+            ε̇yy   .= diff(Vy[2:end-1,:], dims=2)./Δy .- 1.0/3.0*∇V
+            ε̇xy   .= 0.5.*(diff(Vx, dims=2)./Δy .+ diff(Vy, dims=1)./Δx)  
+            # Stresses
+            UpdateStressGeoParams!( τxx, τyy, τxy, ε̇xx, ε̇yy, ε̇xy, τxx0, τyy0, τxy0, MatParam, Δt )
+            if iter==1
+                display(τxx')
+            end
+            # Residuals
+            Rx[2:end-1,:] .= .-diff(Pt, dims=1)./Δx .+ diff(τxx, dims=1)./Δx .+ diff(τxy[2:end-1,:], dims=2)./Δy
+            Ry[:,2:end-1] .= .-diff(Pt, dims=2)./Δy .+ diff(τyy, dims=2)./Δy .+ diff(τxy[:,2:end-1], dims=1)./Δx #.+ av_ya(Rog)
+            Rp            .= .-∇V
+            # PT time step -----------------------------------------------
+            Δτv  .= ρ*min(Δx,Δy)^2 ./ ηv ./ 4.1 * cfl 
+            Δτvx .= (Δτv[:,1:end-1] .+ Δτv[:,2:end]) / 2.
+            Δτvy .= (Δτv[1:end-1,:] .+ Δτv[2:end,:]) / 2.
+            κΔτp .= cfl .* ηc .* Δx ./ Lx
+            # Calculate rate update --------------------------------------
+            dVxdτ          .= (1-ρ) .* dVxdτ .+ Rx
+            dVydτ          .= (1-ρ) .* dVydτ .+ Ry
+            dPdτ           .= Rp
+            # Update velocity and pressure -------------------------------
+            Vx[:,2:end-1]  .+= Δτvx ./ ρ .* dVxdτ
+            Vy[2:end-1,:]  .+= Δτvy ./ ρ .* dVydτ
+            Pt             .+= κΔτp .* dPdτ
+            # convergence check
+            if mod(iter, nout)==0
+                norm_Rx = norm(Rx)/sqrt(length(Rx)); norm_Ry = norm(Ry)/sqrt(length(Ry)); norm_∇V = norm(∇V)/sqrt(length(∇V))
+                err = maximum([norm_Rx, norm_Ry, norm_∇V])
+                push!(err_evo1, err); push!(err_evo2, itg)
+                @printf("it = %03d, iter = %04d, err = %1.3e norm[Rx=%1.3e, Ry=%1.3e, ∇V=%1.3e] \n", it, itg, err, norm_Rx, norm_Ry, norm_∇V)
+            end
+            iter+=1; global itg=iter
+        end
+        t = t + Δt
+        push!(evo_t, t); push!(evo_τxx, maximum(τxx))
+        # Plotting
+        p1 = heatmap(xv, yc, Vx[:,2:end-1]', aspect_ratio=1, xlims=(0, Lx), ylims=(Δy/2, Ly-Δy/2), c=:inferno, title="Vx")
+        p2 = heatmap(xc, yv, Vy[2:end-1,:]', aspect_ratio=1, xlims=(Δx/2, Lx-Δx/2), ylims=(0, Ly), c=:inferno, title="Vy")
+        p3 = heatmap(xc, yc, ε̇xx' , aspect_ratio=1, xlims=(Δx/2, Lx-Δx/2), ylims=(0, Ly), c=:inferno, title="P")
+        p4 = plot(evo_t, evo_τxx , legend=false, xlabel="time", ylabel="max(τxx)", linewiΔth=0, markershape=:circle, framestyle=:box, markersize=3)
+        p4 = plot!(evo_t, 2.0.*εbg.*η0.*(1.0.-exp.(.-evo_t.*G./η0)), linewiΔth=2.0) # analytical solution
+        display(plot(p1, p2, p3, p4))
+    end
+    return
+end
+
+Stokes2D_VE_benchmark()