diode_devsim.py

"""Stage 2: the SAME diode as stage 1, built in DEVSIM.

Same geometry (2 um, junction at 1 um), same doping (NA=1e17 / ND=1e16), same
constant mobilities, same SRH lifetimes, same ni and Vt. DEVSIM's shipped
`simple_physics` package implements the same textbook models we derived in
01_pn_1d/EQUATIONS.md (identical SRH form, identical ohmic-contact BCs), so
after overriding its material defaults the two simulators should agree.
Pass criterion: I-V within a few percent of the from-scratch solver.
"""
import sys
from pathlib import Path

import numpy as np

HERE = Path(__file__).parent
sys.path.insert(0, str(HERE.parent))
from common.bootstrap import import_devsim

ds = import_devsim()
from devsim.python_packages import simple_physics as sp

import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt

# --- identical design to stage 1 ----------------------------------------
L, XJ = 2.0e-4, 1.0e-4
NA, ND = 1.0e17, 1.0e16
Vt = 0.025851997
q = 1.602176634e-19
device, region = "diode1d", "silicon"

ds.create_1d_mesh(mesh="m")
ds.add_1d_mesh_line(mesh="m", pos=0.0, ps=1e-7, tag="anode")
ds.add_1d_mesh_line(mesh="m", pos=XJ, ps=2e-8, tag="mid")     # refine junction
ds.add_1d_mesh_line(mesh="m", pos=L, ps=1e-7, tag="cathode")
ds.add_1d_contact(mesh="m", name="anode", tag="anode", material="metal")
ds.add_1d_contact(mesh="m", name="cathode", tag="cathode", material="metal")
ds.add_1d_region(mesh="m", material="Si", region=region, tag1="anode", tag2="cathode")
ds.finalize_mesh(mesh="m")
ds.create_device(mesh="m", device=device)

# material parameters — overridden to match pn1d.py exactly
sp.SetSiliconParameters(device, region, 300)
for name, value in (("Permittivity", 11.7 * 8.85418782e-14),
                    ("ElectronCharge", q),
                    ("n_i", 1.0e10),
                    ("V_t", Vt),
                    ("kT", Vt * q),
                    ("mu_n", 1400.0),
                    ("mu_p", 450.0),
                    ("taun", 1.0e-7),
                    ("taup", 1.0e-7)):
    ds.set_parameter(device=device, region=region, name=name, value=value)

# doping profile
ds.node_model(device=device, region=region, name="Acceptors",
              equation=f"{NA}*step({XJ}-x)")
ds.node_model(device=device, region=region, name="Donors",
              equation=f"{ND}*step(x-{XJ})")
ds.node_model(device=device, region=region, name="NetDoping",
              equation="Donors-Acceptors")

# equilibrium: Poisson only, then full drift-diffusion
sp.CreateSiliconPotentialOnly(device, region)
for contact in ("anode", "cathode"):
    ds.set_parameter(device=device, name=sp.GetContactBiasName(contact), value=0.0)
    sp.CreateSiliconPotentialOnlyContact(device, region, contact)

ds.solve(type="dc", absolute_error=1e10, relative_error=1e-10, maximum_iterations=50)

from devsim.python_packages.model_create import CreateSolution

for carrier, init in (("Electrons", "IntrinsicElectrons"), ("Holes", "IntrinsicHoles")):
    CreateSolution(device, region, carrier)   # also registers edge values
    ds.set_node_values(device=device, region=region, name=carrier, init_from=init)
sp.CreateSiliconDriftDiffusion(device, region)
for contact in ("anode", "cathode"):
    sp.CreateSiliconDriftDiffusionAtContact(device, region, contact)

ds.solve(type="dc", absolute_error=1e10, relative_error=1e-10, maximum_iterations=50)

x_ds = np.array(ds.get_node_model_values(device=device, region=region, name="x"))
psi_ds = np.array(ds.get_node_model_values(device=device, region=region, name="Potential"))


