pynitevolumes-edu 0.0.3

Description

This library aims at facilitating the coding of algorithms in the Finite Volumes family. These algorithms are used for computing so-called numerical solutions to partial differential equations (PDE). The Pynite Volumes library provide tools to help writing such algorithms.

Usage

At their heart any numerical method for a PDE works by solving a linear system or a sequence of linear systems and as such the numerical methods describes how the matrix and the right-hand side of the linear system are constructed from the data of the problems. Note that the library is concerned with making the linear system (matrix and right-hand side) not solving the linear system which needs to be delegated to dedicated libraries such as Numpy or Scipy.

Let us say that one wants to solve a Laplace (or Poisson) equation with some homogenous Dirichlet boundary conditions

First you would need a mesh of the domain \(\Omega\). Details about the mesh are stored inside an object, let us call it M.

• If the mesh is “good” (we would rather say admissible) and if we want to use the preconstructed tools from the library as black boxes we would do:

  import numpy as np
  from scipy.sparse import csr_array
  from scipy.sparse.linalg import spsolve
  from pynitevolumes.schemes.tpfa.tpfa_basics import (diffusion_matrix, diffusion_rhs,
                                                      build_whole_dirichlet, source)
  
  M = ... #  obtain M from various ways
  
  dof = M.nvol  # the number of unkowns of the system
  
  A = 1.0  #  because the diffusion coefficient in this problem is 1.0
  
  boundary_condition = build_whole_dirichlet(M, 0.0)  #  the value is 0.0 on the boundary
  
  mat_row_d, mat_col_d, mat_dat_d = diffusion_matrix(M, A, boundary_condition)
  rhs_row_d, rhs_dat_d = diffusion_rhs(M, A, boundary_condition)
  rhs_row_source, rhs_dat_source = source(M, 1.0)  #  the right-hand side is the constant function 1.0
  
  mat_row = np.concatenate(mat_row_d)
  mat_col = np.concatenate(mat_col_d)
  mat_dat = np.concatenate(mat_dat_d)
  Mat = csr_array((mat_dat, (mat_row, mat_col)), shape=(dof, dof))  # builds the matrix
  
  # Build the right-hand side
  rhs_row = np.concatenate(rhs_row_d+rhs_row_source)
  rhs_dat = np.concatenate(rhs_dat_d+rhs_dat_source)
  B = csr_array((rhs_dat, (rhs_row, np.zeros(len(rhs_row)))), shape=(dof, 1))  # builds the right-hand side
  
  U = spsolve(Mat, B) #  obtain the numerical solution!

• If the mesh is again “good” and we want to check that we have understood the standard Finite Volume scheme of TPFA for this problem we could try to build the matrix and the right-hand side “by ourselves” in a very detailed way. Then our code could look like that.

  import numpy as np
  from scipy.sparse import csr_array
  from scipy.sparse.linalg import spsolve
  
  M = ... #  obtain M from various ways
  
  dof = M.nvol
  
  A = 1.0  #  because the diffusion coefficient in this problem is 1.0
  
  mat_row = []
  mat_col = []
  mat_dat = []
  
  
  # the next two lines are only there to follow classical mathematical notations for this method
  K = M.diam_i.c[0]
  L = M.diam_i.c[1]
  for s in range(M.nin):  # loop over all inner edges/diamonds because this is how data is stored in our framework
      mat_row.extend([K[s], K[s], L[s], L[s]])
      mat_col.extend([K[s], L[s], K[s], L[s]])
      val_s = M.len_i[s]/M.dKL_i[s]
      mat_dat.extend([val_s, -val_s, -val_s, val_s])
  
  Kbnd = M.diam_b.c
  for s in range(M.nbnd):  # loop now over boundary edges/diamonds
      mat_row.append(Kbnd[s])
      mat_col.append(Kbnd[s])
      mat_dat.append(M.len_b[s]/M.dKs_b[s])
  Mat = csr_array((mat_dat, (mat_row, mat_col)), shape=(dof, dof))  # builds the matrix
  
  B = np.zeros(dof)
  # Build the right-hand side
  for k in range(M.nvol):
      B[k] = M.vol[k]*1.0
  
  U = spsolve(Mat, B) #  obtain the numerical solution!

  (Of course this particular code could be easily improved by taking advantage of Numpy vectorization capabilities)

• Finally more advanced users can dig in the properties of the object M to build more sophisticated examples or even add their own properties to create even more advanced numerical schemes.

In both cases, assuming for simplicity that \(\Omega\) is a simple square that has been subdivided into 4x4 subsquares regularly, one would obtain at the end of the day

>>> print(U)
array([0.0234375, 0.0390625, 0.0390625, 0.0234375, 0.0390625, 0.0703125,
       0.0703125, 0.0390625, 0.0390625, 0.0703125, 0.0703125, 0.0390625,
       0.0234375, 0.0390625, 0.0390625, 0.0234375])

Installation

For the moment the only way to install the library is from the wheel or the sources available at

wheel_url

History behind the project

Documentation

Complete documentation of the library can be found here.