test.py

#!/usr/bin/env python

# Copyright 2018-2020 John T. Foster
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#    http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.
import unittest
import numpy as np
import scipy.sparse
from assignment12 import OneDimReservoir


class TestSolution(unittest.TestCase):

    def setUp(self):

        self.inputs = {
            'conversion factor': 6.33e-3,
            'fluid': {
                'water': {
                    'compressibility': 1e-6, #psi^{-1}
                    'viscosity': 1, #cp
                    'formation volume factor': 1,
                    },
            },
            'reservoir': {
                'permeability': 50, #mD
                'porosity': 0.2,
                'length': 10000, #ft
                'cross sectional area': 200000 #ft^2
            },
            'initial conditions': {
                'pressure': 1000 #psi
            },
            'boundary conditions': {
                'left': {
                    'type': 'prescribed pressure',
                    'value': 2000 #psi
                },
                'right': {
                    'type': 'prescribed flux',
                    'value': 0 #ft^3/day
                }
            },
            'numerical': {
                'solver': 'implicit',
                'number of grids': {'x': 4},
                'time step': 1, #day
                'number of time steps' : 3
            },
            'plots': {
                'frequency': 1
            }
        }

        return

    def test_compute_transmissibility(self):

        problem = OneDimReservoir(self.inputs)

        np.testing.assert_allclose(problem.compute_transmissibility(), 4000.0)

        return

    def test_compute_accumulation(self):

        problem = OneDimReservoir(self.inputs)

        np.testing.assert_allclose(problem.compute_accumulation(), 100.0)

        return

    def test_is_transmissiblity_matrix_sparse(self):

        problem = OneDimReservoir(self.inputs)

        assert scipy.sparse.issparse(problem.T)

        return

    def test_implicit_solve_one_step(self):

        implicit = OneDimReservoir(self.inputs)
        implicit.solve_one_step()
        np.testing.assert_allclose(implicit.get_solution(),
                                   np.array([1295.1463, 1051.1036, 1008.8921, 1001.7998]),
                                   atol=0.5)
        return

    def test_explicit_solve_one_step(self):

        self.inputs['numerical']['solver'] = 'explicit'

        explicit = OneDimReservoir(self.inputs)

        explicit.solve_one_step()

        np.testing.assert_allclose(explicit.get_solution(),
                               np.array([ 1506., 1000.,  1000.,  1000.004]),
                               atol=0.5)
        return

    def test_implicit_solve(self):

        implicit = OneDimReservoir(self.inputs)
        implicit.solve()
        np.testing.assert_allclose(implicit.get_solution(),
                                   np.array([1582.9, 1184.8, 1051.5, 1015.9]),
                                   atol=0.5)
        return

    def test_explicit_solve(self):

        self.inputs['numerical']['solver'] = 'explicit'

        explicit = OneDimReservoir(self.inputs)

        explicit.solve()

        np.testing.assert_allclose(explicit.get_solution(),
                               np.array([1689.8, 1222.3, 1032.4, 1000.0]),
                               atol=0.5)
        return


from copy import deepcopy
from unittest.mock import patch
from verification_cases import case, independent_system


class TestMixedMethod(unittest.TestCase):
    def model(self, theta, dt=.125):
        return OneDimReservoir(case({'mixed method': {'theta': theta}}, dt))

    def test_endpoints_match_existing_solvers(self):
        for theta, solver in [(0., 'implicit'), (1., 'explicit')]:
            with self.subTest(theta=theta):
                mixed = self.model(theta)
                existing = OneDimReservoir(case(solver, .125))
                for _ in range(4):
                    mixed.solve_one_step()
                    existing.solve_one_step()
                np.testing.assert_allclose(mixed.get_solution(), existing.get_solution(), atol=1e-7)

    def test_crank_nicolson(self):
        model = self.model(.5)
        matrix, source = independent_system(model.inputs)
        dt = model.delta_t
        expected = np.linalg.solve(np.eye(model.N) + .5*dt*matrix,
                                   (np.eye(model.N) - .5*dt*matrix) @ model.get_solution() + dt*source)
        model.solve_one_step()
        np.testing.assert_allclose(model.get_solution(), expected, atol=1e-7, rtol=1e-10)

    def test_intermediate_theta_repeated_steps(self):
        model = self.model(.3)
        matrix, source = independent_system(model.inputs)
        dt = model.delta_t
        expected = model.get_solution()
        for _ in range(5):
            expected = np.linalg.solve(np.eye(model.N) + .7*dt*matrix,
                                       (np.eye(model.N) - .3*dt*matrix) @ expected + dt*source)
            model.solve_one_step()
        np.testing.assert_allclose(model.get_solution(), expected, atol=1e-7, rtol=1e-10)

    def test_invalid_theta(self):
        for theta in [-.01, 1.01, float('nan'), float('inf'), '0.5', None, True, [], 1j]:
            with self.subTest(theta=theta):
                model = self.model(theta)
                before = model.get_solution()
                with self.assertRaises(ValueError):
                    model.solve_one_step()
                np.testing.assert_array_equal(model.get_solution(), before)

    def test_invalid_configuration(self):
        for solver in ['unknown', {}, {'mixed method': {}}, {'mixed method': None},
                       {'mixed method': {'theta': .5, 'typo': 1}},
                       {'mixed method': {'theta': .5}, 'other': 1}]:
            with self.subTest(solver=solver):
                model = OneDimReservoir(case(solver))
                with self.assertRaises(ValueError):
                    model.solve_one_step()

    def test_cg_status_and_tolerance(self):
        model = self.model(.5)
        before = model.get_solution()
        with patch('scipy.sparse.linalg.cg', return_value=(np.zeros(model.N), 1)) as cg:
            with self.assertRaises(RuntimeError):
                model.solve_one_step()
        self.assertIsNotNone(cg.call_args, 'Mixed method must call conjugate gradient')
        self.assertEqual(cg.call_args.kwargs, {'atol': 1e-8, 'rtol': 0.0})
        self.assertTrue(scipy.sparse.issparse(cg.call_args.args[0]))
        np.testing.assert_array_equal(model.get_solution(), before)

    def test_yaml_configuration(self):
        import tempfile
        from pathlib import Path
        import yaml
        inputs = case({'mixed method': {'theta': .5}})
        with tempfile.TemporaryDirectory() as folder:
            path = Path(folder) / 'case.yml'
            path.write_text(yaml.safe_dump(inputs))
            model = OneDimReservoir(str(path))
            model.solve_one_step()
            self.assertTrue(np.all(np.isfinite(model.get_solution())))


if __name__ == '__main__':
    unittest.main()