def anode_current():
    In = ds.get_contact_current(device=device, contact="anode", equation=sp.ece_name)
    Ip = ds.get_contact_current(device=device, contact="anode", equation=sp.hce_name)
    return In + Ip


# --- bias sweep, same points as stage 1 ---------------------------------
ref = np.load(HERE.parent / "01_pn_1d" / "pn1d_iv.npz")
iv = {}
for direction in (+1, -1):
    ds.set_parameter(device=device, name="anode_bias", value=0.0)
    ds.solve(type="dc", absolute_error=1e10, relative_error=1e-10, maximum_iterations=50)
    vs = sorted([v for v in ref["V"] if v * direction > 1e-12], key=abs)
    for V in vs:
        ds.set_parameter(device=device, name="anode_bias", value=float(V))
        ds.solve(type="dc", absolute_error=1e10, relative_error=1e-10,
                 maximum_iterations=50)
        iv[round(float(V), 4)] = anode_current()
iv[0.0] = 0.0

V_arr = np.array(sorted(iv))
J_ds = np.array([iv[v] for v in sorted(iv)])

# --- compare against the from-scratch solver ----------------------------
J_ref = ref["J"]
assert np.allclose(V_arr, ref["V"])
fw = V_arr >= 0.1 - 1e-12
ratio = J_ds[fw] / J_ref[fw]
worst = np.max(np.abs(ratio - 1.0))
print("V      J_devsim      J_scratch     ratio")
for v, jd, jr in zip(V_arr[fw], J_ds[fw], J_ref[fw]):
    print(f"{v:5.2f}  {jd:12.4e}  {jr:12.4e}  {jd/jr:6.3f}")
print(f"worst forward-bias deviation (V >= 0.1): {100*worst:.2f}%")

# equilibrium potential profiles (different meshes -> interpolate)
psi_ref_i = np.interp(x_ds, ref["x"], ref["psi0"])
dpsi = np.max(np.abs((psi_ds - psi_ds[0]) - (psi_ref_i - psi_ref_i[0])))
print(f"max equilibrium potential difference: {dpsi*1e3:.2f} mV")

assert worst < 0.05, "I-V disagrees with from-scratch solver by more than 5%"
assert dpsi < 0.005, "equilibrium potential disagrees by more than 5 mV"

# --- plot ---------------------------------------------------------------
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(11, 4.2))
pos = V_arr > 0.01
ax1.semilogy(V_arr[pos], np.abs(J_ds[pos]), "o", ms=5, mfc="none",
             color="tab:blue", label="DEVSIM")
ax1.semilogy(ref["V"][ref["V"] > 0.01], np.abs(J_ref[ref["V"] > 0.01]), "-",
             color="tab:orange", label="from-scratch solver (stage 1)")
neg = V_arr < -0.01
ax1.semilogy(-V_arr[neg], np.abs(J_ds[neg]), "s", ms=4, mfc="none",
             color="tab:blue")
ax1.semilogy(-ref["V"][ref["V"] < -0.01], np.abs(J_ref[ref["V"] < -0.01]), "--",
             color="tab:orange")
ax1.set_xlabel("|V| (V)"); ax1.set_ylabel("|J| (A/cm$^2$)")
ax1.set_title(f"same diode, two simulators (worst dev {100*worst:.1f}%)")
ax1.legend(); ax1.grid(alpha=0.3)

ax2.plot(V_arr[fw], 100 * (ratio - 1.0), "o-", ms=4)
ax2.axhline(0, color="k", lw=0.8)
ax2.set_xlabel("V (V)"); ax2.set_ylabel("deviation (%)")
ax2.set_title("DEVSIM vs from-scratch, point by point")
ax2.grid(alpha=0.3)

fig.suptitle("Stage 2: cross-validation — my solver vs DEVSIM", fontweight="bold")
fig.tight_layout()
out = HERE / "devsim_vs_scratch.png"
fig.savefig(out, dpi=130)
print(f"plot written: {out}")
print("CROSS-VALIDATION PASSED")