Table of Contents

  • 1  Learning goals

  • 2  The problem

    • 2.1  2-D discrete approximation stencil for steady-state diffusion.

      • 2.1.1  2D stencil

    • 2.2  Q1 (2pts) Write your 2d Stencil

    • 2.3  The equations for each gridblock in the mesh.

      • 2.3.1  Row-major numbering

    • 2.4  Q2 Writing equations for each gridblock

      • 2.4.1  Q2 (2pts) entry cell below

    • 2.5  Q3 Put your equations into the matrix

      • 2.5.1  Q3 (4pts) entry cell below

    • 2.6  Q4 (4pts) Solve the system of equations

    • 2.7  Check your answer

  • 3  Reflection

Assignment Finite volumes by hand

Learning goals

  1. Can apply the finite-volume stencil to a 2D steady-state problem, and generate the equations for each unknown in the mesh.

  2. Can place the equations into a system matrix.

  3. Can compute fluxes.

The problem

Consider the two-dimensional diffusion problem below. There are zero-flux boundary conditions on the top and bottom of the domain, and prescribed concentration (more generally known as Dirichlet or first-type boundary conditions)on the to sides. Each gridblock is the same dimension \(\Delta x = \Delta y = 200~m\), and \(\Delta z= 3~m\) (out of the page). The gridblock node is placed at Dirichlet (prescribed concentration) boundaries. The diffusion coefficient in each gridblock is as shown in the figure and the porosity is \(\theta = 0.25\) is the same everywhere.

2-D discrete approximation stencil for steady-state diffusion.

Recall the 1D stencil looks like this.

\begin{align} &\left(D\theta {c_E - c_C \over \Delta x} + D\theta {c_W - c_C \over \Delta x} \right) (\Delta y) (\Delta z) =0 \label{8fvbh1}\\ \end{align}

2D stencil

Change the following 1-D stencil to 2-D by adding the “N” and “S” terms.

\begin{align} &\left(D\theta {c_E - c_C \over \Delta x} + D\theta {c_W - c_C \over \Delta x} \right) (\Delta y) (\Delta z) =0 \label{8fvbh2}\\ \end{align}

You may also want to simplify your stencil for the case of constant diffusion coefficient and constant gridblock dimensions.

Q1 (2pts) Write your 2d Stencil

In the cell below, enter your stencil in markdown.

YOUR ANSWER HERE

The equations for each gridblock in the mesh.

How many unknowns are there? In the 2-d problem above, we have 3 rows and 4 columns so 12 unknowns in total. Problem now our mesh is two dimensional - how can we number our gridblocks and unknowns? In 1-D, we just numbered one gridblock after another: We eventually want to put the system of equations into a matrix that looks like this (for a 5-gridblock 1-D problem with two boundary concentrations \(bc_1\) and \(bc_5\))

\begin{align*} \mathbf{Ac}&=\mathbf{b}\\ {\begin{bmatrix} 1 & 0 & 0 & 0 & 0 \\ 1 & -2 & 1 & 0 & 0\\ 0 & 1 & -2 & 1 & 0\\ 0 & 0 & 1 & -2 & 1\\ 0 & 0 & 0 & 0 & 1\\ \end{bmatrix}} \begin{bmatrix} c_1\\ c_2\\ c_3\\ c_4\\ c_5\\ \end{bmatrix} &=\begin{bmatrix} bc_1\\ 0\\ 0\\ 0\\ bc_5\\ \end{bmatrix} \end{align*}

We will have a similar matrix structure for a 2, 3 or n-dimensional problem. For a multidimensional problem, the key is to number the gridblocks in an organized way. #### Row-major numbering Python is organized so that so-called row-major numbering is the most efficient. That means we should number our gridblocks as shown here:

Where we used the python style and gave the first gridblock the index 0. This numbering can be considered a mapping of the row, column structure to a sequential, 1-D vector of indices. In the 2-D code you were given, we used the following function to map a gridblock index into row and column index that we used to assemble the equations in the algorithm.

[1]:

import numpy as np
def index_to_row_col(ind, nrows, ncol):
    """
    in a 2D array, returns the row and column value
    associated with a 1D index
    Bottom left is index is zero (0-th row, 0-th column)
    while index one is the 0-th row and 1st column
    """
    if ind > nrows * ncol - 1:
        return 0

    row = int(np.floor(ind / ncol))
    col = int(ind - row * ncol)
    return row, col

Let’s see if it returns the proper indices for our situation where nrows = 3 and ncol = 4. According to the figure, the gridblock number 6 should be in python row 1, python column 2 (a python row or column begins with index 0):

[2]:

ind = 6
nrows = 3
ncols = 4
row, column = index_to_row_col(ind, nrows, ncols)
print(f"The gridblock {ind} is in python row {row} and python column {column}")

