usnistgov / fipy

FiPy is a Finite Volume PDE solver written in Python
http://pages.nist.gov/fipy/en/latest
Other
506 stars 148 forks source link

Is there any ways to get upwind face values? #712

Open simulkade opened 4 years ago

simulkade commented 4 years ago

I like to get the face values of a cell variable based on the direction of a velocity vector at the cell faces. Is there any way to do it? I think it is necessary for solving nonlinear equations with an advection term. There is an example here Might also be related to #650

wd15 commented 4 years ago

Hi @simulkade,

This seems to work

from fipy import Grid2D, CellVariable, FaceVariable
import numpy as np

def upwindValues(mesh, field, velocity):
    """Calculate the upwind face values for a field variable

    Note that the mesh.faceNormals point from `id1` to `id2` so if velocity is in the same
    direction as the `faceNormal`s then we take the value from `id1`s and visa-versa.

    Args:
      mesh: a fipy mesh
      field: a fipy cell variable or equivalent numpy array
      velocity: a fipy face variable (rank 1) or equivalent numpy array

    Returns:
      numpy array shaped as a fipy face variable
    """
    # direction is over faces (rank 0)
    direction = np.sum(np.array(mesh.faceNormals * velocity), axis=0)
    # id1, id2 are shaped as faces but contains cell index values
    id1, id2 = mesh._adjacentCellIDs
    return np.where(direction >= 0, field[id1],  field[id2])

mesh = Grid2D(nx=3, ny=3)
print(
    upwindValues(
        mesh,
        np.arange(mesh.numberOfCells),
        2 * np.random.random(size=(2, mesh.numberOfFaces)) - 1
    )
)

This gives you the upwind value as a face variable. I wasn't clear to me whether you wanted the upwind values as a face variable or as a cell to face shaped numpy array (that would be (4, numberOfCells) in this case).

simulkade commented 4 years ago

Thanks a lot, Daniel. I was thinking about something like the arithmeticFaceValue. I'll try your code now and see if I can get what I want do with the Buckley-Leverett equation.

simulkade commented 4 years ago

I wrote about solving this BL problem here, i.e., how to linearize the PDE, use FVM to discretize it, and solve the system of equations. I use my simple code that follows FiPy's approach, but the only thing that I have added is this upwind averaging function and also having ghost cells for easier implementation of the boundary conditions. For this problem, I do need to use upwind because there is a shock front that makes the second order solution unstable. I have also written here why we really need an upwind averaging scheme for a nonlinear advection term that is discretized with the upwind scheme. I tried your suggestion (and it works fine) but I still cannot get FiPy to solve the BL problem. I'm going to try a simpler problem. Maybe something is wrong with my implementation.