Diffusion Chained Tunneling Coupler Modal

dctm.py

# dctm.py - Diffusion Chained Tunneling Coupler Modal
import math

CNST_EP = 2.718281828459045
#___________________________________________________________________________________

def mat_mult_vec(M, v):
    
    return [sum(M[i][j]*v[j] for j in range(len(v))) for i in range(len(M))]
#___________________________________________________________________________________

def zero_matrix(r, c):
    
    return [[0.0]*c for _ in range(r)]
#___________________________________________________________________________________

class DiffusionChainModes:
    
    # Diffusion chain with M-modes per spatial point. Each segment a per-mode
    # D and a lists of len @M(mu_a). The transmission per-interface is setup
    # per-mode or a scalar broadcast. Tunneling provides a coupling matrix K
    # between the 2 interfaces; left bound(prescribed modal phi_left) & right
    # bounded, Dirichlet phi_right, also of length M*... optionals.

    def __init__(self, segments, M, T_between=None, phi_right=None,
                 phi_right_dirichlet=True, tunnel=None, eps_mu=0.0, Tmin=0.0):
        
        self.M, self.eps_mu, self.Tmin = int(M), float(eps_mu), float(Tmin)
        # Validate & copy the segments
        segs = []
        for s in segments:
            L = float(s['L'])
            Dlist = [float(x) for x in s['D_list']]
            mulist = [max(float(x), self.eps_mu) for x in s['mu_list']]
            if len(Dlist) != self.M or len(mulist) != self.M:
                raise ValueError('D_list & mu_list must have a same length M')
            segs.append({'L': L, 'D_list': Dlist, 'mu_list': mulist})
        N = len(segs)
        # Transmissions, list of N-1 entries each either list length M or scalar
        if T_between is None: T_between = [[1.0]*self.M for _ in range(max(0, N-1))]
        else:
            if len(T_between) != max(0, N-1):
                raise ValueError('@T_between length must be N-1')
            tmp = []
            for t in T_between:
                if isinstance(t, (int,float)):
                    tmp.append([max(float(t), self.Tmin) for _ in range(self.M)])
                else:
                    if len(t) != self.M:
                        raise ValueError('each @T_between entry must be list of length M or scalar')
                    tmp.append([max(float(x), self.Tmin) for x in t])
            T_between = tmp
        # Insert tunnel if applied; allows K(MxM) optionals
        if tunnel is not None:
            tparams, pos = tunnel.get('params'), int(tunnel.get('pos', 0))
            if pos < 0 or pos > N:
                raise ValueError('tunnel pos must be in 0..N')
            tL = float(tparams['L'])
            tD = [float(x) for x in tparams['D_list']]
            tmu = [max(float(x), self.eps_mu) for x in tparams['mu_list']]
            if len(tD) != self.M or len(tmu) != self.M:
                raise ValueError('tunnel D_list and mu_list must have a length M')
            tunnel_seg = {'L': tL, 'D_list': tD, 'mu_list': tmu}
            # Build the new segments with tunnel
            segs = segs[:pos]+[tunnel_seg]+segs[pos:]
            # Rebuild T_between so to match the new length; defaults 1.0 per-mode
            newT = [[1.0]*self.M for _ in range(len(segs)-1)]
            # Map original @T_between into newT
            for i in range(len(T_between)):
                if i < pos: newT[i] = T_between[i]
                else: newT[i+1] = T_between[i]
            # From your pos-controller override tunnel-specs provided
            if 'T_to_left' in tunnel and pos>0:
                tleft = tunnel['T_to_left']
                if isinstance(tleft, (int, float)):
                    newT[pos-1] = [max(float(tleft), self.Tmin)]*self.M
                else:
                    if len(tleft)!=self.M:
                        raise ValueError('mismatch T-Left length, check grid format')
                    newT[pos-1] = [max(float(x), self.Tmin) for x in tleft]
            if 'T_to_right' in tunnel and pos < N:
                tright = tunnel['T_to_right']
                if isinstance(tright, (int, float)):
                    newT[pos] = [max(float(tright), self.Tmin)]*self.M
                else:
                    if len(tright)!=self.M:
                        raise ValueError('mismatch T-Right length, check grid format')
                    newT[pos] = [max(float(x), self.Tmin) for x in tright]
            self.tunnel_pos, self.K, T_between = pos, tunnel.get('K'), newT
        else:
            self.tunnel_pos, self.K = None, None
        self.segments, self.N, self.T = segs, len(segs), T_between
        self.phi_right = [0.0]*self.M if phi_right is None else [float(x) for x in phi_right]
        self.phi_right_dirichlet = bool(phi_right_dirichlet)
        # Precompute all per-segment modal alphas
        self.alpha = []
        for seg in self.segments:
            a_list = []
            for m in range(self.M):
                Dm = seg['D_list'][m]; mum = seg['mu_list'][m]
                a = math.sqrt(mum/Dm) if mum > 0 else 0.0
                a_list.append(a)
            self.alpha.append(a_list)
