<commit_before>import numpy as np
import scipy.sparse as sprs
from decimal import Decimal as dc
from openpnm.algorithms import ReactiveTransport
from openpnm.utils import logging
logger = logging.getLogger(__name__)


class TransientReactiveTransport(ReactiveTransport):
    r"""
    A subclass of ReactiveTransport for transient/steady-state simulations

    Parameters
    ----------
    network : OpenPNM Network object
        The Network with which this algorithm is associated.

    project : OpenPNM Project object
        Either a Network or a Project must be specified.

    Notes
    -----

    This subclass performs steady and transient simulations of transport
    phenomena with reactions when source terms are added. It supports 3 time
    discretization schemes; 'steady' to perform a steady-state simulation, and
    'implicit' (fast, 1st order accurate) and 'cranknicolson' (slow, 2nd order
    accurate) both for transient simulations.
    """

    def __init__(self, settings={}, phase=None, **kwargs):
        def_set = {'phase': None,
                   't_initial': 0,
                   't_final': 10,
                   't_step': 0.1,
                   't_output': 1e+08,
                   't_tolerance': 1e-06,
                   'rxn_tolerance': 1e-05,
                   't_precision': 12,
                   't_scheme': 'implicit',
                   'gui': {'setup':        {'phase': None,
                                            'quantity': '',
                                            'conductance': '',
                                            't_initial': None,
                                            't_final': None,
                                            't_step': None,
                                            't_output': None,
                                            't_tolerance': None,
                                            't_precision': None,
                                            't_scheme': ''},
                           'set_IC':       {'values': None},
                           'set_rate_BC':  {'pores': None,
                                            'values': None},
                           'set_value_BC': {'pores': None,
                                            'values': None},
                           'set_source':   {'pores': None,
                                            'propname': ''}
                           }
                   }
        super().__init__(**kwargs)
        self.settings.update(def_set)
        self.settings.update(settings)
        self._A_steady = None  # Initialize the steady sys of eqs A matrix
        if phase is not None:
            self.setup(phase=phase)

    def setup(self, phase=None, quantity='', conductance='',
              t_initial=None, t_final=None, t_step=None, t_output=None,
              t_tolerance=None, t_precision=None, t_scheme='', **kwargs):
        r"""
        This method takes several arguments that are essential to running the
        algorithm and adds them to the settings

        Parameters
        ----------
        phase : OpenPNM Phase object
            The phase on which the algorithm is to be run. If no value is
            given, the existing value is kept.

        quantity : string
            The name of the physical quantity to be calcualted such as
            ``'pore.xxx'``.

        conductance : string
            The name of the pore-scale transport conductance values. These
            are typically calculated by a model attached to a *Physics* object
            associated with the given *Phase*. Example; ``'throat.yyy'``.

        t_initial : scalar, smaller than 't_final'
            The simulation's start time. The default value is 0.

        t_final : scalar, bigger than 't_initial'
            The simulation's end time. The default value is 10.

        t_step : scalar, between 't_initial' and 't_final'
            The simulation's time step. The default value is 0.1.

        t_output : scalar, ND-array, or list
            When 't_output' is a scalar, it is considered as an output interval
            to store transient solutions. The default value is 1e+08. Initial,
            final and steady-state (if reached) fields are always stored. If
            't_output' > 't_final', no transient data is stored. If 't_output'
            is not a multiple of 't_step', 't_output' will be approximated.
            When 't_output' is a list or ND-array, transient solutions
            corresponding to this list or array will be stored.

        output_times : list
            List of output times. The values in the list must be multiples of
            the time step 't_step'.

        t_tolerance : scalar
            Transient solver tolerance. The simulation stops (before reaching
            't_final') when the residual falls below 't_tolerance'. The
            default value is 1e-06. The 'residual' measures the variation from
            one time-step to another in the value of the 'quantity' solved for.

        rxn_tolerance : scalar
            Tolerance to achieve within each time step. The solver passes to
            next time step when 'residual' falls below 'rxn_tolerance'. The
            default value is 1e-05.

        t_precision : integer
            The time precision (number of decimal places).

        t_scheme : string
            The time discretization scheme. Three options available: 'steady'
            to perform a steady-state simulation, and 'implicit' (fast, 1st
            order accurate) and 'cranknicolson' (slow, 2nd order accurate) both
            for transient simulations. The default value is 'implicit'.

        Notes
        -----
        More settings can be adjusted in the presence of a non-linear source
        term such as under-relaxation.
        See the 'ReactiveTransport' class documentation for details.
        """
        if phase:
            self.settings['phase'] = phase.name
        if quantity:
            self.settings['quantity'] = quantity
        if conductance:
            self.settings['conductance'] = conductance
        if t_initial is not None:
            self.settings['t_initial'] = t_initial
        if t_final is not None:
            self.settings['t_final'] = t_final
        if t_step is not None:
            self.settings['t_step'] = t_step
        if t_output is not None:
            self.settings['t_output'] = t_output
        if t_tolerance is not None:
            self.settings['t_tolerance'] = t_tolerance
        if t_precision is not None:
            self.settings['t_precision'] = t_precision
        if t_scheme:
            self.settings['t_scheme'] = t_scheme
        self.settings.update(kwargs)

    def set_IC(self, values):
        r"""
        A method to set simulation initial conditions

        Parameters
        ----------
        values : ND-array or scalar
            Set the initial conditions using an 'Np' long array. 'Np' being
            the number of pores. If a scalar is given, the same value is
            imposed to all pores.
        """
        self[self.settings['quantity']] = values
        converted_array = self[self.settings['quantity']].astype('float64')
        self[self.settings['quantity']] = converted_array

    def _t_update_A(self):
        r"""
        A method to update 'A' matrix at each time step according to 't_scheme'
        """
        network = self.project.network
        Vi = network['pore.volume']
        dt = self.settings['t_step']
        s = self.settings['t_scheme']
        if (s == 'implicit'):
            f1, f2 = 1, 1
        elif (s == 'cranknicolson'):
            f1, f2 = 0.5, 1
        elif (s == 'steady'):
            f1, f2 = 1, 0
        # Compute A (operations involve conversion to 'csr')
        A = ((f2/dt) * sprs.coo_matrix.multiply(
            sprs.coo_matrix(np.reshape(Vi, (self.Np, 1)), shape=(self.Np,)),
            sprs.identity(self.Np, format='coo')) + f1 * self._A_steady)
        # Convert A to 'coo' format to apply BCs
        A = sprs.coo_matrix(A)
        self._A = A
        return A

    def _t_update_b(self):
        r"""
        A method to update 'b' array at each time step according to
        't_scheme' and the source term value
        """
        network = self.project.network
        phase = self.project.phases()[self.settings['phase']]
        Vi = network['pore.volume']
        dt = self.settings['t_step']
        s = self.settings['t_scheme']
        if (s == 'implicit'):
            f1, f2, f3 = 1, 1, 0
        elif (s == 'cranknicolson'):
            f1, f2, f3 = 0.5, 1, 0
        elif (s == 'steady'):
            f1, f2, f3 = 1, 0, 1
        x_old = self[self.settings['quantity']]
        b = (f2*(1-f1)*(-self._A_steady)*x_old +
             f2*(Vi/dt)*x_old +<commit_msg>Fix  sparse matrix ValueError (hopefully)<commit_after>import numpy as np
import scipy.sparse as sprs
from decimal import Decimal as dc
from openpnm.algorithms import ReactiveTransport
from openpnm.utils import logging
logger = logging.getLogger(__name__)


class TransientReactiveTransport(ReactiveTransport):
    r"""
    A subclass of ReactiveTransport for transient/steady-state simulations

    Parameters
    ----------
    network : OpenPNM Network object
        The Network with which this algorithm is associated.

    project : OpenPNM Project object
        Either a Network or a Project must be specified.

    Notes
    -----

    This subclass performs steady and transient simulations of transport
    phenomena with reactions when source terms are added. It supports 3 time
    discretization schemes; 'steady' to perform a steady-state simulation, and
    'implicit' (fast, 1st order accurate) and 'cranknicolson' (slow, 2nd order
    accurate) both for transient simulations.
    """

    def __init__(self, settings={}, phase=None, **kwargs):
        def_set = {'phase': None,
                   't_initial': 0,
                   't_final': 10,
                   't_step': 0.1,
                   't_output': 1e+08,
                   't_tolerance': 1e-06,
                   'rxn_tolerance': 1e-05,
                   't_precision': 12,
                   't_scheme': 'implicit',
                   'gui': {'setup':        {'phase': None,
                                            'quantity': '',
                                            'conductance': '',
                                            't_initial': None,
                                            't_final': None,
                                            't_step': None,
                                            't_output': None,
                                            't_tolerance': None,
                                            't_precision': None,
                                            't_scheme': ''},
                           'set_IC':       {'values': None},
                           'set_rate_BC':  {'pores': None,
                                            'values': None},
                           'set_value_BC': {'pores': None,
                                            'values': None},
                           'set_source':   {'pores': None,
                                            'propname': ''}
                           }
                   }
        super().__init__(**kwargs)
        self.settings.update(def_set)
        self.settings.update(settings)
        self._A_steady = None  # Initialize the steady sys of eqs A matrix
        if phase is not None:
            self.setup(phase=phase)

    def setup(self, phase=None, quantity='', conductance='',
              t_initial=None, t_final=None, t_step=None, t_output=None,
              t_tolerance=None, t_precision=None, t_scheme='', **kwargs):
        r"""
        This method takes several arguments that are essential to running the
        algorithm and adds them to the settings

        Parameters
        ----------
        phase : OpenPNM Phase object
            The phase on which the algorithm is to be run. If no value is
            given, the existing value is kept.

        quantity : string
            The name of the physical quantity to be calcualted such as
            ``'pore.xxx'``.

        conductance : string
            The name of the pore-scale transport conductance values. These
            are typically calculated by a model attached to a *Physics* object
            associated with the given *Phase*. Example; ``'throat.yyy'``.

        t_initial : scalar, smaller than 't_final'
            The simulation's start time. The default value is 0.

        t_final : scalar, bigger than 't_initial'
            The simulation's end time. The default value is 10.

        t_step : scalar, between 't_initial' and 't_final'
            The simulation's time step. The default value is 0.1.

        t_output : scalar, ND-array, or list
            When 't_output' is a scalar, it is considered as an output interval
            to store transient solutions. The default value is 1e+08. Initial,
            final and steady-state (if reached) fields are always stored. If
            't_output' > 't_final', no transient data is stored. If 't_output'
            is not a multiple of 't_step', 't_output' will be approximated.
            When 't_output' is a list or ND-array, transient solutions
            corresponding to this list or array will be stored.

        output_times : list
            List of output times. The values in the list must be multiples of
            the time step 't_step'.

        t_tolerance : scalar
            Transient solver tolerance. The simulation stops (before reaching
            't_final') when the residual falls below 't_tolerance'. The
            default value is 1e-06. The 'residual' measures the variation from
            one time-step to another in the value of the 'quantity' solved for.

        rxn_tolerance : scalar
            Tolerance to achieve within each time step. The solver passes to
            next time step when 'residual' falls below 'rxn_tolerance'. The
            default value is 1e-05.

        t_precision : integer
            The time precision (number of decimal places).

        t_scheme : string
            The time discretization scheme. Three options available: 'steady'
            to perform a steady-state simulation, and 'implicit' (fast, 1st
            order accurate) and 'cranknicolson' (slow, 2nd order accurate) both
            for transient simulations. The default value is 'implicit'.

        Notes
        -----
        More settings can be adjusted in the presence of a non-linear source
        term such as under-relaxation.
        See the 'ReactiveTransport' class documentation for details.
        """
        if phase:
            self.settings['phase'] = phase.name
        if quantity:
            self.settings['quantity'] = quantity
        if conductance:
            self.settings['conductance'] = conductance
        if t_initial is not None:
            self.settings['t_initial'] = t_initial
        if t_final is not None:
            self.settings['t_final'] = t_final
        if t_step is not None:
            self.settings['t_step'] = t_step
        if t_output is not None:
            self.settings['t_output'] = t_output
        if t_tolerance is not None:
            self.settings['t_tolerance'] = t_tolerance
        if t_precision is not None:
            self.settings['t_precision'] = t_precision
        if t_scheme:
            self.settings['t_scheme'] = t_scheme
        self.settings.update(kwargs)

    def set_IC(self, values):
        r"""
        A method to set simulation initial conditions

        Parameters
        ----------
        values : ND-array or scalar
            Set the initial conditions using an 'Np' long array. 'Np' being
            the number of pores. If a scalar is given, the same value is
            imposed to all pores.
        """
        self[self.settings['quantity']] = values
        converted_array = self[self.settings['quantity']].astype('float64')
        self[self.settings['quantity']] = converted_array

    def _t_update_A(self):
        r"""
        A method to update 'A' matrix at each time step according to 't_scheme'
        """
        network = self.project.network
        Vi = network['pore.volume']
        dt = self.settings['t_step']
        s = self.settings['t_scheme']
        if (s == 'implicit'):
            f1, f2 = 1, 1
        elif (s == 'cranknicolson'):
            f1, f2 = 0.5, 1
        elif (s == 'steady'):
            f1, f2 = 1, 0
        # Compute A (operations involve conversion to 'csr')
        temp = np.reshape(Vi, (self.Np, 1))
        A = ((f2/dt) * sprs.coo_matrix.multiply(
            sprs.coo_matrix(temp, shape=(self.Np, 1)),
            sprs.identity(self.Np, format='coo')) + f1 * self._A_steady)
        # Convert A to 'coo' format to apply BCs
        A = sprs.coo_matrix(A)
        self._A = A
        return A

    def _t_update_b(self):
        r"""
        A method to update 'b' array at each time step according to
        't_scheme' and the source term value
        """
        network = self.project.network
        phase = self.project.phases()[self.settings['phase']]
        Vi = network['pore.volume']
        dt = self.settings['t_step']
        s = self.settings['t_scheme']
        if (s == 'implicit'):
            f1, f2, f3 = 1, 1, 0
        elif (s == 'cranknicolson'):
            f1, f2, f3 = 0.5, 1, 0
        elif (s == 'steady'):
            f1, f2, f3 = 1, 0, 1
        x_old = self[self.settings['quantity']]
        b = (f2*(1-f1)*(-self._A_steady)*x_old +
             f2*(Vi/dt)*x_old +
{'fraction_non_alphanumeric': 0.09137542277339346, 'fraction_numerical': 0.008624577226606539, 'mean_word_length': 1.6471202626081767, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 5, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 48, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '3011011', 'n_tokens_mistral': 4811, 'n_tokens_neox': 4512, 'n_words': 1617}