{
"cells": [
{
"outputs": [],
"cell_type": "code",
"source": [
"if isdefined(Main, :is_ci) #hide\n",
" IS_CI = Main.is_ci #hide\n",
"else #hide\n",
" IS_CI = false #hide\n",
"end #hide\n",
"nothing #hide"
],
"metadata": {},
"execution_count": 1
},
{
"cell_type": "markdown",
"source": [
"# Reactive surface\n",
"\n",
"\n",
"\n",
"*Figure 1*: Reactant concentration field of the Gray-Scott model on the unit sphere."
],
"metadata": {}
},
{
"cell_type": "markdown",
"source": [
"## Introduction\n",
"\n",
"This tutorial gives a quick tutorial on how to assemble and solve time-dependent problems\n",
"on embedded surfaces.\n",
"\n",
"For this showcase we use the well known Gray-Scott model, which is a well-known reaction-diffusion\n",
"system to study pattern formation. The strong form is given by\n",
"\n",
"$$\n",
" \\begin{aligned}\n",
" \\partial_t r_1 &= \\nabla \\cdot (D_1 \\nabla r_1) - r_1*r_2^2 + F *(1 - r_1) \\quad \\textbf{x} \\in \\Omega, \\\\\n",
" \\partial_t r_2 &= \\nabla \\cdot (D_2 \\nabla r_2) + r_1*r_2^2 - r_2*(F + k ) \\quad \\textbf{x} \\in \\Omega,\n",
" \\end{aligned}\n",
"$$\n",
"\n",
"where $r_1$ and $r_2$ are the reaction fields, $D_1$ and $D_2$ the diffusion tensors,\n",
"$k$ is the conversion rate, $F$ is the feed rate and $\\Omega$ the domain. Depending on the choice of\n",
"parameters a different pattern can be observed. Please also note that the domain does not have a\n",
"boundary. The corresponding weak form can be derived as usual.\n",
"\n",
"For simplicity we will solve the problem with the Lie-Troter-Godunov operator splitting technique with\n",
"the classical reaction-diffusion split. In this method we split our problem in two problems, i.e. a heat\n",
"problem and a pointwise reaction problem, and solve them alternatingly to advance in time.\n",
"\n",
"## Solver details\n",
"\n",
"The main idea for the Lie-Trotter-Godunov scheme is simple. We can write down the reaction diffusion\n",
"problem in an abstract way as\n",
"$$\n",
" \\partial_t \\mathbf{r} = \\mathcal{D}\\mathbf{r} + R(\\mathbf{r}) \\quad \\textbf{x} \\in \\Omega\n",
"$$\n",
"where $\\mathcal{D}$ is the diffusion operator and $R$ is the reaction operator. Notice that the right\n",
"hand side is just the sum of two operators. Now with our operator splitting scheme we can advance a\n",
"solution $\\mathbf{r}(t_1)$ to $\\mathbf{r}(t_2)$ by first solving a heat problem\n",
"$$\n",
" \\partial_t \\mathbf{r}^{\\mathrm{\\mathrm{A}}} = \\mathcal{D}\\mathbf{r}^{\\mathrm{A}} \\quad \\textbf{x} \\in \\Omega\n",
"$$\n",
"with $\\mathbf{r}^{\\mathrm{A}}(t_1) = \\mathbf{r}(t_1)$ on the time interval $t_1$ to $t_2$ and use\n",
"the solution as the initial condition to solve the reaction problem\n",
"$$\n",
" \\partial_t \\mathbf{r}^{\\mathrm{B}} = R(\\mathbf{r}^{\\mathrm{B}}) \\quad \\textbf{x} \\in \\Omega\n",
"$$\n",
"with $\\mathbf{r}^{\\mathrm{B}}(t_1) = \\mathbf{r}^{\\mathrm{A}}(t_2)$.\n",
"This way we obtain a solution approximation $\\mathbf{r}(t_2) \\approx \\mathbf{r}^{\\mathrm{B}}(t_2)$.\n",
"\n",
"> **Note**\n",
">\n",
"> The operator splitting itself is an approximation, so even if we solve the subproblems analytically\n",
"> we end up with having only a solution approximation. We also do not have a beginner friendly reference\n",
"> for the theory behind operator splitting and can only refer to the original papers for each method."
],
"metadata": {}
},