<commit_before>from SimPEG.EM.Static.DC import Src
from SimPEG import Props
from SimPEG.Utils import sdiag
import scipy.sparse as sp
import numpy as np
from SimPEG.EM.Static.DC import Survey


class StreamingCurrents(Src.BaseSrc):

    L = None
    mesh = None
    # "Hydraulic Head (m)"
    # "Streaming current source (A/m^3)"
    # "Streaming current density (A/m^2)"
    modelType = None

    def __init__(self, rxList, **kwargs):
        Src.BaseSrc.__init__(self, rxList, **kwargs)
        if self.L is None:
            raise Exception("SP source requires cross coupling coefficient L")
        if self.mesh is None:
            raise Exception("SP source requires mesh")
        self.mesh.setCellGradBC("neumann")
        # TODO: remove below when confirmed
        # self.V = sp.block_diag([sdiag(self.mesh.vol)]*3)
        # self.Div = -self.mesh.cellGrad.T

    def eval(self, prob):
        """

            Computing source term using:

            - Hydraulic head: h
            - Cross coupling coefficient: L

            .. math::

                -\nabla \cdot \vec{j}^s = \nabla \cdot L \nabla \phi \\

        """
        if prob._formulation == 'HJ':
            if self.modelType == "Head":
                q = -self.Div*self.MfLiI*prob.Grad*prob.h
            elif self.modelType == "CurrentSource":
                q = prob.q
            elif self.modelType == "CurrentDensity":
                # q = -self.Div*prob.mesh.aveF2CCV.T* self.V *np.r_[prob.jsx, prob.jsy, prob.jsz]
                q = prob.Grad.T*prob.mesh.aveCC2FV*np.r_[prob.jsx, prob.jsy, prob.jsz]
            else:
                raise NotImplementedError()
        elif prob._formulation == 'EB':
            raise NotImplementedError()
        return q

    def evalDeriv(self, prob, v=None, adjoint=False):
        if prob._formulation == 'HJ':
            if adjoint:
                if self.modelType == "Head":
                    srcDeriv = - prob.hDeriv.T * prob.Grad.T * self.MfLiI.T * (self.Div.T * v)
                elif self.modelType == "CurrentSource":
                    srcDeriv = prob.qDeriv.T * v
                elif self.modelType == "CurrentDensity":
                    jsDeriv = sp.vstack((prob.jsxDeriv, prob.jsyDeriv, prob.jszDeriv))
                    # srcDeriv = - jsDeriv.T * prob.mesh.aveF2CCV * self.V * (self.Div.T*v)
                    srcDeriv = jsDeriv.T * prob.mesh.aveCC2FV.T * (prob.Grad*v)
                else:
                    raise NotImplementedError()
            else:
                if self.modelType == "Head":
                    srcDeriv = -self.Div*self.MfLiI*prob.Grad*(prob.hDeriv*v)
                elif self.modelType == "CurrentSource":
                    srcDeriv = prob.qDeriv * v
                elif self.modelType == "CurrentDensity":
                    jsDeriv = sp.vstack((prob.jsxDeriv, prob.jsyDeriv, prob.jszDeriv))
                    # srcDeriv = -self.Div * prob.mesh.aveF2CCV.T* self.V * (jsDeriv*v)
                    srcDeriv = prob.Grad.T * prob.mesh.aveCC2FV*(jsDeriv*v)
                else:
                    raise NotImplementedError()
        elif prob._formulation == 'EB':
            raise NotImplementedError()
        return srcDeriv

    @property
    def MfLi(self):
        """<commit_msg>use `prob.Grad` by letting streaming current, js, at boundary is zero.<commit_after>from SimPEG.EM.Static.DC import Src
from SimPEG import Props
from SimPEG.Utils import sdiag
import scipy.sparse as sp
import numpy as np
from SimPEG.EM.Static.DC import Survey


class StreamingCurrents(Src.BaseSrc):

    L = None
    mesh = None
    # "Hydraulic Head (m)"
    # "Streaming current source (A/m^3)"
    # "Streaming current density (A/m^2)"
    modelType = None

    def __init__(self, rxList, **kwargs):
        Src.BaseSrc.__init__(self, rxList, **kwargs)
        if self.L is None:
            raise Exception("SP source requires cross coupling coefficient L")
        if self.mesh is None:
            raise Exception("SP source requires mesh")
        # TODO: remove below when confirmed
        # self.mesh.setCellGradBC("neumann")
        # self.V = sp.block_diag([sdiag(self.mesh.vol)]*3)
        # self.Div = -self.mesh.cellGrad.T

    def eval(self, prob):
        """

            Computing source term using:

            - Hydraulic head: h
            - Cross coupling coefficient: L

            .. math::

                -\nabla \cdot \vec{j}^s = \nabla \cdot L \nabla \phi \\

        """
        if prob._formulation == 'HJ':
            if self.modelType == "Head":
                q = prob.Grad.T*self.MfLiI*prob.Grad*prob.h
            elif self.modelType == "CurrentSource":
                q = prob.q
            elif self.modelType == "CurrentDensity":
                # q = -self.Div*prob.mesh.aveF2CCV.T* self.V *np.r_[prob.jsx, prob.jsy, prob.jsz]
                q = prob.Grad.T*prob.mesh.aveCC2FV*np.r_[prob.jsx, prob.jsy, prob.jsz]
            else:
                raise NotImplementedError()
        elif prob._formulation == 'EB':
            raise NotImplementedError()
        return q

    def evalDeriv(self, prob, v=None, adjoint=False):
        if prob._formulation == 'HJ':
            if adjoint:
                if self.modelType == "Head":
                    srcDeriv = prob.hDeriv.T * prob.Grad.T * self.MfLiI.T * (prob.Grad * v)
                elif self.modelType == "CurrentSource":
                    srcDeriv = prob.qDeriv.T * v
                elif self.modelType == "CurrentDensity":
                    jsDeriv = sp.vstack((prob.jsxDeriv, prob.jsyDeriv, prob.jszDeriv))
                    # srcDeriv = - jsDeriv.T * prob.mesh.aveF2CCV * self.V * (self.Div.T*v)
                    srcDeriv = jsDeriv.T * prob.mesh.aveCC2FV.T * (prob.Grad*v)
                else:
                    raise NotImplementedError()
            else:
                if self.modelType == "Head":
                    srcDeriv = prob.Grad.T*self.MfLiI*prob.Grad*(prob.hDeriv*v)
                elif self.modelType == "CurrentSource":
                    srcDeriv = prob.qDeriv * v
                elif self.modelType == "CurrentDensity":
                    jsDeriv = sp.vstack((prob.jsxDeriv, prob.jsyDeriv, prob.jszDeriv))
                    # srcDeriv = -self.Div * prob.mesh.aveF2CCV.T* self.V * (jsDeriv*v)
                    srcDeriv = prob.Grad.T * prob.mesh.aveCC2FV*(jsDeriv*v)
                else:
                    raise NotImplementedError()
        elif prob._formulation == 'EB':
            raise NotImplementedError()
        return srcDeriv

    @property
    def MfLi(self):
        """
{'fraction_non_alphanumeric': 0.11706948640483383, 'fraction_numerical': 0.002719033232628399, 'mean_word_length': 1.6838265099310903, 'pattern_counts': {'":': 18, '<': 3, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '94231', 'n_tokens_mistral': 2111, 'n_tokens_neox': 2018, 'n_words': 458}