Cahn-Hillard with Gridap, Julia

Implementation

using Gridap
using Random
using Printf
using LinearAlgebra

function run_ch()
    kappa        = 1e-3
    dt           = 1e-4
    n_steps      = 500
    plot_interval = 1

    domain    = (0.0, 1.0, 0.0, 1.0)
    partition = (2^8, 2^8)
    model = CartesianDiscreteModel(domain, partition)

    order = 1
    reffe = ReferenceFE(lagrangian, Float64, order)

    V = TestFESpace(model, reffe; conformity = :H1)
    U = TrialFESpace(V)

    Y = MultiFieldFESpace([V, V])  # test : (v, q)
    X = MultiFieldFESpace([U, U])  # trial: (c, μ)

    Ω  = Triangulation(model)
    deg = 2 * order
    dΩ = Measure(Ω, deg)

    # ============================
    # 初期条件
    # ============================
    Random.seed!(0)
    N = num_free_dofs(U)

    mean_c0 = 0.0
    c_vec = mean_c0 .+ 0.01 .* randn(N)
    c_vec .-= mean(c_vec) .- mean_c0

    mkpath("ch_output")
    vtufiles = String[]
    times = Float64[]

    c0_h = FEFunction(U, c_vec)
    writevtk(Ω, "ch_output/ch_0000", cellfields = ["c" => c0_h])
    push!(vtufiles, "ch_0000.vtu")
    push!(times, 0.0)

    a((c, μ), (v, q)) = ∫(
          c * v
        + dt * (∇(μ) ⋅ ∇(v))
        - kappa * (∇(c) ⋅ ∇(q))
        + μ * q
    ) * dΩ
    A = assemble_matrix(a, X, Y)
    factA = lu(A)

    # ============================
    # Time Evolution Loop
    # ============================
    for step in 1:n_steps

        c_old_h = FEFunction(U, c_vec)

        l((v, q)) = ∫(
              c_old_h * v
            + (c_old_h * c_old_h * c_old_h - c_old_h) * q
        ) * dΩ

        b = assemble_vector(l, Y)
        x_vec = factA \ b
        xh = FEFunction(X, x_vec)

        c_h, μ_h = xh

        c_vec = get_free_dof_values(c_h)

        if (step % plot_interval == 0) || (step == n_steps)
            fname = @sprintf("ch_output/ch_%04d", step)
            writevtk(Ω, fname, cellfields = ["c" => c_h])
            push!(vtufiles, @sprintf("ch_%04d.vtu", step))
            push!(times, step * dt)
            @info "t-step = $step : wrote $fname.vtu"
        end
    end
    write_pvd_collection("ch_output/ch_series.pvd", vtufiles, times)
end

function write_pvd_collection(path, files, times)
    open(path, "w") do io
        println(io, """<?xml version="1.0"?>""")
        println(io, """<VTKFile type="Collection" version="0.1" byte_order="LittleEndian">""")
        println(io, "  <Collection>")
        for (f, t) in zip(files, times)
            println(io, @sprintf("""    <DataSet timestep="%.6f" group="" part="0" file="%s"/>""", t, f))
        end
        println(io, "  </Collection>")
        println(io, "</VTKFile>")
    end
end

run_ch()