#___________________________________________________________________________________

    # Basis functions per-segment/per-mode, same @ scalar cases with ext M optional
    
    def _phi_basis_val(self, i: int, m: int, s: float, AB: tuple) -> float:
        a = self.alpha[i][m]
        A, B = AB
        if a > 0.0:
            return A*(CNST_EP**(a*s))+B*(CNST_EP**(-a*s))
        else:
            return A+B*s

    def _dphi_basis_val(self, i: int, m: int, s: float, AB: tuple) -> float:
        
        a = self.alpha[i][m]
        A, B = AB
        if a > 0.0:
            return a*(A*(CNST_EP**(a*s))-B*(CNST_EP**(-a*s)))
        else:
            return B
#___________________________________________________________________________________

    def solve(self, phi_left: list) -> list:
        
        # @phi_left; list length M-modal values @ known left boundary
        
        # Solves linear system for all modal coefficients. Returns an
        # coefficents: Per-segment lists of (A_m,B_m) for m in 0..M-1
        
        M, N = self.M, self.N
        nvar = 2*N*M
        # Build a dense matrix & rhs
        A, b = [[0.0]*nvar for _ in range(nvar)], [0.0]*nvar
        
        def var_index(isA: bool, seg_i: int, mode_m: int) -> int:
            return ((seg_i*M+mode_m)*2)+(0 if isA else 1)

        # Left boundary, for each mode m, phi_0(0)_m == phi_left[m]
        for m in range(M):
            row = m
            # phi_0(0) with per-basis AB coefficients
            A[row][var_index(True, 0, m)] = self._phi_basis_val(0, m, 0.0, (1.0, 0.0))
            A[row][var_index(False, 0, m)] = self._phi_basis_val(0, m, 0.0, (0.0, 1.0))
            b[row] = float(phi_left[m])
        
        # Fill interface equations for each interface and each mode
        # *right boundary rows reserved at end, rows nvar-M..nvar-1
        row = M
        for i in range(N-1):
            L_i, T_im = float(self.segments[i]['L']), self.T[i]
            for m in range(M):
                # phi_i(L_i)_m*T_im[m]-phi_{i+1}(0)_m = 0
                Ai, Bi = var_index(True, i, m), var_index(False, i, m)
                Aj, Bj = var_index(True, i+1, m), var_index(False, i+1, m)
                A[row][Ai] = T_im[m]*self._phi_basis_val(i, m, L_i, (1.0, 0.0))
                A[row][Bi] = T_im[m]*self._phi_basis_val(i, m, L_i, (0.0, 1.0))
                A[row][Aj] = -self._phi_basis_val(i+1, m, 0.0, (1.0, 0.0))
                A[row][Bj] = -self._phi_basis_val(i+1, m, 0.0, (0.0, 1.0))
                b[row] = 0.0
                row+=1
                # And current continuity: D_i, m*dphi_i(L)-D_{i+1, m}*dphi_{i+1}(0) = 0
                A[row][Ai] = float(self.segments[i]['D_list'][m])*self._dphi_basis_val(i, m, L_i, (1.0, 0.0))
                A[row][Bi] = float(self.segments[i]['D_list'][m])*self._dphi_basis_val(i, m, L_i, (0.0, 1.0))
                A[row][Aj] = -float(self.segments[i+1]['D_list'][m])*self._dphi_basis_val(i+1, m, 0.0, (1.0, 0.0))
                A[row][Bj] = -float(self.segments[i+1]['D_list'][m])*self._dphi_basis_val(i+1, m, 0.0, (0.0, 1.0))
                b[row] = 0.0
                row+=1

        # If tunnel coupling K is provided, replace the continuity equation @
        # the tunnel interfaces with the coupled versions. Assume @tunnel_pos
        # is an index where tunnel segment was inserted.
        if self.tunnel_pos is not None and self.K is not None:
            p = self.tunnel_pos
            # Continuity eq between p-1 & p(left interface of tunnel) already set
            # above for each m. Need to replace that set with phi_{p-1}(L)*T_left
            # equals sum_n K_mn*phi_p(0)_n; find starting row corresponding to an
            # interface @i=p-1, compute where rows are. Each interface contribute
            # 2*M rows optional and/or interfaces are then in order i=0..N-2. Row
            # index for interface i optional will begin @ M+i*(2*M):
            start_row = M+(p-1)*(2*M)
            # Overwrite the optional 2*M rows starting at @start_row
            for rm in range(2*M):
                for col in range(nvar): A[start_row+rm][col] = 0.0
                b[start_row+rm] = 0.0
            # Fill coupled continuity(first M-rows) --> for each m in 0..M-1
            T_left = self.T[p-1]
            for m in range(M):
                r0 = start_row+m
                # phi_{p-1}(L)*T_left[m]
                Ai, Bi = var_index(True, p-1, m), var_index(False, p-1, m)
                A[r0][Ai] = T_left[m]*self._phi_basis_val(p-1, m, self.segments[p-1]['L'], (1.0, 0.0))
                A[r0][Bi] = T_left[m]*self._phi_basis_val(p-1, m, self.segments[p-1]['L'], (0.0, 1.0))
                # right sum_n K_mn*phi_p(0)_n
                for n in range(M):
                    # *right here Job'e(not job) does love you, semi answered->rO|《k_mn,I+[n^rO],..M-n
                    Aj, Bj = var_index(True, p, n), var_index(False, p, n)
                    k_mn = float(self.K[m][n])
                    A[r0][Aj]+=-k_mn*self._phi_basis_val(p, n, 0.0, (1.0, 0.0))
                    A[r0][Bj]+=-k_mn*self._phi_basis_val(p, n, 0.0, (0.0, 1.0))
            # And current continuity; next M rows, keep @ scalar per-mode continuity
            # Eliminates cross-coupling with currents from any optional M*? squaring
            for m in range(M):
                r1 = start_row+M+m
                Ai, Bi = var_index(True, p-1, m), var_index(False, p-1, m)
                Aj, Bj = var_index(True, p, m), var_index(False, p, m)
                A[r1][Ai] = float(self.segments[p-1]['D_list'][m])*self._dphi_basis_val(p-1, m, self.segments[p-1]['L'], (1.0, 0.0))
                A[r1][Bi] = float(self.segments[p-1]['D_list'][m])*self._dphi_basis_val(p-1, m, self.segments[p-1]['L'], (0.0, 1.0))
                A[r1][Aj] = -float(self.segments[p]['D_list'][m])*self._dphi_basis_val(p, m, 0.0, (1.0, 0.0))
                A[r1][Bj] = -float(self.segments[p]['D_list'][m])*self._dphi_basis_val(p, m, 0.0, (0.0, 1.0))
                b[r1] = 0.0
        # The current R-boundary rows, place Dirichlet | Neumann per-mode @ interface
        for m in range(M):
            r = nvar-M+m
            Ai, Bi = var_index(True, N-1, m), var_index(False, N-1, m)
            if self.phi_right_dirichlet:
                A[r][Ai] = self._phi_basis_val(N-1, m, self.segments[N-1]['L'], (1.0, 0.0))
                A[r][Bi] = self._phi_basis_val(N-1, m, self.segments[N-1]['L'], (0.0, 1.0))
                b[r] = float(self.phi_right[m])
            else:
                A[r][Ai] = float(self.segments[N-1]['D_list'][m])*self._dphi_basis_val(N-1, m, self.segments[N-1]['L'], (1.0, 0.0))
                A[r][Bi] = float(self.segments[N-1]['D_list'][m])*self._dphi_basis_val(N-1, m, self.segments[N-1]['L'], (0.0, 1.0))
                b[r] = 0.0
        x = self._solve_linear_system(A, b)
        # First, unpack the coefficients
        coeffs = []
        for i in range(N):
            seg_coeffs = []
            for m in range(M):
                seg_coeffs.append((x[var_index(True, i, m)], x[var_index(False, i, m)]))
            coeffs.append(seg_coeffs)
        self.coeffs = coeffs
        # Returns the modal current density into y: j_m = -D_{N-1, m}*dphi/ds @ s=L
        j_modal = []
        for m in range(M):
            AB = coeffs[N-1][m]
            dphi_at_L = self._dphi_basis_val(N-1, m, self.segments[N-1]['L'], AB)
            j = -float(self.segments[N-1]['D_list'][m])*dphi_at_L
            j_modal.append(j)
        self.j_modal = j_modal
        return coeffs
