Merge variational BCs (#14)

* Add new scripts

* modify

* Add data preparation script

* WIP

* WIP

* WIP

* Remove tmp script

* Remove deps

* WIP

* Remove commented code

* Remove unused file

* Remove old files

* Add files

* fix function name

* Bugfix!

* Fix style

* Limit to only a single "phase"

* Fix residual Stokes

* Fix plot axes

* Fix volume fraction symmetry

* Work with new TinyKernels and union for levelset.

* Add nonlinear viscosity (needs cutoff)

* Remove comments

* WIP support for field extension

* Add new dependencies

* Implement new `TK` API

* Add fields API

* Add logging API

* Change scripts

* Remove unneeded argument

* Mixup duplicate

* Add 2D TM draft

* Fix code structure and figure resolution

* Add thermodynamics

* WIP

* Add 3D TM

* Fix typo

* Fix typos

* WIP 3D TM

* WIP 3D TM

* Ignore pictures and add script

* 3D TM WIP

* Update test cases

* Add "dream script"

---------

Co-authored-by: Ludovic Raess <[email protected]>

.gitignore

@@ -27,8 +27,18 @@ docs/site/
 Manifest.toml
 LocalPreferences.toml
 GeoData/Manifest.toml
+
+# OSX-specific files
 .DS_Store
+
+# VSCode-specific files
 .vscode
+
+# Project-specific files
 out_visu/
 gitignore/
 data/
+*.h5
+*.xdmf3
+*.png
+*.mp4

Project.toml

@@ -4,23 +4,29 @@ authors = ["Ludovic Raess <[email protected]>, Ivan Utkin and contributors"
 version = "0.1.0"
 
 [deps]
-AMDGPU = "21141c5a-9bdb-4563-92ae-f87d6854732e"
 BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
 CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba"
-CellArrays = "d35fcfd7-7af4-4c67-b1aa-d78070614af4"
+CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0"
 Colors = "5ae59095-9a9b-59fe-a467-6f913c188581"
 ElasticArrays = "fdbdab4c-e67f-52f5-8c3f-e7b388dad3d4"
+FileIO = "5789e2e9-d7fb-5bc7-8068-2c6fae9b9549"
+ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
+GLMakie = "e9467ef8-e4e7-5192-8a1a-b1aee30e663a"
 GPUCompiler = "61eb1bfa-7361-4325-ad38-22787b887f55"
 GeometryBasics = "5c1252a2-5f33-56bf-86c9-59e7332b4326"
 HDF5 = "f67ccb44-e63f-5c2f-98bd-6dc0ccc4ba2f"
 ImplicitGlobalGrid = "4d7a3746-15be-11ea-1130-334b0c4f5fa0"
 Interpolations = "a98d9a8b-a2ab-59e6-89dd-64a1c18fca59"
+JLD2 = "033835bb-8acc-5ee8-8aae-3f567f8a3819"
 LightXML = "9c8b4983-aa76-5018-a973-4c85ecc9e179"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+Logging = "56ddb016-857b-54e1-b83d-db4d58db5568"
 MPI = "da04e1cc-30fd-572f-bb4f-1f8673147195"
 MPIPreferences = "3da0fdf6-3ccc-4f1b-acd9-58baa6c99267"
+Metal = "dde4c033-4e86-420c-a63e-0dd931031962"
 NetCDF = "30363a11-5582-574a-97bb-aa9a979735b9"
-ParallelStencil = "94395366-693c-11ea-3b26-d9b7aac5d958"
-Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
+Preferences = "21216c6a-2e73-6563-6e65-726566657250"
 Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
+TinyKernels = "f7cbc414-f748-44bf-86e6-e44e9a55e39d"
+UnicodePlots = "b8865327-cd53-5732-bb35-84acbb429228"

scripts2D_variational/app_inclusion2D.jl

@@ -0,0 +1,164 @@
+using FastIce
+using TinyKernels
+using CairoMakie
+using ElasticArrays
+using Printf
+
+include("bcs.jl")
+include("helpers_tmp.jl")
+include("level_sets.jl")
+include("stokes.jl")
+include("volume_fractions.jl")
+
+@views av1(A) = 0.5 .* (A[1:end-1] .+ A[2:end])
+@views inn_x(A) = A[2:end-1,:]
+@views inn_y(A) = A[:,2:end-1]
+@views inn(A) = A[2:end-1,2:end-1]
+
+@views function runsim(::Type{DAT}; nx=127) where {DAT}
+    # physics
+    # lx, ly   = 2.0, 1.0
+    lx, ly   = 1.0, 1.0
+    ox, oy   = -0.5lx, -0.5ly
+    # xb1, yb1 = ox + 0.5lx, oy + 0.0ly
+    # xb2, yb2 = ox + 0.5lx, oy + 3.0ly
+    xb1, yb1 = ox + 0.5lx, oy + 0.5ly
+    rinc     = 0.14ly
+    # rair     = 2.3ly
+    ηs0      = 1.0
+    ebg      = 1.0
+    ρg0      = 0.0
+    α        = 0.0
+    npow     = 1
+    # numerics
+    ny       = ceil(Int, (nx + 1) * ly / lx) - 1
+    maxiter  = 40nx
+    ncheck   = 2nx
+    ϵtol     = (1e-6, 1e-6, 1e-6)
+    nt       = 1
+    χ        = 1.0
+    # preprocessing
+    dx, dy   = lx / nx, ly / ny
+    xv, yv   = LinRange(ox, ox + lx, nx + 1), LinRange(oy, oy + ly, ny + 1)
+    xc, yc   = av1(xv), av1(yv)
+    mc1      = to_device(make_marker_chain_circle(Point(xb1, yb1), rinc, min(dx, dy)))
+    # mc2      = to_device(make_marker_chain_circle(Point(xb2, yb2), rair, min(dx, dy)))
+    ρg       = (x=ρg0 .* sin(α), y=ρg0 .* cos(α))
+    mpow     = -(1 - 1 / npow)
+    # PT parameters
+    r        = 0.7
+    re_mech  = 6π
+    lτ       = min(lx, ly)
+    vdτ      = min(dx, dy) / sqrt(2.1)
+    θ_dτ     = lτ * (r + 4 / 3) / (re_mech * vdτ)
+    nudτ     = vdτ * lτ / re_mech
+    dτ_r     = 1.0 / (θ_dτ + 1.0)
+    # level set
+    Ψ  = (
+        not_air = field_array(DAT, nx + 1, ny + 1),  # liquid
+    )
+    wt = (
+        not_solid = volfrac_field(DAT, nx, ny), # fluid
+        not_air   = volfrac_field(DAT, nx, ny), # liquid
+    )
+    # mechanics
+    Pr = scalar_field(DAT, nx, ny)
+    τ  = tensor_field(DAT, nx, ny)
+    ε  = tensor_field(DAT, nx, ny)
+    V  = vector_field(DAT, nx, ny)
+    ηs = scalar_field(DAT, nx, ny)
+    τII = scalar_field(DAT, nx, ny)
+    εII = scalar_field(DAT, nx, ny)
+    # residuals
+    Res = (
+        Pr = scalar_field(DAT, nx    , ny    ),
+        V  = vector_field(DAT, nx - 2, ny - 2),
+    )
+    # visualisation
+    Vmag = field_array(DAT, nx - 2, ny - 2)
+    Ψav = (
+        not_solid = field_array(DAT, nx - 2, ny - 2),
+        not_air   = field_array(DAT, nx - 2, ny - 2),
+    )
+    # initial and boundary conditions
+    @info "computing the level set for the inclusion"
+    for comp in eachindex(Ψ) fill!(Ψ[comp], 1.0) end
+    init!(Pr, τ, V, ηs, ebg, ηs0, xv, yv)
+    fill!(τII, 0.0)
+    fill!(εII, 0.0)
+    Ψ.not_air .= Inf # needs init now
+    compute_levelset!(Ψ.not_air, xv, yv, mc1)
+    # compute_levelset!(Ψ.not_air, xv, yv, mc2)
+    Ψ.not_air .= .-Ψ.not_air
+
+    @info "computing volume fractions from level sets"
+    compute_volume_fractions_from_level_set!(wt.not_air, Ψ.not_air, dx, dy)
+    for comp in eachindex(wt.not_solid) fill!(wt.not_solid[comp], 1.0) end
+
+    update_vis!(Vmag, Ψav, V, Ψ)
+    # convergence history
+    iter_evo = Float64[]
+    errs_evo = ElasticArray{Float64}(undef, length(ϵtol), 0)
+    # figures
+    fig = Figure(resolution=(2500, 1600), fontsize=32)
+    ax = (
+        Pr  =Axis(fig[1, 1][1, 1]; aspect=DataAspect(), title="p"),
+        τII =Axis(fig[1, 2][1, 1]; aspect=DataAspect(), title="τII"),
+        Vmag=Axis(fig[2, 1][1, 1]; aspect=DataAspect(), title="|v|"),
+        εII =Axis(fig[2, 2][1, 1]; aspect=DataAspect(), title="εII"),
+        wt  =Axis(fig[1, 3][1, 1]; aspect=DataAspect(), title="Volume fraction"),
+        errs=Axis(fig[2, 3]      ; yscale=log10, title="Convergence", xlabel="#iter/ny", ylabel="error"),
+    )
+    plt = (
+        fields=(
+            Pr  =heatmap!(ax.Pr  , xc, yc, to_host(Pr  ); colormap=:turbo),
+            τII =heatmap!(ax.τII , xc, yc, to_host(τII ); colormap=:turbo),
+            Vmag=heatmap!(ax.Vmag, xc, yc, to_host(Vmag); colormap=:turbo),
+            εII =heatmap!(ax.εII , xc, yc, to_host(εII ); colormap=:turbo),
+            wt  =heatmap!(ax.wt  , xc, yc, to_host(wt.not_air.c); colormap=Reverse(:grays)),
+        ),
+        errs=[scatterlines!(ax.errs, Point2.(iter_evo, errs_evo[ir, :])) for ir in eachindex(ϵtol)],
+    )
+    Colorbar(fig[1, 1][1, 2], plt.fields.Pr  )
+    Colorbar(fig[1, 2][1, 2], plt.fields.τII )
+    Colorbar(fig[2, 1][1, 2], plt.fields.Vmag)
+    Colorbar(fig[2, 2][1, 2], plt.fields.εII )
+    Colorbar(fig[1, 3][1, 2], plt.fields.wt  )
+    display(fig)
+
+    @info "running simulation 🚀"
+    for it in 1:nt
+        @printf "it # %d\n" it
+        # iteration loop
+        empty!(iter_evo); resize!(errs_evo, length(ϵtol), 0)
+        iter = 0; errs = 2.0 .* ϵtol
+        while any(errs .>= ϵtol) && (iter += 1) <= maxiter
+            update_σ!(Pr, ε, τ, V, ηs, wt, r, θ_dτ, dτ_r, dx, dy)
+            compute_invariants!(εII, τII, ε, τ, ηs, χ, mpow)
+            update_V!(V, Pr, τ, ηs, wt, nudτ, ρg, dx, dy)
+            if iter % ncheck == 0
+                compute_residual!(Res, Pr, V, τ, wt, ρg, dx, dy)
+                errs = (maximum(abs.(Res.V.x)), maximum(abs.(Res.V.y)), maximum(abs.(Res.Pr)))
+                @printf "  iter/nx # %2.1f, errs: [ Vx = %1.3e, Vy = %1.3e, Pr = %1.3e ]\n" iter / nx errs...
+                push!(iter_evo, iter / nx); append!(errs_evo, errs)
+                # visu
+                for ir in eachindex(plt.errs)
+                    plt.errs[ir][1] = Point2.(iter_evo, errs_evo[ir, :])
+                end
+                autolimits!(ax.errs)
+                update_vis!(Vmag, Ψav, V, Ψ)
+                plt.fields[1][3] = to_host(to_host(Pr))
+                plt.fields[2][3] = to_host(to_host(τII))
+                plt.fields[3][3] = to_host(to_host(Vmag))
+                # plt.fields[4][3] = to_host(to_host(εII))
+                plt.fields[4][3] = to_host(to_host(log10.(ηs)))
+                plt.fields[5][3] = to_host(to_host(wt.not_air.c))
+                # plt.fields[4][3] = to_host(to_host(Ψ.not_air))
+                display(fig)
+            end
+        end
+    end
+    return
+end
+
+runsim(Float64, nx=127)

scripts2D_variational/bc_kernels.jl

@@ -0,0 +1,43 @@
+@tiny function _kernel_bc_x_dirichlet!(val,arrays...)
+    iy, = @indices
+    for A in arrays
+        if iy ∈ axes(A,2)
+            @inbounds A[1  ,iy] = val
+            @inbounds A[end,iy] = val
+        end
+    end
+    return
+end
+
+@tiny function _kernel_bc_y_dirichlet!(val, arrays...)
+    ix, = @indices
+    for A in arrays
+        if ix ∈ axes(A,1)
+            @inbounds A[ix,1  ] = val
+            @inbounds A[ix,end] = val
+        end
+    end
+    return
+end
+
+@tiny function _kernel_bc_x_neumann!(val, arrays...)
+    iy, = @indices
+    for A in arrays
+        if iy ∈ axes(A,2)
+            @inbounds A[1  ,iy] = A[2    ,iy] + val
+            @inbounds A[end,iy] = A[end-1,iy] + val
+        end
+    end
+    return
+end
+
+@tiny function _kernel_bc_y_neumann!(val, arrays...)
+    ix, = @indices
+    for A in arrays
+        if ix ∈ axes(A,1)
+            @inbounds A[ix,1  ] = A[ix,2    ] + val
+            @inbounds A[ix,end] = A[ix,end-1] + val
+        end
+    end
+    return
+end

scripts2D_variational/bcs.jl

@@ -0,0 +1,33 @@
+include("bc_kernels.jl")
+
+const _bc_x_dirichlet! = _kernel_bc_x_dirichlet!(get_device())
+const _bc_y_dirichlet! = _kernel_bc_y_dirichlet!(get_device())
+
+const _bc_x_neumann! = _kernel_bc_x_neumann!(get_device())
+const _bc_y_neumann! = _kernel_bc_y_neumann!(get_device())
+
+for fname in (:bx_x_dirichlet!, :bc_x_neumann!)
+    @eval begin
+        function $fname(val, arrays...)
+            ax = axes(arrays[1], 2)
+            for A in arrays[2:end]
+                ax = union.(ax, axes(A, 2))
+            end
+            wait($(Symbol(:_, fname))(val, arrays...; ndrange=ax))
+            return
+        end
+    end
+end
+
+for fname in (:bx_y_dirichlet!, :bc_y_neumann!)
+    @eval begin
+        function $fname(val, arrays...)
+            ax = axes(arrays[1], 1)
+            for A in arrays[2:end]
+                ax = union.(ax, axes(A, 1))
+            end
+            wait($(Symbol(:_, fname))(val, arrays...; ndrange=ax))
+            return
+        end
+    end
+end

scripts2D_variational/geometry.jl

@@ -0,0 +1,24 @@
+using LinearAlgebra, GeometryBasics
+
+function make_marker_chain_circle(rc, rad, hmax)
+    np = ceil(Int, 2π * rad / hmax)
+    return [rc + rad .* Point2(reverse(sincospi(2 * (i - 1) / np))...) for i in 1:np]
+end
+
+function signed_distance(p::Point2{T}, poly::AbstractVector{Point2{T}}) where {T}
+    d = dot(p - poly[1], p - poly[1])
+    s = 1.0
+    j = length(poly)
+    for i in eachindex(poly)
+        e = poly[j] - poly[i]
+        w = p - poly[i]
+        b = w - e .* clamp(dot(w, e) / dot(e, e), 0.0, 1.0)
+        d = min(d, dot(b, b))
+        c = p[2] >= poly[i][2], p[2] < poly[j][2], e[1] * w[2] > e[2] * w[1]
+        if all(c) || all(.!c)
+            s = -s
+        end
+        j = i
+    end
+    return s * sqrt(d)
+end

scripts2D_variational/helpers_tmp.jl

@@ -0,0 +1,67 @@
+@inline scalar_field(::Type{T}, nx, ny) where {T}  = field_array(T, nx, ny)
+@inline vector_field(::Type{T}, nx, ny) where {T}  =  (x = field_array(T, nx + 1, ny    ),
+                                                       y = field_array(T, nx    , ny + 1))
+@inline tensor_field(::Type{T}, nx, ny) where {T}  = (xx = field_array(T, nx    , ny    ),
+                                                      yy = field_array(T, nx    , ny    ),
+                                                      xy = field_array(T, nx - 1, ny - 1))
+@inline volfrac_field(::Type{T}, nx, ny) where {T} = (c  = field_array(T, nx    , ny    ),
+                                                      x  = field_array(T, nx + 1, ny    ),
+                                                      y  = field_array(T, nx    , ny + 1),
+                                                      xy = field_array(T, nx - 1, ny - 1))
+
+@tiny function _kernel_init!(Pr, τ, V, ηs, ebg, ηs0, xv, yv)
+    ix, iy = @indices()
+    @inbounds if ix ∈ axes(Pr, 1) && iy ∈ axes(Pr, 2)
+        Pr[ix, iy]   = 0.0
+        τ.xx[ix, iy] = 0.0
+        τ.yy[ix, iy] = 0.0
+        ηs[ix, iy]   = ηs0
+    end
+    if ix ∈ axes(τ.xy, 1) && iy ∈ axes(τ.xy, 2)
+        @inbounds τ.xy[ix, iy] = 0.0
+    end
+    if ix ∈ axes(V.x, 1) && iy ∈ axes(V.x, 2)
+        @inbounds V.x[ix, iy] = -xv[ix] * ebg
+    end
+    if ix ∈ axes(V.y, 1) && iy ∈ axes(V.y, 2)
+        @inbounds V.y[ix, iy] = yv[iy] * ebg
+    end
+end
+
+@tiny function _kernel_update_vis_fields!(Vmag, Ψav, V, Ψ)
+    ix, iy = @indices
+    @inline isin(A) = checkbounds(Bool, A, ix, iy)
+    @inbounds if isin(Ψ.not_air)
+        pav = 0.0
+        for idy = 1:2, idx = 1:2
+            pav += Ψ.not_air[ix+idx, iy+idy]
+        end
+        Ψav.not_air[ix, iy] = pav / 8
+    end
+    # @inbounds if isin(Ψ.not_solid)
+    #     pav = 0.0
+    #     for idy = 1:2, idx = 1:2
+    #         pav += Ψ.not_solid[ix+idx, iy+idy]
+    #     end
+    #     Ψav.not_solid[ix, iy] = pav / 8
+    # end
+    @inbounds if isin(Vmag)
+        vxc = 0.5 * (V.x[ix+1, iy+1] + V.x[ix+2, iy+1])
+        vyc = 0.5 * (V.y[ix+1, iy+1] + V.y[ix+1, iy+2])
+        Vmag[ix, iy] = sqrt(vxc^2 + vyc^2)
+    end
+    return
+end
+
+const _init! = _kernel_init!(get_device())
+const _update_vis! = _kernel_update_vis_fields!(get_device())
+
+function init!(Pr, τ, V, ηs, ebg, ηs0, xv, yv)
+    wait(_init!(Pr, τ, V, ηs, ebg, ηs0, xv, yv; ndrange=size(Pr) .+ 1))
+    return
+end
+
+function update_vis!(Vmag, Ψav, V, Ψ)
+    wait(_update_vis!(Vmag, Ψav, V, Ψ; ndrange=axes(Vmag)))
+    return
+end