# Mean Curvature Flow (revisited)¶

As discussed before we can simulate the shrinking of a sphere under mean curvature flow using a finite element approach based on the following time discrete approximation: \begin{align} \int_{\Gamma^n} \big( U^{n+1} - {\rm id}\big) \cdot \varphi + \tau \int_{\Gamma^n} \big( \theta\nabla_{\Gamma^n} U^{n+1} + (1-\theta) I \big) \colon\nabla_{\Gamma^n}\varphi =0~. \end{align} Here $$U^n$$ parametrizes $$\Gamma(t^{n+1})$$ over $$\Gamma^n:=\Gamma(t^{n})$$, $$I$$ is the identity matrix, $$\tau$$ is the time step and $$\theta\in[0,1]$$ is a discretization parameter.

If the initial surface $$\Gamma^0$$ is a sphere of radius $$R_0$$, the surface remains sphere and we have an exact formula for the evolution of the radius of the surface

$R(t) = \sqrt{R_0^2 - 4t}.$

To compare the accuracy of the surface approximation we compute an average radius of the discrete surface in each time step $$t^n$$ using

$R_h^n = \frac{ \int_{\Gamma^n} |x| }{ |\Gamma^n| }.$

Computing $$R_h^n$$ requires a grid traversal and a number of calls to interface methods on each element. Doing this on the Python side has a potential performance impact which we investigate here by comparing a pure python implementation with a hybrid approach where computing $$R_h^n$$ is implemented in C++ using the dune.generator.algorithm functionality.

:

import math, time
import pickle
import numpy as np

from ufl import *
import dune.ufl
from dune.generator import algorithm
import dune.geometry as geometry
import dune.fem as fem
from mcf_cmp_plot import plot

# polynomial order of surface approximation
order = 2

R0 = 2.

# end time
endTime = 0.1


Main function for calculating the mean curvature flow of a given surface. If first argument is True the radius of the computed surface is computed using an algorithm implemented in C++ otherwise the computation is done in Python.

:

def calcRadius(surface):
R,vol = 0, 0
for e in surface.elements:
for p in rule:
geo = e.geometry
weight = geo.volume * p.weight
R   += geo.toGlobal(p.position).two_norm * weight
vol += weight
return R/vol



Timings for a number of different grid refinements is dumped to disk

:

from dune.fem.view import geometryGridView as geoGridView
from dune.fem.space import lagrange as solutionSpace
from dune.fem.scheme import galerkin as solutionScheme
def calculate(use_cpp, gridView):
# space on Gamma_0 to describe position of Gamma(t)
space = solutionSpace(gridView, dimRange=gridView.dimWorld, order=order)
u = TrialFunction(space)
v = TestFunction(space)
x = SpatialCoordinate(space.cell())
positions = space.interpolate(x, name="position")

# space for discrete solution on Gamma(t)
surface = geoGridView(positions)
space = solutionSpace(surface, dimRange=surface.dimWorld, order=order)
solution  = space.interpolate(x, name="solution")

# set up model using theta scheme
theta = 0.5   # Crank-Nicholson

I = Identity(3)
dt = dune.ufl.Constant(0,"dt")

a = (inner(u - x, v) + dt * inner(theta*grad(u)
+ (1 - theta)*I, grad(v))) * dx

scheme = solutionScheme(a == 0, space, solver="cg")

Rexact = lambda t: math.sqrt(R0*R0 - 4.*t)

scheme.model.dt = 0.02

numberOfLoops = 3
times = np.zeros(numberOfLoops)
errors = np.zeros(numberOfLoops)
totalIterations = np.zeros(numberOfLoops, np.dtype(np.uint32))
gridSizes = np.zeros(numberOfLoops, np.dtype(np.uint32))
for i in range(numberOfLoops):
positions.interpolate(x * (R0/sqrt(dot(x,x))))
solution.interpolate(x)
t = 0.
iterations = 0
start = time.time()
while t < endTime:
info = scheme.solve(target=solution)
# move the surface
positions.dofVector.assign(solution.dofVector)
# store some information about the solution process
iterations += int( info["linear_iterations"] )
t          += scheme.model.dt
print("time used:", time.time() - start)
times[i] = time.time() - start
errors[i] = error
totalIterations[i] = iterations
gridSizes[i] = gridView.size(2)
if i < numberOfLoops - 1:
gridView.hierarchicalGrid.globalRefine(1)
scheme.model.dt /= 2
eocs = np.log(errors[0:][:numberOfLoops-1] / errors[1:]) / math.log(math.sqrt(2))
try:
import pandas as pd
keys = {'size': gridSizes, 'error': errors, "eoc": np.insert(eocs, 0, None), 'iterations': totalIterations}
table = pd.DataFrame(keys, index=range(numberOfLoops),columns=['size', 'error', 'eoc', 'iterations'])
print(table)
except ImportError:
print("pandas could not be used to show table with results")
pass
return gridSizes, times


Compute the mean curvature flow evolution of a spherical surface. Compare computational time of a pure Python implementation and using a C++ algorithm to compute the radius of the surface for verifying the algorithm.

:

# set up reference domain Gamma_0
results = []
from dune.alugrid import aluConformGrid as leafGridView
gridView = leafGridView("sphere.dgf", dimgrid=2, dimworld=3)
results += [calculate(True, gridView)]

gridView = leafGridView("sphere.dgf", dimgrid=2, dimworld=3)
results += [calculate(False, gridView)]

time used: 0.3605988025665283
time used: 2.0574331283569336
time used: 8.871084690093994
size     error       eoc  iterations
0   318  0.001060       NaN          94
1   766  0.000605  1.619154         271
2  1745  0.000275  2.275142         594
time used: 0.4666929244995117
time used: 2.5975165367126465
time used: 11.073150396347046
size     error       eoc  iterations
0   318  0.001060       NaN          94
1   766  0.000605  1.619154         271
2  1745  0.000275  2.275142         594


Compare the hybrid and pure Python versions

:

plot(results,results) 