Source code for pyfvtool.pdesolver

import numpy as np

from scipy.sparse import csr_array
from scipy.sparse.linalg import spsolve
from scipy.sparse.linalg import use_solver
use_solver(useUmfpack=False) 
# For reproducibility, do not automatically use any installed `scikits.umfpack`
# solver. Always use the built-in SuperLU by default. In PyFVTool, the
# `scikits.umfpack.spsolve` solver (if installed) should be supplied via
# the `externalsolver` keyword argument of `solvePDE`, as is the case for
# the `pypardiso.spsolve` solver (and any other external solvers).

from .mesh import MeshStructure
from .cell import CellVariable
from .utilities import TrackedArray





def solvePDE(phi: CellVariable, eqnterms: list, 
              externalsolver = None) -> CellVariable:
    """
    Solve a PDE using the finite volume method
    
    This solver routine constructs the matrix equation from the FVM-discretized
    terms. These terms are provided by the user as a list each term being the 
    output of a prior call to the appropriate xxxTerm() function. 
    
    The terms provided by the user should not include any terms related to the 
    boundary conditions, as these are handled automatically by solvePDE.
    
    The CellVariable is updated with the solution values, and a
    reference to this one and the same variable is returned. If the 'old' 
    input CellVariable is to be conserved, it should be copied beforehand to
    a separate instance.

    Parameters
    ----------
    phi : CellVariable
        Solution variable subject to the boundary conditions represented by
        bcterm.
    eqnterms : list
        List of matrix equation terms generated by the xxxTerm functions.
        The elements of the list can either be a tuple (M, RHS), a simple 
        csr_array matrix (M) or a 1D np.ndarray (RHS). These will be added
        up to create both the M and RHS of the matrix equation to be solved.
    externalsolver : function, optional
        If specified, use an external sparse solver via a function call
        having the same interface as the default solver
        `scipy.sparse.linalg.spsolve`. If not specified, this default built-in
        SuperLU solver will be used.

    Raises
    ------
    TypeError
        If any of the provided terms is not conform to expectations.

    Returns
    -------
    CellVariable
        The updated CellVariable.

    """
    if externalsolver is None:
        solver = spsolve
    else:
        solver = externalsolver

    if phi.BCs.modified or phi.value.modified:
        phi.apply_BCs()
    
    # Construct BCs Term
    # Mbc, RHSbc = boundaryConditionsTerm(phi.BCs)
    # Retrieve pre-constructed BCs Term
    Mbc, RHSbc = phi._BCsTerm
    
    # Initialize overall cumulative matrix and right-hand side for 
    # matrix equation
    M = Mbc.copy() # need to copy, so that original 'bcterm' is protected
    RHS = RHSbc.copy() # need to copy, so that original 'bcterm' is protected
    
    # Cumulate all equation terms
    for term in eqnterms:
        if isinstance(term, tuple):
            Mterm, RHSterm = term
            if (Mterm.ndim != 2) or (RHSterm.ndim != 1):
                raise TypeError('Unknown term')
            M += Mterm
            RHS += RHSterm
        elif term.ndim == 1:
            RHS += term
        elif term.ndim == 2:
            M += term
        else:
            raise TypeError('Unknown term')

    # Solve!
    phi_new_values = solver(M, RHS)
    
    # Update phi
    phi._value = TrackedArray(np.reshape(phi_new_values, phi.domain.dims+2))
    phi.apply_BCs()
    
    return phi






def solveMatrixPDE(m: MeshStructure, M:csr_array, RHS: np.ndarray,
             externalsolver = None) -> CellVariable:
    """
    Solve the PDE discretized according to the finite volume method    
        
    This 'expert-level' solver routine uses the M matrix and RHS right-hand side 
    vector directly. These matrices should be constructed beforehand by 
    combining the different M matrices and RHS vectors generated using the
    xxxTerm routines, including those related to the boundary conditions.
    
    Returns a new CellVariable without changing the input ('old') CellVariable.
    
    Application of boundary condition terms and updating the related ghost
    cells should be handled entirely by the (expert-level) user.
    
    Parameters
    ----------
    m: MeshStructure
        Mesh structure
    M: csr_array
        Matrix of the linear system
    RHS: np.ndarray
        Right hand side of the linear system
    externalsolver: function (optional)
        If provided, use an external sparse solver via a function call
        having the same interface as the default built-in SuperLU solver
        scipy.sparse.linalg.spsolve.
    
    Returns
    -------
    phi: CellVariable
        Solution of the PDE (newly created CellVariable instance)
    
    """

    if externalsolver is None:
        solver = spsolve
    else:
        solver = externalsolver
    phi = solver(M, RHS)
    return CellVariable(m, np.reshape(phi, m.dims+2)) # BCs handled by user






def solveExplicitPDE(phi_old: CellVariable, 
                     dt: float, 
                     RHS: np.ndarray) -> CellVariable:
    """
    Solve the PDE using an explicit finite volume method.

    Parameters
    ----------
    phi_old: CellVariable
        Solution of the previous time step
    dt: float
        Time step
    RHS: np.ndarray
        Right hand side of the linear system

    
    Returns
    -------
    phi: CellVariable
        Solution of the PDE
    
    """
    if phi_old.BCs.modified or phi_old.value.modified:
        phi_old.apply_BCs()
    
    x = phi_old._value + dt*RHS.reshape(phi_old._value.shape)
    phi = CellVariable(phi_old.domain, 0.0, phi_old.BCs, 
                       BCsTerm_precalc = False)
    phi._value = TrackedArray(x)
    phi.apply_BCs()
    return phi