tcad-devsim_triac/mos_transfer_single_temp.py

# mos_transfer_single_temp.py
# 2D MOSFET Ids-Vgs Transfer Curve Sweep at a Specific Temperature
# Sweeps Vgs from 0.0V to 10.0V under constant Vds = 1.0V

import os
import sys
import glob
import gc
import time
import argparse
import pickle
import numpy as np

# Limit thread count to avoid system resource starvation
os.environ["OMP_NUM_THREADS"] = "4"
os.environ["MKL_NUM_THREADS"] = "4"
os.environ["TBB_NUM_THREADS"] = "4"
os.environ["OPENBLAS_NUM_THREADS"] = "4"
os.environ["MKL_DISABLE_FAST_MM"] = "1"

# Load Intel MKL runtime library if available
mkl_libs = glob.glob(os.path.join(os.path.dirname(sys.executable), "../lib/libmkl_rt.so.*"))
if mkl_libs:
    os.environ["DEVSIM_MATH_LIBS"] = mkl_libs[0]

import devsim

# Use MKL PARDISO solver
devsim.set_parameter(name="solver_type", value="pardiso")

# Parse arguments
parser = argparse.ArgumentParser(description="Run transfer characteristics sweep at a single temperature.")
parser.add_argument("--temp", type=float, required=True, help="Temperature in Celsius.")
parser.add_argument("--out_dir", type=str, required=True, help="Output directory path.")
parser.add_argument("--vds", type=float, default=1.0, help="Drain-Source Voltage Vds.")
args = parser.parse_args()

temp_C = args.temp
T_kelvin = temp_C + 273.15
OUT_DIR = os.path.join(args.out_dir, "")
os.makedirs(OUT_DIR, exist_ok=True)

# Default device directory to LDMOS
DEV_DIR = os.environ.get("DEV_DIR", "devices/LDMOS")
sys.path.insert(0, os.path.abspath(DEV_DIR))

# Import geometry parameters and physics creators
from device_config import *
from physics.model_create import *
from physics.new_physics import *

device = "device_2d"

# --- 1. Load Mesh ---
mesh_file = os.path.join(DEV_DIR, "device_2d.msh")
print(f"[{temp_C} C] Loading mesh: {mesh_file}...")
devsim.create_gmsh_mesh(mesh=device, file=mesh_file)
devsim.add_gmsh_region(mesh=device, gmsh_name="Silicon", region="Silicon", material="Silicon")
devsim.add_gmsh_region(mesh=device, gmsh_name="Oxide", region="Oxide", material="Oxide")
devsim.add_gmsh_region(mesh=device, gmsh_name="Molding", region="Molding", material="Molding")

# Add contacts
devsim.add_gmsh_contact(mesh=device, gmsh_name="contR_Si", name="contR_Si", region="Silicon", material="metal")
devsim.add_gmsh_contact(mesh=device, gmsh_name="contL_Si", name="contL_Si", region="Silicon", material="metal")
devsim.add_gmsh_contact(mesh=device, gmsh_name="gate_Ox", name="gate_Ox", region="Oxide", material="metal")

devsim.add_gmsh_contact(mesh=device, gmsh_name="gate_Mold", name="gate_Mold", region="Molding", material="metal")
devsim.add_gmsh_contact(mesh=device, gmsh_name="contR_Mold", name="contR_Mold", region="Molding", material="metal")
devsim.add_gmsh_contact(mesh=device, gmsh_name="contL_Mold", name="contL_Mold", region="Molding", material="metal")

# Interfaces
devsim.add_gmsh_interface(mesh=device, gmsh_name="Si_Ox_Interface", name="Si_Ox", region0="Silicon", region1="Oxide")
devsim.add_gmsh_interface(mesh=device, gmsh_name="Ox_Mold_Interface", name="Ox_Mold", region0="Oxide", region1="Molding")

devsim.finalize_mesh(mesh=device)
devsim.create_device(mesh=device, device=device)

# --- 2. Setup Doping (PCAD) ---
try:
    from device_pcad_config import apply_pcad_doping_2d
    apply_pcad_doping_2d(device, region="Silicon")
    print(f"[{temp_C} C] Applied PCAD doping successfully.")
