Pressure Poisson equation giving wrong results

[saitta-s]

Hello, I'm trying to recover the pressure field from a known fluid velocity field using the pressure Poisson equation.
I'm able to compute spatial derivatives of the velocity vector components. I know these are correct: I have ground truth derivatives.

I have assembled the right hand side of the pressure Poisson equation as:

bx = - rho * (u_t + rho * (u * u_x + v * u_y + w * u_z)) + mu * (u_xx + u_yy + u_zz)
by = - rho * (v_t + rho * (u * v_x + v * v_y + w * v_z)) + mu * (v_xx + v_yy + v_zz)
bz = - rho * (w_t + rho * (u * w_x + v * w_y + w * w_z)) + mu * (w_xx + w_yy + w_zz)
rhs = Wx.dot(bx) + Wy.dot(by) + Wz.dot(bz)

with subscripts indicating derivatives, rho and mu are fluid density and viscosity.
Then, I assembled the RBF weight matrix as: Wppe = weight_matrix(nodes[interiorIDs], nodes, n=80, [[2,0,0], [0,2,0], [0,0,2]], phi='phs5', order=5), and applied boundary conditions to Wppe and rhs.
However, when I solve the linear system, I obtain wrong pressure fields, both in patterns and values (15 orders of magnitude higher!).

I've tried lots of combinations for assembling the rhs and the coefficient matrix with no success.

Is there something wrong in the way I'm using the weight_matrix function? Or could the problem be related to my rhs?

Any advice or help would be greatly appreciated.
Thanks!

[treverhines]

Your usage looks correct; you are creating a matrix that maps the pressure at nodes to the Laplacian of the pressure at nodes[InteriorIDs]. What are your boundary conditions and how are you applying them?

[saitta-s]

Thanks for the reply.
interiorIDs is the list of node indices where I don't want to apply a zero pressure condition. On the other hand, outletIDs is the list of node indices where I want to apply zero pressure.
This is how I am enforcing the zero pressure boundary condition:

A_boundary  = weight_matrix(nodes[outletIDs], nodes, n=1, diffs=[0, 0, 0], phi=phi, order=0)
WLap        = expand_rows(Wppe, interiorIDs, nNodes)
WLap        += expand_rows(A_boundary, outletIDs, nNodes)

rhs[outletIDs] = 0.0

p = spsolve(WLap, rhs)

I believe this is how b.c. are enforced in the examples provided with the repo. Any idea of what I might be doing wrong?

Many thanks!

[treverhines]

That all looks good to me. If you set rhs[InteriorIDs] = 1 before you solve, do you get something that looks like the solution to Poisson's equation? If not, then maybe the issue is with the node locations. Do the outletIDs nodes completely encompass the domain?

[saitta-s]

Setting rhs[interiorIDs] = 1 and rhs[outletIDs] = 0 before solving, gives a typical Poisson's equation solution; a high pressure area that smoothly goes to zeros on outletIDs nodes. The outletIDs nodes are just the nodes on the zero pressure boundary.

I think I'm doing something wrong when assembling the right hand side from the known velocity field. Do you know where I could find a minimal code example of how to use RBF for solving laminar, steady state Navier-Stokes? I'm working on 3D data, but even a 2D basic channel flow example would be great.

Thanks again for the help!

[treverhines]

I noticed that the units were not matching up in your rhs. Try creating it as:

bx = - rho * u_t - rho * (u * u_x + v * u_y + w * u_z) + mu * (u_xx + u_yy + u_zz)
by = - rho * v_t - rho * (u * v_x + v * v_y + w * v_z) + mu * (v_xx + v_yy + v_zz)
bz = - rho * w_t - rho * (u * w_x + v * w_y + w * w_z) + mu * (w_xx + w_yy + w_zz)
rhs = Wx.dot(bx) + Wy.dot(by) + Wz.dot(bz)

[saitta-s]

