# How to compute the potential of a 2D vector field? **URL:** \ **Category:** variational formulation\ **Created:** [April 7, 2021, 11:34am UTC](https://fenicsproject.discourse.group/t/how-to-compute-the-potential-of-a-2d-vector-field/5522 "2021-04-07T11:34:11Z")\ **Posts on this page:** 4\ **Page:** 1
**Author:** ![Jose\_Galaz](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/jose_galaz/32/2309_2.png) [@Jose\_Galaz](https://fenicsproject.discourse.group/u/Jose_Galaz)\ **Post date:** [April 7, 2021, 11:34am UTC](https://fenicsproject.discourse.group/t/how-to-compute-the-potential-of-a-2d-vector-field/5522/1 "2021-04-07T11:34:11Z")
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! ```auto 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) ``` ````auto --------------------------------------------------------------------------- RuntimeError Traceback (most recent call last) in 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 *** ------------------------------------------------------------------------- ​``` ```` ---
**Author:** ![kamensky](https://avatars.discourse-cdn.com/v4/letter/k/e95f7d/32.png) [@kamensky](https://fenicsproject.discourse.group/u/kamensky)\ **Post date:** [April 7, 2021, 1:58pm UTC](https://fenicsproject.discourse.group/t/how-to-compute-the-potential-of-a-2d-vector-field/5522/2 "2021-04-07T13:58:39Z")
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: ```python 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] **Author:** ![Jose\_Galaz](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/jose_galaz/32/2309_2.png) [@Jose\_Galaz](https://fenicsproject.discourse.group/u/Jose_Galaz)\ **Post date:** [April 7, 2021, 2:40pm UTC](https://fenicsproject.discourse.group/t/how-to-compute-the-potential-of-a-2d-vector-field/5522/3 "2021-04-07T14:40:59Z") 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!! ---
**Author:** ![kamensky](https://avatars.discourse-cdn.com/v4/letter/k/e95f7d/32.png) [@kamensky](https://fenicsproject.discourse.group/u/kamensky)\ **Post date:** [April 7, 2021, 4:11pm UTC](https://fenicsproject.discourse.group/t/how-to-compute-the-potential-of-a-2d-vector-field/5522/4 "2021-04-07T16:11:16Z")
> 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](https://en.wikipedia.org/wiki/Babu%C5%A1ka%E2%80%93Lax%E2%80%93Milgram_theorem) 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.