except ImportError:
    print(f"[{temp_C} C] Error: device_pcad_config not found. Make sure DEV_DIR points to the correct LDMOS folder.")
    sys.exit(1)

# --- 3. Initialize Electrostatic Solution (Poisson) ---
CreateSolution(device, "Silicon", "Potential")
CreateSiliconPotentialOnly(device, "Silicon")

# Set temperature T on both device level and region level to prevent SetSiliconParameters overriding it back to 300K
devsim.set_parameter(device=device, name="T", value=T_kelvin)
devsim.set_parameter(device=device, region="Silicon", name="T", value=T_kelvin)


def CreateOxidePotentialOnly(device, region):
    if not InNodeModelList(device, region, "Potential"):
        CreateSolution(device, region, "Potential")
    devsim.set_parameter(device=device, region=region, name="Permittivity", value=3.9 * 8.85e-14)
    efield = "(Potential@n0 - Potential@n1)*EdgeInverseLength"
    CreateEdgeModel(device, region, "EField", efield)
    CreateEdgeModelDerivatives(device, region, "EField", efield, "Potential")
    dfield = "Permittivity*EField"
    CreateEdgeModel(device, region, "PotentialEdgeFlux", dfield)
    CreateEdgeModelDerivatives(device, region, "PotentialEdgeFlux", dfield, "Potential")
    devsim.equation(device=device, region=region, name="PotentialEquation", variable_name="Potential",
                    edge_model="PotentialEdgeFlux", variable_update="default", min_error=1e-3)

CreateOxidePotentialOnly(device, "Oxide")

def CreateMoldingPotentialOnly(device, region):
    if not InNodeModelList(device, region, "Potential"):
        CreateSolution(device, region, "Potential")
    devsim.set_parameter(device=device, region=region, name="Permittivity", value=4.0 * 8.85e-14)
    efield = "(Potential@n0 - Potential@n1)*EdgeInverseLength"
    CreateEdgeModel(device, region, "EField", efield)
    CreateEdgeModelDerivatives(device, region, "EField", efield, "Potential")
    dfield = "Permittivity*EField"
    CreateEdgeModel(device, region, "PotentialEdgeFlux", dfield)
    CreateEdgeModelDerivatives(device, region, "PotentialEdgeFlux", dfield, "Potential")
    devsim.equation(device=device, region=region, name="PotentialEquation", variable_name="Potential",
                    edge_model="PotentialEdgeFlux", variable_update="default", min_error=1e-3)

CreateMoldingPotentialOnly(device, "Molding")

# Interfaces continuous potential
def CreateContinuousPotentialInterface(device, interface):
    model_name = CreateContinuousInterfaceModel(device, interface, "Potential")
    devsim.interface_equation(device=device, interface=interface, name="PotentialEquation",
                              interface_model=model_name, type="continuous")

CreateContinuousPotentialInterface(device, "Si_Ox")
CreateContinuousPotentialInterface(device, "Ox_Mold")

# Set contact potential boundary equations
silicon_contacts = ["contL_Si", "contR_Si"]
for c in silicon_contacts:
    devsim.set_parameter(device=device, name=GetContactBiasName(c), value=0.0)
    CreateSiliconPotentialOnlyContact(device, "Silicon", c)

def CreateOxidePotentialOnlyContact(device, region, contact):
    contact_bias = GetContactBiasName(contact)
    contact_model = f"Potential - {contact_bias}"
    contact_model_name = f"{contact}_bc"
    CreateContactNodeModel(device, contact, contact_model_name, contact_model)
    CreateContactNodeModelDerivative(device, contact, contact_model_name, contact_model, "Potential")
    devsim.contact_equation(device=device, contact=contact, name="PotentialEquation",
                            node_model=contact_model_name, edge_charge_model="PotentialEdgeFlux")

for c in ["gate_Ox"]:
    devsim.set_parameter(device=device, name=GetContactBiasName(c), value=0.0)
    CreateOxidePotentialOnlyContact(device, "Oxide", c)

