import abc
from cosapp.base import System, Port


class Voltage(Port):
    def setup(self):
        self.add_variable('V', unit='V')


class Intensity(Port):
    def setup(self):
        self.add_variable('I', unit='A')


class Dipole(System):
    """Abstract directional dipole model computing
    current from end values of electric potential.
    """
    def setup(self):
        self.add_input(Voltage, 'V_in')
        self.add_input(Voltage, 'V_out')

        self.add_output(Intensity, 'I')
        self.add_outward('deltaV', unit='V')

    def compute(self):
        self.deltaV = self.V_in.V - self.V_out.V
        self.compute_I()

    @abc.abstractmethod
    def compute_I(self) -> None:
        pass

Next, we define the following elements:

• Resistors

• Diodes

• Electric nodes

• Source

• Ground

from __future__ import annotations
from cosapp.base import System
import math


class Resistor(Dipole):

    def setup(self, R=1.0):
        super().setup()
        self.add_inward('R', R, desc="Resistance in Ohms")

    def compute_I(self):
        self.I.I = self.deltaV / self.R


class Diode(Dipole):
    """Regularized diode model

    The current intensity flowing through the diode is calculated based on

    $ I = I_s \\exp \\left( \\dfrac{V_{in} - V_{out}}{V_t} - 1 \\right) $
    """
    def setup(self):
        super().setup()
        self.add_inward('Is', 1e-15, desc='Saturation current in Amps')
        self.add_inward('Vt', 0.025875, desc='Thermal voltage in Volts')

    def compute_I(self):
        """Regularized diode model"""
        self.I.I = self.Is * math.exp(self.deltaV / self.Vt - 1.)


class Node(System):
    """Electric node model with `n_in` incoming and `n_out` outgoing currents.
    """
    def setup(self, n_in=1, n_out=1):
        self.add_property('n_in', max(1, int(n_in)))
        self.add_property('n_out', max(1, int(n_out)))

        incoming = tuple(
            self.add_input(Intensity, f"I_in{i}")
            for i in range(self.n_in)
        )
        outgoing = tuple(
            self.add_input(Intensity, f"I_out{i}")
            for i in range(self.n_out)
        )
        self.add_property('incoming_currents', incoming)
        self.add_property('outgoing_currents', outgoing)

        self.add_inward('V', 1.0, unit='V')
        self.add_outward('sum_I_in', 0., unit='A', desc='Sum of all incoming currents')
        self.add_outward('sum_I_out', 0., unit='A', desc='Sum of all outgoing currents')

        self.add_unknown('V')
        self.add_equation('sum_I_in == sum_I_out', name='current balance')

    def compute(self):
        self.sum_I_in = sum(current.I for current in self.incoming_currents)
        self.sum_I_out = sum(current.I for current in self.outgoing_currents)

    @classmethod
    def make(
        cls,
        name: str,
        parent: System,
        incoming: list[Dipole]=[],
        outgoing: list[Dipole]=[],
        pulling=None,
    ) -> Node:
        """Factory creating new node within `parent`, with
        appropriate connections with incoming and outgoing dipoles.
        """
        node = cls(name, n_in=len(incoming), n_out=len(outgoing))
        parent.add_child(node, pulling=pulling)

        for dipole, current in zip(incoming, node.incoming_currents):
            # print(dipole.name, type(dipole), dipole.I)
            parent.connect(dipole.I, current)
            parent.connect(dipole.V_out, node.inwards, 'V')

        for dipole, current in zip(outgoing, node.outgoing_currents):
            parent.connect(dipole.I, current)
            parent.connect(dipole.V_in, node.inwards, 'V')

        return node


class Source(System):

    def setup(self, I=0.1):
        self.add_inward('I', I, unit='A')
        self.add_output(Intensity, 'I_out', {'I': I})

    def compute(self):
        self.I_out.I = self.I


class Ground(System):

    def setup(self, V=0.):
        self.add_inward('V', V, unit='V')
        self.add_output(Voltage, 'V_out', {'V': V})

    def compute(self):
        self.V_out.V = self.V

Building the circuit

The circuit is built by reusing Node and Resistor elements customized thanks to keywords arguments that will be processed in the setup method of those classes.

class Circuit(System):

    def setup(self):
        R1 = self.add_child(Resistor('R1', R=100.), pulling={'V_out': 'Vg'})
        R2 = self.add_child(Resistor('R2', R=10000.))
        D1 = self.add_child(Diode('D1'), pulling={'V_out': 'Vg'})

        # Define nodes
        Node.make('n1',
            parent=self,
            pulling={'I_in0': 'I_in'},
            outgoing=[R1, R2],
        )
        Node.make('n2',
            parent=self,
            incoming=[R2],
            outgoing=[D1],
        )


p = System('model')
p.add_child(Source('source', I=0.1))
p.add_child(Ground('ground', V=0.0))
p.add_child(Circuit('circuit'))

p.connect(p.source.I_out, p.circuit.I_in)
p.connect(p.ground.V_out, p.circuit.Vg)

Solving the problem

from cosapp.drivers import NonLinearSolver

p.add_driver(NonLinearSolver('solver'))
p.run_drivers()

print(
    f"{p.circuit.n1.V = }",
    f"{p.circuit.n2.V = }",
    f"{p.circuit.R1.I = }",
    f"{p.circuit.R2.I = }",
    f"{p.circuit.D1.I = }",
    sep="\n",
)

print(f"\nSanity check: {p.circuit.R1.I.I + p.circuit.D1.I.I = } ({p.source.I} expected)")

p.circuit.n1.V = 9.908302820430912
p.circuit.n2.V = 0.7385848635220567
p.circuit.R1.I = Intensity: {'I': 0.09908302820430911}
p.circuit.R2.I = Intensity: {'I': 0.0009169717956908855}
p.circuit.D1.I = Intensity: {'I': 0.0009169717956908876}

Sanity check: p.circuit.R1.I.I + p.circuit.D1.I.I = 0.1 (0.1 expected)

Show mathematical problem:

p.drivers['solver'].problem

Unknowns [2]
  circuit.n1_V = 9.908302820430912
  circuit.n2_V = 0.7385848635220567
Equations [2]
  circuit.n1: current balance := 0.0
  circuit.n2: current balance := -2.0599841277224584e-18