HeatSteadyState.py

def init(params: Parameters, T0: Temperature, q: HeatFlow, ts:float) -> np.array:
    """
    The time independent solution to heat equation.

    :param params: The model parameters.
    :param T0:     The surface temperature [°C].
    :param q:      The mantle heat flow density [mW/m^2].
    :param ts:     The time step in years.
    :returns:      The temperature field [°C].
    """
    H, k, n, d, _, _ = params.values()
    q = - abs(q) # Always negative => "blbuvzdorný"
                 
    # Setup the domain.
    dx = d / (n - 1)          # The node spacing.
    dt = ts * YEAR_IN_SECONDS # The time step in seconds.

    # The coeficient matrix + condistions.
    A = spdiags([np.ones(n), -2 * np.ones(n), np.ones(n)], [-1, 0, 1], n, n, 'csr')
    
    A[0, :2] = [1, 0]
    A[-1, -2:] = [2, -2]

    # The vector of constant terms + conditions.
    b = np.ones(n) * (-H * dx**2) / k
    b[0] = T0
    b[-1] = (q * 2 * dx - H * dx**2) / k
    
    return spsolve(A, b)