Pathological Asymmetry in Projected Function

[zaccandels]

First thing, I’m using FEniCS 2019.1.0. In the MWE pasted below, I initialize a function phi_n by projecting the right-hand side of the equation

\phi(x,y) = 0.5 - 0.5 \tanh\biggl( 2 \frac{ \sqrt{ (x - x_c)^2 + (y- y_c)^2 ) } -R }{\epsilon} \biggr)

onto a function space, V. Then I initialize a function mu_n by determining the representation of

\mu = 4\beta\phi(\phi - 1)(\phi - 1/2) - \kappa\nabla^2\phi

in the space V - that is to say, I use L2 projection to get a finite-dimensional approximation to the above.

Now, mu_n should be perfectly symmetric. Given two points, (x_0 - dx, y_0) and (x_0 + dx, y_0), I should see that \mu(x_0 - dx, y_0) - \mu(x_0 + dx, y_0)= \mathcal{O}(1\times10^{-16}). However, this is not what I see at all.

After solving the requisite linear system for mu_n, I find that (for example) mu_n(2.5, 0.1) = 0.002777… and mu_n(7.5, 0.1) = 0.003796…. What accounts for this discrepancy?

I’ve tried different solvers, changing the mesh resolution, changing the degree of the function space… the list goes on. Nothing dispenses with this disparity in mu_n.

import fenics as fe
import os
import numpy as np

initDropDiam = 5
L_x, L_y = 2*initDropDiam, 2*initDropDiam
nx, ny = 200, 200
h = min(L_x/nx, L_y/ny)

# Relevant constants
sigma = 0.005 

Cn = 0.05

eps = Cn * initDropDiam

beta = 12*sigma/eps

kappa = 3*sigma*eps/2 


# Droplet center
center_init_x, center_init_y = L_x/2, 0

# Set up domain.

mesh = fe.RectangleMesh(fe.Point(0, 0), fe.Point(L_x, L_y), nx, ny)

# Set periodic boundary conditions at left and right endpoints


class PeriodicBoundaryX(fe.SubDomain):
    def inside(self, point, on_boundary):
        return fe.near(point[0], 0.0) and on_boundary

    def map(self, right_bdy, left_bdy):
        # Map left boundary to the right
        left_bdy[0] = right_bdy[0] - L_x
        left_bdy[1] = right_bdy[1]


pbc = PeriodicBoundaryX()


V = fe.FunctionSpace(mesh, "P", 1, constrained_domain=pbc)


# Define trial and test functions, as well as
# finite element functions at previous timesteps

f_trial = fe.TrialFunction(V)
phi_n = fe.Function(V)
mu_n = fe.Function(V)

v = fe.TestFunction(V)

    
# Initialize \phi
phi_init_expr = fe.Expression(
    "0.5 - 0.5 * tanh( 2.0 * (sqrt(pow(x[0]-xc,2) + pow(x[1]-yc,2)) - R) / eps )",
    degree=1,  # polynomial degree used for interpolation
    xc=center_init_x,
    yc=center_init_y,
    R=initDropDiam/2,
    eps=eps
)

phi_n = fe.project(phi_init_expr, V)



# Define variational problems

bilin_form_mu = f_trial * v * fe.dx

lin_form_mu = 4*beta*(phi_n - 1)*(phi_n - 0)*(phi_n - 0.5)*v*fe.dx\
    + kappa*fe.dot(fe.grad(phi_n),fe.grad(v))*fe.dx 

mu_mat = fe.assemble(bilin_form_mu)
mu_solver = fe.KrylovSolver("cg", "ilu")
mu_solver.set_operator(mu_mat)


rhs_mu = fe.assemble(lin_form_mu)

mu_solver.solve(mu_n.vector(), rhs_mu)

[dokken]

There are several issues here:

1. Mesh dependence: Your mesh does not use crossed diagonal, i.e. all triangles are skewed in one direction. You can fix this by using diagonal="crossed" when creating the mesh.

This will make the code give a difference of

x0 = 5
dx = 2.5
y0 = 0.1
print(mu_n((x0 -dx, y0)), mu_n((x0+dx, y0)),mu_n((x0 -dx, y0))-mu_n((x0+dx, y0)))

0.0018872926763350327 0.0018871372270226824 1.554493123502934e-07

2. Usage of an iterative solver. Use LUSolver with 1., and you get

   mesh = fe.RectangleMesh(fe.Point(0, 0), fe.Point(L_x, L_y), nx, ny, diagonal="crossed")
   # Code the same here:
   # ....
   mu_solver = fe.LUSolver("mumps")
   

   0.0018872024997417568 0.0018872024997478832 -6.126392795846591e-15
   

3. Usage of projection into a lower order space.
   As you are doing a projection anyhow, I see no point of first projecting phi_init_expr into phi_n. I would simply use:

   try: 
      from ufl import tanh
   except ImportError:
       from ufl_legacy import tanh
   x = fe.SpatialCoordinate(mesh)
   R = initDropDiam / 2
   phi_n = 0.5 - 0.5 * tanh( 2.0 * (fe.sqrt((x[0]-center_init_x)**2 + (x[1]-center_init_y)**2) - R) / eps )
   

   As this gives a way better result if visualizing it:

   with fe.XDMFFile("mu_field.xdmf") as xdmf:
       xdmf.write(mesh)
       xdmf.write(mu_n, 0.0)
   

New approach:

Old approach:

Zoomed in:

[zaccandels]

Many thanks for your response! Would you (or anyone else reading this) mind referring me to a paper in which the notion of a crossed-diagonal FE mesh is developed? I’ve never once encountered this.

I’m very surprised that the mere asymmetry of the mesh produces such a drastic asymmetry in \nabla \phi_n, and by extension, \mu.

Surely this must be a well-documented shortcoming of uniform triangulations. To your knowledge, have there been other methods developed to contend with this?