The gridblock 6 is in python row 1 and python column 2

Q2 Writing equations for each gridblock

Now, write your equations using your stencil, the gridblock dimensions, diffusion coefficient and boundary conditions. I’ll write a few (but not all) of the simple equations at the Dirichlet (specified concentration) boundaries:

\begin{align*} &c_0 = 64\\ &c_4 = 64\\ &c_3= 125\\ &c_7 = 125\\ \end{align*}

Use your stencil and write the other equations here in markdown. Zero-flux boundaries

Hint: the 2-D stencil is basically just \(J_{WC}+J_{EC}+J_{SC} + J_{NC}=0\). When the C gridblock is adjacent to a zero-flux boundary, one of the terms in the stencil is set to zero. For example in gridblock 1 above, we have this mapping of the stencil \(\rightarrow\) gridblock index: C \(\rightarrow\) 1, W \(\rightarrow\) 0, N \(\rightarrow\) 5, but to the S we have a zero-flux boundary. Accordingly, in the stencil for gridblock with index 1, \(J_{SC}=0\) and the stencil simplifies to \(J_{WC}+J_{EC}+ J_{NC}=0\).

Q2 (2pts) entry cell below

YOUR ANSWER HERE

Q3 Put your equations into the matrix

I’ve put some of the easy ones in. You do the rest by replacing the zeros in the markdown array below with the correct coefficients from your equations. Double click on this cell to get the raw latex formatting, then copy it below into the answer box and edit.

\begin{align*} \mathbf{Ac} &= \mathbf{b}\\ {\begin{bmatrix} 1 & 0 & 0 & 0 & 0 &0 &0 &0 &0 &0 &0 &0\\ 0 & 0 & 0 & 0 & 0 &0 &0 &0 &0 &0 &0 &0 \\ 0 & 0 & 0 & 0 & 0 &0 &0 &0 &0 &0 &0 &0 \\ 0 & 0 & 0 & 1 & 0 &0 &0 &0 &0 &0 &0 &0 \\ 0 & 0 & 0 & 0 & 1 &0 &0 &0 &0 &0 &0 &0 \\ 0 & 0 & 0 & 0 & 0 &0 &0 &0 &0 &0 &0 &0 \\ 0 & 0 & 0 & 0 & 0 &0 &0 &0 &0 &0 &0 &0 \\ 0 & 0 & 0 & 0 & 0 &0 &0 &1 &0 &0 &0 &0 \\ 0 & 0 & 0 & 0 & 0 &0 &0 &0 &0 &0 &0 &0 \\ 0 & 0 & 0 & 0 & 0 &0 &0 &0 &0 &0 &0 &0 \\ 0 & 0 & 0 & 0 & 0 &0 &0 &0 &0 &0 &0 &0 \\ 0 & 0 & 0 & 0 & 0 &0 &0 &0 &0 &0 &0 &0 \\ \end{bmatrix}} \begin{bmatrix} c_0\\ c_1\\ c_2\\ c_3\\ c_4\\ c_5\\ c_6\\ c_7\\ c_8\\ c_9\\ c_{10}\\ c_{11}\\ \end{bmatrix} &= \begin{bmatrix} 64\\ 0\\ 0\\ 125\\ 64\\ 0\\ 0\\ 125\\ 0\\ 0\\ 0\\ 0\\ \end{bmatrix} \end{align*}

Q3 (4pts) entry cell below

YOUR ANSWER HERE

Q4 (4pts) Solve the system of equations

You have a few options here: 1. Guess the answer and check if the equations are satisfied. 2. Use python. See how the code uses the statement c = np.linalg.solve(A, b). The trick is to put your equations into numpy arrays A and b. To check that you have solved the equations, we ask that you report the concentrations in gridblocks 5 and 6 for autograding in the cell below.

[3]:

# replace -9999 with the values of the concentration in gridblocks 5 and 6 below
# c5 = -9999
# c6 = -9999
# YOUR CODE HERE
raise NotImplementedError()

---------------------------------------------------------------------------
NotImplementedError                       Traceback (most recent call last)
<ipython-input-3-ef678f5605af> in <module>
      3 # c6 = -9999
      4 # YOUR CODE HERE
----> 5 raise NotImplementedError()

NotImplementedError:

[ ]:

# Hidden test follows here.

Check your answer

Take the equation for any gridblock and substitute in the concentrations you computed. Is the equation satisfied?

Reflection

  1. Each row in the coefficient matrix has only a few non-zero coefficients. Why?

  2. What would this system matrix look like for a 3-d discretization?

  3. How would the python function index_to_row_col change for a 3-d discretization, where we now have row, columns and layers?