pchang718
689 lines
32 KiB
Python
689 lines
32 KiB
Python
# mos_sweep.py
|
|
# 2D MOSFET Ids-Vds Curve Sweep for LDMOS
|
|
# Sweeps Vds under multiple constant Vgs values
|
|
|
|
import os
|
|
import sys
|
|
import glob
|
|
import gc
|
|
import time
|
|
import pickle
|
|
import numpy as np
|
|
import matplotlib.pyplot as plt
|
|
|
|
# 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")
|
|
|
|
# Default device directory to LDMOS
|
|
DEV_DIR = os.environ.get("DEV_DIR", "devices/LDMOS")
|
|
sys.path.insert(0, os.path.abspath(DEV_DIR))
|
|
|
|
# Respect OUT_DIR env variable, defaulting to devices/LDMOS/output_yymmdd_03 based on current date
|
|
default_out_dirname = f"output_{time.strftime('%y%m%d')}_03"
|
|
OUT_DIR = os.path.join(os.environ.get("OUT_DIR", os.path.join(DEV_DIR, default_out_dirname)), "")
|
|
os.makedirs(OUT_DIR, exist_ok=True)
|
|
|
|
# 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"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("Applied PCAD doping successfully.")
|
|
except ImportError:
|
|
print("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
|
|
temp_val = os.environ.get("TEMP", "300")
|
|
devsim.set_parameter(device=device, name="T", value=temp_val)
|
|
devsim.set_parameter(device=device, region="Silicon", name="T", value=temp_val)
|
|
|
|
|
|
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")
|
|
|
|
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")
|
|
|
|
# Contact Potential helper functions
|
|
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")
|
|
|
|
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")
|
|
|
|
# Potential boundary conditions on contacts
|
|
silicon_contacts = ["contL_Si", "contR_Si"]
|
|
oxide_contacts = ["gate_Ox"]
|
|
molding_contacts = ["gate_Mold", "contL_Mold", "contR_Mold"]
|
|
|
|
for c in silicon_contacts:
|
|
devsim.set_parameter(device=device, name=GetContactBiasName(c), value=0.0)
|
|
CreateSiliconPotentialOnlyContact(device, "Silicon", c)
|
|
|
|
for c in oxide_contacts:
|
|
devsim.set_parameter(device=device, name=GetContactBiasName(c), value=0.0)
|
|
CreateOxidePotentialOnlyContact(device, "Oxide", c)
|
|
|
|
for c in molding_contacts:
|
|
devsim.set_parameter(device=device, name=GetContactBiasName(c), value=0.0)
|
|
CreateMoldingPotentialOnlyContact(device, "Molding", c)
|
|
|
|
print("Solving zero-bias Poisson...")
|
|
devsim.solve(type="dc", absolute_error=1.0, relative_error=1e-10, maximum_iterations=50)
|
|
print("Poisson converged.")
|
|
|
|
# --- 4. Setup Drift-Diffusion (DD) ---
|
|
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")
|
|
|
|
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="IntrinsicCharge", equation="Holes - Electrons + NetDoping")
|
|
devsim.node_model(device=device, region="Silicon", name="IntrinsicCharge:Potential", equation="0")
|
|
devsim.node_model(device=device, region="Silicon", name="IntrinsicCharge:Electrons", equation="-1")
|
|
devsim.node_model(device=device, region="Silicon", name="IntrinsicCharge: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'])
|
|
|
|
print("Solving zero-bias Drift-Diffusion...")
|
|
devsim.solve(type="dc", absolute_error=1e10, relative_error=1e30, charge_error=1e12, maximum_iterations=50)
|
|
print("Initial Drift-Diffusion converged.")
|
|
|
|
# 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="default", min_error=1e-3)
|
|
|
|
|
|
# --- 5. State Saving/Restoration and Sweeping Logic ---
|
|
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"])
|
|
|
|
# Keep the original zero-bias state as reference
|
|
zero_bias_state = save_state(device)
|
|
|
|
|
|
def ramp_gate_bias(device, target_vgs):
|
|
"""
|
|
Ramps gate voltage from 0.0V to target_vgs under Vds = 0.0V.
|
|
"""
|
|
# Sync bias parameters with zero_bias_state before ramping
|
|
for c in ["gate_Ox", "gate_Mold", "contR_Si", "contR_Mold", "contL_Si", "contL_Mold"]:
|
|
devsim.set_parameter(device=device, name=f"{c}_bias", value=0.0)
|
|
|
|
current_vgs = 0.0
|
|
if abs(current_vgs - target_vgs) < 1e-5:
|
|
return True
|
|
|
|
print(f"Ramping gate bias from {current_vgs:.2f}V to {target_vgs:.2f}V...")
|
|
|
|
g_step = 0.5 if target_vgs > current_vgs else -0.5
|
|
v_curr = current_vgs
|
|
|
|
while abs(v_curr - target_vgs) > 1e-5:
|
|
v_next = v_curr + g_step
|
|
if g_step > 0 and v_next > target_vgs:
|
|
v_next = target_vgs
|
|
elif g_step < 0 and v_next < target_vgs:
|
|
v_next = target_vgs
|
|
|
|
for c in ["gate_Ox", "gate_Mold"]:
|
|
devsim.set_parameter(device=device, name=f"{c}_bias", value=v_next)
|
|
|
|
try:
|
|
devsim.solve(type="dc", absolute_error=1e10, relative_error=3e-3, charge_error=1e12, maximum_iterations=20)
|
|
v_curr = v_next
|
|
except devsim.error:
|
|
# Revert to last converged bias and reduce step size
|
|
for c in ["gate_Ox", "gate_Mold"]:
|
|
devsim.set_parameter(device=device, name=f"{c}_bias", value=v_curr)
|
|
g_step *= 0.5
|
|
if abs(g_step) < 1e-3:
|
|
print("Error: Gate bias ramping failed to converge.")
|
|
return False
|
|
print(f"Gate bias successfully ramped to {target_vgs:.2f}V.")
|
|
return True
|
|
|
|
|
|
# Whitelist of variables to output for visualization.
|
|
# Unwanted variables can be commented out below to reduce file size without losing the ability to enable them later.
|
|
visual_variables = {
|
|
"x", "y",
|
|
"Potential", "Electrons", "Holes",
|
|
"NetDoping", "Donors", "Acceptors", "LogNetDoping",
|
|
"Emag", "logEmag",
|
|
"Jmag", "logJmag",
|
|
"AvalancheGeneration_mag", "logAvalancheGen",
|
|
|
|
# --- Vector components (commented out by default, uncomment if needed) ---
|
|
# "EField_x", "EField_y",
|
|
# "Jn_x", "Jn_y",
|
|
# "Jp_x", "Jp_y",
|
|
|
|
# --- Geometry/Volume metrics (commented out by default, uncomment if needed) ---
|
|
# "NodeVolume", "ElementNodeVolume", "ElementEdgeCouple",
|
|
|
|
# --- Derivative/Jacobian variables (commented out by default to save space/speed up loading) ---
|
|
# "Potential:Potential@n0", "Potential:Potential@n1",
|
|
# "Electrons:Electrons@n0", "Electrons:Electrons@n1",
|
|
# "Holes:Holes@n0", "Holes:Holes@n1",
|
|
# "AvalancheGeneration:Electrons@n0", "AvalancheGeneration:Electrons@n1",
|
|
# "AvalancheGeneration:Holes@n0", "AvalancheGeneration:Holes@n1",
|
|
# "AvalancheGeneration:Potential@n0", "AvalancheGeneration:Potential@n1",
|
|
}
|
|
|
|
|
|
def save_visualization(device, vgs, vds, prefix):
|
|
# Calculate current density components on elements for visualization
|
|
for reg in ["Silicon"]:
|
|
jn_name = opts['Jn']
|
|
jp_name = opts['Jp']
|
|
devsim.element_from_edge_model(edge_model=jn_name, device=device, region=reg)
|
|
devsim.element_from_edge_model(edge_model=jp_name, device=device, region=reg)
|
|
devsim.element_model(device=device, region=reg, name="Jn_x", equation=f"{jn_name}_x")
|
|
devsim.element_model(device=device, region=reg, name="Jn_y", equation=f"{jn_name}_y")
|
|
devsim.element_model(device=device, region=reg, name="Jp_x", equation=f"{jp_name}_x")
|
|
devsim.element_model(device=device, region=reg, name="Jp_y", equation=f"{jp_name}_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)")
|
|
|
|
devsim.element_from_edge_model(edge_model="AvalancheGeneration", device=device, region=reg)
|
|
devsim.element_model(device=device, region=reg, name="AvalancheGeneration_raw",
|
|
equation="(AvalancheGeneration_x^2 + AvalancheGeneration_y^2)^(0.5)")
|
|
devsim.element_model(device=device, region=reg, name="AvalancheGeneration_mag",
|
|
equation="(AvalancheGeneration_raw > 1e-20) * AvalancheGeneration_raw")
|
|
devsim.element_model(device=device, region=reg, name="logAvalancheGen",
|
|
equation="log(AvalancheGeneration_mag + 1e-40) / log(10.0)")
|
|
|
|
# Define placeholder variables in non-conducting regions to prevent Tecplot column mismatch in ParaView
|
|
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")
|
|
devsim.element_model(device=device, region=reg, name="AvalancheGeneration_raw", equation="0.0")
|
|
devsim.element_model(device=device, region=reg, name="AvalancheGeneration_mag", equation="0.0")
|
|
devsim.element_model(device=device, region=reg, name="logAvalancheGen", equation="-40.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)")
|
|
|
|
# Save tecplot file
|
|
tec_file = os.path.join(OUT_DIR, f"mos_sweep_vgs_{vgs:.1f}V_vds_{prefix}.tec")
|
|
devsim.write_devices(file=tec_file, type="tecplot", include_test=lambda x: x in visual_variables)
|
|
print(f"Saved 2D field visualization to {tec_file}")
|
|
|
|
# Save VTK XML files (.vtm/.vtu) (commented out to keep output folder clean)
|
|
# orig_cwd = os.getcwd()
|
|
# try:
|
|
# os.chdir(OUT_DIR)
|
|
# vtk_file = f"mos_sweep_vgs_{vgs:.1f}V_vds_{prefix}"
|
|
# devsim.write_devices(file=vtk_file, type="vtk", include_test=lambda x: x in visual_variables)
|
|
# print(f"Saved 2D VTK field visualization to {os.path.join(OUT_DIR, vtk_file)}.vtm")
|
|
# finally:
|
|
# os.chdir(orig_cwd)
|
|
|
|
|
|
def sweep_vds(device, vgs, vds_max):
|
|
"""
|
|
Sweeps Vds from 0.0V to vds_max while keeping Vgs constant.
|
|
"""
|
|
print(f"Starting Vds sweep up to Vds_max = {vds_max:.2f}V for Vgs = {vgs:.2f}V...")
|
|
|
|
log_path = os.path.join(OUT_DIR, "simulation_time.log")
|
|
file_exists = os.path.exists(log_path)
|
|
time_log = open(log_path, "a", buffering=1)
|
|
if not file_exists:
|
|
time_log.write("Time\tVoltage(V)\tStep(V)\tCurrent(A)\tIav_eval(A)\tIterations\tTimeTaken(s)\tMemStage1\tMemStage2\n")
|
|
|
|
# Check if there is an existing seed checkpoint to resume from (checked in descending Vds order)
|
|
resume_data = None
|
|
for checkpoint in [60.0, 30.0, 10.0]:
|
|
seed_path = os.path.join(OUT_DIR, f"seed_vgs_{vgs:.1f}V_vds_{int(checkpoint)}V.pkl")
|
|
if os.path.exists(seed_path):
|
|
try:
|
|
with open(seed_path, "rb") as f:
|
|
resume_data = pickle.load(f)
|
|
print(f"Found seed checkpoint file: {seed_path}")
|
|
print(f"Successfully loaded state and resuming sweep from Vds = {checkpoint:.1f}V...")
|
|
break
|
|
except Exception as e:
|
|
print(f"Warning: Failed to load seed file {seed_path}: {e}")
|
|
|
|
if resume_data:
|
|
v_curr = resume_data["voltage"]
|
|
step_size = resume_data["step_size"]
|
|
vds_list = list(resume_data["vds_list"])
|
|
ids_list = list(resume_data["ids_list"])
|
|
iav_list = list(resume_data["iav_list"])
|
|
# Restore node values
|
|
restore_state(device, resume_data["state"])
|
|
# Restore contact bias parameters in DEVSIM
|
|
for c in ["contR_Si", "contR_Mold"]:
|
|
devsim.set_parameter(device=device, name=f"{c}_bias", value=v_curr)
|
|
for c in ["gate_Ox", "gate_Mold"]:
|
|
devsim.set_parameter(device=device, name=f"{c}_bias", value=vgs)
|
|
step_state = save_state(device)
|
|
else:
|
|
vds_list = [0.0]
|
|
# Compute initial current at Vds = 0
|
|
i_n_si = devsim.get_contact_current(device=device, contact="contR_Si", equation="ElectronContinuityEquation")
|
|
i_p_si = devsim.get_contact_current(device=device, contact="contR_Si", equation="HoleContinuityEquation")
|
|
ids_list = [i_n_si + i_p_si]
|
|
|
|
g_av_init = devsim.get_edge_model_values(device=device, region="Silicon", name="AvalancheGeneration")
|
|
edge_vol_init = devsim.get_edge_model_values(device=device, region="Silicon", name="EdgeNodeVolume")
|
|
init_av = 2.0 * sum(g * v for g, v in zip(g_av_init, edge_vol_init))
|
|
iav_list = [init_av]
|
|
|
|
v_curr = 0.0
|
|
step_size = 0.05 # Start with small step size to capture the linear region smoothly
|
|
step_state = save_state(device)
|
|
|
|
max_step = max(1.0, min(50.0, vds_max / 10.0)) # Dynamic limit for high voltage sweeps
|
|
min_step = 1e-4
|
|
|
|
consecutive_fails = 0
|
|
step_ctl_reduce = 12
|
|
step_ctl_enlarge = 8
|
|
iter_history = []
|
|
sn = 0
|
|
|
|
while v_curr < vds_max:
|
|
v_next = v_curr + step_size
|
|
# Force exact hits on checkpoint targets (e.g. 10V, 30V, 40V, 60V)
|
|
for target in [10.0, 30.0, 40.0, 60.0]:
|
|
if v_curr < target <= v_next:
|
|
v_next = target
|
|
break
|
|
v_next = min(v_next, vds_max)
|
|
|
|
# Apply new Vds to drain contacts
|
|
for c in ["contR_Si", "contR_Mold"]:
|
|
devsim.set_parameter(device=device, name=f"{c}_bias", value=v_next)
|
|
|
|
t_start = time.time()
|
|
try:
|
|
# Stage 1: Pre-conditioning
|
|
iters1 = 10
|
|
try:
|
|
res1 = devsim.solve(type="dc", absolute_error=1e10, relative_error=1e-1, charge_error=1e12, maximum_iterations=10, info=True)
|
|
iters1 = len(res1.get("iterations", []))
|
|
except devsim.error:
|
|
pass
|
|
|
|
import psutil
|
|
mem_stage1 = psutil.Process(os.getpid()).memory_info().rss / (1024**3)
|
|
# Stage 2: Precision Newton
|
|
res = devsim.solve(type="dc", absolute_error=1e10, relative_error=3e-3, charge_error=1e12, maximum_iterations=15, info=True)
|
|
iters2 = len(res.get("iterations", []))
|
|
total_iters = iters2
|
|
|
|
mem_stage2 = psutil.Process(os.getpid()).memory_info().rss / (1024**3)
|
|
|
|
if not res.get("converged", False):
|
|
raise devsim.error("Convergence failure")
|
|
|
|
# Convergence succeeded! Extract currents
|
|
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 ionization current (avalanche onset evaluation)
|
|
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))
|
|
|
|
v_curr = v_next
|
|
step_state = save_state(device)
|
|
|
|
vds_list.append(v_curr)
|
|
ids_list.append(total_curr)
|
|
iav_list.append(total_av_curr)
|
|
|
|
t_taken = time.time() - t_start
|
|
print(f" Vds = {v_curr:.3f}V, Ids = {total_curr:.4e}A, I_av_eval = {total_av_curr:.4e}A (step: {step_size:.3f}V, iters: {iters1}+{iters2})")
|
|
time_log.write(f"{time.strftime('%X')}\t{v_curr:.4f}\t{step_size:.4f}\t{total_curr:.4e}\t{total_av_curr:.4e}\t{iters1}+{iters2}\t{t_taken:.2f}\t{mem_stage1:.1f}GB\t{mem_stage2:.1f}GB\n")
|
|
|
|
# Save visual checkpoints at exactly 40.0V and 60.0V
|
|
if abs(v_curr - 40.0) < 1e-5:
|
|
save_visualization(device, vgs, v_curr, "40.0V")
|
|
elif abs(v_curr - 60.0) < 1e-5:
|
|
save_visualization(device, vgs, v_curr, "60.0V")
|
|
|
|
# Save state seed files (.pkl) at exactly 10.0V, 30.0V, and 60.0V
|
|
for cp in [10.0, 30.0, 60.0]:
|
|
if abs(v_curr - cp) < 1e-5:
|
|
seed_path = os.path.join(OUT_DIR, f"seed_vgs_{vgs:.1f}V_vds_{int(cp)}V.pkl")
|
|
seed_data = {
|
|
"voltage": v_curr,
|
|
"step_size": step_size,
|
|
"vds_list": vds_list,
|
|
"ids_list": ids_list,
|
|
"iav_list": iav_list,
|
|
"state": save_state(device)
|
|
}
|
|
with open(seed_path, "wb") as f:
|
|
pickle.dump(seed_data, f)
|
|
print(f"Saved state seed checkpoint at Vds = {v_curr:.1f}V to {seed_path}")
|
|
|
|
# Adaptive step control
|
|
consecutive_fails = 0
|
|
iter_history.append(total_iters)
|
|
if len(iter_history) > 3:
|
|
iter_history.pop(0)
|
|
avg_iters = sum(iter_history) / len(iter_history)
|
|
|
|
eff_iters = max(0.0, avg_iters - sn)
|
|
if avg_iters > step_ctl_reduce:
|
|
step_size = max(step_size * 0.5, min_step)
|
|
sn = 0
|
|
elif eff_iters < step_ctl_enlarge:
|
|
step_size = min(step_size * 1.5, max_step)
|
|
sn = 0
|
|
else:
|
|
sn += 1.25
|
|
|
|
gc.collect()
|
|
|
|
except devsim.error:
|
|
consecutive_fails += 1
|
|
attempted_step = step_size
|
|
shrink_factor = 0.3 if consecutive_fails == 1 else 0.1
|
|
step_size *= shrink_factor
|
|
t_taken = time.time() - t_start
|
|
print(f" Convergence failure at Vds = {v_next:.3f}V. Reverting state and cutting step size to {step_size:.4f}V")
|
|
time_log.write(f"{time.strftime('%X')}\t{v_next:.4f}\t{attempted_step:.4f}\tFAILED\t-\t{t_taken:.2f}\n")
|
|
|
|
restore_state(device, step_state)
|
|
iter_history = [step_ctl_reduce, step_ctl_reduce, step_ctl_reduce]
|
|
sn = 0
|
|
|
|
if step_size < min_step:
|
|
print(" Step size fell below minimum limit. Aborting sweep for this Vgs.")
|
|
break
|
|
|
|
time_log.close()
|
|
return np.array(vds_list), np.array(ids_list), np.array(iav_list)
|
|
|
|
|
|
# --- 6. Run Sweeps ---
|
|
# Configurable Vgs list and max Vds
|
|
#vgs_list = [0.0, 1.0, 2.0, 5.0, 10.0, 15.0]
|
|
vgs_list = [0.0, 2.0, 5.0, 10.0]
|
|
vds_max = 60.0
|
|
|
|
print(f"MOS Sweep Configuration:")
|
|
print(f" Vgs List: {vgs_list}")
|
|
print(f" Vds Sweep Range: 0.0V to {vds_max}V")
|
|
print(f" Output Directory: {OUT_DIR}")
|
|
|
|
sweep_results_active = {}
|
|
|
|
for vgs in vgs_list:
|
|
print(f"\n========================================")
|
|
print(f"Processing Vgs = {vgs:.2f}V")
|
|
print(f"========================================")
|
|
|
|
# Restore to clean zero-bias state
|
|
restore_state(device, zero_bias_state)
|
|
|
|
# Ramp gate to target Vgs
|
|
success = ramp_gate_bias(device, vgs)
|
|
if not success:
|
|
print(f"Failed to ramp gate bias to Vgs = {vgs}V. Skipping this sweep.")
|
|
continue
|
|
|
|
# Sweep Vds
|
|
vds_vals, ids_vals, iav_vals = sweep_vds(device, vgs, vds_max)
|
|
sweep_results_active[vgs] = (vds_vals, ids_vals, iav_vals)
|
|
|
|
# Save CSV for this Vgs
|
|
csv_file = os.path.join(OUT_DIR, f"mos_sweep_vgs_{vgs:.1f}V.csv")
|
|
np.savetxt(csv_file, np.column_stack((vds_vals, ids_vals, iav_vals)),
|
|
header="Vds(V),Ids(A),Iav_eval(A)", delimiter=",", comments="")
|
|
print(f"Saved Vds-Ids-Iav data for Vgs = {vgs:.1f}V to {csv_file}")
|
|
|
|
# Calculate current density components on elements for visualization
|
|
# Save final visual outputs using our global function
|
|
save_visualization(device, vgs, vds_vals[-1], "final_vds")
|
|
|
|
|
|
# --- 7. Plotting and Global Data Save ---
|
|
# Populate final sweep_results by combining active run and disk files
|
|
all_vgs_targets = sorted(list(set(vgs_list + [0.0, 1.0, 2.0, 5.0, 10.0, 15.0])))
|
|
sweep_results = {}
|
|
for v in all_vgs_targets:
|
|
if v in sweep_results_active:
|
|
sweep_results[v] = sweep_results_active[v]
|
|
else:
|
|
csv_file = os.path.join(OUT_DIR, f"mos_sweep_vgs_{v:.1f}V.csv")
|
|
if os.path.exists(csv_file):
|
|
try:
|
|
data = np.loadtxt(csv_file, delimiter=",", skiprows=1)
|
|
if data.ndim == 2:
|
|
if data.shape[1] == 3:
|
|
sweep_results[v] = (data[:, 0], data[:, 1], data[:, 2])
|
|
elif data.shape[1] == 2:
|
|
sweep_results[v] = (data[:, 0], data[:, 1], np.zeros(len(data[:, 0])))
|
|
except Exception as e:
|
|
print(f"Warning: Could not load {csv_file}: {e}")
|
|
|
|
if not sweep_results:
|
|
print("Warning: No sweep results found to plot or combine.")
|
|
else:
|
|
# Combine all curves to a single combined CSV file
|
|
combined_csv_file = os.path.join(OUT_DIR, "mos_sweep_all_curves.csv")
|
|
with open(combined_csv_file, "w") as f_combined:
|
|
# Write header
|
|
headers = []
|
|
for vgs in sorted(sweep_results.keys()):
|
|
headers.append(f"Vds_Vgs_{vgs:.1f}V(V)")
|
|
headers.append(f"Ids_Vgs_{vgs:.1f}V(A)")
|
|
headers.append(f"Iav_eval_Vgs_{vgs:.1f}V(A)")
|
|
f_combined.write(",".join(headers) + "\n")
|
|
|
|
# Find max length to align columns
|
|
max_len = max(len(v[0]) for v in sweep_results.values())
|
|
for idx in range(max_len):
|
|
row_cells = []
|
|
for vgs in sorted(sweep_results.keys()):
|
|
vds_vals, ids_vals, iav_vals = sweep_results[vgs]
|
|
if idx < len(vds_vals):
|
|
row_cells.append(f"{vds_vals[idx]:.6e}")
|
|
row_cells.append(f"{ids_vals[idx]:.6e}")
|
|
row_cells.append(f"{iav_vals[idx]:.6e}")
|
|
else:
|
|
row_cells.append("")
|
|
row_cells.append("")
|
|
row_cells.append("")
|
|
f_combined.write(",".join(row_cells) + "\n")
|
|
|
|
print(f"\nSaved combined curve data to {combined_csv_file}")
|
|
|
|
# Plot 1: Standard Linear Scale Plot (Ids in mA)
|
|
plt.figure(figsize=(10, 7))
|
|
for vgs in sorted(sweep_results.keys()):
|
|
vds_vals, ids_vals, _ = sweep_results[vgs]
|
|
plt.plot(vds_vals, ids_vals * 1e3, 'o-', markersize=3, label=f"Vgs = {vgs:.1f} V")
|
|
plt.grid(True, which="both", linestyle="--", alpha=0.5)
|
|
plt.xlabel("Drain-Source Voltage Vds (V)", fontsize=12)
|
|
plt.ylabel("Drain-Source Current Ids (mA)", fontsize=12)
|
|
plt.title("LDMOS 2D Ids vs Vds Output Characteristics", fontsize=14, fontweight="bold")
|
|
plt.legend(loc="upper left", fontsize=10)
|
|
plt.tight_layout()
|
|
plot_path = os.path.join(OUT_DIR, "mos_sweep_iv_curves.png")
|
|
plt.savefig(plot_path, dpi=300)
|
|
plt.close()
|
|
print(f"Saved linear visualization plot to {plot_path}")
|
|
|
|
# Plot 2: Logarithmic Scale Plot (Ids and Iav_eval in Amperes)
|
|
plt.figure(figsize=(10, 7))
|
|
for vgs in sorted(sweep_results.keys()):
|
|
vds_vals, ids_vals, iav_vals = sweep_results[vgs]
|
|
line, = plt.plot(vds_vals, np.abs(ids_vals), '-', linewidth=1.5, label=f"Vgs = {vgs:.1f} V")
|
|
plt.plot(vds_vals, np.abs(iav_vals), '--', color=line.get_color(), linewidth=1.2)
|
|
|
|
# Add dummy entries for line style legend
|
|
plt.plot([], [], 'k-', label='Total Ids (Solid)')
|
|
plt.plot([], [], 'k--', label='Iav_eval (Dashed)')
|
|
|
|
plt.yscale('log')
|
|
plt.ylim(1e-14, 1.0) # Limit range to see leakage and active region clearly
|
|
plt.grid(True, which="both", linestyle="--", alpha=0.5)
|
|
plt.xlabel("Drain-Source Voltage Vds (V)", fontsize=12)
|
|
plt.ylabel("Current Magnitude (A)", fontsize=12)
|
|
plt.title("LDMOS 2D Ids & Integrated Ionization Current (Log Scale)", fontsize=14, fontweight="bold")
|
|
plt.legend(loc="lower right", fontsize=10)
|
|
plt.tight_layout()
|
|
plot_path_log = os.path.join(OUT_DIR, "mos_sweep_iv_curves_log.png")
|
|
plt.savefig(plot_path_log, dpi=300)
|
|
plt.close()
|
|
print(f"Saved log visualization plot to {plot_path_log}")
|
|
|
|
# Run speed plotting script
|
|
log_path = os.path.join(OUT_DIR, "simulation_time.log")
|
|
speed_png = os.path.join(OUT_DIR, "simulation_speed.png")
|
|
if os.path.exists(log_path):
|
|
print(f"Plotting simulation speed to {speed_png}...")
|
|
import os; os.system(f'{sys.executable} plot_speed.py {log_path} {speed_png}')
|
|
|
|
print("\nMOS Sweep finished successfully!")
|