def CreateMoldingPotentialOnlyContact(device, region, contact):
    contact_bias = GetContactBiasName(contact)
    contact_model = f"Potential - {contact_bias}"
    contact_model_name = f"{contact}_bc"
    CreateContactNodeModel(device, contact, contact_model_name, contact_model)
    CreateContactNodeModelDerivative(device, contact, contact_model_name, contact_model, "Potential")
    devsim.contact_equation(device=device, contact=contact, name="PotentialEquation",
                            node_model=contact_model_name, edge_charge_model="PotentialEdgeFlux")

for c in ["gate_Mold", "contR_Mold", "contL_Mold"]:
    devsim.set_parameter(device=device, name=GetContactBiasName(c), value=0.0)
    CreateMoldingPotentialOnlyContact(device, "Molding", c)

# Solve initial zero-bias Poisson
print(f"[{temp_C} C] Solving initial Poisson at thermal equilibrium...")
devsim.solve(type="dc", absolute_error=1.0, relative_error=1e-10, maximum_iterations=50)
print(f"[{temp_C} C] Initial Poisson converged.")

# --- 4. Setup Drift-Diffusion (Poisson + Continuity) ---
CreateSolution(device, "Silicon", "Electrons")
CreateSolution(device, "Silicon", "Holes")
devsim.set_node_values(device=device, region="Silicon", name="Electrons", init_from="IntrinsicElectrons")
devsim.set_node_values(device=device, region="Silicon", name="Holes", init_from="IntrinsicHoles")

print(f"[{temp_C} C] Redefining equilibrium models to prevent high-bias overflow...")
devsim.node_model(device=device, region="Silicon", name="IntrinsicElectrons", equation="Electrons")
devsim.node_model(device=device, region="Silicon", name="IntrinsicElectrons:Potential", equation="0")
devsim.node_model(device=device, region="Silicon", name="IntrinsicElectrons:Electrons", equation="1")
devsim.node_model(device=device, region="Silicon", name="IntrinsicElectrons:Holes", equation="0")

devsim.node_model(device=device, region="Silicon", name="IntrinsicHoles", equation="Holes")
devsim.node_model(device=device, region="Silicon", name="IntrinsicHoles:Potential", equation="0")
devsim.node_model(device=device, region="Silicon", name="IntrinsicHoles:Electrons", equation="0")
devsim.node_model(device=device, region="Silicon", name="IntrinsicHoles:Holes", equation="1")

devsim.node_model(device=device, region="Silicon", name="PotentialIntrinsicCharge", equation="0")
devsim.node_model(device=device, region="Silicon", name="PotentialIntrinsicCharge:Potential", equation="0")

opts = CreateAroraMobilityLF(device, "Silicon")
CreateSiliconDriftDiffusion(device, "Silicon", **opts)
CreateAvalancheGeneration(device, "Silicon", opts['Jn'], opts['Jp'])

for c in silicon_contacts:
    CreateSiliconDriftDiffusionContact(device, "Silicon", c, opts['Jn'], opts['Jp'])

# Set up positive update continuity equations for stable sweeping
devsim.equation(device=device, region="Silicon", name="ElectronContinuityEquation", variable_name="Electrons",
                time_node_model="NCharge", edge_model=opts['Jn'], edge_volume_model="", 
                variable_update="positive", node_model="ElectronGeneration", min_error=1e5)
devsim.equation(device=device, region="Silicon", name="HoleContinuityEquation", variable_name="Holes",
                time_node_model="PCharge", edge_model=opts['Jp'], edge_volume_model="", 
                variable_update="positive", node_model="HoleGeneration", min_error=1e5)
devsim.equation(device=device, region="Silicon", name="PotentialEquation", variable_name="Potential",
                node_model="PotentialNodeCharge", edge_model="DField", variable_update="log_damp", min_error=3e-3)


# Helper functions to save and restore simulator state
def save_state(device):
    state = {}
    for region in ["Silicon", "Oxide", "Molding"]:
        state[region] = {
            "Potential": list(devsim.get_node_model_values(device=device, region=region, name="Potential"))
        }
    state["Silicon"]["Electrons"] = list(devsim.get_node_model_values(device=device, region="Silicon", name="Electrons"))
    state["Silicon"]["Holes"] = list(devsim.get_node_model_values(device=device, region="Silicon", name="Holes"))
    return state