It's exactly the same, isn't it? I get identical results. :(

[treverhines]

It would be the same if rho=1. Are Wx, Wy, and Wz also RBF-FD matrices? How are you making those?

[saitta-s]

rho = 1060 which is my fluid density.
This is how I create the RBF-FD matrices:

K = 65
phi = 'phs7'
q = 5
Wx = weight_matrix(nodes, nodes, K, [1, 0, 0], phi=phi, order=q)
Wy = weight_matrix(nodes, nodes, K, [0, 1, 0], phi=phi, order=q)
Wz = weight_matrix(nodes, nodes, K, [0, 0, 1], phi=phi, order=q)

Wxx = weight_matrix(nodes, nodes, K, [2, 0, 0], phi=phi, order=q)
Wyy = weight_matrix(nodes, nodes, K, [0, 2, 0], phi=phi, order=q)
Wzz = weight_matrix(nodes, nodes, K, [0, 0, 2], phi=phi, order=q)

And then I compute, for instance, the derivative of u w.r.t. x as u_x = Wx.dot(u).

Are these correct?

Thanks for the help, very much appreciated.

[treverhines]

that looks right. Do the results look like the issue is numerical instability, or just solving the wrong equation? If there are numerical instability issues, maybe try setting q=2 and phi="phs3"

[saitta-s]

I tried lots of combinations for q and phi with no significant difference. It just looks like I'm solving the wrong equation.

[treverhines]

What are you using for the time derivatives u_t, v_t, and w_t? Are they zero or is there a time integration component as well?

[saitta-s]

I'm working with velocity data from a transient CFD simulation. A curved pipe with one inlet and one outlet. A constant velocity is applied at the inlet, therefore the velocity fields at to consecutive timesteps are almost identical. The time derivatives are in fact close to zero. I tried removing them, but there's no significant difference.

[treverhines]

I am not too familiar with fluid mechanics, but wouldn't you also need some sort of boundary condition like u=0 at the sides of the pipe?

[saitta-s]

That's correct, but in my case, I'm not actually solving for velocity. I want to recover the pressure field from the pre-computed velocity field. So I assumed I should not worry about the no slip condition when solving for pressure or when I compute spatial derivatives.

[treverhines]

I see your point that the no slip condition would not be helpful for your problem. Are there any other boundary conditions that would make sense at the sides of the pipe? From the perspective of solving Poisson's equation, my understanding is that you should have some combination of Dirichlet and Neumann boundary conditions around the entire domain in order for the problem to be well-posed. What is the condition number for your LHS? Would it make sense to enforce that dp/dn = 0 at the boundary, where p is pressure and n is the normal vector to the boundary?

[saitta-s]

np.linalg.cond(WLap.toarray()) = 5.986725487273041e+17. Could this be an issue?
Something I noticed is that if I calculate p = spsolve(WLap, rhs) or p = np.linalg.solve(WLap.toarray(), rhs) I obtain very different results. They both look like Poisson's equation results but very different distributions and values.

I guess I could try to enforce dp/dn = 0 at the wall. I should use ghost nodes to do that, right? I was able to add one layer of ghost nodes along each wall node's normal, but I'm not sure how distant to the wall I should place them. If I do this, how would I implement the dp/dn = 0 condition?

Thanks!

[treverhines]

Yup, that would be the issue. The matrix is very close to singular.

You could first try it without ghost nodes. The solution would be much less accurate without the ghost nodes, but you would be able to see if the dp/dn = 0 boundary conditions get you closer to the pressures you are expecting.

Implementing the dp/dn = 0 boundary conditions without ghost nodes is easy as long as you have the outward normal vectors at the boundaries. Creating the matrix that maps the pressures at nodes to dp/dn at boundary_nodes, which have normal vectors boundary_normals, would be something like:

 weight_matrix(
    boundary_nodes,
    nodes,
    n=n, 
    diffs=[[1, 0], [0, 1]], 
    coeffs=[boundary_normals[:,0], boundary_normals[:, 1]], 
    phi=phi, 
    order=order)

Once you have that matrix, implementing the dp/dn=0 boundary condition is just like implementing the fixed boundary condition.

Ghost nodes require a little more effort to implement. If you look at the second example down here https://rbf.readthedocs.io/en/latest/fd.html, you can see an example of solving Poisson's equation with a mix of fixed and free boundary conditions, where the free boundaries have ghost nodes. When I generate ghost nodes, each ghost node's spacing from the boundary is equal to half the nearest neighbor distance for its corresponding boundary node. However, that was a pretty arbitrary decision.

[saitta-s]

I think the use of coeffs in irregular 3D geometries has not been tested. When trying to create the matrix as you suggested, I get a ValueError: operands could not be broadcast together with remapped shapes [original->remapped]: (2092,) and requested shape (1000,). This is because my boundary_normals is an array of shape (2092 x 3), with 2092 being the number of boundary_nodes.
Is there a way to avoid using coeffs?

Thanks

[treverhines]

Thanks for pointing this out. It is due to a bug I recently introduced. Can you try creating the matrix with the argument chunk_size=None? If you have too many boundary nodes, then it may cause a memory error. Otherwise, you can create the matrix for each term in the differential operator, multiply the resulting matrix by scipy.sparse.diags(boundary_normals[:, i]), then add the matrices together.

[treverhines]

The bug is now fixed in master

[saitta-s]

Alright, now I obtain something much closer to the expected solution, although still one order of magnitude larger. I think this is due to the error introduced by the rbf derivative operators, which I noticed being especially large near the walls. I might give ghost nodes a try.
Just to be sure, this is how I implemented the dp/dn = 0 boundary condition:

Wxxyyzz       = weight_matrix(nodes[interiorIDs], nodes, K, LaplacianDiff, phi=phi, order=q)
A_boundary    = weight_matrix(nodes[outletIDs], nodes[0:nNodes], n=1, diffs=[0,0,0], phi=phi, order=0)
Wpn           = weight_matrix(nodes[wall_idx], nodes, K, diffs=[[1, 0, 0], [0, 1, 0], [0, 0, 1]],
                              coeffs=[wall_normals[:, 0], wall_normals[:, 1], wall_normals[:, 2]],
                              phi=phi, order=2, chunk_size=None)
WLap         = expand_rows(Wxxyyzz, interiorIDs, nNodes)
WLap         += expand_rows(A_boundary, outletIDs, nNodes)
WLap         += expand_rows(Wpn, wall_idx, nNodes)

and then before solving:

rhs[outletIDs] = 0.0
rhs[wall_idx] = 0.0

I also noticed that the error increases if I use higher degree phs or higher order monomials, which is not what I was expecting.

Anyway, thanks a lot.

[treverhines]

That looks right. The increased error at the edges with higher order monomials and phs is due to Runge's phenomenon. You may get better results at the edges using phi="phs2" and order=1. Adding ghost nodes would be the best option though.

[saitta-s]

No luck.
I implemented ghost nodes and tried multiple combinations of phi, order and n. Results are always different and not what I expect.

[treverhines]

Are the results at least consistent? Or are you getting significant variations in the results for the different parameter choices?

[saitta-s]

Solving the linear systems gives the same results, no matter what solver I use. However, if I use a greater number of stencil knots, I get very different results (worse, compared to the desired solution). Same thing for the choice of phi; if I change it to phs2 or anything that is not phs3, my results drastically change.

[treverhines]

That is odd and it has me thinking that there may be an issue with your nodes/domain. I have several question: Are your nodes the vertices of the mesh used for the CFD simulation? Can you show me what your domain looks like? Can you show me what some of these incorrect solutions look like? What is the condition number of your LHS matrix? Does your rhs also look significantly different for different parameter choices?

[saitta-s]

First is a visualization of the mesh nodes and my geometry. Second, a pressure solution obtained with RBF-FD (left) and the desired solution (right).
The condition number is around 10^8.
The rhs doesn't change significantly, is just a sum of spatial derivatives.

[treverhines]

Thanks! The domain does not look like an issue. It looks like the problem is still ill-conditioned. Are you missing a BC on the left end? I realize it is cheating, but if you fix the left end at 20,000, would that give you the right answer? Does enforcing dp/dn=0 at the left end give you a better solution?

[saitta-s]

I don't think I'm missing any conditions. I tried to fix the solution on one node at the left end of my domain; the solution changes only slightly. The distribution you see in the figure was calculated with the dp/dn = 0 condtion. If I don't do that, I obtain much higher values (up to 10^20) and nonsense distributions.

[treverhines]

Would you be able to send me your mesh file and velocities? The files may be small enough to upload here.

[saitta-s]

ts_1_0_0.zip

Everything is in the .vtu file.
This is the piece of code I use to obtain outlet nodes indices where p = 0:

outletIDs = np.arange(0, nNodes)[abs(nodes[:, 0] - 0.1000) < 1e-5].tolist()
interiorIDs = np.arange(0, nNodes)
ii = [i for i, e in enumerate(interiorIDs.tolist()) if e not in outletIDs]
interiorIDs = interiorIDs[ii].tolist()

[treverhines]

I think I see the problem. The solution is mostly being controlled by the large velocity derivatives at the left end of the pipe. It also looks like the node spacing is too coarse to accurately resolve the velocity derivatives.

Consequently, the RHS is very sensitive to the choice of phi and order. For example, here is the RHS with phi="phs3" and order=2

The corresponding pressures are

and here is the RHS with phi="phs5" and order=3

the corresponding pressures are

You can see that the sign of the pressure solution flips.

I am not sure what the best solution would be to mitigate this issue. Do you care about what happens at the left end? or are the large velocity gradients just a boundary artifact to you?

Here is the code that I used to generate these results

import numpy as np
import meshio
from rbf.pde.geometry import simplex_outward_normals
from rbf.pde.fd import weight_matrix
from rbf.sputils import expand_rows
from scipy.spatial import cKDTree

from scipy.sparse.linalg import spsolve

phi = 'phs5'
order = 3
n = 50
rho = 1060.0
mu = 1.0 # not sure what this should be

data = meshio.read('ts_1_0_0.vtu')
nodes = data.points
N = len(nodes)
u, v, w = data.point_data['velocity'].T

# compute the RHS at `nodes`. Note I have removed the time derivative term
Wx = weight_matrix(nodes, nodes, n, [1, 0, 0], phi=phi, order=order)
Wy = weight_matrix(nodes, nodes, n, [0, 1, 0], phi=phi, order=order)
Wz = weight_matrix(nodes, nodes, n, [0, 0, 1], phi=phi, order=order)
Wxx = weight_matrix(nodes, nodes, n, [2, 0, 0], phi=phi, order=order)
Wyy = weight_matrix(nodes, nodes, n, [0, 2, 0], phi=phi, order=order)
Wzz = weight_matrix(nodes, nodes, n, [0, 0, 2], phi=phi, order=order)
bx = -rho*(u*Wx.dot(u) + v*Wy.dot(u) + w*Wz.dot(u)) + mu*(Wxx.dot(u) + Wyy.dot(u) + Wzz.dot(u))
by = -rho*(u*Wx.dot(v) + v*Wy.dot(v) + w*Wz.dot(v)) + mu*(Wxx.dot(v) + Wyy.dot(v) + Wzz.dot(v))
bz = -rho*(u*Wx.dot(w) + v*Wy.dot(w) + w*Wz.dot(w)) + mu*(Wxx.dot(w) + Wyy.dot(w) + Wzz.dot(w))
rhs = Wx.dot(bx) + Wy.dot(by) + Wz.dot(bz)

# break the tets up into faces. If a face is unique, then it is part of the
# boundary
cells = data.cells[0].data
faces = cells[:, [[0, 1, 2], [0, 1, 3], [0, 2, 3], [1, 2, 3]]]
faces = faces.reshape(-1, 3)
unique_faces, count = np.unique(np.sort(faces, axis=1), axis=0, return_counts=True)
boundary_faces = unique_faces[count == 1]

# break the boundary faces up into faces that get fixed/free BCs
fixed_faces, free_faces = [], []
for face in boundary_faces:
    if np.all(np.abs(nodes[face, 0] - 0.1) < 1e-5):
        fixed_faces.append(face)
    else:
        free_faces.append(face)

# get the node indices making up the different boundaries
fixed_ids = np.unique(fixed_faces)
free_ids = np.array([i for i in np.unique(free_faces) if i not in fixed_ids])
interior_ids = np.array([i for i in range(N) if ((i not in fixed_ids) & (i not in free_ids))])

# Get the normal vectors for each free node. The normal at each node will be
# the average of the normals for each face containing the node
free_normals = np.zeros((len(free_ids), 3))
face_normals = simplex_outward_normals(nodes, free_faces)
for idx, node_id in enumerate(free_ids):
    for face, nrm in zip(free_faces, face_normals):
        if node_id in face:
            free_normals[idx] += nrm

free_normals /= np.linalg.norm(free_normals, axis=1)[:, None]

# create and append the ghost nodes. dx is the distance from the free nodes to
# their nearest neighbors
dx = cKDTree(nodes).query(nodes[free_ids], 2)[0][:, 1]
ghost_nodes = nodes[free_ids] + 0.5*dx[:, None]*free_normals

nodes = np.vstack((nodes, ghost_nodes))
ghost_ids = np.arange(N, N + len(ghost_nodes))
N = len(nodes)

# build the LHS
A_interior = weight_matrix(
    x=nodes[interior_ids],
    p=nodes,
    n=n,
    diffs=[[2, 0, 0], [0, 2, 0], [0, 0, 2]],
    phi=phi,
    order=order
    )

A_ghost = weight_matrix(
    x=nodes[free_ids],
    p=nodes,
    n=n,
    diffs=[[2, 0, 0], [0, 2, 0], [0, 0, 2]],
    phi=phi,
    order=order
    )

A_free = weight_matrix(
    x=nodes[free_ids],
    p=nodes,
    n=n,
    diffs=[[1, 0, 0], [0, 1, 0], [0, 0, 1]],
    coeffs=free_normals.T,
    phi=phi,
    order=order
    )

A_fixed = weight_matrix(
    x=nodes[fixed_ids],
    p=nodes,
    n=1,
    diffs=[0, 0, 0],
    phi=phi,
    order=0
    )

A = expand_rows(A_interior, interior_ids, N)
A += expand_rows(A_ghost, ghost_ids, N)
A += expand_rows(A_free, free_ids, N)
A += expand_rows(A_fixed, fixed_ids, N)
print('condition number: %.2e' % np.linalg.cond(A.A))

d = np.zeros((N,))
d[interior_ids] = rhs[interior_ids]
d[ghost_ids] = rhs[free_ids]

soln = spsolve(A, d)
# remove the solution at ghost nodes
soln = soln[:len(data.points)]

data.point_data['rbf-fd-pressure'] = soln
data.point_data['rhs'] = rhs
data.write('soln.vtu')

[saitta-s]

I see. I'll try to run another test with a finer mesh to check if the coarse node spacing is the issue.
Anyway, thank you very much for all your help.