Finite Volume Implementation

Question

I am trying to implement a simple finite volume method solver. I had a class on FVM a while back, but am still aware of the principal concepts. But implementing the FVM for non-cartesian or 1D meshes is not straightforward to me.

E.g. I have a mesh representing my control volumes. I compute the cell centroids to be used as my DOF later on. How would I compute the connectivity between cell centroids? I can compute this by searching for cells that share a cell surface, but this is very naive and hence computationally very inefficient. Further, I just don't see sophisticated softwares doing so.

I am looking for anything that might help me understand some of the tricks involved in implementing a FVM solver efficiently! Things that are not directly involved with the underlying mathematical/numerical schemes, but rather the set-up.

Answer 1

I have 1D finite-volume code written in python for a cell-centred mesh,

First generate a sequence which is the location of the faces, $ \{ x_{j-1/2} \}$. For example, for uniform spacing over the domain [0,1] this is a simple as,

a = 0
b = 1
J = 50 
faces = numpy.linspace(a, b, J}

After you have the faces the rest is easy because the faces uniquely determine the cell centres,

$$ x_j = \frac{1}{2}\left(x_{j-1/2} + x_{j+1/2} \right) $$

and the cell volumes (only a true volume in 3D but we stick with the wording),

$$ h_j = x_{j+1/2} - x_{j-1/2}. $$

For connivence you can define the centroid spacing, this is useful when interpolating values to/from centres/faces. Following Hundsorfer's approach,

$$ h_{-} = x_j - x_{j-1} = \frac{1}{2}\left( h_{j-1} + h_j \right) \\ h_{+} = x_{j+1} - x_j = \frac{1}{2}\left( h_{j} + h_{j+1} \right) $$

The following project on github might be useful, https://github.com/danieljfarrell/FVM to see how this can be implemented in practice. I use a Mesh and CellVariable objects (inspired by Fipy). You can instantiate the mesh and the query it for centroid position, mesh volume etc. The CellVariable class using the mesh to do transparent linear interpolation of solution variable.

Here is the mesh class and an example of how to use it. It might be a useful basis for your implementation,

class Mesh(object):
    """A 1D cell centered mesh defined by faces for the finite volume method."""

    def __init__(self, faces):
        super(Mesh, self).__init__()

        # Check for duplicated points
        if len(faces) != len(set(faces)):
            raise ValueError("The faces array contains duplicated positions. No cell can have zero volume so please update with unique face positions.")
        self.faces = np.array(faces)
        self.cells = 0.5 * (self.faces[0:-1] + self.faces[1:])
        self.J = len(self.cells)
        self.cell_widths = (self.faces[1:] - self.faces[0:-1])

    def h(self, i):
        """Returns the width of the cell at the specified index."""
        return self.cell_widths[i]

    def hm(self, i):
        """Distance between centroids in the backwards direction."""
        if not check_index_within_bounds(i,1,self.J-1):
            raise ValueError("hm index runs out of bounds")
        return (self.cells[i] - self.cells[i-1])

    def hp(self, i):
        """Distance between centroids in the forward direction."""
        if not check_index_within_bounds(i,0,self.J-2):
            raise ValueError("hp index runs out of bounds")
        return (self.cells[i+1] - self.cells[i])

import numpy
faces = numpy.linspace(0,1,6)
mesh = Mesh(faces)
mesh.cells # i.e. cell centres
# array([ 0.1,  0.3,  0.5,  0.7,  0.9])
mesh.cell_widths # i.e. cell volumes
# array([ 0.2,  0.2,  0.2,  0.2,  0.2])