
{
"cell_type": "markdown",
"source": [
"## Commented Program\n",
"\n",
"Now we solve the problem in Ferrite. What follows is a program spliced with comments.\n",
"\n",
"First we load Ferrite, and some other packages we need"
],
"metadata": {}
},
{
"outputs": [],
"cell_type": "code",
"source": [
"using Ferrite, FerriteGmsh\n",
"using BlockArrays, SparseArrays, LinearAlgebra, WriteVTK"
],
"metadata": {},
"execution_count": 2
},
{
"cell_type": "markdown",
"source": [
"### Assembly routines\n",
"Before we head into the assembly, we define a helper struct to control the dispatches."
],
"metadata": {}
},
{
"outputs": [],
"cell_type": "code",
"source": [
"struct GrayScottMaterial{T}\n",
" D₁::T\n",
" D₂::T\n",
" F::T\n",
" k::T\n",
"end;"
],
"metadata": {},
"execution_count": 3
},
{
"cell_type": "markdown",
"source": [
"The following assembly routines are written analogue to these found in previous tutorials."
],
"metadata": {}
},
{
"outputs": [],
"cell_type": "code",
"source": [
"function assemble_element_mass!(Me::Matrix, cellvalues::CellValues)\n",
" n_basefuncs = getnbasefunctions(cellvalues)\n",
" # The mass matrices between the reactions are not coupled, so we get a blocked-strided matrix.\n",
" num_reactants = 2\n",
" r₁range = 1:num_reactants:(num_reactants * n_basefuncs)\n",
" r₂range = 2:num_reactants:(num_reactants * n_basefuncs)\n",
" Me₁ = @view Me[r₁range, r₁range]\n",
" Me₂ = @view Me[r₂range, r₂range]\n",
" # Reset to 0\n",
" fill!(Me, 0)\n",
" # Loop over quadrature points\n",
" for q_point in 1:getnquadpoints(cellvalues)\n",
" # Get the quadrature weight\n",
" dΩ = getdetJdV(cellvalues, q_point)\n",
" # Loop over test shape functions\n",
" for i in 1:n_basefuncs\n",
" δuᵢ = shape_value(cellvalues, q_point, i)\n",
" # Loop over trial shape functions\n",
" for j in 1:n_basefuncs\n",
" δuⱼ = shape_value(cellvalues, q_point, j)\n",
" # Add contribution to Ke\n",
" Me₁[i, j] += (δuᵢ * δuⱼ) * dΩ\n",
" Me₂[i, j] += (δuᵢ * δuⱼ) * dΩ\n",
" end\n",
" end\n",
" end\n",
" return nothing\n",
"end\n",
"\n",
"function assemble_element_diffusion!(De::Matrix, cellvalues::CellValues, material::GrayScottMaterial)\n",
" n_basefuncs = getnbasefunctions(cellvalues)\n",
" D₁ = material.D₁\n",
" D₂ = material.D₂\n",
" # The diffusion between the reactions is not coupled, so we get a blocked-strided matrix.\n",
" num_reactants = 2\n",
" r₁range = 1:num_reactants:(num_reactants * n_basefuncs)\n",
" r₂range = 2:num_reactants:(num_reactants * n_basefuncs)\n",
" De₁ = @view De[r₁range, r₁range]\n",
" De₂ = @view De[r₂range, r₂range]\n",
" # Reset to 0\n",
" fill!(De, 0)\n",
" # Loop over quadrature points\n",
" for q_point in 1:getnquadpoints(cellvalues)\n",
" # Get the quadrature weight\n",
" dΩ = getdetJdV(cellvalues, q_point)\n",
" # Loop over test shape functions\n",
" for i in 1:n_basefuncs\n",
" ∇δuᵢ = shape_gradient(cellvalues, q_point, i)\n",
" # Loop over trial shape functions\n",
" for j in 1:n_basefuncs\n",
" ∇δuⱼ = shape_gradient(cellvalues, q_point, j)\n",
" # Add contribution to Ke\n",
" De₁[i, j] += D₁ * (∇δuᵢ ⋅ ∇δuⱼ) * dΩ\n",
" De₂[i, j] += D₂ * (∇δuᵢ ⋅ ∇δuⱼ) * dΩ\n",
" end\n",
" end\n",
" end\n",
" return nothing\n",
"end\n",
"\n",
"function assemble_matrices!(M::SparseMatrixCSC, D::SparseMatrixCSC, cellvalues::CellValues, dh::DofHandler, material::GrayScottMaterial)\n",
" n_basefuncs = getnbasefunctions(cellvalues)\n",
"\n",
" # Allocate the element stiffness matrix and element force vector\n",
" Me = zeros(2 * n_basefuncs, 2 * n_basefuncs)\n",
" De = zeros(2 * n_basefuncs, 2 * n_basefuncs)\n",
"\n",
" # Create an assembler\n",
" M_assembler = start_assemble(M)\n",
" D_assembler = start_assemble(D)\n",
" # Loop over all cels\n",
" for cell in CellIterator(dh)\n",
" # Reinitialize cellvalues for this cell\n",
" reinit!(cellvalues, cell)\n",
" # Compute element contribution\n",
" assemble_element_mass!(Me, cellvalues)\n",
" assemble!(M_assembler, celldofs(cell), Me)\n",
"\n",
" assemble_element_diffusion!(De, cellvalues, material)\n",
" assemble!(D_assembler, celldofs(cell), De)\n",
" end\n",
" return nothing\n",
"end;"
],
"metadata": {},
"execution_count": 4
},