def restore_state(device, state):
    for region in ["Silicon", "Oxide", "Molding"]:
        devsim.set_node_values(device=device, region=region, name="Potential", values=state[region]["Potential"])
    devsim.set_node_values(device=device, region="Silicon", name="Electrons", values=state["Silicon"]["Electrons"])
    devsim.set_node_values(device=device, region="Silicon", name="Holes", values=state["Silicon"]["Holes"])

def robust_solve(max_iters=40):
    # Stage 1: Pre-conditioning solve (loose tolerance)
    try:
        devsim.solve(type="dc", absolute_error=1e13, relative_error=1e-1, charge_error=1e13, maximum_iterations=10)
    except devsim.error:
        pass
    # Stage 2: Precision Newton solve
    devsim.solve(type="dc", absolute_error=1e12, relative_error=3e-2, charge_error=1e13, maximum_iterations=max_iters)

# Solve zero-bias Drift-Diffusion state
print(f"[{temp_C} C] Solving initial zero-bias Drift-Diffusion...")
robust_solve(max_iters=40)
print(f"[{temp_C} C] Initial Drift-Diffusion converged.")

zero_bias_state = save_state(device)

# --- 5. Phase 1: Ramp Vds to target voltage at Vgs = 0.0V ---
vds_target = args.vds
print(f"[{temp_C} C] Ramping Vds from 0.0V to {vds_target:.3f}V (Vgs = 0.0V)...")
vds_curr = 0.0
vds_step = 0.1

# Sync gate contacts at 0V
for c in ["gate_Ox", "gate_Mold"]:
    devsim.set_parameter(device=device, name=f"{c}_bias", value=0.0)

state_vds = save_state(device)

while vds_curr < vds_target:
    vds_next = min(vds_curr + vds_step, vds_target)
    for c in ["contR_Si", "contR_Mold"]:
        devsim.set_parameter(device=device, name=f"{c}_bias", value=vds_next)
    
    try:
        robust_solve(max_iters=40)
        vds_curr = vds_next
        state_vds = save_state(device)
        vds_step = min(vds_step * 1.5, 0.2)
    except devsim.error:
        # Revert and shrink step
        for c in ["contR_Si", "contR_Mold"]:
            devsim.set_parameter(device=device, name=f"{c}_bias", value=vds_curr)
        restore_state(device, state_vds)
        vds_step *= 0.5
        if vds_step < 1e-4:
            print(f"[{temp_C} C] Error: Vds ramping failed to converge at {vds_next:.3f}V. Aborting.")
            sys.exit(1)

print(f"[{temp_C} C] Vds successfully ramped to 1.0V.")

# --- 6. Phase 2: Sweep Vgs from 0.0V to 10.0V at constant Vds = 1.0V ---
print(f"[{temp_C} C] Sweeping Vgs from 0.0V to 10.0V at constant Vds = 1.0V...")
vgs_target = 10.0
vgs_curr = 0.0
vgs_step = 0.05  # Start with fine step to capture the subthreshold region precisely

vgs_list = [vgs_curr]
# Extract current at Vgs = 0
i_n = devsim.get_contact_current(device=device, contact="contR_Si", equation="ElectronContinuityEquation")
i_p = devsim.get_contact_current(device=device, contact="contR_Si", equation="HoleContinuityEquation")
ids_list = [i_n + i_p]

g_av = devsim.get_edge_model_values(device=device, region="Silicon", name="AvalancheGeneration")
edge_vol = devsim.get_edge_model_values(device=device, region="Silicon", name="EdgeNodeVolume")
iav_list = [2.0 * sum(g * v for g, v in zip(g_av, edge_vol))]

state_vgs = save_state(device)

