How to compute the potential of a 2D vector field?

variational formulation
1

I’m trying to compute the potential \Phi of a 2D velocity field (U,W) but I’m unable to represent the source correctly and I’m getting a RuntimeError that says Error: Unable to successfully call PETSc function 'KSPSolve'..

If I understand correctly the variational formulation should be

\int \nabla \Phi \cdot v dx= \int (U,W)^T \cdot v dx

with v a vector test function.

Below are the python code and error trace. I would appreciate any hints on how to improve this code!

I’m running this from this docker image quay.io/fenicsproject/stable:current in jupyter.

Thanks!

import fenics
import matplotlib.pyplot as plt
%matplotlib inline

mesh = fenics.UnitSquareMesh(12,12)

V = fenics.FunctionSpace(mesh, 'P', 1)
V_vec = fenics.VectorFunctionSpace(mesh, 'P', 1, dim=2)

class VelocityFieldExpression(fenics.UserExpression):
    def eval(self, values, x):
        values[0] = x[0]
        values[1] = x[1]
    def value_shape(self):
        return (2,)

u = fenics.TrialFunction(V)
v = fenics.TestFunction(V_vec)
velocity = VelocityFieldExpression()

a = fenics.dot(fenics.grad(u), v)*fenics.dx
L = fenics.dot(velocity, v)*fenics.dx
phi = fenics.Function(V)

fenics.solve(a == L, phi)

---------------------------------------------------------------------------
RuntimeError                              Traceback (most recent call last)
<ipython-input-7-c19cd53db90f> in <module>
     23 phi = fenics.Function(V)
     24 
---> 25 fenics.solve(a == L, phi)

/usr/local/lib/python3.6/dist-packages/dolfin/fem/solving.py in solve(*args, **kwargs)
    218     # tolerance)
    219     elif isinstance(args[0], ufl.classes.Equation):
--> 220         _solve_varproblem(*args, **kwargs)
    221 
    222     # Default case, just call the wrapped C++ solve function

/usr/local/lib/python3.6/dist-packages/dolfin/fem/solving.py in _solve_varproblem(*args, **kwargs)
    245         solver = LinearVariationalSolver(problem)
    246         solver.parameters.update(solver_parameters)
--> 247         solver.solve()
    248 
    249     # Solve nonlinear variational problem

RuntimeError: 

*** -------------------------------------------------------------------------
*** DOLFIN encountered an error. If you are not able to resolve this issue
*** using the information listed below, you can ask for help at
***
***     fenics-support@googlegroups.com
***
*** Remember to include the error message listed below and, if possible,
*** include a *minimal* running example to reproduce the error.
***
*** -------------------------------------------------------------------------
*** Error:   Unable to successfully call PETSc function 'KSPSolve'.
*** Reason:  PETSc error code is: 62 (Invalid argument).
*** Where:   This error was encountered inside /tmp/dolfin/dolfin/la/PETScKrylovSolver.cpp.
*** Process: 0
*** 
*** DOLFIN version: 2019.1.0
*** Git changeset:  74d7efe1e84d65e9433fd96c50f1d278fa3e3f3f
*** -------------------------------------------------------------------------
​```
2

The problem is that V and V_vec have different numbers of degrees of freedom, and thus the bilinear form a assembles to a non-square matrix. You can fix this by projecting in an H^1 sense: Find \phi\in H^1(\Omega) such that for all v\in H^1(\Omega)

\int_\Omega\nabla\phi\cdot\nabla v\,d\Omega = \int_\Omega \mathbf{u}\cdot\nabla v\,d\Omega\text{ ,}

where \mathbf{u} is the velocity field. This is still ill-posed, because \phi is determined only up to a constant, but you can apply a Dirichlet BC to one degree of freedom in the computation.

Here is the modified code:

import fenics
import matplotlib.pyplot as plt

N = 12
mesh = fenics.UnitSquareMesh(N,N)

V = fenics.FunctionSpace(mesh, 'P', 1)

class VelocityFieldExpression(fenics.UserExpression):
    def eval(self, values, x):
        values[0] = x[0]
        values[1] = x[1]
    def value_shape(self):
        return (2,)

u = fenics.TrialFunction(V)
v = fenics.TestFunction(V)
velocity = VelocityFieldExpression()
n = fenics.FacetNormal(mesh)

# Use H^1 projection:
a = fenics.dot(fenics.grad(u),fenics.grad(v))*fenics.dx
L = fenics.dot(velocity,fenics.grad(v))*fenics.dx

# Pin one corner to eliminate free constant in potential and have a unique
# solution.
bc = fenics.DirichletBC(V,fenics.Constant(0),
                        "x[0]<DOLFIN_EPS && x[1]<DOLFIN_EPS",
                        "pointwise")
phi = fenics.Function(V)
fenics.solve(a == L, phi, bc)

# Check H^1 error in potential; converges at expected rate.
import math
grad_e = fenics.grad(phi) - velocity
print(math.sqrt(fenics.assemble(fenics.dot(grad_e,grad_e)*fenics.dx)))
3

Thanks! your explanation is very clear and works perfectly!

I have two questions:

Is “projecting in another sense” the same as changing the inner product used for the projection? If yes, then, is there a rule to “safely” choose a different inner product?

Also, does this choice come from a mathematical or an implementation constraint? At least on the infinite-dimension space, I don’t understand why one couldn’t formulate the problem on a vector test space.

Thanks!!

4

Is “projecting in another sense” the same as changing the inner product used for the projection? If yes, then, is there a rule to “safely” choose a different inner product?

For my formulation to be understood as a projection, you need to assume that \exists \phi_\text{exact}\in H^1 such that \mathbf{u} := \nabla\phi_\text{exact}. If you already had \phi_\text{exact}, you could of course project it onto a finite element space with respect to whatever inner product you want, but, if you only have access to \mathbf{u}, that limits the choices somewhat. It doesn’t constrain you to a unique choice, though; you could still easily formulate, e.g., a weighted H^1 inner product projection.

Also, does this choice come from a mathematical or an implementation constraint? At least on the infinite-dimension space, I don’t understand why one couldn’t formulate the problem on a vector test space.

You would still need to satisfy inf-sup stability in the infinite-dimensional setting. If the vector-valued test space were too large, it could include functions orthogonal to all gradients of scalar functions in the trial space, which would not satisfy the non-degeneracy condition required of the bilinear form. You might formally use a vector-valued test space constrained to, say, an irrotational subspace of \mathbf{H}(\text{curl}), but then the easiest way to construct a finite-dimensional subspace of that would be to define it as the space of gradients of functions in a finite-dimensional scalar space, which could be implemented as above.

1 Like