Core Functionality

pypde.solvers.pde_solver(Q0, tf, L, F=None, B=None, S=None, boundaryTypes='transitive', cfl=0.9, order=2, ndt=100, flux='rusanov', stiff=True, nThreads=-1)

Solves PDEs of the following form:

\[\frac{\partial\mathbf{Q}}{\partial t}+\frac{\partial}{\partial x_{1}}\mathbf{F}_{1}\left(\mathbf{Q},\frac{\partial\mathbf{Q}}{\partial x_{1}},\ldots,\frac{\partial\mathbf{Q}}{\partial x_{n}}\right)+\cdots+\frac{\partial}{\partial x_{n}}\mathbf{F}_{n}\left(\mathbf{Q},\frac{\partial\mathbf{Q}}{\partial x_{1}},\ldots,\frac{\partial\mathbf{Q}}{\partial x_{n}}\right)+B_{1}\left(\mathbf{Q}\right)\frac{\partial\mathbf{Q}}{\partial x_{1}}+\cdots+B_{n}\left(\mathbf{Q}\right)\frac{\partial\mathbf{Q}}{\partial x_{n}}=\mathbf{S}\left(\mathbf{Q}\right)\]

or, more succinctly:

\[\frac{\partial\mathbf{Q}}{\partial t}+\nabla\mathbf{F}\left(\mathbf{Q},\nabla\mathbf{Q}\right)+B\left(\mathbf{Q}\right)\cdot\nabla\mathbf{Q}=\mathbf{S}\left(\mathbf{Q}\right)\]

where \(\mathbf{Q},\mathbf{F}_i,\mathbf{S}\) are vectors of \(n_{var}\) variables and \(B_i\) are matrices of shape \(\left(n_{var}\times n_{var}\right)\). Of course, \(\mathbf{F},\mathbf{B},\mathbf{S}\) can be 0.

Define \(\Omega=\left[0,L_1\right]\times\cdots\times\left[0,L_n\right]\). Given \(\mathbf{Q}\left(\mathbf{x},0\right)\) for \(\mathbf{x}\in\Omega\), pde_solver finds \(\mathbf{Q}\left(\mathbf{x},t\right)\) for any \(t>0\).

Taking integers \(m_1,\ldots,m_n>0\) we split \(\Omega\) into \(\left(m_1\times\cdots\times m_n\right)\) cells with volume \(dx_1 dx_2 \ldots dx_n\) where \(dx_i=\frac{L_i}{m_i}\).

Parameters:

• Q0 (ndarray) –

  An array with shape \(\left(m_1,\ldots m_n,n_{var}\right)\).

  Q0[i1, i2, ..., in, j] is equal to:

  \[\mathbf{Q}_{j+1}\left(\frac{\left(i_1+0.5\right) dx_1}{m_1},\ldots,\frac{\left(i_n+0.5\right) dx_n}{m_n}, 0\right)\]
• tf (double) – The final time at which to return the value of \(\mathbf{Q}\left(\mathbf{x},t\right)\).
• L (ndarray or list) – An array of length \(n\), L[i] is equal to \(L_{i+1}\).
• F (callable, optional) –

  The flux terms, with signature F(Q, DQ, d) -> ndarray, corresponding to \(\mathbf{F}_{d+1}\left(\mathbf{Q},\nabla\mathbf{Q}\right)\).

  If \(\mathbf{F}\) has no \(\nabla\mathbf{Q}\) dependence, the signature may be F(Q, d) -> ndarray, corresponding to \(\mathbf{F}_{d+1}\left(\mathbf{Q}\right)\).

  If \(n=1\) and \(\mathbf{F}\) has no \(\nabla\mathbf{Q}\) dependence, the signature may be F(Q) -> ndarray, corresponding to \(\mathbf{F}_1\left(\mathbf{Q}\right)\).

  • Q is a 1-D array with shape \(\left(n_{var},\right)\)
  • DQ is an 2-D array with shape \(\left(n,n_{var}\right)\) (as DQ[i] is equal to \(\frac{\partial\mathbf{Q}}{\partial x_{i}}\))
  • d is an integer (ranging from \(0\) to \(n-1\))
  • the returned array has shape \(\left(n_{var},\right)\)

• B (callable, optional) –

  The non-conservative terms, with signature B(Q, d) -> ndarray, corresponding to \(\mathbf{B}_{d+1}\left(\mathbf{Q}\right)\).

  If \(n=1\), the signature may be B(Q) -> ndarray, corresponding to \(\mathbf{B}_1\left(\mathbf{Q}\right)\).

  • Q is a 1-D array with shape \(\left(n_{var},\right)\)
  • d is an integer (ranging from \(0\) to \(n-1\))
  • the returned array has shape \(\left(n_{var},n_{var}\right)\)

• S (callable, optional) –

  The source terms, with signature B(Q, d) -> ndarray, corresponding to \(\mathbf{S}\left(\mathbf{Q}\right)\).

  • Q is a 1-D array with shape \(\left(n_{var},\right)\)
  • the returned array has shape \(\left(n_{var},\right)\)

• boundaryTypes (string or list, optional) –
  • If a string, must be one of 'transitive', 'periodic'. In this case, all boundaries will take the stated form.
  • If a list, must have length \(n\), containing strings taken from 'transitive', 'periodic'. In this case, the first element of the list describes the boundaries at \(x_1=0,L_1\), the second describes the boundaries at \(x_2=0,L_2\), etc.
• cfl (double, optional) – The CFL number: \(0<CFL<1\) (default 0.9).
• order (int, optional) – The order of the polynomial reconstructions used in the solver, must be an integer greater than 0 (default 2)
• ndt (int, optional) – The number of timesteps at which to return the value of the grid (default 100) e.g. if tf=5, and ndt=4 then the value of the grid will be returned at t=1.25,2.5,3.75,5.
• flux (string, optional) – The kind of flux to use at cell boundaries (default 'rusanov') Must be one of 'rusanov', 'roe', 'osher'.
• stiff (bool, optional) – Whether to use a stiff solver for the Discontinuous Galerkin step (default True) If the equations are stiff (i.e. the source terms are much larger than the other terms), then this is probably required to make the method converge. Otherwise, it can be turned off to improve speed.
• nThreads (int, optional) – The number of threads to use in the solver (default -1). If less than 1, the thread count will default to (number of cores) - 1.

Returns:

out – An array with shape \(\left(ndt,m_1,\ldots m_n,n_{var}\right)\).

Defining \(dt=\frac{tf}{ndt}\), then out[i, i1, i2, ..., in, j] is equal to:

\[\mathbf{Q}_{j+1}\left(\frac{\left(i_1+0.5\right) dx_1}{m_1},\ldots,\frac{\left(i_n+0.5\right) dx_n}{m_n}, \left(i+1\right) dt\right)\]

Return type:

ndarray