while vgs_curr < vgs_target:
    vgs_next = min(vgs_curr + vgs_step, vgs_target)
    for c in ["gate_Ox", "gate_Mold"]:
        devsim.set_parameter(device=device, name=f"{c}_bias", value=vgs_next)
        
    try:
        robust_solve(max_iters=40)
        vgs_curr = vgs_next
        state_vgs = save_state(device)
        
        # Extract current
        i_n = devsim.get_contact_current(device=device, contact="contR_Si", equation="ElectronContinuityEquation")
        i_p = devsim.get_contact_current(device=device, contact="contR_Si", equation="HoleContinuityEquation")
        total_curr = i_n + i_p
        
        # Integrated avalanche current
        g_av = devsim.get_edge_model_values(device=device, region="Silicon", name="AvalancheGeneration")
        edge_vol = devsim.get_edge_model_values(device=device, region="Silicon", name="EdgeNodeVolume")
        total_av_curr = 2.0 * sum(g * v for g, v in zip(g_av, edge_vol))
        
        vgs_list.append(vgs_curr)
        ids_list.append(total_curr)
        iav_list.append(total_av_curr)
        
        print(f"  Vgs = {vgs_curr:.3f}V, Ids = {total_curr:.4e}A, I_av = {total_av_curr:.4e}A (step: {vgs_step:.3f}V)")
        
        # Adaptive step size
        vgs_step = min(vgs_step * 1.5, 0.2)
        gc.collect()
        
    except devsim.error:
        # Revert and shrink step
        for c in ["gate_Ox", "gate_Mold"]:
            devsim.set_parameter(device=device, name=f"{c}_bias", value=vgs_curr)
        restore_state(device, state_vgs)
        vgs_step *= 0.3
        if vgs_step < 1e-4:
            print(f"[{temp_C} C] Error: Vgs sweep failed to converge at Vgs = {vgs_next:.3f}V. Aborting.")
            break

# Save CSV output for this temperature
csv_file = os.path.join(OUT_DIR, f"transfer_vgs_sweep_{temp_C:+.1f}C.csv")
np.savetxt(csv_file, np.column_stack((vgs_list, ids_list, iav_list)), 
           header="Vgs(V),Ids(A),Iav_eval(A)", delimiter=",", comments="")
print(f"[{temp_C} C] Saved transfer data to {csv_file}")

# Save final field visualization for ParaView
visual_variables = {"x", "y", "Potential", "Electrons", "Holes", "NetDoping", "Emag", "logEmag", "Jmag", "logJmag"}
for reg in ["Silicon"]:
    devsim.element_from_edge_model(edge_model=opts['Jn'], device=device, region=reg)
    devsim.element_from_edge_model(edge_model=opts['Jp'], device=device, region=reg)
    devsim.element_model(device=device, region=reg, name="Jn_x", equation=f"{opts['Jn']}_x")
    devsim.element_model(device=device, region=reg, name="Jn_y", equation=f"{opts['Jn']}_y")
    devsim.element_model(device=device, region=reg, name="Jp_x", equation=f"{opts['Jp']}_x")
    devsim.element_model(device=device, region=reg, name="Jp_y", equation=f"{opts['Jp']}_y")
    devsim.element_model(device=device, region=reg, name="Jmag", equation="((Jn_x + Jp_x)^2 + (Jn_y + Jp_y)^2)^(0.5)")
    devsim.element_model(device=device, region=reg, name="logJmag", equation="log(Jmag + 1e-20) / log(10.0)")

for reg in ["Oxide", "Molding"]:
    devsim.element_model(device=device, region=reg, name="Jn_x", equation="0.0")
    devsim.element_model(device=device, region=reg, name="Jn_y", equation="0.0")
    devsim.element_model(device=device, region=reg, name="Jp_x", equation="0.0")
    devsim.element_model(device=device, region=reg, name="Jp_y", equation="0.0")
    devsim.element_model(device=device, region=reg, name="Jmag", equation="0.0")
    devsim.element_model(device=device, region=reg, name="logJmag", equation="-20.0")

for reg in ["Silicon", "Oxide", "Molding"]:
    devsim.element_from_edge_model(edge_model="EField", device=device, region=reg)
    devsim.element_model(device=device, region=reg, name="Emag", equation="(EField_x^2 + EField_y^2)^(0.5)")
    devsim.element_model(device=device, region=reg, name="logEmag", equation="log(Emag + 1e-20) / log(10.0)")

tec_file = os.path.join(OUT_DIR, f"transfer_sweep_{temp_C:+.1f}C_final.tec")
devsim.write_devices(file=tec_file, type="tecplot", include_test=lambda x: x in visual_variables)
print(f"[{temp_C} C] Saved 2D field visualization to {tec_file}")

print(f"[{temp_C} C] Transfer sweep finished successfully!")