From 2bcb3e5e98e2f99a30bff69cf442bd936a44a36b Mon Sep 17 00:00:00 2001 From: tduretz Date: Sat, 7 Dec 2024 13:46:07 +0100 Subject: [PATCH] update visualisation and stash Cosserat --- .../Main_MatrixCheckSolve_MD7.jl | 7 +- .../ChristmasTreeAnisotropy.jl | 21 +- .../RiftingAnisotropy/ResolutionTest.jl | 346 ++---------- .../RiftingAnisotropy/ResolutionTestStress.jl | 515 +++++++++++++++++ .../RiftingAnisotropy/TimeEvolution.jl | 516 ++++++++++++++++++ .../ViscousAnisotropySimpleV.jl | 45 +- .../ViscousAnisotropySimpleVEP.jl | 4 +- 7 files changed, 1133 insertions(+), 321 deletions(-) create mode 100644 JuliaVisualisation/_SpecificStudies/RiftingAnisotropy/ResolutionTestStress.jl create mode 100644 JuliaVisualisation/_SpecificStudies/RiftingAnisotropy/TimeEvolution.jl diff --git a/JuliaVisualisation/Main_MatrixCheckSolve_MD7.jl b/JuliaVisualisation/Main_MatrixCheckSolve_MD7.jl index de09c7d2..b03df53e 100644 --- a/JuliaVisualisation/Main_MatrixCheckSolve_MD7.jl +++ b/JuliaVisualisation/Main_MatrixCheckSolve_MD7.jl @@ -42,7 +42,7 @@ function main() # File path ="/Users/tduretz/REPO/MDOODZ7.0/MDLIB/" - filename = string(path, "Stokes_01cpu_step1429_iter00.gzip.h5") + filename = string(path, "Stokes_01cpu_step01_iter00.gzip.h5") # path ="/Users/tduretz/REPO/MDOODZ7.0/RUNS/RiverTom/" # filename = string(path, "Stokes_01cpu_step229_iter00.gzip.h5") # Matrix block A @@ -75,6 +75,11 @@ function main() model = ExtractData(filename,"/model/params"); nx, nz, Δx, Δz = model + @show size(MA) + @show size(MB) + @show size(MC) + @show size(MD) + # Spies Msc = MA - MB*(MD*MC); @show norm(MA-MA') diff --git a/JuliaVisualisation/_SpecificStudies/RiftingAnisotropy/ChristmasTreeAnisotropy.jl b/JuliaVisualisation/_SpecificStudies/RiftingAnisotropy/ChristmasTreeAnisotropy.jl index 4a4b616d..96223161 100644 --- a/JuliaVisualisation/_SpecificStudies/RiftingAnisotropy/ChristmasTreeAnisotropy.jl +++ b/JuliaVisualisation/_SpecificStudies/RiftingAnisotropy/ChristmasTreeAnisotropy.jl @@ -49,10 +49,10 @@ import Statistics:mean # Plasticity C = 5e7 fric = 30*π/180 - dil = 10*π/180 - ηvp = 0.0 + dil = 0*π/180 + ηvp = 0*1e21 # Pressure in extension - X = 1 .+ (δ .- 1).*sind.(2*θ) + X = 1 / (sqrt((δ^2-1)*cos(2*θ)^2 + 1)/δ) ϕeff = asin(sin(fric) / X) @show ϕeff*180/π σH = -Pl*(1-sin(ϕeff))/(1+sin(ϕeff)) @@ -187,15 +187,10 @@ import Statistics:mean 2.0*ani*d1 2.0-2.0*ani*d1 2.0*ani*d2; -2.0*ani*d2 2.0*ani*d2 1.0 + 2.0*ani*(d1 - 0.5);] - - E = [ε̇xxd; ε̇yyd; ε̇xyd] τ0 = [τxx0/ηe; τyy0/ηe; τxy0/ηe] - Eeff = E + A\τ0 - - @show Eeff - - + b = A\τ0 + Eeff = E + b./[1; 1; 2] Da11 = 2.0 - 2.0*ani*d1; Da12 = 2.0*ani*d1; @@ -252,8 +247,8 @@ import Statistics:mean ##################### - τii_yield = P*sin(fric) .+ C#*cos(fric) - τii_yield_min = τii_yield * (1 / (1 + (δ-1)*sin(2*θ[i]))) + τii_yield = P*sin(fric) .+ C*cos(fric) + τii_yield_min = τii_yield*sqrt((δ^2-1)*cos(2*θ[i])^2 + 1)/δ ax3 = Axis(f[1, 2], title = L"$$Stress profile ($\delta = 3$)", xlabel = L"$τ_\mathrm{II}$ [GPa]") lines!(ax3, τii_rot1./1e9, yc./1e3, label=L"$\sqrt{Y_2^{'}}$", color=:blue ) @@ -268,6 +263,8 @@ import Statistics:mean axislegend(position=:rb) display(f) + save("/Users/tduretz/PowerFolders/_manuscripts/RiftingAnisotropy/Figures/ChristmasTreeAnisotropy.png", f, px_per_unit = 4) + @show τii_rot1./τii_cart1 end diff --git a/JuliaVisualisation/_SpecificStudies/RiftingAnisotropy/ResolutionTest.jl b/JuliaVisualisation/_SpecificStudies/RiftingAnisotropy/ResolutionTest.jl index 3bfb002c..5ccf4250 100644 --- a/JuliaVisualisation/_SpecificStudies/RiftingAnisotropy/ResolutionTest.jl +++ b/JuliaVisualisation/_SpecificStudies/RiftingAnisotropy/ResolutionTest.jl @@ -20,6 +20,11 @@ function ExtractField(filename, field, size, mask_air, mask) end function AddCountourQuivers!(PlotOnTop, ax1, coords, V, T, ϕ, σ1, ε̇1, Fab, height, Lc, cm_y, group_phases, Δ) + if PlotOnTop.quiver_origin + xc2D = coords.c.x .+ 0*coords.c.z' + zc2D = 0*coords.c.x .+ coords.c.z' + scatter!(ax1, xc2D[1:V.step:end,1:V.step:end][:]./Lc, zc2D[1:V.step:end,1:V.step:end][:]./Lc) + end if PlotOnTop.T_contours contour!(ax1, coords.c.x./Lc, coords.c.z./Lc, T, levels=0:200:1400, linewidth = 4, color=:black ) end @@ -157,7 +162,7 @@ function ReadFile(path, step, scales, options, PlotOnTop) end - return (tMy=tMy,length_unit=length_unit, Lx=Lx, Lz=Lz, xc=xc, zc=zc, ε̇II=ε̇II, + return (tMy=tMy,length_unit=length_unit, Lx=Lx, Lz=Lz, xc=xc, zc=zc, ε̇II=ε̇II, τII=τII, coords=coords, V=V, T=T, ϕ=ϕ, σ1=σ1, ε̇1=ε̇1, Fab=Fab, height=height, group_phases=group_phases, Δ=Δ) end @@ -171,8 +176,8 @@ end # File numbers - step = [120; 200; 500] - # step = [500; 1000; 2000] + step = [120; 160; 360] + step = [500; 1100; 2220] # Select field to visualise # field = :Phases # field = :Cohesion @@ -196,13 +201,13 @@ end # field = :EffectiveFrictionTime # field = :ChristmasTree - # Define Tuple for enlargment window - zoom = ( - xmin = 75, - xmax = 200, - zmin = -30, - zmax = 0.e3, - ) + # # Define Tuple for enlargment window + # zoom = ( + # xmin = 125, + # xmax = 150, + # zmin = -15, + # zmax = 0.e3, + # ) # Switches PlotOnTop = ( @@ -210,13 +215,14 @@ end fabric = false, # add fabric quiver (normal to director) T_contours = false, # add temperature contours topo = false, - σ1_axis = true, - ε̇1_axis = true, + quiver_origin = false, + σ1_axis = false, + ε̇1_axis = false, vel_vec = false, - ϕ_contours = true, + ϕ_contours = false, ) options = ( - printfig = false, # print figures to disk + printfig = true, # print figures to disk printvid = false, framerate = 6, α_heatmap = 0.5, # transparency of heatmap @@ -264,7 +270,7 @@ end @unpack Lc, tc, Vc, τc = scales model = ReadFile(path_LR, step[1], scales, options, PlotOnTop) - @unpack tMy, length_unit, Lx, Lz, xc, zc, ε̇II, coords, V, T, ϕ, σ1, Fab, height, group_phases, Δ = model + @unpack tMy, length_unit, Lx, Lz, xc, zc, ε̇II, τII, coords, V, T, ϕ, σ1, Fab, height, group_phases, Δ = model if @isdefined zoom Lx = (zoom.xmax - zoom.xmin)./scales.Lc @@ -283,80 +289,6 @@ end ftsz = 12*options.resol/500 f = Figure(size = (1.1*Lx/4/Lz*options.resol*1.2, options.resol), fontsize=ftsz) - # if field==:Phases - # ax1 = Axis(f[1, 1], title = L"Phases at $t$ = %$(tMy) Ma", xlabel = L"$x$ [m]", ylabel = L"$y$ [m]") - # hm = heatmap!(ax1, xc_hr./Lc, zc_hr./Lc, ph_hr, colormap = phase_colors) - # # hm = heatmap!(ax1, xc_hr./Lc, zc_hr./Lc, ph_hr, colormap = :turbo) - # hm = heatmap!(ax1, xc_hr./Lc, zc_hr./Lc, ph_dual_hr, colormap = :turbo) - # AddCountourQuivers!(PlotOnTop, ax1, coords, V, T, ϕ, σ1, Fab, height, Lc, cm_y, group_phases, Δ) - # colsize!(f.layout, 1, Aspect(1, Lx/Lz)) - # Mak.Colorbar(f[1, 2], hm, label = "Phases", width = 20, labelsize = ftsz, ticklabelsize = ftsz ) - # Mak.colgap!(f.layout, 20) - # xlims!(ax1, window.xmin, window.xmax) - # ylims!(ax1, window.zmin, window.zmax) - # if printfig Print2Disk( f, path, string(field), istep) end - # end - - # if field==:Viscosity - # ax1 = Axis(f[1, 1], title = L"$\eta$ at $t$ = %$(tMy) Ma", xlabel = L"$x$ [m]", ylabel = L"$y$ [m]") - # hm = heatmap!(ax1, xc./Lc, zc./Lc, log10.(ηc), colormap = (:turbo, α_heatmap)) - # AddCountourQuivers!(PlotOnTop, ax1, coords, V, T, ϕ, σ1, Fab, height, Lc, cm_y, group_phases, Δ) - # colsize!(f.layout, 1, Aspect(1, Lx/Lz)) - # Mak.Colorbar(f[1, 2], hm, label = L"$\eta$ [Pa.s]", width = 20, labelsize = ftsz, ticklabelsize = ftsz ) - # Mak.colgap!(f.layout, 20) - # xlims!(ax1, window.xmin, window.xmax) - # ylims!(ax1, window.zmin, window.zmax) - # if printfig Print2Disk( f, path, string(field), istep) end - # end - - # if field==:Density - # ax1 = Axis(f[1, 1], title = L"$\rho$ at $t$ = %$(tMy) Ma", xlabel = L"$x$ [km]", ylabel = L"$y$ [km]") - # hm = heatmap!(ax1, xc./Lc, zc./Lc, ρc, colormap = (:turbo, α_heatmap), colorrange=(2600, 3000)) - # AddCountourQuivers!(PlotOnTop, ax1, coords, V, T, ϕ, σ1, Fab, height, Lc, cm_y, group_phases, Δ) - # colsize!(f.layout, 1, Aspect(1, Lx/Lz)) - # Mak.Colorbar(f[1, 2], hm, label = L"$\rho$ [kg.m$^{-3}$]", width = 20, labelsize = ftsz, ticklabelsize = ftsz ) - # Mak.colgap!(f.layout, 20) - # xlims!(ax1, window.xmin, window.xmax) - # ylims!(ax1, window.zmin, window.zmax) - # if printfig Print2Disk( f, path, string(field), istep) end - # end - - # if field==:Stress - # ax1 = Axis(f[1, 1], title = L"$\tau_\textrm{II}$ at $t$ = %$(tMy) Ma", xlabel = L"$x$ [km]", ylabel = L"$y$ [km]") - # hm = heatmap!(ax1, xc./Lc, zc./Lc, τII./1e6, colormap = (:turbo, α_heatmap)) - # AddCountourQuivers!(PlotOnTop, ax1, coords, V, T, ϕ, σ1, Fab, height, Lc, cm_y, group_phases, Δ) - # colsize!(f.layout, 1, Aspect(1, Lx/Lz)) - # Mak.Colorbar(f[1, 2], hm, label = L"$\tau_\textrm{II}$ [MPa]", width = 20, labelsize = ftsz, ticklabelsize = ftsz ) - # Mak.colgap!(f.layout, 20) - # xlims!(ax1, window.xmin, window.xmax) - # ylims!(ax1, window.zmin, window.zmax) - # if printfig Print2Disk( f, path, string(field), istep) end - # end - - # if field==:Pressure - # ax1 = Axis(f[1, 1], title = L"$P$ at $t$ = %$(tMy) Ma", xlabel = L"$x$ [km]", ylabel = L"$y$ [km]") - # hm = heatmap!(ax1, xc./Lc, zc./Lc, P./1e9, colormap = (:turbo, α_heatmap)) #, colorrange=(1,1.2)1e4*365*24*3600 - # AddCountourQuivers!(PlotOnTop, ax1, coords, V, T, ϕ, σ1, Fab, height, Lc, cm_y, group_phases, Δ) - # colsize!(f.layout, 1, Aspect(1, Lx/Lz)) - # Mak.Colorbar(f[1, 2], hm, label = L"$P$ [GPa]", width = 20, labelsize = ftsz, ticklabelsize = ftsz ) - # Mak.colgap!(f.layout, 20) - # xlims!(ax1, window.xmin, window.xmax) - # ylims!(ax1, window.zmin, window.zmax) - # if printfig Print2Disk( f, path, string(field), istep) end - # end - - # if field==:Divergence - # ax1 = Axis(f[1, 1], title = L"∇⋅V at $t$ = %$(tMy) Ma", xlabel = L"$x$ [km]", ylabel = L"$y$ [km]") - # hm = heatmap!(ax1, xc./Lc, zc./Lc, divu, colormap = (:turbo, α_heatmap)) - # AddCountourQuivers!(PlotOnTop, ax1, coords, V, T, ϕ, σ1, Fab, height, Lc, cm_y, group_phases, Δ) - # colsize!(f.layout, 1, Aspect(1, Lx/Lz)) - # Mak.Colorbar(f[1, 2], hm, label = L"∇⋅V [s$^{-1}$]", width = 20, labelsize = ftsz, ticklabelsize = ftsz ) - # Mak.colgap!(f.layout, 20) - # xlims!(ax1, window.xmin, window.xmax) - # ylims!(ax1, window.zmin, window.zmax) - # if printfig Print2Disk( f, path, string(field), istep) end - # end - if field==:StrainRate options = ( @@ -366,7 +298,7 @@ end α_heatmap = 0.5, # transparency of heatmap vel_arrow = 5, vel_scale = 10000, - vel_step = 5, + vel_step = 2, nap = 0.1, # pause for animation resol = 1000, Lx = 1.0, @@ -376,7 +308,7 @@ end ) model = ReadFile(path_LR, step[1], scales, options, PlotOnTop) - @unpack tMy, length_unit, Lx, Lz, xc, zc, ε̇II, coords, V, T, ϕ, σ1, ε̇1, Fab, height, group_phases, Δ = model + @unpack tMy, length_unit, Lx, Lz, xc, zc, ε̇II, τII, coords, V, T, ϕ, σ1, ε̇1, Fab, height, group_phases, Δ = model ax1 = Axis(f[1, 1], title = L"$\dot{\varepsilon}_\textrm{II}$ at $t$ = %$(tMy) Ma - 300 $\times$ 200 cells", ylabel = L"$y$ [%$(length_unit)]") hm = heatmap!(ax1, xc./Lc, zc./Lc, log10.(ε̇II), colormap = (:turbo, options.α_heatmap), colorrange=(-17, -13)) @@ -386,6 +318,11 @@ end Mak.colgap!(f.layout, 20) xlims!(ax1, window.xmin, window.xmax) ylims!(ax1, window.zmin, window.zmax) + + # lines!(ax1, xc./Lc, log10.(ε̇II[:, 170]) ) + # lines!(ax1, xc./Lc, atand.(ε̇1.z[:, 170]) ./ ε̇1.x[:, 170] ) + # xlims!(ax1, window.xmin, window.xmax) + # ylims!(ax1, -180, 180) ############################################################## @@ -397,7 +334,7 @@ end α_heatmap = 0.5, # transparency of heatmap vel_arrow = 5, vel_scale = 10000, - vel_step = 10, + vel_step = 4, nap = 0.1, # pause for animation resol = 1000, Lx = 1.0, @@ -407,32 +344,34 @@ end ) model = ReadFile(path_MR, step[2], scales, options, PlotOnTop) - @unpack tMy, length_unit, Lx, Lz, xc, zc, ε̇II, coords, V, T, ϕ, σ1, ε̇1, Fab, height, group_phases, Δ = model + @unpack tMy, length_unit, Lx, Lz, xc, zc, ε̇II, τII, coords, V, T, ϕ, σ1, ε̇1, Fab, height, group_phases, Δ = model ax1 = Axis(f[2, 1], title = L"$\dot{\varepsilon}_\textrm{II}$ at $t$ = %$(tMy) Ma - 600 $\times$ 400 cells", ylabel = L"$y$ [%$(length_unit)]") - # hm = heatmap!(ax1, xc./Lc, zc./Lc, log10.(ε̇II), colormap = (:turbo, options.α_heatmap), colorrange=(-17, -13)) - # AddCountourQuivers!(PlotOnTop, ax1, coords, V, T, ϕ, σ1, ε̇1, Fab, height, Lc, cm_y, group_phases, Δ) - # colsize!(f.layout, 1, Aspect(1, Lx/Lz)) - # Mak.Colorbar(f[2, 2], hm, label = L"$\dot{\varepsilon}_\textrm{II}$ [s$^{-1}$]", width = 20, labelsize = ftsz, ticklabelsize = ftsz ) - # Mak.colgap!(f.layout, 20) - # xlims!(ax1, window.xmin, window.xmax) - # ylims!(ax1, window.zmin, window.zmax) - # lines!(ax1, xc./Lc, log10.(ε̇II[:, 350]) ) - lines!(ax1, xc./Lc, atand.(σ1.z[:, 350]) ./ σ1.x[:, 350] ) - + hm = heatmap!(ax1, xc./Lc, zc./Lc, log10.(ε̇II), colormap = (:turbo, options.α_heatmap), colorrange=(-17, -13)) + AddCountourQuivers!(PlotOnTop, ax1, coords, V, T, ϕ, σ1, ε̇1, Fab, height, Lc, cm_y, group_phases, Δ) + colsize!(f.layout, 1, Aspect(1, Lx/Lz)) + Mak.Colorbar(f[2, 2], hm, label = L"$\dot{\varepsilon}_\textrm{II}$ [s$^{-1}$]", width = 20, labelsize = ftsz, ticklabelsize = ftsz ) + Mak.colgap!(f.layout, 20) xlims!(ax1, window.xmin, window.xmax) + ylims!(ax1, window.zmin, window.zmax) + + + # lines!(ax1, xc./Lc, log10.(ε̇II[:, 350]) ) + # lines!(ax1, xc./Lc, atand.(ε̇1.z[:, 350]) ./ ε̇1.x[:, 350] ) + # xlims!(ax1, window.xmin, window.xmax) + # ylims!(ax1, -180, 180) # ############################################################## options = ( - printfig = false, # print figures to disk + printfig = true, # print figures to disk printvid = false, framerate = 6, α_heatmap = 0.5, # transparency of heatmap vel_arrow = 5, vel_scale = 10000, - vel_step = 20, + vel_step = 8, nap = 0.1, # pause for animation resol = 1000, Lx = 1.0, @@ -453,194 +392,18 @@ end xlims!(ax1, window.xmin, window.xmax) ylims!(ax1, window.zmin, window.zmax) - ############################################################## + # lines!(ax1, xc./Lc, log10.(ε̇II[:, 700]) ) + # lines!(ax1, xc./Lc, atand.(ε̇1.z[:, 700]) ./ ε̇1.x[:, 700] ) + # xlims!(ax1, window.xmin, window.xmax) + # ylims!(ax1, -180, 180) + - if options.printfig Print2Disk( f, path, string(field), istep) end #if printfig Print2Disk( f, path, string(field),ε̇BG) end + ############################################################## + + if options.printfig Print2Disk( f, path_HR, string(field), step[3]) end #if printfig Print2Disk( f, path, string(field),ε̇BG) end end - # if field==:PlasticStrainrate - # ε̇pl[ε̇pl.==0.0] .= 1e-30 - # ax1 = Axis(f[1, 1], title = L"$\dot{\varepsilon}_\textrm{II}^\textrm{pl}$ at $t$ = %$(tMy) Ma", xlabel = L"$x$ [km]", ylabel = L"$y$ [km]") - # hm = heatmap!(ax1, xc./Lc, zc./Lc, log10.(ε̇pl), colormap = (:turbo, α_heatmap)) - # AddCountourQuivers!(PlotOnTop, ax1, coords, V, T, ϕ, σ1, Fab, height, Lc, cm_y, group_phases, Δ) - # colsize!(f.layout, 1, Aspect(1, Lx/Lz)) - # Mak.Colorbar(f[1, 2], hm, label = L"$\dot{\varepsilon}_\textrm{II}^\textrm{pl}$ [s$^{-1}$]", width = 20, labelsize = ftsz, ticklabelsize = ftsz ) - # Mak.colgap!(f.layout, 20) - # xlims!(ax1, window.xmin, window.xmax) - # ylims!(ax1, window.zmin, window.zmax) - # if printfig Print2Disk( f, path, string(field), istep) end - # end - - # if field==:Velocity - # ax1 = Axis(f[1, 1], title = L"$V$ at $t$ = %$(tMy) Ma", xlabel = L"$x$ [km]", ylabel = L"$y$ [km]") - # Vx_BG = 0*xc .- 2*zc' - # V = sqrt.( (Vxc .- 0.0*Vx_BG).^2 + (Vzc).^2) - # hm = heatmap!(ax1, xc./Lc, zc./Lc, V, colormap = (:jet, α_heatmap))#, colorrange=(0., 0.6) - # AddCountourQuivers!(PlotOnTop, ax1, coords, V, T, ϕ, σ1, Fab, height, Lc, cm_y, group_phases, Δ) - # # xlims!(ax1, 0., 3.e-3) - # # ylims!(ax1, 0., 3.e-3) - # colsize!(f.layout, 1, Aspect(1, Lx/Lz)) - # Mak.Colorbar(f[1, 2], hm, label = L"$V$ [m.s$^{-1}$]", width = 20, labelsize = ftsz, ticklabelsize = ftsz ) - # Mak.colgap!(f.layout, 20) - # xlims!(ax1, window.xmin, window.xmax) - # ylims!(ax1, window.zmin, window.zmax) - # if printfig Print2Disk( f, path, string(field), istep) end - # end - - # if field==:Velocity_x - # ax1 = Axis(f[1, 1], title = L"$Vx$ at $t$ = %$(tMy) Ma", xlabel = L"$x$ [km]", ylabel = L"$y$ [km]") - # hm = heatmap!(ax1, xv./Lc, zvx[2:end-1]./Lc, Vx[2:end-1,:]*cm_yr, colormap = (:jet, α_heatmap))#, colorrange=(0., 0.6) - # AddCountourQuivers!(PlotOnTop, ax1, coords, V, T, ϕ, σ1, Fab, height, Lc, cm_y, group_phases, Δ) - # colsize!(f.layout, 1, Aspect(1, Lx/Lz)) - # Mak.Colorbar(f[1, 2], hm, label = L"$Vx$ [cm.yr$^{-1}$]", width = 20, labelsize = ftsz, ticklabelsize = ftsz ) - # Mak.colgap!(f.layout, 20) - # xlims!(ax1, window.xmin, window.xmax) - # ylims!(ax1, window.zmin, window.zmax) - # if printfig Print2Disk( f, path, string(field), istep) end - # end - - # if field==:Velocity_z - # ax1 = Axis(f[1, 1], title = L"$Vz$ at $t$ = %$(tMy) Ma", xlabel = L"$x$ [km]", ylabel = L"$y$ [km]") - # hm = heatmap!(ax1, xvz[2:end-1]./Lc, zv./Lc, Vz[:,2:end-1]*cm_yr, colormap = (:jet, α_heatmap))#, colorrange=(0., 0.6) - # AddCountourQuivers!(PlotOnTop, ax1, coords, V, T, ϕ, σ1, Fab, height, Lc, cm_y, group_phases, Δ) - # colsize!(f.layout, 1, Aspect(1, Lx/Lz)) - # Mak.Colorbar(f[1, 2], hm, label = L"$Vz$ [cm.yr$^{-1}$]", width = 20, labelsize = ftsz, ticklabelsize = ftsz ) - # Mak.colgap!(f.layout, 20) - # xlims!(ax1, window.xmin, window.xmax) - # ylims!(ax1, window.zmin, window.zmax) - # if printfig Print2Disk( f, path, string(field), istep) end - # end - - # if field==:GrainSize - # ax1 = Axis(f[1, 1], title = L"$d$ at $t$ = %$(tMy) Ma", xlabel = L"$x$ [km]", ylabel = L"$y$ [km]") - # hm = heatmap!(ax1, xc./Lc, zc./Lc, log10.(d.*1e6), colormap = (:turbo, α_heatmap), colorrange=(1, 3)) - # AddCountourQuivers!(PlotOnTop, ax1, coords, V, T, ϕ, σ1, Fab, height, Lc, cm_y, group_phases, Δ) - # xminz, xmaxz = -0.4, 0.4 - # zminz, zmaxz = -0.17, 0.17 - # Lx = xmaxz - xminz - # Lz = zmaxz - zminz - # xlims!(ax1, -0.4, 0.4) - # ylims!(ax1, -0.17, 0.17) - # colsize!(f.layout, 1, Aspect(1, Lx/Lz)) - # Mak.Colorbar(f[1, 2], hm, label = "d", width = 20, labelsize = ftsz, ticklabelsize = ftsz ) - # Mak.colgap!(f.layout, 20) - # xlims!(ax1, window.xmin, window.xmax) - # ylims!(ax1, window.zmin, window.zmax) - # if printfig Print2Disk( f, path, string(field), istep) end - # end - - # if field==:AnisotropyFactor - # ax1 = Axis(f[1, 1], title = L"$δ_\textrm{ani}$ at $t$ = %$(tMy) Ma", xlabel = L"$x$ [km]", ylabel = L"$y$ [km]") - # hm = heatmap!(ax1, xc./Lc, zc./Lc, δani, colormap = (:bilbao, α_heatmap)) - # AddCountourQuivers!(PlotOnTop, ax1, coords, V, T, ϕ, σ1, Fab, height, Lc, cm_y, group_phases, Δ) - # # xminz, xmaxz = -0.4, 0.4 - # # zminz, zmaxz = -0.17, 0.17 - # # Lx = xmaxz - xminz - # # Lz = zmaxz - zminz - # # xlims!(ax1, -0.4, 0.4) - # # ylims!(ax1, -0.17, 0.17) - # Mak.Colorbar(f[1, 2], hm, label = L"$δ_\textrm{ani}$", width = 20, labelsize = ftsz, ticklabelsize = ftsz ) - # Mak.colgap!(f.layout, 20) - # colsize!(f.layout, 1, Aspect(1, Lx/Lz)) - # xlims!(ax1, window.xmin, window.xmax) - # ylims!(ax1, window.zmin, window.zmax) - # if printfig Print2Disk( f, path, string(field), istep) end - # end - - # if field==:MeltFraction - # ax1 = Axis(f[1, 1], title = L"ϕ at $t$ = %$(tMy) Ma", xlabel = L"$x$ [km]", ylabel = L"$y$ [km]") - # hm = heatmap!(ax1, xc./Lc, zc./Lc, ϕ, colormap = (:bilbao, α_heatmap)) - # AddCountourQuivers!(PlotOnTop, ax1, coords, V, T, ϕ, σ1, Fab, height, Lc, cm_y, group_phases, Δ) - # Mak.Colorbar(f[1, 2], hm, label = L"$ϕ$", width = 20, labelsize = ftsz, ticklabelsize = ftsz ) - # Mak.colgap!(f.layout, 20) - # colsize!(f.layout, 1, Aspect(1, Lx/Lz)) - # xlims!(ax1, window.xmin, window.xmax) - # ylims!(ax1, window.zmin, window.zmax) - # if printfig Print2Disk( f, path, string(field), istep) end - # end - - # if field==:Cohesion - # ax1 = Axis(f[1, 1], title = L"$C$ [MPa] at $t$ = %$(tMy) Ma", xlabel = L"$x$ [km]", ylabel = L"$y$ [km]") - # hm = heatmap!(ax1, xc./Lc, zc./Lc, C./1e6, colormap = (Reverse(:bilbao), α_heatmap)) - # AddCountourQuivers!(PlotOnTop, ax1, coords, V, T, ϕ, σ1, Fab, height, Lc, cm_y, group_phases, Δ) - # Mak.Colorbar(f[1, 2], hm, label = L"$C$ [MPa]", width = 20, labelsize = ftsz, ticklabelsize = ftsz ) - # Mak.colgap!(f.layout, 20) - # colsize!(f.layout, 1, Aspect(1, Lx/Lz)) - # xlims!(ax1, window.xmin, window.xmax) - # ylims!(ax1, window.zmin, window.zmax) - # if printfig Print2Disk( f, path, string(field), istep) end - # end - - # if field==:Temperature - # ax1 = Axis(f[1, 1], title = L"$T$ [C] at $t$ = %$(tMy) Ma", xlabel = L"$x$ [km]", ylabel = L"$y$ [km]") - # hm = heatmap!(ax1, xc./Lc, zc./Lc, T, colormap = Reverse(:bilbao)) - # AddCountourQuivers!(PlotOnTop, ax1, coords, V, T, ϕ, σ1, Fab, height, Lc, cm_y, group_phases, Δ) - # Mak.Colorbar(f[1, 2], hm, label = L"$T$ [C]", width = 20, labelsize = ftsz, ticklabelsize = ftsz ) - # Mak.colgap!(f.layout, 20) - # colsize!(f.layout, 1, Aspect(1, Lx/Lz)) - # xlims!(ax1, window.xmin, window.xmax) - # ylims!(ax1, window.zmin, window.zmax) - # if printfig Print2Disk( f, path, string(field), istep) end - # end - - # if field==:Topography - # height = Float64.(ExtractData( filename, "/Topo/z_grid")); - # Vx_grid = Float64.(ExtractData( filename, "/Topo/Vx_grid")); - # Vz_grid = Float64.(ExtractData( filename, "/Topo/Vz_grid")); - # Vx_mark = Float64.(ExtractData( filename, "/Topo/Vx_mark")); - # Vz_mark = Float64.(ExtractData( filename, "/Topo/Vz_mark")); - # x_mark = Float64.(ExtractData( filename, "/Topo/x_mark")); - # z_mark = Float64.(ExtractData( filename, "/Topo/z_mark")); - - # ax1 = Axis(f[1, 1], title = L"Topography at $t$ = %$(tMy) Ma", xlabel = L"$x$ [km]", ylabel = L"$h$ [km]") - # @show mean(height) - # @show mean(z_mark) - # lines!(ax1, xv./Lc, height./Lc) - # scatter!(ax1, x_mark./Lc, z_mark./Lc) - # xlims!(ax1, window.xmin, window.xmax) - - # ax2 = Axis(f[2, 1], xlabel = L"$x$ [km]", ylabel = L"$Vx$ [km]") - # lines!(ax2, xv./Lc, Vx_grid./Vc) - # scatter!(ax2, x_mark./Lc, Vx_mark./Vc) - # xlims!(ax1, window.xmin, window.xmax) - - # ax3 = Axis(f[3, 1], xlabel = L"$x$ [km]", ylabel = L"$Vz$ [km]") - # lines!(ax3, xc./Lc, Vz_grid[2:end-1]./Vc) - # scatter!(ax3, x_mark/Lc, Vz_mark./Vc) - # xlims!(ax1, window.xmin, window.xmax) - # end - - # if field==:TimeSeries - # ax1 = Axis(f[1, 1], title = L"$τ_{xz}$", xlabel = L"$t$", ylabel = L"$\tau_{xz}$") - # τxz_t[1] = 0. - # @show .-τxz_t - # lines!(ax1, t_t, .-τxz_t./(.-P_t.+τzz_t)) - # end - - # if field==:EffectiveFrictionTime - # ϕ_eff = -mean(τxzc[:,end] ./ (τzz[:,end] .- P[:,end]) ) - # if istep==0 ϕ_eff = 0. end - # push!(probe.ϕeff, ϕ_eff) - # push!(probe.t, t) - # if istep==file_end - # ax1 = Axis(f[1, 1], title = L"$ϕ_\mathrm{eff}$", xlabel = L"$t$", ylabel = L"$ϕ_\mathrm{eff}$") - # lines!(ax1, Float64.(probe.t), Float64.(probe.ϕeff)) - # end - # end - - # if field==:ChristmasTree - # ax1 = Axis(f[1, 1], title = L"Stress profile at $t$ = %$(tMy) Ma", xlabel = L"$τII$ [MPa]", ylabel = L"$z$ [km]") - # lines!(ax1, mean(τII, dims=1)[:]/τc, coords.c.z./Lc/1e3, ) - # lines!(ax1, mean(T, dims=1)[:], coords.c.z./Lc/1e3, ) - # ax2 = Axis(f[1, 2], title = L"Temperature profile at $t$ = %$(tMy) Ma", xlabel = L"$T$ [C]", ylabel = L"$h$ [km]") - # lines!(ax2, mean(T, dims=1)[:], coords.c.z./Lc/1e3, ) - # end - - # if field!=:EffectiveFrictionTime || istep!=file_end - # DataInspector(f) - display(f) - # sleep(nap) - # end + display(f) end @@ -668,7 +431,8 @@ end function Print2Disk( f, path, field, istep; res=4) path1 = path*"/_$field/" mkpath(path1) - save(path1*"$field"*@sprintf("%05d", istep)*".png", f, px_per_unit = res) + @show name = path1*"$field"*@sprintf("%05d", istep)*".png" + save(name, f, px_per_unit = res) end main() \ No newline at end of file diff --git a/JuliaVisualisation/_SpecificStudies/RiftingAnisotropy/ResolutionTestStress.jl b/JuliaVisualisation/_SpecificStudies/RiftingAnisotropy/ResolutionTestStress.jl new file mode 100644 index 00000000..1f5c7fec --- /dev/null +++ b/JuliaVisualisation/_SpecificStudies/RiftingAnisotropy/ResolutionTestStress.jl @@ -0,0 +1,515 @@ +import Pkg +Pkg.activate(normpath(joinpath(@__DIR__, "../.."))) +using HDF5, Printf, Colors, ColorSchemes, MathTeXEngine, LinearAlgebra, FFMPEG, Statistics, UnPack +using CairoMakie#, GLMakie +Mak = CairoMakie +Makie.update_theme!( fonts = (regular = texfont(), bold = texfont(:bold), italic = texfont(:italic))) +# fontsize_theme = Theme(fontsize=200) +# set_theme!(fontsize_theme) +const y = 365*24*3600 +const My = 1e6*y +const cm_y = y*100. + +function ExtractField(filename, field, size, mask_air, mask) + field = try (Float64.(reshape(ExtractData( filename, field), size...))) + catch + @warn "$field not found" + end + mask_air ? field[mask] .= NaN : nothing + return field +end + +function AddCountourQuivers!(PlotOnTop, ax1, coords, V, T, ϕ, σ1, ε̇1, Fab, height, Lc, cm_y, group_phases, Δ) + if PlotOnTop.quiver_origin + xc2D = coords.c.x .+ 0*coords.c.z' + zc2D = 0*coords.c.x .+ coords.c.z' + scatter!(ax1, xc2D[1:V.step:end,1:V.step:end][:]./Lc, zc2D[1:V.step:end,1:V.step:end][:]./Lc) + end + if PlotOnTop.T_contours + contour!(ax1, coords.c.x./Lc, coords.c.z./Lc, T, levels=0:200:1400, linewidth = 4, color=:black ) + end + if PlotOnTop.ph_contours + contour!(ax1, coords.c_hr.x./Lc, coords.c_hr.z./Lc, group_phases, levels=-1:1:maximum(group_phases), linewidth = 4, color=:white ) + end + if PlotOnTop.fabric + arrows!(ax1, coords.c.x[1:V.step:end]./Lc, coords.c.z[1:V.step:end]./Lc, Fab.x[1:V.step:end,1:V.step:end], Fab.z[1:V.step:end,1:V.step:end], arrowsize = 0, lengthscale=10Δ/1.5) + end + if PlotOnTop.σ1_axis + arrows!(ax1, coords.c.x[1:V.step:end]./Lc, coords.c.z[1:V.step:end]./Lc, σ1.x[1:V.step:end,1:V.step:end], σ1.z[1:V.step:end,1:V.step:end], arrowsize = 0, lengthscale=10Δ/1.5, color=:red) + end + if PlotOnTop.ε̇1_axis + arrows!(ax1, coords.c.x[1:V.step:end]./Lc, coords.c.z[1:V.step:end]./Lc, ε̇1.x[1:V.step:end,1:V.step:end], ε̇1.z[1:V.step:end,1:V.step:end], arrowsize = 0, lengthscale=10Δ/1.5, color=:black) + end + if PlotOnTop.vel_vec + arrows!(ax1, coords.c.x[1:V.step:end]./Lc, coords.c.z[1:V.step:end]./Lc, V.x[1:V.step:end,1:V.step:end]*cm_y, V.z[1:V.step:end,1:V.step:end]*cm_y, arrowsize = V.arrow, lengthscale = V.scale) + end + if PlotOnTop.topo + lines!(ax1, xv./Lc, height./Lc) + end + if PlotOnTop.ϕ_contours + contour!(ax1, coords.c.x./Lc, coords.c.z./Lc, ϕ, levels=0:1:0.1, linewidth = 6, color=:white, linestyle= :dashdot ) + end +end + +function ReadFile(path, step, scales, options, PlotOnTop) + + name = @sprintf("Output%05d.gzip.h5", step) + @info "Reading $(name)" + filename = string(path, name) + model = ExtractData( filename, "/Model/Params") + xc = ExtractData( filename, "/Model/xc_coord") + zc = ExtractData( filename, "/Model/zc_coord") + xv = ExtractData( filename, "/Model/xg_coord") + zv = ExtractData( filename, "/Model/zg_coord") + xvz = ExtractData( filename, "/Model/xvz_coord") + zvx = ExtractData( filename, "/Model/zvx_coord") + xv_hr = ExtractData( filename, "/VizGrid/xviz_hr") + zv_hr = ExtractData( filename, "/VizGrid/zviz_hr") + τzz_t = ExtractData( filename, "TimeSeries/szzd_mean_time") + P_t = ExtractData( filename, "TimeSeries/P_mean_time") + τxz_t = ExtractData( filename, "TimeSeries/sxz_mean_time") + t_t = ExtractData( filename, "TimeSeries/Time_time") + + xc_hr = 0.5.*(xv_hr[1:end-1] .+ xv_hr[2:end]) + zc_hr = 0.5.*(zv_hr[1:end-1] .+ zv_hr[2:end]) + ncx_hr, ncz_hr = length(xc_hr), length(zc_hr) + coords = ( c=(x=xc, z=zc), c_hr=(x=xc_hr, z=zc_hr), v=(x=xv, z=zv)) + + t = model[1] + tMy = round(t/scales.tc, digits=6) + nvx = Int(model[4]) + nvz = Int(model[5]) + ncx, ncz = nvx-1, nvz-1 + xmin, xmax = xv[1], xv[end] + zmin, zmax = zv[1], zv[end] + Lx, Lz = (xmax-xmin)/scales.Lc, (zv[end]-zv[1])/scales.Lc + Δx, Δz, Δ = Lx/ncx, Lz/ncz, sqrt( (Lx/ncx)^2 + (Lz/ncz)^2) + (xmax-xmin)>1e3 ? length_unit="km" : length_unit="m" + + @info "Model info" + @show "Model apect ratio" Lx/Lz + @show "Model time" t/My + @show length_unit + centroids, vertices = (ncx, ncz), (nvx, nvz) + + ph = ExtractField(filename, "/VizGrid/compo", centroids, false, 0) + mask_air = ph .== -1.00 + ph_hr = ExtractField(filename, "/VizGrid/compo_hr", (ncx_hr, ncz_hr), false, 0) + ph_dual_hr = ExtractField(filename, "/VizGrid/compo_dual_hr", (ncx_hr, ncz_hr), false, 0) + group_phases = copy(ph_hr); ph_hr[ph_hr.==-1.00] .= NaN; ph_dual_hr[ph_dual_hr.==-1.00] .= NaN; + ηc = ExtractField(filename, "/Centers/eta_n", centroids, true, mask_air) + ρc = Float64.(reshape(ExtractData( filename, "/Centers/rho_n"), ncx, ncz)); ρc[mask_air] .= NaN + P = Float64.(reshape(ExtractData( filename, "/Centers/P"), ncx, ncz)); P[mask_air] .= NaN + T = Float64.(reshape(ExtractData( filename, "/Centers/T"), ncx, ncz)) .- 273.15; T[mask_air] .= NaN + d = Float64.(reshape(ExtractData( filename, "/Centers/d"), ncx, ncz)); d[mask_air] .= NaN + ε̇pl = Float64.(reshape(ExtractData( filename, "/Centers/eII_pl"), ncx, ncz)); ε̇pl[mask_air] .= NaN + Vx = Float64.(reshape(ExtractData( filename, "/VxNodes/Vx"), (ncx+1), (ncz+2))) + Vz = Float64.(reshape(ExtractData( filename, "/VzNodes/Vz"), (ncx+2), (ncz+1))) + τxx = Float64.(reshape(ExtractData( filename, "/Centers/sxxd"), ncx, ncz)) + τzz = Float64.(reshape(ExtractData( filename, "/Centers/szzd"), ncx, ncz)) + τyy = -(τzz .+ τxx) + τxz = Float64.(reshape(ExtractData( filename, "/Vertices/sxz"), nvx, nvz)) + ε̇xx = Float64.(reshape(ExtractData( filename, "/Centers/exxd"), ncx, ncz)) + ε̇zz = Float64.(reshape(ExtractData( filename, "/Centers/ezzd"), ncx, ncz)) + ε̇yy = -(ε̇xx .+ ε̇zz) + ε̇xz = Float64.(reshape(ExtractData( filename, "/Vertices/exz"), nvx, nvz)) + τII = sqrt.( 0.5*(τxx.^2 .+ τyy.^2 .+ τzz.^2 .+ 0.5*(τxz[1:end-1,1:end-1].^2 .+ τxz[2:end,1:end-1].^2 .+ τxz[1:end-1,2:end].^2 .+ τxz[2:end,2:end].^2 ) ) ); τII[mask_air] .= NaN + ε̇II = sqrt.( 0.5*(ε̇xx.^2 .+ ε̇yy.^2 .+ ε̇zz.^2 .+ 0.5*(ε̇xz[1:end-1,1:end-1].^2 .+ ε̇xz[2:end,1:end-1].^2 .+ ε̇xz[1:end-1,2:end].^2 .+ ε̇xz[2:end,2:end].^2 ) ) ); ε̇II[mask_air] .= NaN + τxzc = 0.25*(τxz[1:end-1,1:end-1] .+ τxz[2:end,1:end-1] .+ τxz[1:end-1,2:end] .+ τxz[2:end,2:end]) + C = Float64.(reshape(ExtractData( filename, "/Centers/cohesion"), ncx, ncz)) + ϕ = ExtractField(filename, "/Centers/phi", centroids, false, 0) + divu = ExtractField(filename, "/Centers/divu", centroids, false, 0) + T_hr = zeros(size(ph_hr)); T_hr[1:2:end-1,1:2:end-1] .= T; T_hr[2:2:end-0,2:2:end-0] .= T + if options.LAB_color + ph_hr[(ph_hr.==2 .|| ph_hr.==3) .&& T_hr.options.LAB_T] .= 3 + ph_dual_hr[(ph_dual_hr.==2 .|| ph_dual_hr.==3) .&& T_hr.options.LAB_T] .= 3 + ph_dual_hr[(ph_dual_hr.==2 .|| ph_dual_hr.==6) .&& T_hr.options.LAB_T] .= 3 + end + + Fab = 0. + if PlotOnTop.fabric + δani = ExtractField(filename, "/Centers/ani_fac", centroids, false, 0) + Nx = Float64.(reshape(ExtractData( filename, "/Centers/nx"), ncx, ncz)) + Nz = Float64.(reshape(ExtractData( filename, "/Centers/nz"), ncx, ncz)) + Fab = (x=-Nz./Nx, z=ones(size(Nz))) + nrm = sqrt.(Fab.x.^2 .+ Fab.z.^2) + Fab.x ./= nrm + Fab.z ./= nrm + end + height = 0. + if PlotOnTop.topo + height = Float64.(ExtractData( filename, "/Topo/z_grid")); + Vx_grid = Float64.(ExtractData( filename, "/Topo/Vx_grid")); + Vz_grid = Float64.(ExtractData( filename, "/Topo/Vz_grid")); + Vx_mark = Float64.(ExtractData( filename, "/Topo/Vx_mark")); + Vz_mark = Float64.(ExtractData( filename, "/Topo/Vz_mark")); + x_mark = Float64.(ExtractData( filename, "/Topo/x_mark")); + z_mark = Float64.(ExtractData( filename, "/Topo/z_mark")); + end + Vxc = 0.5 .* (Vx[1:end-1,2:end-1] .+ Vx[2:end-0,2:end-1]) + Vzc = 0.5 .* (Vz[2:end-1,1:end-1] .+ Vz[2:end-1,2:end-0]) + V = (x=Vxc, z=Vzc, arrow=options.vel_arrow, step=options.vel_step, scale=options.vel_scale) + σ1, ε̇1 = 0., 0. + if PlotOnTop.σ1_axis + σ1 = PrincipalStress(τxx, τzz, τxz, P) + end + if PlotOnTop.ε̇1_axis + ε̇1 = PrincipalStress(ε̇xx, ε̇zz, ε̇xz, -1/3*divu) + @show extrema(ε̇1.x) + + end + + return (tMy=tMy,length_unit=length_unit, Lx=Lx, Lz=Lz, xc=xc, zc=zc, ε̇II=ε̇II, τII=τII, + coords=coords, V=V, T=T, ϕ=ϕ, σ1=σ1, ε̇1=ε̇1, Fab=Fab, height=height, group_phases=group_phases, Δ=Δ) +end + + +@views function main() + + # Set the path to your files + path_LR ="/Users/tduretz/REPO/MDOODZ7.0/RUNS/RiftingAnisotropy/ref_d4_LR/" + path_MR ="/Users/tduretz/REPO/MDOODZ7.0/RUNS/RiftingAnisotropy/ref_d4_MR/" + path_HR ="/Users/tduretz/REPO/MDOODZ7.0/RUNS/RiftingAnisotropy/ref_d4_HR/" + + + # File numbers + step = [120; 160; 360] + step = [500; 1100; 2220] + # Select field to visualise + # field = :Phases + # field = :Cohesion + # field = :Density + # field = :Viscosity + # field = :PlasticStrainrate + # field = :Stress + field = :StrainRate + # field = :Pressure + # field = :Divergence + # field = :Temperature + # field = :Velocity_x + # field = :Velocity_z + # field = :Velocity + # field = :GrainSize + # field = :Topography + # field = :TimeSeries + # field = :AnisotropyFactor + # field = :MeltFraction + # field = :TimeSeries + # field = :EffectiveFrictionTime + # field = :ChristmasTree + + # # Define Tuple for enlargment window + # zoom = ( + # xmin = 125, + # xmax = 150, + # zmin = -15, + # zmax = 0.e3, + # ) + + # Switches + PlotOnTop = ( + ph_contours = true, # add phase contours + fabric = true, # add fabric quiver (normal to director) + T_contours = false, # add temperature contours + topo = false, + quiver_origin = false, + σ1_axis = false, + ε̇1_axis = false, + vel_vec = false, + ϕ_contours = false, + ) + options = ( + printfig = true, # print figures to disk + printvid = false, + framerate = 6, + α_heatmap = 0.5, # transparency of heatmap + vel_arrow = 5, + vel_scale = 10000, + vel_step = 20, + nap = 0.1, # pause for animation + resol = 1000, + Lx = 1.0, + Lz = 1.0, + LAB_color = true, + LAB_T = 1250, + ) + + # Scaling + # Lc = 1000. + # tc = My + # Vc = 1e-9 + + scales = ( + Lc = 1e3, + tc = My, + Vc = 1.0, + τc = 1e6, + ) + + probe = (ϕeff = Float64.([]), t = Float64.([])) + cm_yr = 100.0*3600.0*24.0*365.25 + + ##################################### + + # Color palette for phase map + cmap = zeros(RGB{Float64}, 7) + cmap[1] = RGBA{Float64}(210/255, 218/255, 205/255, 1.) + cmap[2] = RGBA{Float64}(207/255, 089/255, 087/255, 1.) + cmap[3] = RGBA{Float64}(117/255, 164/255, 148/255, 1.) + cmap[4] = RGBA{Float64}(213/255, 213/255, 209/255, 1.) + cmap[5] = RGBA{Float64}(117/255, 099/255, 097/255, 1.) + cmap[6] = RGBA{Float64}(244/255, 218/255, 205/255, 1.) + cmap[7] = RGBA{Float64}(223/255, 233/255, 219/255, 1.) + phase_colors = cgrad(cmap, length(cmap), categorical=true, rev=false) + + ##################################### + + @unpack Lc, tc, Vc, τc = scales + + model = ReadFile(path_LR, step[1], scales, options, PlotOnTop) + @unpack tMy, length_unit, Lx, Lz, xc, zc, ε̇II, τII, coords, V, T, ϕ, σ1, Fab, height, group_phases, Δ = model + + if @isdefined zoom + Lx = (zoom.xmax - zoom.xmin)./scales.Lc + Lz = (zoom.zmax - zoom.zmin)./scales.Lc + window = zoom + else + window = ( + xmin = minimum(xc)/scales.Lc, + xmax = maximum(xc)/scales.Lc, + zmin = minimum(zc)/scales.Lc, + zmax = maximum(zc)/scales.Lc, + ) + end + + ##################################### + ftsz = 12*options.resol/500 + f = Figure(size = (1.1*Lx/4/Lz*options.resol*1.2, options.resol), fontsize=ftsz) + + # if field==:Phases + # ax1 = Axis(f[1, 1], title = L"Phases at $t$ = %$(tMy) Ma", xlabel = L"$x$ [m]", ylabel = L"$y$ [m]") + # hm = heatmap!(ax1, xc_hr./Lc, zc_hr./Lc, ph_hr, colormap = phase_colors) + # # hm = heatmap!(ax1, xc_hr./Lc, zc_hr./Lc, ph_hr, colormap = :turbo) + # hm = heatmap!(ax1, xc_hr./Lc, zc_hr./Lc, ph_dual_hr, colormap = :turbo) + # AddCountourQuivers!(PlotOnTop, ax1, coords, V, T, ϕ, σ1, Fab, height, Lc, cm_y, group_phases, Δ) + # colsize!(f.layout, 1, Aspect(1, Lx/Lz)) + # Mak.Colorbar(f[1, 2], hm, label = "Phases", width = 20, labelsize = ftsz, ticklabelsize = ftsz ) + # Mak.colgap!(f.layout, 20) + # xlims!(ax1, window.xmin, window.xmax) + # ylims!(ax1, window.zmin, window.zmax) + # if printfig Print2Disk( f, path, string(field), istep) end + # end + + # if field==:Viscosity + # ax1 = Axis(f[1, 1], title = L"$\eta$ at $t$ = %$(tMy) Ma", xlabel = L"$x$ [m]", ylabel = L"$y$ [m]") + # hm = heatmap!(ax1, xc./Lc, zc./Lc, log10.(ηc), colormap = (:turbo, α_heatmap)) + # AddCountourQuivers!(PlotOnTop, ax1, coords, V, T, ϕ, σ1, Fab, height, Lc, cm_y, group_phases, Δ) + # colsize!(f.layout, 1, Aspect(1, Lx/Lz)) + # Mak.Colorbar(f[1, 2], hm, label = L"$\eta$ [Pa.s]", width = 20, labelsize = ftsz, ticklabelsize = ftsz ) + # Mak.colgap!(f.layout, 20) + # xlims!(ax1, window.xmin, window.xmax) + # ylims!(ax1, window.zmin, window.zmax) + # if printfig Print2Disk( f, path, string(field), istep) end + # end + + # if field==:Density + # ax1 = Axis(f[1, 1], title = L"$\rho$ at $t$ = %$(tMy) Ma", xlabel = L"$x$ [km]", ylabel = L"$y$ [km]") + # hm = heatmap!(ax1, xc./Lc, zc./Lc, ρc, colormap = (:turbo, α_heatmap), colorrange=(2600, 3000)) + # AddCountourQuivers!(PlotOnTop, ax1, coords, V, T, ϕ, σ1, Fab, height, Lc, cm_y, group_phases, Δ) + # colsize!(f.layout, 1, Aspect(1, Lx/Lz)) + # Mak.Colorbar(f[1, 2], hm, label = L"$\rho$ [kg.m$^{-3}$]", width = 20, labelsize = ftsz, ticklabelsize = ftsz ) + # Mak.colgap!(f.layout, 20) + # xlims!(ax1, window.xmin, window.xmax) + # ylims!(ax1, window.zmin, window.zmax) + # if printfig Print2Disk( f, path, string(field), istep) end + # end + + # if field==:Stress + # ax1 = Axis(f[1, 1], title = L"$\tau_\textrm{II}$ at $t$ = %$(tMy) Ma", xlabel = L"$x$ [km]", ylabel = L"$y$ [km]") + # hm = heatmap!(ax1, xc./Lc, zc./Lc, τII./1e6, colormap = (:turbo, α_heatmap)) + # AddCountourQuivers!(PlotOnTop, ax1, coords, V, T, ϕ, σ1, Fab, height, Lc, cm_y, group_phases, Δ) + # colsize!(f.layout, 1, Aspect(1, Lx/Lz)) + # Mak.Colorbar(f[1, 2], hm, label = L"$\tau_\textrm{II}$ [MPa]", width = 20, labelsize = ftsz, ticklabelsize = ftsz ) + # Mak.colgap!(f.layout, 20) + # xlims!(ax1, window.xmin, window.xmax) + # ylims!(ax1, window.zmin, window.zmax) + # if printfig Print2Disk( f, path, string(field), istep) end + # end + + # if field==:Pressure + # ax1 = Axis(f[1, 1], title = L"$P$ at $t$ = %$(tMy) Ma", xlabel = L"$x$ [km]", ylabel = L"$y$ [km]") + # hm = heatmap!(ax1, xc./Lc, zc./Lc, P./1e9, colormap = (:turbo, α_heatmap)) #, colorrange=(1,1.2)1e4*365*24*3600 + # AddCountourQuivers!(PlotOnTop, ax1, coords, V, T, ϕ, σ1, Fab, height, Lc, cm_y, group_phases, Δ) + # colsize!(f.layout, 1, Aspect(1, Lx/Lz)) + # Mak.Colorbar(f[1, 2], hm, label = L"$P$ [GPa]", width = 20, labelsize = ftsz, ticklabelsize = ftsz ) + # Mak.colgap!(f.layout, 20) + # xlims!(ax1, window.xmin, window.xmax) + # ylims!(ax1, window.zmin, window.zmax) + # if printfig Print2Disk( f, path, string(field), istep) end + # end + + # if field==:Divergence + # ax1 = Axis(f[1, 1], title = L"∇⋅V at $t$ = %$(tMy) Ma", xlabel = L"$x$ [km]", ylabel = L"$y$ [km]") + # hm = heatmap!(ax1, xc./Lc, zc./Lc, divu, colormap = (:turbo, α_heatmap)) + # AddCountourQuivers!(PlotOnTop, ax1, coords, V, T, ϕ, σ1, Fab, height, Lc, cm_y, group_phases, Δ) + # colsize!(f.layout, 1, Aspect(1, Lx/Lz)) + # Mak.Colorbar(f[1, 2], hm, label = L"∇⋅V [s$^{-1}$]", width = 20, labelsize = ftsz, ticklabelsize = ftsz ) + # Mak.colgap!(f.layout, 20) + # xlims!(ax1, window.xmin, window.xmax) + # ylims!(ax1, window.zmin, window.zmax) + # if printfig Print2Disk( f, path, string(field), istep) end + # end + + if field==:StrainRate + + cmap = :vik + + options = ( + printfig = false, # print figures to disk + printvid = false, + framerate = 6, + α_heatmap = 0.75, # transparency of heatmap + vel_arrow = 5, + vel_scale = 10000, + vel_step = 5, + nap = 0.1, # pause for animation + resol = 1000, + Lx = 1.0, + Lz = 1.0, + LAB_color = true, + LAB_T = 1250, + ) + + model = ReadFile(path_LR, step[1], scales, options, PlotOnTop) + @unpack tMy, length_unit, Lx, Lz, xc, zc, ε̇II, τII, coords, V, T, ϕ, σ1, ε̇1, Fab, height, group_phases, Δ = model + + tMy_string = @sprintf("%1.2lf", tMy) + ax1 = Axis(f[1, 1], title = L"A) $\tau_\textrm{II}$ at $t$ = %$(tMy_string) Ma - 300 $\times$ 200 cells", ylabel = L"$y$ [%$(length_unit)]") + hm = heatmap!(ax1, xc./Lc, zc./Lc, log10.(τII), colormap = (cmap, options.α_heatmap), colorrange=(6.69, 8.69)) + AddCountourQuivers!(PlotOnTop, ax1, coords, V, T, ϕ, σ1, ε̇1, Fab, height, Lc, cm_y, group_phases, Δ/2) + colsize!(f.layout, 1, Aspect(1, Lx/Lz)) + Mak.Colorbar(f[1, 2], hm, label = L"$\tau_\textrm{II}$ [s$^{-1}$]", width = 20, labelsize = ftsz, ticklabelsize = ftsz ) + Mak.colgap!(f.layout, 20) + xlims!(ax1, window.xmin, window.xmax) + ylims!(ax1, window.zmin, window.zmax) + + # lines!(ax1, xc./Lc, log10.(ε̇II[:, 170]) ) + # lines!(ax1, xc./Lc, atand.(ε̇1.z[:, 170]) ./ ε̇1.x[:, 170] ) + # xlims!(ax1, window.xmin, window.xmax) + # ylims!(ax1, -180, 180) + + ############################################################## + + + options = ( + printfig = false, # print figures to disk + printvid = false, + framerate = 6, + α_heatmap = 0.75, # transparency of heatmap + vel_arrow = 5, + vel_scale = 10000, + vel_step = 10, + nap = 0.1, # pause for animation + resol = 1000, + Lx = 1.0, + Lz = 1.0, + LAB_color = true, + LAB_T = 1250, + ) + + model = ReadFile(path_MR, step[2], scales, options, PlotOnTop) + @unpack tMy, length_unit, Lx, Lz, xc, zc, ε̇II, τII, coords, V, T, ϕ, σ1, ε̇1, Fab, height, group_phases, Δ = model + + tMy_string = @sprintf("%1.2lf", tMy) + ax1 = Axis(f[2, 1], title = L"B) $\tau_\textrm{II}$ at $t$ = %$(tMy_string) Ma - 600 $\times$ 400 cells", ylabel = L"$y$ [%$(length_unit)]") + hm = heatmap!(ax1, xc./Lc, zc./Lc, log10.(τII), colormap = (cmap, options.α_heatmap), colorrange=(6.69, 8.69)) + AddCountourQuivers!(PlotOnTop, ax1, coords, V, T, ϕ, σ1, ε̇1, Fab, height, Lc, cm_y, group_phases, Δ) + colsize!(f.layout, 1, Aspect(1, Lx/Lz)) + Mak.Colorbar(f[2, 2], hm, label = L"$\tau_\textrm{II}$ [s$^{-1}$]", width = 20, labelsize = ftsz, ticklabelsize = ftsz ) + Mak.colgap!(f.layout, 20) + xlims!(ax1, window.xmin, window.xmax) + ylims!(ax1, window.zmin, window.zmax) + + + # lines!(ax1, xc./Lc, log10.(ε̇II[:, 350]) ) + # lines!(ax1, xc./Lc, atand.(ε̇1.z[:, 350]) ./ ε̇1.x[:, 350] ) + # xlims!(ax1, window.xmin, window.xmax) + # ylims!(ax1, -180, 180) + + + # ############################################################## + + options = ( + printfig = true, # print figures to disk + printvid = false, + framerate = 6, + α_heatmap = 0.75, # transparency of heatmap + vel_arrow = 5, + vel_scale = 10000, + vel_step = 20, + nap = 0.1, # pause for animation + resol = 1000, + Lx = 1.0, + Lz = 1.0, + LAB_color = true, + LAB_T = 1250, + ) + + model = ReadFile(path_HR, step[3], scales, options, PlotOnTop) + @unpack tMy, length_unit, Lx, Lz, xc, zc, ε̇II, coords, V, T, ϕ, σ1, ε̇1, Fab, height, group_phases, Δ = model + + tMy_string = @sprintf("%1.2lf", tMy) + ax1 = Axis(f[3, 1], title = L"C) $\tau_\textrm{II}$ at $t$ = %$(tMy_string) Ma - 1200 $\times$ 800 cells", xlabel = L"$x$ [%$(length_unit)]", ylabel = L"$y$ [%$(length_unit)]") + hm = heatmap!(ax1, xc./Lc, zc./Lc, log10.(τII), colormap = (cmap, options.α_heatmap), colorrange=(6.69, 8.69)) + AddCountourQuivers!(PlotOnTop, ax1, coords, V, T, ϕ, σ1, ε̇1, Fab, height, Lc, cm_y, group_phases, Δ*2) + colsize!(f.layout, 1, Aspect(1, Lx/Lz)) + Mak.Colorbar(f[3, 2], hm, label = L"$\tau_\textrm{II}$ [s$^{-1}$]", width = 20, labelsize = ftsz, ticklabelsize = ftsz ) + Mak.colgap!(f.layout, 20) + xlims!(ax1, window.xmin, window.xmax) + ylims!(ax1, window.zmin, window.zmax) + + # lines!(ax1, xc./Lc, log10.(ε̇II[:, 700]) ) + # lines!(ax1, xc./Lc, atand.(ε̇1.z[:, 700]) ./ ε̇1.x[:, 700] ) + # xlims!(ax1, window.xmin, window.xmax) + # ylims!(ax1, -180, 180) + + + ############################################################## + + save("/Users/tduretz/PowerFolders/_manuscripts/RiftingAnisotropy/Figures/ResolutionTest.png", f, px_per_unit = 4) end + + display(f) +end + +function PrincipalStress(τxx, τzz, τxz, P) + σ1 = (x=zeros(size(τxx)), z=zeros(size(τxx)) ) + τxzc = 0.25*(τxz[1:end-1,1:end-1] .+ τxz[2:end-0,1:end-1] .+ τxz[1:end-1,2:end-0] .+ τxz[2:end-0,2:end-0]) + for i in eachindex(τxzc) + if P[i]>-1e-13 + σ = [-P[i]+τxx[i] τxzc[i]; τxzc[i] -P[i]+τzz[i]] + v = eigvecs(σ) + σ1.x[i] = v[1,1] + σ1.z[i] = v[2,1] + end + end + return σ1 +end + +function ExtractData( file_path, data_path) + data = h5open(file_path, "r") do file + read(file, data_path) + end + return data +end + +function Print2Disk( f, path, field, istep; res=4) + path1 = path*"/_$field/" + mkpath(path1) + @show name = path1*"$field"*@sprintf("%05d", istep)*".png" + save(name, f, px_per_unit = res) +end + +main() \ No newline at end of file diff --git a/JuliaVisualisation/_SpecificStudies/RiftingAnisotropy/TimeEvolution.jl b/JuliaVisualisation/_SpecificStudies/RiftingAnisotropy/TimeEvolution.jl new file mode 100644 index 00000000..d67e9f22 --- /dev/null +++ b/JuliaVisualisation/_SpecificStudies/RiftingAnisotropy/TimeEvolution.jl @@ -0,0 +1,516 @@ +import Pkg +Pkg.activate(normpath(joinpath(@__DIR__, "../.."))) +using HDF5, Printf, Colors, ColorSchemes, MathTeXEngine, LinearAlgebra, FFMPEG, Statistics, UnPack +using CairoMakie#, GLMakie +Mak = CairoMakie +Makie.update_theme!( fonts = (regular = texfont(), bold = texfont(:bold), italic = texfont(:italic))) +# fontsize_theme = Theme(fontsize=200) +# set_theme!(fontsize_theme) +const y = 365*24*3600 +const My = 1e6*y +const cm_y = y*100. + +function ExtractField(filename, field, size, mask_air, mask) + field = try (Float64.(reshape(ExtractData( filename, field), size...))) + catch + @warn "$field not found" + end + mask_air ? field[mask] .= NaN : nothing + return field +end + +function AddCountourQuivers!(PlotOnTop, ax1, coords, V, T, ϕ, σ1, ε̇1, Fab, height, Lc, cm_y, group_phases, Δ) + if PlotOnTop.quiver_origin + xc2D = coords.c.x .+ 0*coords.c.z' + zc2D = 0*coords.c.x .+ coords.c.z' + scatter!(ax1, xc2D[1:V.step:end,1:V.step:end][:]./Lc, zc2D[1:V.step:end,1:V.step:end][:]./Lc) + end + if PlotOnTop.T_contours + contour!(ax1, coords.c.x./Lc, coords.c.z./Lc, T, levels=0:200:1400, linewidth = 4, color=:black ) + end + if PlotOnTop.ph_contours + contour!(ax1, coords.c_hr.x./Lc, coords.c_hr.z./Lc, group_phases, levels=-1:1:maximum(group_phases), linewidth = 4, color=:white ) + end + if PlotOnTop.fabric + arrows!(ax1, coords.c.x[1:V.step:end]./Lc, coords.c.z[1:V.step:end]./Lc, Fab.x[1:V.step:end,1:V.step:end], Fab.z[1:V.step:end,1:V.step:end], arrowsize = 0, lengthscale=10Δ/1.5) + end + if PlotOnTop.σ1_axis + arrows!(ax1, coords.c.x[1:V.step:end]./Lc, coords.c.z[1:V.step:end]./Lc, σ1.x[1:V.step:end,1:V.step:end], σ1.z[1:V.step:end,1:V.step:end], arrowsize = 0, lengthscale=10Δ/1.5, color=:red) + end + if PlotOnTop.ε̇1_axis + arrows!(ax1, coords.c.x[1:V.step:end]./Lc, coords.c.z[1:V.step:end]./Lc, ε̇1.x[1:V.step:end,1:V.step:end], ε̇1.z[1:V.step:end,1:V.step:end], arrowsize = 0, lengthscale=10Δ/1.5, color=:black) + end + if PlotOnTop.vel_vec + arrows!(ax1, coords.c.x[1:V.step:end]./Lc, coords.c.z[1:V.step:end]./Lc, V.x[1:V.step:end,1:V.step:end]*cm_y, V.z[1:V.step:end,1:V.step:end]*cm_y, arrowsize = V.arrow, lengthscale = V.scale) + end + if PlotOnTop.topo + lines!(ax1, xv./Lc, height./Lc) + end + if PlotOnTop.ϕ_contours + contour!(ax1, coords.c.x./Lc, coords.c.z./Lc, ϕ, levels=0:1:0.1, linewidth = 6, color=:white, linestyle= :dashdot ) + end +end + +function ReadFile(path, step, scales, options, PlotOnTop) + + name = @sprintf("Output%05d.gzip.h5", step) + @info "Reading $(name)" + filename = string(path, name) + model = ExtractData( filename, "/Model/Params") + xc = ExtractData( filename, "/Model/xc_coord") + zc = ExtractData( filename, "/Model/zc_coord") + xv = ExtractData( filename, "/Model/xg_coord") + zv = ExtractData( filename, "/Model/zg_coord") + xvz = ExtractData( filename, "/Model/xvz_coord") + zvx = ExtractData( filename, "/Model/zvx_coord") + xv_hr = ExtractData( filename, "/VizGrid/xviz_hr") + zv_hr = ExtractData( filename, "/VizGrid/zviz_hr") + τzz_t = ExtractData( filename, "TimeSeries/szzd_mean_time") + P_t = ExtractData( filename, "TimeSeries/P_mean_time") + τxz_t = ExtractData( filename, "TimeSeries/sxz_mean_time") + t_t = ExtractData( filename, "TimeSeries/Time_time") + + xc_hr = 0.5.*(xv_hr[1:end-1] .+ xv_hr[2:end]) + zc_hr = 0.5.*(zv_hr[1:end-1] .+ zv_hr[2:end]) + ncx_hr, ncz_hr = length(xc_hr), length(zc_hr) + coords = ( c=(x=xc, z=zc), c_hr=(x=xc_hr, z=zc_hr), v=(x=xv, z=zv)) + + t = model[1] + tMy = round(t/scales.tc, digits=6) + nvx = Int(model[4]) + nvz = Int(model[5]) + ncx, ncz = nvx-1, nvz-1 + xmin, xmax = xv[1], xv[end] + zmin, zmax = zv[1], zv[end] + Lx, Lz = (xmax-xmin)/scales.Lc, (zv[end]-zv[1])/scales.Lc + Δx, Δz, Δ = Lx/ncx, Lz/ncz, sqrt( (Lx/ncx)^2 + (Lz/ncz)^2) + (xmax-xmin)>1e3 ? length_unit="km" : length_unit="m" + + @info "Model info" + @show "Model apect ratio" Lx/Lz + @show "Model time" t/My + @show length_unit + centroids, vertices = (ncx, ncz), (nvx, nvz) + + ph = ExtractField(filename, "/VizGrid/compo", centroids, false, 0) + mask_air = ph .== -1.00 + ph_hr = ExtractField(filename, "/VizGrid/compo_hr", (ncx_hr, ncz_hr), false, 0) + ph_dual_hr = ExtractField(filename, "/VizGrid/compo_dual_hr", (ncx_hr, ncz_hr), false, 0) + group_phases = copy(ph_hr); ph_hr[ph_hr.==-1.00] .= NaN; ph_dual_hr[ph_dual_hr.==-1.00] .= NaN; + ηc = ExtractField(filename, "/Centers/eta_n", centroids, true, mask_air) + ρc = Float64.(reshape(ExtractData( filename, "/Centers/rho_n"), ncx, ncz)); ρc[mask_air] .= NaN + P = Float64.(reshape(ExtractData( filename, "/Centers/P"), ncx, ncz)); P[mask_air] .= NaN + T = Float64.(reshape(ExtractData( filename, "/Centers/T"), ncx, ncz)) .- 273.15; T[mask_air] .= NaN + d = Float64.(reshape(ExtractData( filename, "/Centers/d"), ncx, ncz)); d[mask_air] .= NaN + ε̇pl = Float64.(reshape(ExtractData( filename, "/Centers/eII_pl"), ncx, ncz)); ε̇pl[mask_air] .= NaN + Vx = Float64.(reshape(ExtractData( filename, "/VxNodes/Vx"), (ncx+1), (ncz+2))) + Vz = Float64.(reshape(ExtractData( filename, "/VzNodes/Vz"), (ncx+2), (ncz+1))) + τxx = Float64.(reshape(ExtractData( filename, "/Centers/sxxd"), ncx, ncz)) + τzz = Float64.(reshape(ExtractData( filename, "/Centers/szzd"), ncx, ncz)) + τyy = -(τzz .+ τxx) + τxz = Float64.(reshape(ExtractData( filename, "/Vertices/sxz"), nvx, nvz)) + ε̇xx = Float64.(reshape(ExtractData( filename, "/Centers/exxd"), ncx, ncz)) + ε̇zz = Float64.(reshape(ExtractData( filename, "/Centers/ezzd"), ncx, ncz)) + ε̇yy = -(ε̇xx .+ ε̇zz) + ε̇xz = Float64.(reshape(ExtractData( filename, "/Vertices/exz"), nvx, nvz)) + τII = sqrt.( 0.5*(τxx.^2 .+ τyy.^2 .+ τzz.^2 .+ 0.5*(τxz[1:end-1,1:end-1].^2 .+ τxz[2:end,1:end-1].^2 .+ τxz[1:end-1,2:end].^2 .+ τxz[2:end,2:end].^2 ) ) ); τII[mask_air] .= NaN + ε̇II = sqrt.( 0.5*(ε̇xx.^2 .+ ε̇yy.^2 .+ ε̇zz.^2 .+ 0.5*(ε̇xz[1:end-1,1:end-1].^2 .+ ε̇xz[2:end,1:end-1].^2 .+ ε̇xz[1:end-1,2:end].^2 .+ ε̇xz[2:end,2:end].^2 ) ) ); ε̇II[mask_air] .= NaN + τxzc = 0.25*(τxz[1:end-1,1:end-1] .+ τxz[2:end,1:end-1] .+ τxz[1:end-1,2:end] .+ τxz[2:end,2:end]) + C = Float64.(reshape(ExtractData( filename, "/Centers/cohesion"), ncx, ncz)) + ϕ = ExtractField(filename, "/Centers/phi", centroids, false, 0) + divu = ExtractField(filename, "/Centers/divu", centroids, false, 0) + T_hr = zeros(size(ph_hr)); T_hr[1:2:end-1,1:2:end-1] .= T; T_hr[2:2:end-0,2:2:end-0] .= T + if options.LAB_color + ph_hr[(ph_hr.==2 .|| ph_hr.==3) .&& T_hr.options.LAB_T] .= 3 + ph_dual_hr[(ph_dual_hr.==2 .|| ph_dual_hr.==3) .&& T_hr.options.LAB_T] .= 3 + ph_dual_hr[(ph_dual_hr.==2 .|| ph_dual_hr.==6) .&& T_hr.options.LAB_T] .= 3 + end + + Fab = 0. + if PlotOnTop.fabric + δani = ExtractField(filename, "/Centers/ani_fac", centroids, false, 0) + Nx = Float64.(reshape(ExtractData( filename, "/Centers/nx"), ncx, ncz)) + Nz = Float64.(reshape(ExtractData( filename, "/Centers/nz"), ncx, ncz)) + Fab = (x=-Nz./Nx, z=ones(size(Nz))) + nrm = sqrt.(Fab.x.^2 .+ Fab.z.^2) + Fab.x ./= nrm + Fab.z ./= nrm + end + height = 0. + if PlotOnTop.topo + height = Float64.(ExtractData( filename, "/Topo/z_grid")); + Vx_grid = Float64.(ExtractData( filename, "/Topo/Vx_grid")); + Vz_grid = Float64.(ExtractData( filename, "/Topo/Vz_grid")); + Vx_mark = Float64.(ExtractData( filename, "/Topo/Vx_mark")); + Vz_mark = Float64.(ExtractData( filename, "/Topo/Vz_mark")); + x_mark = Float64.(ExtractData( filename, "/Topo/x_mark")); + z_mark = Float64.(ExtractData( filename, "/Topo/z_mark")); + end + Vxc = 0.5 .* (Vx[1:end-1,2:end-1] .+ Vx[2:end-0,2:end-1]) + Vzc = 0.5 .* (Vz[2:end-1,1:end-1] .+ Vz[2:end-1,2:end-0]) + V = (x=Vxc, z=Vzc, arrow=options.vel_arrow, step=options.vel_step, scale=options.vel_scale) + σ1, ε̇1 = 0., 0. + if PlotOnTop.σ1_axis + σ1 = PrincipalStress(τxx, τzz, τxz, P) + end + if PlotOnTop.ε̇1_axis + ε̇1 = PrincipalStress(ε̇xx, ε̇zz, ε̇xz, -1/3*divu) + @show extrema(ε̇1.x) + + end + + return (tMy=tMy,length_unit=length_unit, Lx=Lx, Lz=Lz, xc=xc, zc=zc, ε̇II=ε̇II, τII=τII, + coords=coords, V=V, T=T, ϕ=ϕ, σ1=σ1, ε̇1=ε̇1, Fab=Fab, height=height, group_phases=group_phases, Δ=Δ) +end + + +@views function main() + + # Set the path to your files + path_LR ="/Users/tduretz/REPO/MDOODZ7.0/RUNS/RiftingAnisotropy/ref_d4_LR/" + path_MR ="/Users/tduretz/REPO/MDOODZ7.0/RUNS/RiftingAnisotropy/ref_d4_MR/" + path_HR ="/Users/tduretz/REPO/MDOODZ7.0/RUNS/RiftingAnisotropy/ref_d4_HR/" + + + # File numbers + step = [120; 160; 360] + step = [500; 1100; 2220] + # Select field to visualise + # field = :Phases + # field = :Cohesion + # field = :Density + # field = :Viscosity + # field = :PlasticStrainrate + # field = :Stress + field = :StrainRate + # field = :Pressure + # field = :Divergence + # field = :Temperature + # field = :Velocity_x + # field = :Velocity_z + # field = :Velocity + # field = :GrainSize + # field = :Topography + # field = :TimeSeries + # field = :AnisotropyFactor + # field = :MeltFraction + # field = :TimeSeries + # field = :EffectiveFrictionTime + # field = :ChristmasTree + + # # Define Tuple for enlargment window + # zoom = ( + # xmin = 125, + # xmax = 150, + # zmin = -15, + # zmax = 0.e3, + # ) + + # Switches + PlotOnTop = ( + ph_contours = true, # add phase contours + fabric = true, # add fabric quiver (normal to director) + T_contours = false, # add temperature contours + topo = false, + quiver_origin = false, + σ1_axis = false, + ε̇1_axis = false, + vel_vec = false, + ϕ_contours = false, + ) + options = ( + printfig = true, # print figures to disk + printvid = false, + framerate = 6, + α_heatmap = 0.5, # transparency of heatmap + vel_arrow = 5, + vel_scale = 10000, + vel_step = 20, + nap = 0.1, # pause for animation + resol = 1000, + Lx = 1.0, + Lz = 1.0, + LAB_color = true, + LAB_T = 1250, + ) + + # Scaling + # Lc = 1000. + # tc = My + # Vc = 1e-9 + + scales = ( + Lc = 1e3, + tc = My, + Vc = 1.0, + τc = 1e6, + ) + + probe = (ϕeff = Float64.([]), t = Float64.([])) + cm_yr = 100.0*3600.0*24.0*365.25 + + ##################################### + + # Color palette for phase map + cmap = zeros(RGB{Float64}, 7) + cmap[1] = RGBA{Float64}(210/255, 218/255, 205/255, 1.) + cmap[2] = RGBA{Float64}(207/255, 089/255, 087/255, 1.) + cmap[3] = RGBA{Float64}(117/255, 164/255, 148/255, 1.) + cmap[4] = RGBA{Float64}(213/255, 213/255, 209/255, 1.) + cmap[5] = RGBA{Float64}(117/255, 099/255, 097/255, 1.) + cmap[6] = RGBA{Float64}(244/255, 218/255, 205/255, 1.) + cmap[7] = RGBA{Float64}(223/255, 233/255, 219/255, 1.) + phase_colors = cgrad(cmap, length(cmap), categorical=true, rev=false) + + ##################################### + + @unpack Lc, tc, Vc, τc = scales + + model = ReadFile(path_LR, step[1], scales, options, PlotOnTop) + @unpack tMy, length_unit, Lx, Lz, xc, zc, ε̇II, τII, coords, V, T, ϕ, σ1, Fab, height, group_phases, Δ = model + + if @isdefined zoom + Lx = (zoom.xmax - zoom.xmin)./scales.Lc + Lz = (zoom.zmax - zoom.zmin)./scales.Lc + window = zoom + else + window = ( + xmin = minimum(xc)/scales.Lc, + xmax = maximum(xc)/scales.Lc, + zmin = minimum(zc)/scales.Lc, + zmax = maximum(zc)/scales.Lc, + ) + end + + ##################################### + ftsz = 12*options.resol/500 + f = Figure(size = (1.1*Lx/4/Lz*options.resol*1.2, options.resol), fontsize=ftsz) + + # if field==:Phases + # ax1 = Axis(f[1, 1], title = L"Phases at $t$ = %$(tMy) Ma", xlabel = L"$x$ [m]", ylabel = L"$y$ [m]") + # hm = heatmap!(ax1, xc_hr./Lc, zc_hr./Lc, ph_hr, colormap = phase_colors) + # # hm = heatmap!(ax1, xc_hr./Lc, zc_hr./Lc, ph_hr, colormap = :turbo) + # hm = heatmap!(ax1, xc_hr./Lc, zc_hr./Lc, ph_dual_hr, colormap = :turbo) + # AddCountourQuivers!(PlotOnTop, ax1, coords, V, T, ϕ, σ1, Fab, height, Lc, cm_y, group_phases, Δ) + # colsize!(f.layout, 1, Aspect(1, Lx/Lz)) + # Mak.Colorbar(f[1, 2], hm, label = "Phases", width = 20, labelsize = ftsz, ticklabelsize = ftsz ) + # Mak.colgap!(f.layout, 20) + # xlims!(ax1, window.xmin, window.xmax) + # ylims!(ax1, window.zmin, window.zmax) + # if printfig Print2Disk( f, path, string(field), istep) end + # end + + # if field==:Viscosity + # ax1 = Axis(f[1, 1], title = L"$\eta$ at $t$ = %$(tMy) Ma", xlabel = L"$x$ [m]", ylabel = L"$y$ [m]") + # hm = heatmap!(ax1, xc./Lc, zc./Lc, log10.(ηc), colormap = (:turbo, α_heatmap)) + # AddCountourQuivers!(PlotOnTop, ax1, coords, V, T, ϕ, σ1, Fab, height, Lc, cm_y, group_phases, Δ) + # colsize!(f.layout, 1, Aspect(1, Lx/Lz)) + # Mak.Colorbar(f[1, 2], hm, label = L"$\eta$ [Pa.s]", width = 20, labelsize = ftsz, ticklabelsize = ftsz ) + # Mak.colgap!(f.layout, 20) + # xlims!(ax1, window.xmin, window.xmax) + # ylims!(ax1, window.zmin, window.zmax) + # if printfig Print2Disk( f, path, string(field), istep) end + # end + + # if field==:Density + # ax1 = Axis(f[1, 1], title = L"$\rho$ at $t$ = %$(tMy) Ma", xlabel = L"$x$ [km]", ylabel = L"$y$ [km]") + # hm = heatmap!(ax1, xc./Lc, zc./Lc, ρc, colormap = (:turbo, α_heatmap), colorrange=(2600, 3000)) + # AddCountourQuivers!(PlotOnTop, ax1, coords, V, T, ϕ, σ1, Fab, height, Lc, cm_y, group_phases, Δ) + # colsize!(f.layout, 1, Aspect(1, Lx/Lz)) + # Mak.Colorbar(f[1, 2], hm, label = L"$\rho$ [kg.m$^{-3}$]", width = 20, labelsize = ftsz, ticklabelsize = ftsz ) + # Mak.colgap!(f.layout, 20) + # xlims!(ax1, window.xmin, window.xmax) + # ylims!(ax1, window.zmin, window.zmax) + # if printfig Print2Disk( f, path, string(field), istep) end + # end + + # if field==:Stress + # ax1 = Axis(f[1, 1], title = L"$\tau_\textrm{II}$ at $t$ = %$(tMy) Ma", xlabel = L"$x$ [km]", ylabel = L"$y$ [km]") + # hm = heatmap!(ax1, xc./Lc, zc./Lc, τII./1e6, colormap = (:turbo, α_heatmap)) + # AddCountourQuivers!(PlotOnTop, ax1, coords, V, T, ϕ, σ1, Fab, height, Lc, cm_y, group_phases, Δ) + # colsize!(f.layout, 1, Aspect(1, Lx/Lz)) + # Mak.Colorbar(f[1, 2], hm, label = L"$\tau_\textrm{II}$ [MPa]", width = 20, labelsize = ftsz, ticklabelsize = ftsz ) + # Mak.colgap!(f.layout, 20) + # xlims!(ax1, window.xmin, window.xmax) + # ylims!(ax1, window.zmin, window.zmax) + # if printfig Print2Disk( f, path, string(field), istep) end + # end + + # if field==:Pressure + # ax1 = Axis(f[1, 1], title = L"$P$ at $t$ = %$(tMy) Ma", xlabel = L"$x$ [km]", ylabel = L"$y$ [km]") + # hm = heatmap!(ax1, xc./Lc, zc./Lc, P./1e9, colormap = (:turbo, α_heatmap)) #, colorrange=(1,1.2)1e4*365*24*3600 + # AddCountourQuivers!(PlotOnTop, ax1, coords, V, T, ϕ, σ1, Fab, height, Lc, cm_y, group_phases, Δ) + # colsize!(f.layout, 1, Aspect(1, Lx/Lz)) + # Mak.Colorbar(f[1, 2], hm, label = L"$P$ [GPa]", width = 20, labelsize = ftsz, ticklabelsize = ftsz ) + # Mak.colgap!(f.layout, 20) + # xlims!(ax1, window.xmin, window.xmax) + # ylims!(ax1, window.zmin, window.zmax) + # if printfig Print2Disk( f, path, string(field), istep) end + # end + + # if field==:Divergence + # ax1 = Axis(f[1, 1], title = L"∇⋅V at $t$ = %$(tMy) Ma", xlabel = L"$x$ [km]", ylabel = L"$y$ [km]") + # hm = heatmap!(ax1, xc./Lc, zc./Lc, divu, colormap = (:turbo, α_heatmap)) + # AddCountourQuivers!(PlotOnTop, ax1, coords, V, T, ϕ, σ1, Fab, height, Lc, cm_y, group_phases, Δ) + # colsize!(f.layout, 1, Aspect(1, Lx/Lz)) + # Mak.Colorbar(f[1, 2], hm, label = L"∇⋅V [s$^{-1}$]", width = 20, labelsize = ftsz, ticklabelsize = ftsz ) + # Mak.colgap!(f.layout, 20) + # xlims!(ax1, window.xmin, window.xmax) + # ylims!(ax1, window.zmin, window.zmax) + # if printfig Print2Disk( f, path, string(field), istep) end + # end + + if field==:StrainRate + + cmap = :vik + + options = ( + printfig = false, # print figures to disk + printvid = false, + framerate = 6, + α_heatmap = 0.75, # transparency of heatmap + vel_arrow = 5, + vel_scale = 10000, + vel_step = 5, + nap = 0.1, # pause for animation + resol = 1000, + Lx = 1.0, + Lz = 1.0, + LAB_color = true, + LAB_T = 1250, + ) + + model = ReadFile(path_LR, step[1], scales, options, PlotOnTop) + @unpack tMy, length_unit, Lx, Lz, xc, zc, ε̇II, τII, coords, V, T, ϕ, σ1, ε̇1, Fab, height, group_phases, Δ = model + + tMy_string = @sprintf("%1.2lf", tMy) + ax1 = Axis(f[1, 1], title = L"A) $\tau_\textrm{II}$ at $t$ = %$(tMy_string) Ma - 300 $\times$ 200 cells", ylabel = L"$y$ [%$(length_unit)]") + hm = heatmap!(ax1, xc./Lc, zc./Lc, log10.(τII), colormap = (cmap, options.α_heatmap), colorrange=(6.69, 8.69)) + AddCountourQuivers!(PlotOnTop, ax1, coords, V, T, ϕ, σ1, ε̇1, Fab, height, Lc, cm_y, group_phases, Δ/2) + colsize!(f.layout, 1, Aspect(1, Lx/Lz)) + Mak.Colorbar(f[1, 2], hm, label = L"$\tau_\textrm{II}$ [Pa]", width = 20, labelsize = ftsz, ticklabelsize = ftsz ) + Mak.colgap!(f.layout, 20) + xlims!(ax1, window.xmin, window.xmax) + ylims!(ax1, window.zmin, window.zmax) + + # lines!(ax1, xc./Lc, log10.(ε̇II[:, 170]) ) + # lines!(ax1, xc./Lc, atand.(ε̇1.z[:, 170]) ./ ε̇1.x[:, 170] ) + # xlims!(ax1, window.xmin, window.xmax) + # ylims!(ax1, -180, 180) + + ############################################################## + + + options = ( + printfig = false, # print figures to disk + printvid = false, + framerate = 6, + α_heatmap = 0.75, # transparency of heatmap + vel_arrow = 5, + vel_scale = 10000, + vel_step = 10, + nap = 0.1, # pause for animation + resol = 1000, + Lx = 1.0, + Lz = 1.0, + LAB_color = true, + LAB_T = 1250, + ) + + model = ReadFile(path_MR, step[2], scales, options, PlotOnTop) + @unpack tMy, length_unit, Lx, Lz, xc, zc, ε̇II, τII, coords, V, T, ϕ, σ1, ε̇1, Fab, height, group_phases, Δ = model + + tMy_string = @sprintf("%1.2lf", tMy) + ax1 = Axis(f[2, 1], title = L"B) $\tau_\textrm{II}$ at $t$ = %$(tMy_string) Ma - 600 $\times$ 400 cells", ylabel = L"$y$ [%$(length_unit)]") + hm = heatmap!(ax1, xc./Lc, zc./Lc, log10.(τII), colormap = (cmap, options.α_heatmap), colorrange=(6.69, 8.69)) + AddCountourQuivers!(PlotOnTop, ax1, coords, V, T, ϕ, σ1, ε̇1, Fab, height, Lc, cm_y, group_phases, Δ) + colsize!(f.layout, 1, Aspect(1, Lx/Lz)) + Mak.Colorbar(f[2, 2], hm, label = L"$\tau_\textrm{II}$ [Pa]", width = 20, labelsize = ftsz, ticklabelsize = ftsz ) + Mak.colgap!(f.layout, 20) + xlims!(ax1, window.xmin, window.xmax) + ylims!(ax1, window.zmin, window.zmax) + + + # lines!(ax1, xc./Lc, log10.(ε̇II[:, 350]) ) + # lines!(ax1, xc./Lc, atand.(ε̇1.z[:, 350]) ./ ε̇1.x[:, 350] ) + # xlims!(ax1, window.xmin, window.xmax) + # ylims!(ax1, -180, 180) + + + # ############################################################## + + options = ( + printfig = true, # print figures to disk + printvid = false, + framerate = 6, + α_heatmap = 0.75, # transparency of heatmap + vel_arrow = 5, + vel_scale = 10000, + vel_step = 20, + nap = 0.1, # pause for animation + resol = 1000, + Lx = 1.0, + Lz = 1.0, + LAB_color = true, + LAB_T = 1250, + ) + + model = ReadFile(path_HR, step[3], scales, options, PlotOnTop) + @unpack tMy, length_unit, Lx, Lz, xc, zc, ε̇II, coords, V, T, ϕ, σ1, ε̇1, Fab, height, group_phases, Δ = model + + tMy_string = @sprintf("%1.2lf", tMy) + ax1 = Axis(f[3, 1], title = L"C) $\tau_\textrm{II}$ at $t$ = %$(tMy_string) Ma - 1200 $\times$ 800 cells", xlabel = L"$x$ [%$(length_unit)]", ylabel = L"$y$ [%$(length_unit)]") + hm = heatmap!(ax1, xc./Lc, zc./Lc, log10.(τII), colormap = (cmap, options.α_heatmap), colorrange=(6.69, 8.69)) + AddCountourQuivers!(PlotOnTop, ax1, coords, V, T, ϕ, σ1, ε̇1, Fab, height, Lc, cm_y, group_phases, Δ*2) + colsize!(f.layout, 1, Aspect(1, Lx/Lz)) + Mak.Colorbar(f[3, 2], hm, label = L"$\tau_\textrm{II}$ [Pa]", width = 20, labelsize = ftsz, ticklabelsize = ftsz ) + Mak.colgap!(f.layout, 20) + xlims!(ax1, window.xmin, window.xmax) + ylims!(ax1, window.zmin, window.zmax) + + # lines!(ax1, xc./Lc, log10.(ε̇II[:, 700]) ) + # lines!(ax1, xc./Lc, atand.(ε̇1.z[:, 700]) ./ ε̇1.x[:, 700] ) + # xlims!(ax1, window.xmin, window.xmax) + # ylims!(ax1, -180, 180) + + + ############################################################## + + if options.printfig Print2Disk( f, path_HR, "ResolutionTest", step[3]) end #if printfig Print2Disk( f, path, string(field),ε̇BG) end + end + + display(f) +end + +function PrincipalStress(τxx, τzz, τxz, P) + σ1 = (x=zeros(size(τxx)), z=zeros(size(τxx)) ) + τxzc = 0.25*(τxz[1:end-1,1:end-1] .+ τxz[2:end-0,1:end-1] .+ τxz[1:end-1,2:end-0] .+ τxz[2:end-0,2:end-0]) + for i in eachindex(τxzc) + if P[i]>-1e-13 + σ = [-P[i]+τxx[i] τxzc[i]; τxzc[i] -P[i]+τzz[i]] + v = eigvecs(σ) + σ1.x[i] = v[1,1] + σ1.z[i] = v[2,1] + end + end + return σ1 +end + +function ExtractData( file_path, data_path) + data = h5open(file_path, "r") do file + read(file, data_path) + end + return data +end + +function Print2Disk( f, path, field, istep; res=4) + path1 = path*"/_$field/" + mkpath(path1) + @show name = path1*"$field"*@sprintf("%05d", istep)*".png" + save(name, f, px_per_unit = res) +end + +main() \ No newline at end of file diff --git a/JuliaVisualisation/_SpecificStudies/RiftingAnisotropy/ViscousAnisotropySimpleV.jl b/JuliaVisualisation/_SpecificStudies/RiftingAnisotropy/ViscousAnisotropySimpleV.jl index 85c9bfeb..25965095 100644 --- a/JuliaVisualisation/_SpecificStudies/RiftingAnisotropy/ViscousAnisotropySimpleV.jl +++ b/JuliaVisualisation/_SpecificStudies/RiftingAnisotropy/ViscousAnisotropySimpleV.jl @@ -3,17 +3,16 @@ import Statistics:mean function main_simple_ani_vis() # Material parameters - δ = 3.0 - a = sqrt(δ) + δ = 2.0 ηv = 1.0 # normal viscosity # Kinematics pure_shear = 1 - ε̇xxd = pure_shear*5.0 + ε̇xxd = pure_shear*.5 ε̇yyd = -ε̇xxd ε̇xyd = ((1.0 - pure_shear)*5.0) ε̇bg = sqrt(ε̇xxd^2 + ε̇xyd^2) # Arrays - θ = LinRange( 0.0, π/2, 51 ) + θ = LinRange( 0.0, π, 51 ) τii_cart = zero(θ) τii_cart_MD7 = zero(θ) τii_rot1 = zero(θ) @@ -105,18 +104,34 @@ function main_simple_ani_vis() end f = Figure(size = (500, 500), fontsize=18) - ax1 = Axis(f[1, 1], title="Stress invariant", xlabel="θ", ylabel="τᵢᵢ") - lines!(ax1, θ*180/π, τii_cart ) - scatter!(ax1, θ*180/π, τii_cart_MD7 ) - lines!(ax1, θ*180/π, τii_rot1 ) - lines!(ax1, θ*180/π, τii_rot2 ) - lines!(ax1, θ*180/π, τii_cart2 ) - # scatter!(ax1, θ*180/π, τii_cart, marker=:xcross, markersize=10 ) + ax1 = Axis(f[1, 1], title=L"$$A) Stress invariant", xlabel=L"$\theta$ [$^\circ$]", ylabel=L"$\tau_{II}$, $\tau_{II}'$ [-]") + lines!(ax1, θ*180/π, τii_cart, label=L"$\tau_{II}$ $(\delta$=2)$") + # scatter!(ax1, θ*180/π, τii_cart_MD7 ) + lines!(ax1, θ*180/π, τii_rot1, label=L"$\tau_{II}'$ $(\delta$=2)$") + text!(ax1, 50, 0.6; text=L"$\tau_{II} = \tau_{II}' \frac{\sqrt{d}}{\delta} $", fontsize=16) + + tii = 1.0*sqrt.((δ^2-1)*cos.(2*θ).^2 .+ 1)/δ + # scatter!(ax1, θ*180/π, tii, ) + - ax3 = Axis(f[2, 1:2], title="Flow enveloppe", xlabel=L"$τ_{xx}$", ylabel=L"$τ_{xy}$", aspect=DataAspect()) - lines!(ax3, τxx_cart, τxy_cart, marker=:xcross ) - scatter!(ax3, τxx_rot, τxy_rot, marker=:xcross ) - scatter!(ax3, τxx_cart2, τxy_cart2, marker=:xcross ) + scatter!(ax1, θ[1:5:end]*180/π, τii_rot2[1:5:end], label=L"$\tau_{II}$ $(\delta$=1)$" ) + # lines!(ax1, θ*180/π, τii_cart2 ) + # scatter!(ax1, θ*180/π, τii_cart, marker=:xcross, markersize=10 ) + axislegend(position=:rb) + + ax3 = Axis(f[2, 1], title=L"$$B) Flow enveloppe", xlabel=L"$\tau_{xx}'$ [-]", ylabel=L"$\tau_{xy}'$ [-]", aspect=DataAspect()) + lines!(ax3, LinRange(0,1,10), zeros(10), color=:black) + text!(ax3, 0.3, 0.05; text=L"$\tau_{xx}' = \tau_{II}'$", fontsize=16) + lines!(ax3, zeros(10), LinRange(0,0.5,10), color=:black) + text!(ax3, -0.05, 0.0; text=L"$\tau_{xy}' = \frac{\tau_{II}'}{\delta}$", rotation=π/2, fontsize=16) + + a = LinRange(0, 2π, 100) + txx = 1*cos.(a) + txy = 1*sin.(a)/δ + + # lines!(ax3, τxx_cart, τxy_cart, marker=:xcross ) + lines!(ax3, τxx_rot, τxy_rot, label=L"$\delta$=2$" ) + # scatter!(ax3, txx, txy, marker=:xcross ) display(f) diff --git a/JuliaVisualisation/_SpecificStudies/RiftingAnisotropy/ViscousAnisotropySimpleVEP.jl b/JuliaVisualisation/_SpecificStudies/RiftingAnisotropy/ViscousAnisotropySimpleVEP.jl index 37430555..51b440ee 100644 --- a/JuliaVisualisation/_SpecificStudies/RiftingAnisotropy/ViscousAnisotropySimpleVEP.jl +++ b/JuliaVisualisation/_SpecificStudies/RiftingAnisotropy/ViscousAnisotropySimpleVEP.jl @@ -11,8 +11,8 @@ function main_simple_ani_vis() ηe = G*Δt ηve = (1/ηv + 1/ηe)^-1 C = 5 - nt = 10 - + nt = 100 + # Kinematics pure_shear = 1 ε̇xxd = pure_shear*5.0