Fix finite strain aniso (#119)

* fix bugs and redundancies in computation of ani_fac

* fix leak

* add new default preconditionner (-1) for anisotropy

* addres Roman's comments

* try to remove anisotropic contribution to reconditionner

* fix CI

* fix memory leak and fix matrix symmetry

* revise anisotropic shear heating

* Update InputOutput.c

* fix initialisation of aniso_fac

* add a fix for jumps in periodic anisotropic models
add a benchmark for oceanic cooling
add a double subduction setup (slab extraction)
add a simple coesite-quartz rim generation model

* add Roman setup

* add my viz

* not sure why I have to commit this file

* add new Plasticitylayers setup

* add new Plasticitylayers setup

* add layers to PlasticityLayers

* add safety in case reaction phases are not set up

JuliaVisualisation/Project.toml

@@ -3,6 +3,7 @@ CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0"
 ColorSchemes = "35d6a980-a343-548e-a6ea-1d62b119f2f4"
 Colors = "5ae59095-9a9b-59fe-a467-6f913c188581"
 DelimitedFiles = "8bb1440f-4735-579b-a4ab-409b98df4dab"
+FFMPEG = "c87230d0-a227-11e9-1b43-d7ebe4e7570a"
 FileIO = "5789e2e9-d7fb-5bc7-8068-2c6fae9b9549"
 GLMakie = "e9467ef8-e4e7-5192-8a1a-b1aee30e663a"
 GeometryTypes = "4d00f742-c7ba-57c2-abde-4428a4b178cb"
@@ -17,4 +18,5 @@ Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
 Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
 Revise = "295af30f-e4ad-537b-8983-00126c2a3abe"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
+SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"

JuliaVisualisation/_InitialConditionsStandAloneScripts/IsostasyMDoodz.jl

@@ -0,0 +1,91 @@
+import Pkg
+Pkg.activate(normpath(joinpath(@__DIR__, ".")))
+using HDF5, GLMakie, Printf#, Colors, ColorSchemes, MathTeXEngine, LinearAlgebra, FFMPEG
+Makie.update_theme!(fonts = (regular = texfont(), bold = texfont(:bold), italic = texfont(:italic)))
+
+y    = 365*24*3600
+My   = 1e6*y
+cm_y = y*100.
+
+function main()
+    x  = (min=-1e6, max=1e6)
+    z  = (min=-6e5, max=0)
+    L  = (x=x.max-x.min, z=z.max-z.min)
+    nc = (x=500, z=300)
+    Δ  = (x=L.x/nc.x, z=L.z/nc.z)
+    xc = LinRange(x.min, x.max, nc.x) 
+    zc = LinRange(z.min, z.max, nc.z) 
+    ρc = 3300.0.*ones(nc...)
+    # ρc[(0.0*xc .+ zc') .> -1.5e5] .= 2700
+    hc = 400.0.*ones(nc.x)
+    # ρc[abs.(xc .+ 0.0*zc') .< 2e5 .&& (0.0*xc .+ zc') .> -1.5e5] .= 2800
+
+    for j=1:nc.z, i=1:nc.x
+        # Introduce a circle
+        z0 = -250e3 
+        r  = 25e3
+        x  = xc[i]
+        z  = zc[j]
+        if (x-0.0)^2 + (z-z0)^2 < r^2
+            ρc[i,j] = 3450.0 
+        end
+
+        # Introduce a margin
+        Lm  = 100e3 # margin center
+        L0  = -2e5   # margin position
+        HW  = 1e5   # lithosphere thickness
+        HE  = 2e4   # lithosphere thickness
+        # 2 points  
+        p1  = (x=L0-Lm/2, z=-HW)
+        p2  = (x=L0+Lm/2, z=-HE)
+        a   = (p2.z - p1.z) / (p2.x - p1.x)
+        b   = p2.z - a*p2.x
+        zm  = a*x +b
+        if (abs(x-L0)<Lm && z>zm && z<-HE)
+            ρc[i,j] = 2700.0
+        end
+
+        # Introduce a margin
+        Lm  = 100e3 # margin center
+        L0  = 2e5   # margin position
+        HW  = 1e5   # lithosphere thickness
+        HE  = 2e4   # lithosphere thickness
+        # 2 points  
+        p1  = (x=L0-Lm/2, z=-HE)
+        p2  = (x=L0+Lm/2, z=-HW)
+        a   = (p2.z - p1.z) / (p2.x - p1.x)
+        b   = p2.z - a*p2.x
+        zm  = a*x +b
+        if (abs(x-L0)<Lm && z>zm && z<-HE)
+            ρc[i,j] = 2700.0
+        end
+
+
+    end
+
+    # Compute west column load in units of P/g
+    for i=2:nc.x
+        PW     = sum(ρc[1,:]*Δ.z)
+        PE     = sum(ρc[i,:]*Δ.z)
+        Δh     = -(PE - PW) / ρc[i,1]
+        hc[i] += Δh
+    end
+
+    ##################################################################
+    resol       = 1000
+    f = Figure(resolution = (L.x/L.z*resol, resol), fontsize=25)
+
+    ax1 = Axis(f[1, 1], title = L"$h$", xlabel = L"$x$ [km]", ylabel = L"$h$ [m]")
+    lines!(ax1, xc/1e3, hc)
+
+    ax2 = Axis(f[2, 1], title = L"$ρ$", xlabel = L"$x$ [km]", ylabel = L"$y$ [km]")
+    hm = heatmap!(ax2, xc/1e3, zc/1e3, ρc, colormap = (:turbo, 1.0))           
+    colsize!(f.layout, 1, Aspect(1, L.x/L.z))
+    GLMakie.Colorbar(f[2, 2], hm, label = L"$ρ$ [kg.m$^3$]", width = 20, labelsize = 25, ticklabelsize = 14 )
+    GLMakie.colgap!(f.layout, 20)
+    display(f)
+    DataInspector(f)
+
+end
+
+main()