verification_cases.py

"""Provided small dense benchmarks assembled independently from input data."""
from pathlib import Path
from copy import deepcopy
import numpy as np
from scipy.linalg import expm
import yaml


def case(solver='implicit', dt=0.125):
    inputs = yaml.safe_load((Path(__file__).parent / 'inputs.yml').read_text())
    inputs.pop('plots', None)
    inputs['numerical'].update(solver=solver, **{'time step': dt})
    return inputs


def independent_system(inputs):
    """Return dense T, diagonal storage, and Q without using the model."""
    n = inputs['numerical']['number of grids']['x']
    r, w = inputs['reservoir'], inputs['fluid']['water']
    def vector(v):
        return np.broadcast_to(np.asarray(v, dtype=float), (n,)).copy()
    k, area, phi = [vector(r[key]) for key in
                    ['permeability', 'cross sectional area', 'porosity']]
    width = vector(inputs['numerical'].get('delta x', r['length'] / n))
    # Full-block resistances; interior faces join half of each resistance.
    resistance = w['viscosity'] * w['formation volume factor'] * width / (k * area)
    factor = inputs.get('conversion factor', 1.0)
    t, q = np.zeros((n, n)), np.zeros(n)
    for i in range(n - 1):
        edge = np.zeros(n); edge[i] = 1; edge[i + 1] = -1
        t += (2 * factor / (resistance[i] + resistance[i + 1])) * np.outer(edge, edge)
    for side, i in [('left', 0), ('right', n - 1)]:
        bc = inputs['boundary conditions'][side]
        if bc['type'] == 'prescribed pressure':
            face = 2 * factor / resistance[i]
            t[i, i] += face
            q[i] += face * bc['value']
        elif bc['type'] != 'prescribed flux' or bc['value'] != 0:
            raise ValueError('only pressure or zero-flux boundaries are supported')
    b = area * width * phi * w['compressibility'] / w['formation volume factor']
    return t, b, q


def exact_pressure(inputs, time):
    """Exact-in-time discrete solution, not a continuum exact solution."""
    t, b, q = independent_system(inputs)
    n = len(b)
    augmented = np.zeros((n + 1, n + 1))
    augmented[:n, :n] = -t / b[:, None]
    augmented[:n, n] = q / b
    initial = np.r_[np.full(n, inputs['initial conditions']['pressure']), 1.0]
    return (expm(time * augmented) @ initial)[:n]