
{
"cell_type": "markdown",
"source": [
"### Initial condition setup\n",
"Time-dependent problems always need an initial condition from which the time evolution starts.\n",
"In this tutorial we set the concentration of reactant 1 to $1$ and the concentration of reactant\n",
"2 to $0$ for all nodal dof with associated coordinate $z \\leq 0.9$ on the sphere. Since the\n",
"simulation would be pretty boring with a steady-state initial condition, we introduce some\n",
"heterogeneity by setting the dofs associated to top part of the sphere (i.e. dofs with $z > 0.9$\n",
"to store the reactant concentrations of $0.5$ and $0.25$ for the reactants 1 and 2 respectively."
],
"metadata": {}
},
{
"outputs": [],
"cell_type": "code",
"source": [
"function setup_initial_conditions!(u₀::Vector, cellvalues::CellValues, dh::DofHandler)\n",
" u₀ .= ones(ndofs(dh))\n",
" u₀[2:2:end] .= 0.0\n",
"\n",
" n_basefuncs = getnbasefunctions(cellvalues)\n",
"\n",
" for cell in CellIterator(dh)\n",
" reinit!(cellvalues, cell)\n",
"\n",
" coords = getcoordinates(cell)\n",
" dofs = celldofs(cell)\n",
" uₑ = @view u₀[dofs]\n",
" rv₀ₑ = reshape(uₑ, (2, n_basefuncs))\n",
"\n",
" for i in 1:n_basefuncs\n",
" if coords[i][3] > 0.9\n",
" rv₀ₑ[1, i] = 0.5\n",
" rv₀ₑ[2, i] = 0.25\n",
" end\n",
" end\n",
" end\n",
"\n",
" u₀ .+= 0.01 * rand(ndofs(dh))\n",
" return\n",
"end;"
],
"metadata": {},
"execution_count": 5
},
{
"cell_type": "markdown",
"source": [
"### Mesh generation\n",
"In this section we define a routine to create a surface mesh with the help of FerriteGmsh.jl."
],
"metadata": {}
},
{
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": "create_embedded_sphere (generic function with 1 method)"
},
"metadata": {},
"execution_count": 6
}
],
"cell_type": "code",
"source": [
"function create_embedded_sphere(refinements)\n",
" gmsh.initialize()\n",
"\n",
" # Add a unit sphere in 3D space\n",
" gmsh.model.occ.addSphere(0.0, 0.0, 0.0, 1.0)\n",
" gmsh.model.occ.synchronize()\n",
"\n",
" # Generate nodes and surface elements only, hence we need to pass 2 into generate\n",
" gmsh.model.mesh.generate(2)\n",
"\n",
" # To get good solution quality refine the elements several times\n",
" for _ in 1:refinements\n",
" gmsh.model.mesh.refine()\n",
" end\n",
"\n",
" # Now we create a Ferrite grid out of it. Note that we also call toelements\n",
" # with our surface element dimension to obtain these.\n",
" nodes = tonodes()\n",
" elements, _ = toelements(2)\n",
" gmsh.finalize()\n",
" return Grid(elements, nodes)\n",
"end"
],
"metadata": {},
"execution_count": 6
},