#___________________________________________________________________________________

    def _solve_linear_system(self, A: list, b: list) -> list:
        
        n = len(A)
        M, rhs = [row[:] for row in A], b[:]
        for k in range(n):
            piv, maxv = k, abs(M[k][k])
            for i in range(k+1, n):
                if abs(M[i][k]) > maxv:
                    maxv, piv = abs(M[i][k]), i
            if maxv == 0.0:
                raise RuntimeError('singular matrix')
            if piv != k:
                M[k], M[piv] = M[piv], M[k]
                rhs[k], rhs[piv] = rhs[piv], rhs[k]
            pv = M[k][k]
            for j in range(k, n): M[k][j]/=pv
            rhs[k]/=pv
            for i in range(k+1, n):
                fac = M[i][k]
                if fac == 0.0:
                    continue
                for j in range(k, n): M[i][j]-=fac*M[k][j]
                rhs[i]-=fac*rhs[k]
        x = [0.0]*n
        for i in range(n-1, -1, -1):
            s = rhs[i]
            for j in range(i+1, n): s-=M[i][j]*x[j]
            x[i] = s/M[i][i] if M[i][i] != 0 else s
        return x
#___________________________________________________________________________________
    
    def _apply_timed_tunnel(self, t, tau_func=None, omega_list=None, r_min=0.0, r_local=None, s_shape=8.0):
        
        # *If used to build the phased K&T --> done so before calling solve()
        
        # @t; an cascade progress in [0, 1]
        # @tau_func; function t->tau(default tau(t)=t)
        # @omega_list; list length M of the angular rates(rad per-unit tau)
        # @r_min; cutoff distance: if r_local < r_min then coupling suppressed
        # @r_local; list length N-1 giving local refraction distances per interface
        # @s_shape; steepness of smooth gate g(r) = sigmoid((r-r_min)*s_shape)
        
        # *Builds self._T_eff(list of per-interface lists length M or 2M if phases used)
        # *Builds self._K_eff(M×M real matrix or 2M×2M if phases used) tunneling
        # *If phases used, expands self.M -> 2*M temporarily for the solve
        
        if tau_func is None: tau_func = lambda tt: tt
        if omega_list is None: omega_list = [0.0]*self.M
        if r_local is None: r_local = [1e6]*(max(0, self.N-1))
        
        def gate(r):
            z = (r-r_min)*s_shape
            return 1.0/(1.0+CNST_EP**(-z))
            
        # Blend function for transmissions/tunnel's strength
        blend = lambda a, b: (1.0-t)*a+t*b
        self._T_eff = []
        for i in range(len(self.T)):
            # The list length @ M
            base = self.T[i]
            # If tunnel present & interface is at tunnel pos, use tunnel-specific target if given
            target = base[:]
            if self.tunnel_pos is not None and (i == self.tunnel_pos-1 or i == self.tunnel_pos):
                # *self.T already contains tunnel entries; assume target is final grids stored
                target = self.T[i]
            # Apply bridged-gate based on @r_local
            g = gate(r_local[i]) if i < len(r_local) else 1.0
            self._T_eff.append([g*blend(base[m], target[m]) for m in range(self.M)])
        # Merge/build @K_eff if tunnel & K is current present
        self._K_eff = None
        if self.tunnel_pos is not None and self.K is not None:
            g = gate(r_local[self.tunnel_pos-1] if self.tunnel_pos-1 < len(r_local) else 1.0)
            K_final = self.K
            K_blend = [[(1.0-t)*(1.0 if i==j else 0.0)+t*K_final[i][j] for j in range(self.M)] for i in range(self.M)]
            K_gated = [[g*K_blend[i][j] for j in range(self.M)] for i in range(self.M)]
            # Apply phase, diagonal(exp(i -> omega -> tau))
            tau, use_phase = tau_func(t), any(abs(w) > 0.0 for w in omega_list)
            if not use_phase: self._K_eff = K_gated
            else:
                # Expand to a real 2M×2M, K_complex = diagonal(e^{iωτ}) @ K_gated
                # Complex multiplication on it's vector [Re;Im] by real block matrix
                # If phase p_m= cos+i->sin -> [[Re(Kp), -Im(Kp)], [Im(Kp), Re(Kp)]]
                import math
                P = [complex(math.cos(w*tau), math.sin(w*tau)) for w in omega_list]
                # Compute Kp = K_gated*diag(P), apply phase on the 2nd factor
                Kp = [[ K_gated[i][j]*P[j] for j in range(self.M)] for i in range(self.M)]
                # Build the 2M×2M real matrix
                M2 = 2*self.M
                Keff = [[0.0]*M2 for _ in range(M2)]
                for i in range(self.M):
                    for j in range(self.M):
                        re, im = Kp[i][j].real, Kp[i][j].imag
                        Keff[i][j] = re
                        Keff[i][j+self.M] = -im
                        Keff[i+self.M][j] = im
                        Keff[i+self.M][j+self.M] = re
                self._K_eff = Keff
                # Need to expand @T_eff & an modal counts; set flag for solver to use
                # doubled M, expand @T_eff per-interface -> each modal scalar becomes
                # two identical entries yet with phase mixing applied in K only
                self._T_eff = [sum([[tval, tval] for tval in row ], []) for row in self._T_eff ]
                self._use_phase_expansion = True
                self._omega_list = omega_list
        else: self._use_phase_expansion = False
        # If no phase expansion, solver will use direct @T_eff & @K_eff
        self.T_work = getattr(self, '_T_eff', self.T)
        self.K_work = getattr(self, '_K_eff', self.K)
        return
