Reactive surface¶

function create_embedded_sphere(refinements)
    gmsh.initialize()

    # Add a unit sphere in 3D space
    gmsh.model.occ.addSphere(0.0, 0.0, 0.0, 1.0)
    gmsh.model.occ.synchronize()

    # Generate nodes and surface elements only, hence we need to pass 2 into generate
    gmsh.model.mesh.generate(2)

    # To get good solution quality refine the elements several times
    for _ in 1:refinements
        gmsh.model.mesh.refine()
    end

    # Now we create a Ferrite grid out of it. Note that we also call toelements
    # with our surface element dimension to obtain these.
    nodes = tonodes()
    elements, _ = toelements(2)
    gmsh.finalize()
    return Grid(elements, nodes)
end

function gray_scott_on_sphere(material::GrayScottMaterial, Δt::Real, T::Real, refinements::Integer)
    # We start by setting up grid, dof handler and the matrices for the heat problem.
    grid = create_embedded_sphere(refinements)

    # Next we are creating our element assembly helper for surface elements.
    # The only change which we need to introduce here is to pass in a geometrical
    # interpolation with the same dimension as the physical space into which our
    # elements are embedded into, which is in this example 3.
    ip = Lagrange{RefTriangle, 1}()
    qr = QuadratureRule{RefTriangle}(2)
    cellvalues = CellValues(qr, ip, ip^3)

    # We have two options to add the reactants to the dof handler, which will give us slightly
    # different resulting dof distributions:
    # A) We can add a scalar-valued interpolation for each reactant.
    # B) We can add one vectorized interpolation whose dimension is the number of reactants
    #    number of reactants.
    # In this tutorial we opt for B, because the dofs are distributed per cell entity -- or
    # to be specific for this tutorial, we use an isoparametric concept such that the nodes
    # of our grid and the nodes of our solution approximation coincide. This way a reaction
    # we can create simply reshape the solution vector u to a matrix where the inner index
    # corresponds to the index of the reactant. Note that we will still use the scalar
    # interpolation for the assembly procedure.
    dh = DofHandler(grid)
    add!(dh, :reactants, ip^2)
    close!(dh)

    # We can save some memory by telling the sparsity pattern that the matrices are not coupled.
    M = allocate_matrix(dh; coupling = [true false; false true])
    D = allocate_matrix(dh; coupling = [true false; false true])

    # Since the heat problem is linear and has no time dependent parameters, we precompute the
    # decomposition of the system matrix to speed up the linear system solver.
    assemble_matrices!(M, D, cellvalues, dh, material)
    A = M + Δt .* D
    cholA = cholesky(A)

    # Now we setup buffers for the time dependent solution and fill the initial condition.
    uₜ = zeros(ndofs(dh))
    uₜ₋₁ = ones(ndofs(dh))
    setup_initial_conditions!(uₜ₋₁, cellvalues, dh)

    # And prepare output for visualization.
    pvd = paraview_collection("reactive-surface")
    VTKGridFile("reactive-surface-0", dh) do vtk
        write_solution(vtk, dh, uₜ₋₁)
        pvd[0.0] = vtk
    end

    # This is now the main solve loop.
    F = material.F
    k = material.k
    for (iₜ, t) in enumerate(Δt:Δt:T)
        # First we solve the heat problem
        uₜ .= cholA \ (M * uₜ₋₁)

        # Then we solve the point-wise reaction problem with the solution of
        # the heat problem as initial guess. 2 is the number of reactants.
        num_individual_reaction_dofs = ndofs(dh) ÷ 2
        rvₜ = reshape(uₜ, (2, num_individual_reaction_dofs))
        for i in 1:num_individual_reaction_dofs
            r₁ = rvₜ[1, i]
            r₂ = rvₜ[2, i]
            rvₜ[1, i] += Δt * (-r₁ * r₂^2 + F * (1 - r₁))
            rvₜ[2, i] += Δt * (r₁ * r₂^2 - r₂ * (F + k))
        end

        # The solution is then stored every 10th step to vtk files for
        # later visualization purposes.
        if (iₜ % 10) == 0
            VTKGridFile("reactive-surface-$(iₜ)", dh) do vtk
                write_solution(vtk, dh, uₜ₋₁)
                pvd[t] = vtk
            end
        end

        # Finally we totate the solution to initialize the next timestep.
        uₜ₋₁ .= uₜ
    end
    vtk_save(pvd)
    return
end

# This parametrization gives the spot pattern shown in the gif above.
material = GrayScottMaterial(0.00016, 0.00008, 0.06, 0.062)
    gray_scott_on_sphere(material, 10.0, 32000.0, 3)