
{
"cell_type": "markdown",
"source": [
"### Simulation routines\n",
"Now we define a function to setup and solve the problem with given feed and conversion rates\n",
"$F$ and $k$, as well as the time step length and for how long we want to solve the model."
],
"metadata": {}
},
{
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Info : Meshing 1D...\n",
"Info : [ 40%] Meshing curve 2 (Circle)\n",
"Info : Done meshing 1D (Wall 0.000107039s, CPU 0.000107s)\n",
"Info : Meshing 2D...\n",
"Info : Meshing surface 1 (Sphere, Frontal-Delaunay)\n",
"Info : Done meshing 2D (Wall 0.0094344s, CPU 0.009429s)\n",
"Info : 162 nodes 332 elements\n",
"Info : Refining mesh...\n",
"Info : Meshing order 2 (curvilinear on)...\n",
"Info : [ 0%] Meshing curve 1 order 2\n",
"Info : [ 30%] Meshing curve 2 order 2\n",
"Info : [ 50%] Meshing curve 3 order 2\n",
"Info : [ 70%] Meshing surface 1 order 2\n",
"Info : [ 90%] Meshing volume 1 order 2\n",
"Info : Surface mesh: worst distortion = 0.968739 (0 elements in ]0, 0.2]); worst gamma = 0.533321\n",
"Info : Done meshing order 2 (Wall 0.00112749s, CPU 0.001126s)\n",
"Info : Done refining mesh (Wall 0.00146639s, CPU 0.001465s)\n",
"Info : Refining mesh...\n",
"Info : Meshing order 2 (curvilinear on)...\n",
"Info : [ 0%] Meshing curve 1 order 2\n",
"Info : [ 30%] Meshing curve 2 order 2\n",
"Info : [ 50%] Meshing curve 3 order 2\n",
"Info : [ 70%] Meshing surface 1 order 2\n",
"Info : [ 90%] Meshing volume 1 order 2\n",
"Info : Surface mesh: worst distortion = 0.991679 (0 elements in ]0, 0.2]); worst gamma = 0.514495\n",
"Info : Done meshing order 2 (Wall 0.00431626s, CPU 0.004315s)\n",
"Info : Done refining mesh (Wall 0.00559021s, CPU 0.00559s)\n",
"Info : Refining mesh...\n",
"Info : Meshing order 2 (curvilinear on)...\n",
"Info : [ 0%] Meshing curve 1 order 2\n",
"Info : [ 30%] Meshing curve 2 order 2\n",
"Info : [ 50%] Meshing curve 3 order 2\n",
"Info : [ 70%] Meshing surface 1 order 2\n",
"Info : [ 90%] Meshing volume 1 order 2\n",
"Info : Surface mesh: worst distortion = 0.997887 (0 elements in ]0, 0.2]); worst gamma = 0.509717\n",
"Info : Done meshing order 2 (Wall 0.0163413s, CPU 0.01634s)\n",
"Info : Done refining mesh (Wall 0.0213992s, CPU 0.021399s)\n"
]
}
],
"cell_type": "code",
"source": [
"function gray_scott_on_sphere(material::GrayScottMaterial, Δt::Real, T::Real, refinements::Integer)\n",
" # We start by setting up grid, dof handler and the matrices for the heat problem.\n",
" grid = create_embedded_sphere(refinements)\n",
"\n",
" # Next we are creating our element assembly helper for surface elements.\n",
" # The only change which we need to introduce here is to pass in a geometrical\n",
" # interpolation with the same dimension as the physical space into which our\n",
" # elements are embedded into, which is in this example 3.\n",
" ip = Lagrange{RefTriangle, 1}()\n",
" qr = QuadratureRule{RefTriangle}(2)\n",
" cellvalues = CellValues(qr, ip, ip^3)\n",
"\n",
" # We have two options to add the reactants to the dof handler, which will give us slightly\n",
" # different resulting dof distributions:\n",
" # A) We can add a scalar-valued interpolation for each reactant.\n",
" # B) We can add one vectorized interpolation whose dimension is the number of reactants\n",
" # number of reactants.\n",
" # In this tutorial we opt for B, because the dofs are distributed per cell entity -- or\n",
" # to be specific for this tutorial, we use an isoparametric concept such that the nodes\n",
" # of our grid and the nodes of our solution approximation coincide. This way a reaction\n",
" # we can create simply reshape the solution vector u to a matrix where the inner index\n",
" # corresponds to the index of the reactant. Note that we will still use the scalar\n",
" # interpolation for the assembly procedure.\n",
" dh = DofHandler(grid)\n",
" add!(dh, :reactants, ip^2)\n",
" close!(dh)\n",
"\n",
" # We can save some memory by telling the sparsity pattern that the matrices are not coupled.\n",
" M = allocate_matrix(dh; coupling = [true false; false true])\n",
" D = allocate_matrix(dh; coupling = [true false; false true])\n",
"\n",
" # Since the heat problem is linear and has no time dependent parameters, we precompute the\n",
" # decomposition of the system matrix to speed up the linear system solver.\n",
" assemble_matrices!(M, D, cellvalues, dh, material)\n",
" A = M + Δt .* D\n",
" cholA = cholesky(A)\n",
"\n",
" # Now we setup buffers for the time dependent solution and fill the initial condition.\n",
" uₜ = zeros(ndofs(dh))\n",
" uₜ₋₁ = ones(ndofs(dh))\n",
" setup_initial_conditions!(uₜ₋₁, cellvalues, dh)\n",
"\n",
" # And prepare output for visualization.\n",
" pvd = paraview_collection(\"reactive-surface\")\n",
" VTKGridFile(\"reactive-surface-0\", dh) do vtk\n",
" write_solution(vtk, dh, uₜ₋₁)\n",
" pvd[0.0] = vtk\n",
" end\n",
"\n",
" # This is now the main solve loop.\n",
" F = material.F\n",
" k = material.k\n",
" for (iₜ, t) in enumerate(Δt:Δt:T)\n",
" # First we solve the heat problem\n",
" uₜ .= cholA \\ (M * uₜ₋₁)\n",
"\n",
" # Then we solve the point-wise reaction problem with the solution of\n",
" # the heat problem as initial guess. 2 is the number of reactants.\n",
" num_individual_reaction_dofs = ndofs(dh) ÷ 2\n",
" rvₜ = reshape(uₜ, (2, num_individual_reaction_dofs))\n",
" for i in 1:num_individual_reaction_dofs\n",
" r₁ = rvₜ[1, i]\n",
" r₂ = rvₜ[2, i]\n",
" rvₜ[1, i] += Δt * (-r₁ * r₂^2 + F * (1 - r₁))\n",
" rvₜ[2, i] += Δt * (r₁ * r₂^2 - r₂ * (F + k))\n",
" end\n",
"\n",
" # The solution is then stored every 10th step to vtk files for\n",
" # later visualization purposes.\n",
" if (iₜ % 10) == 0\n",
" VTKGridFile(\"reactive-surface-$(iₜ)\", dh) do vtk\n",
" write_solution(vtk, dh, uₜ₋₁)\n",
" pvd[t] = vtk\n",
" end\n",
" end\n",
"\n",
" # Finally we totate the solution to initialize the next timestep.\n",
" uₜ₋₁ .= uₜ\n",
" end\n",
" vtk_save(pvd)\n",
" return\n",
"end\n",
"\n",
"# This parametrization gives the spot pattern shown in the gif above.\n",
"material = GrayScottMaterial(0.00016, 0.00008, 0.06, 0.062)\n",
" gray_scott_on_sphere(material, 10.0, 32000.0, 3)"
],
"metadata": {},
"execution_count": 7
},
{
"cell_type": "markdown",
"source": [
"---\n",
"\n",
"*This notebook was generated using [Literate.jl](https://github.com/fredrikekre/Literate.jl).*"
],
"metadata": {}
}
],
"nbformat_minor": 3,
"metadata": {
"language_info": {
"file_extension": ".jl",
"mimetype": "application/julia",
"name": "julia",
"version": "1.11.3"
},
"kernelspec": {
"name": "julia-1.11",
"display_name": "Julia 1.11.3",
"language": "julia"
}
},
"nbformat": 4
}