#___________________________________________________________________________________

def test():
    
    # Two-mode chain: left seg, right seg
    segs = [{'L':1.0, 'D_list':[0.5, 0.3], 'mu_list':[1e-6, 0.02]},
            {'L':1.2, 'D_list':[0.4, 0.25], 'mu_list':[0.01, 0.03]}]
            
    # Insert tunnel between @segs pos=1 with mixed K
    tunnel = {'params': {'L':0.05, 'D_list':[0.2, 0.2], 'mu_list':[1e-4, 1e-4]},
              'pos': 1,
              'K': [[0.9, 0.1], [0.2, 0.8]],
              'T_to_left': [0.8, 0.7],
              'T_to_right': [0.7, 0.6]}
              
    M1 = 2
    M2 = M1*2
         
    chain = DiffusionChainModes(segs, M=M1, T_between=None, phi_right=[0.0, 0.0],
                                phi_right_dirichlet=True, tunnel=tunnel, eps_mu=1e-8,
                                Tmin=1e-9)
    # Tunnel params:
    t = 0.7
    tau_func = lambda tt: tt*0.8
    omega = [2.75, 4.75]   # per-mode angular rates
    r_local = [0.58, 0.72] # interface distance(one)
    r_min = 1.21
    s_shape = 42.4
    chain._apply_timed_tunnel(t, tau_func, omega, r_min, r_local, s_shape)
    
    chain.T = chain.T_work
    chain.K = chain.K_work
    coeffs = chain.solve(phi_left=[1.45, 0.45])
    print(f'\nModal currents into Y: {chain.j_modal}')
    
    # *To use in phase expansion. Above temporary, segments duplicated in
    #  modal dimension -> @D_list & mu_list then T_between = self.T_work,
    #  tunnel K = self.K_work(real 2M×2M), whereof phi_left duplicated as
    #  [Re,Im]=[phi_left, 0] from using the same class repeated solves.
test()