wisetop/tcad-devsim_triac
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity
Files
main
tcad-devsim_triac/resume_run.py
T

979 lines
51 KiB
Python
Raw Permalink Blame History

import pickle
import os
import sys
import glob
import gc
DEV_DIR = os.environ.get("DEV_DIR", "devices/Triac_rp")
sys.path.insert(0, os.path.abspath(DEV_DIR))
OUT_DIR = os.path.join(os.environ.get("OUT_DIR", os.path.join(DEV_DIR, "output_this_run")), "")
# Limit the thread count for parallel solvers to prevent WSL from resource starvation/disconnecting
os.environ["OMP_NUM_THREADS"] = "6"
os.environ["MKL_NUM_THREADS"] = "6"
os.environ["TBB_NUM_THREADS"] = "6"
os.environ["OPENBLAS_NUM_THREADS"] = "6"
# Disable MKL internal memory manager to prevent memory accumulation (OOM) across solve calls
os.environ["MKL_DISABLE_FAST_MM"] = "1"
# Auto-detect and load Intel MKL before devsim is imported
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
import numpy as np
is_avalanche_enabled = os.environ.get("AVALANCHE", "false").lower() == "true"
is_btbt_enabled = os.environ.get("BTBT", "false").lower() == "true"
refine_v_step = float(os.environ.get("REFINE_V_STEP", "50.0"))
is_refine_enabled = os.environ.get("REFINE", "false").lower() == "true"
if refine_v_step < 1.0:
is_refine_enabled = False
print(f"refine_v_step < 1.0 V detected ({refine_v_step} V): Dynamic refinement is completely DISABLED.")
print(f"Option: AVALANCHE={is_avalanche_enabled}, BTBT={is_btbt_enabled}, REFINE={is_refine_enabled}, REFINE_V_STEP={refine_v_step}V")
import matplotlib.pyplot as plt
import time
# Enable Intel MKL PARDISO multi-threaded sparse solver
devsim.set_parameter(name="solver_type", value="pardiso")
# ----------------- LOAD CHECKPOINT (EARLY DETECT MESH) -----------------
# 讀取指定引數
specified_checkpoint = sys.argv[1] if len(sys.argv) > 1 else None
valid_checkpoints = []
if specified_checkpoint:
# 1. 優先使用使用者指定的 Checkpoint 檔案
if os.path.exists(specified_checkpoint):
try:
with open(specified_checkpoint, "rb") as f:
data = pickle.load(f)
# 驗證資料結構完整性
if isinstance(data, dict) and ("v_current" in data or "voltage" in data) and "state" in data:
mtime = os.path.getmtime(specified_checkpoint)
v_val = data.get("v_current", data.get("voltage", 0.0))
valid_checkpoints.append((mtime, v_val, specified_checkpoint, data))
else:
print(f"Error: Specified checkpoint {specified_checkpoint} format is invalid.")
except Exception as e:
print(f"Error reading specified checkpoint {specified_checkpoint}: {e}")
else:
print(f"Error: Specified checkpoint path '{specified_checkpoint}' does not exist.")
else:
# 2. 自動搜尋邏輯（優先搜尋目前裝置資料夾內的輸出）
search_dirs = [OUT_DIR]
for d in glob.glob(os.path.join(DEV_DIR, "output_*")):
if os.path.isdir(d):
search_dirs.append(d + "/")
search_dirs.append(os.path.join(DEV_DIR, "")) # 搜尋裝置主目錄
search_dirs.append("./") # 最後搜尋根目錄
for d in search_dirs:
# 搜尋 WSL Checkpoints
for name in ["wsl_recovery_checkpoint_A.pkl", "wsl_recovery_checkpoint_B.pkl"]:
cp_path = os.path.join(d, name)
if os.path.exists(cp_path):
try:
with open(cp_path, "rb") as f:
data = pickle.load(f)
if isinstance(data, dict) and ("v_current" in data or "voltage" in data) and "state" in data:
mtime = os.path.getmtime(cp_path)
v_val = data.get("v_current", data.get("voltage", 0.0))
valid_checkpoints.append((mtime, v_val, cp_path, data))
except Exception:
pass
# 搜尋 Recon Seeds (做為備份恢復源)
for seed_path in glob.glob(os.path.join(d, "seed_*.pkl")):
if os.path.exists(seed_path):
try:
with open(seed_path, "rb") as f:
data = pickle.load(f)
if isinstance(data, dict) and ("v_current" in data or "voltage" in data) and "state" in data:
mtime = os.path.getmtime(seed_path)
v_val = data.get("v_current", data.get("voltage", 0.0))
valid_checkpoints.append((mtime, v_val, seed_path, data))
except Exception:
pass
if not valid_checkpoints:
raise RuntimeError("No valid rolling checkpoint or seed file found. Cannot resume run.")
# 我們優先選擇「模擬電壓最大」的合法檔案作為起點；若電壓相同，則選修改時間最新者
valid_checkpoints.sort(key=lambda x: (x[1], x[0]), reverse=True)
latest_mtime, latest_voltage, latest_cp, cp_data = valid_checkpoints[0]
cp_dir = os.path.dirname(latest_cp)
if cp_dir:
OUT_DIR = os.path.join(cp_dir, "")
print(f"Resuming outputs to directory: {OUT_DIR}")
os.makedirs(OUT_DIR, exist_ok=True)
import shutil
try:
mesh_orig = os.path.join(DEV_DIR, "device_2d.msh")
if os.path.exists(mesh_orig):
shutil.copy2(mesh_orig, os.path.join(OUT_DIR, "device_2d.msh"))
config_orig = os.path.join(DEV_DIR, "device_config.py")
if os.path.exists(config_orig):
shutil.copy2(config_orig, os.path.join(OUT_DIR, "device_config.py"))
if __file__:
shutil.copy2(__file__, os.path.join(OUT_DIR, os.path.basename(__file__)))
except Exception as e:
print(f"Warning: Failed to back up reproduction archive files: {e}")
print(f"Loading latest valid restore source: {latest_cp} (voltage: {latest_voltage} V, mtime: {time.ctime(latest_mtime)})")
v_current = cp_data.get("v_current", cp_data.get("voltage"))
step_size = cp_data.get("step_size", 0.1)
step_count = cp_data.get("step_count", 0)
voltage_list = cp_data.get("voltage_list", [v_current])
current_list = cp_data.get("current_list", [0.0])
next_recon_v = cp_data.get("next_recon_v", (v_current // refine_v_step + 1) * refine_v_step)
state = cp_data["state"]
# Check if checkpoint is refined by comparing Silicon node count
is_refined = False
silicon_state_len = len(cp_data["state"]["Silicon"]["Potential"])
try:
# Load coarse mesh temporarily to count Silicon nodes
devsim.create_gmsh_mesh(mesh="temp_check", file=os.path.join(DEV_DIR, "device_2d.msh"))
devsim.add_gmsh_region(mesh="temp_check", gmsh_name="Silicon", region="Silicon", material="Silicon")
devsim.add_gmsh_region(mesh="temp_check", gmsh_name="Oxide", region="Oxide", material="Oxide")
devsim.add_gmsh_region(mesh="temp_check", gmsh_name="Molding", region="Molding", material="Molding")
devsim.finalize_mesh(mesh="temp_check")
devsim.create_device(mesh="temp_check", device="temp_check")
coarse_silicon_nodes = len(devsim.get_node_model_values(device="temp_check", region="Silicon", name="x"))
devsim.delete_device(device="temp_check")
print(f"Coarse mesh Silicon nodes: {coarse_silicon_nodes}, Checkpoint Silicon nodes: {silicon_state_len}")
if silicon_state_len > coarse_silicon_nodes:
is_refined = True
except Exception as e:
print(f"Warning: Failed to load coarse mesh for node count verification: {e}")
# Fallback to old heuristic
if silicon_state_len > 70000:
is_refined = True
if is_refined:
step_size = max(step_size / 2.0, 0.1) # 載入細網格 Checkpoint 時，縮小步長至 1/2 以防首步發散
if is_refined:
# Look for refined mesh in the same directory as the checkpoint or in output_this_run/
cp_dir = os.path.dirname(latest_cp)
mesh_candidates = [
os.path.join(cp_dir, "device_2d_refined.msh") if cp_dir else None,
os.path.join(OUT_DIR, "device_2d_refined.msh")
]
mesh_file = None
for cand in mesh_candidates:
if cand and os.path.exists(cand):
mesh_file = cand
break
if not mesh_file:
import glob
cands = glob.glob("**/device_2d_refined.msh", recursive=True)
if cands:
mesh_file = cands[0]
if not mesh_file:
raise RuntimeError("Refined checkpoint requires device_2d_refined.msh, but it was not found.")
print(f"Refined checkpoint detected. Loading refined mesh: {mesh_file}")
else:
mesh_file = os.path.join(DEV_DIR, "device_2d.msh")
print(f"Standard checkpoint detected. Loading base mesh: {mesh_file}")
# sys.path already configured at top
from device_config import *
from physics.model_create import *
from physics.new_physics import *
device = "device_2d"
# 1. Load the mesh
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 for Silicon region
devsim.add_gmsh_contact(mesh=device, gmsh_name="MT1_Si", name="MT1_Si", region="Silicon", material="metal")
devsim.add_gmsh_contact(mesh=device, gmsh_name="MT2_Si", name="MT2_Si", region="Silicon", material="metal")
devsim.add_gmsh_contact(mesh=device, gmsh_name="MT1_P12_Si", name="MT1_P12_Si", region="Silicon", material="metal")
devsim.add_gmsh_contact(mesh=device, gmsh_name="MT2_P12_Si", name="MT2_P12_Si", region="Silicon", material="metal")
# Add contacts for Oxide region
devsim.add_gmsh_contact(mesh=device, gmsh_name="MT1_Ox", name="MT1_Ox", region="Oxide", material="metal")
devsim.add_gmsh_contact(mesh=device, gmsh_name="MT2_Ox", name="MT2_Ox", region="Oxide", material="metal")
# Add contacts for Molding region
devsim.add_gmsh_contact(mesh=device, gmsh_name="MT1_Mold", name="MT1_Mold", region="Molding", material="metal")
devsim.add_gmsh_contact(mesh=device, gmsh_name="MT2_Mold", name="MT2_Mold", region="Molding", material="metal")
# Add 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.add_gmsh_interface(mesh=device, gmsh_name="Si_Mold_Interface", name="Si_Mold", region0="Silicon", region1="Molding")
devsim.finalize_mesh(mesh=device)
devsim.create_device(mesh=device, device=device)
# 2. Set up doping in Silicon region
use_pcad = os.environ.get("USE_PCAD", "false").lower() == "true"
if use_pcad:
try:
from device_pcad_config import apply_pcad_doping_2d
apply_pcad_doping_2d(device, region="Silicon")
print("Applied PCAD doping configuration.")
except (ImportError, ModuleNotFoundError):
print("Warning: device_pcad_config not found. Falling back to analytical doping.")
use_pcad = False
if not use_pcad:
devsim.node_model(device=device, region="Silicon", name="nD_sub", equation=f"{N_SUB}")
def get_erfc_expr(peak, x1, x2, hdiff, vdiff):
return f"{peak} * erfc(y / {vdiff}) * 0.5 * (erf((x - ({x1})) / {hdiff}) - erf((x - ({x2})) / {hdiff}))"
p11_left_expr = get_erfc_expr(P11_PEAK, -P11_X2, -P11_X1, P_WELL_HDDIFF, P_WELL_VDDIFF)
p11_right_expr = get_erfc_expr(P11_PEAK, P11_X1, P11_X2, P_WELL_HDDIFF, P_WELL_VDDIFF)
devsim.node_model(device=device, region="Silicon", name="nA_p11_l", equation=p11_left_expr)
devsim.node_model(device=device, region="Silicon", name="nA_p11_r", equation=p11_right_expr)
p12_left_expr = get_erfc_expr(P12_PEAK, -P12_X2, -P12_X1, P_WELL_HDDIFF, P_WELL_VDDIFF)
p12_right_expr = get_erfc_expr(P12_PEAK, P12_X1, P12_X2, P_WELL_HDDIFF, P_WELL_VDDIFF)
devsim.node_model(device=device, region="Silicon", name="nA_p12_l", equation=p12_left_expr)
devsim.node_model(device=device, region="Silicon", name="nA_p12_r", equation=p12_right_expr)
p13_left_expr = get_erfc_expr(P13_PEAK, -P13_X2, -P13_X1, P_WELL_HDDIFF, P_WELL_VDDIFF)
p13_right_expr = get_erfc_expr(P13_PEAK, P13_X1, P13_X2, P_WELL_HDDIFF, P_WELL_VDDIFF)
devsim.node_model(device=device, region="Silicon", name="nA_p13_l", equation=p13_left_expr)
devsim.node_model(device=device, region="Silicon", name="nA_p13_r", equation=p13_right_expr)
nplus_left_expr = get_erfc_expr(NPLUS_PEAK, -NPLUS_X2, -NPLUS_X1, NPLUS_HDDIFF, NPLUS_VDDIFF)
nplus_right_expr = get_erfc_expr(NPLUS_PEAK, NPLUS_X1, NPLUS_X2, NPLUS_HDDIFF, NPLUS_VDDIFF)
devsim.node_model(device=device, region="Silicon", name="nD_nplus_l", equation=nplus_left_expr)
devsim.node_model(device=device, region="Silicon", name="nD_nplus_r", equation=nplus_right_expr)
mring_l_expr = get_erfc_expr(NPLUS_PEAK, -W_DEVICE, -MRING_X1, NPLUS_HDDIFF, NPLUS_VDDIFF)
mring_r_expr = get_erfc_expr(NPLUS_PEAK, MRING_X1, W_DEVICE, NPLUS_HDDIFF, NPLUS_VDDIFF)
devsim.node_model(device=device, region="Silicon", name="nD_mring_l", equation=mring_l_expr)
devsim.node_model(device=device, region="Silicon", name="nD_mring_r", equation=mring_r_expr)
devsim.node_model(device=device, region="Silicon", name="Donors",
equation="nD_sub + nD_nplus_l + nD_nplus_r + nD_mring_l + nD_mring_r")
devsim.node_model(device=device, region="Silicon", name="Acceptors",
equation="1e10 + nA_p11_l + nA_p11_r + nA_p12_l + nA_p12_r + nA_p13_l + nA_p13_r")
devsim.node_model(device=device, region="Silicon", name="NetDoping", equation="Donors - Acceptors")
devsim.node_model(device=device, region="Silicon", name="LogNetDoping", equation="asinh(NetDoping / 2.0) / log(10.0)")
# 3. Initialize electrostatic potential simulation (Poisson only)
CreateSolution(device, "Silicon", "Potential")
temp_val = os.environ.get("TEMP", "300")
devsim.set_parameter(device=device, name="T", value=temp_val)
CreateSiliconPotentialOnly(device, "Silicon")
# Oxide potential equations
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")
# Molding potential equations
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")
CreateContinuousPotentialInterface(device, "Si_Mold")
# Potential contacts setup
silicon_contacts = ["MT1_Si", "MT2_Si", "MT1_P12_Si", "MT2_P12_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")
oxide_contacts = ["MT1_Ox", "MT2_Ox"]
for c in oxide_contacts:
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")
molding_contacts = ["MT1_Mold", "MT2_Mold"]
for c in molding_contacts:
devsim.set_parameter(device=device, name=GetContactBiasName(c), value=0.0)
CreateMoldingPotentialOnlyContact(device, "Molding", c)
# Skip solving initial zero-bias Poisson since we will restore from checkpoint
print("Skipping initial Poisson solve, preparing to restore from checkpoint...")
# 4. Set up carrier solutions for Silicon Drift-Diffusion
# Compute initial guess for Electrons and Holes based on Potential
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")
# Redefine IntrinsicElectrons, IntrinsicHoles, and related models to avoid potential exponential overflow at high bias.
print("Redefining equilibrium models to prevent high-bias exponential 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="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")
# Mobility and drift diffusion equations
opts = CreateAroraMobilityLF(device, "Silicon")
# Bypassing HFMobility to prevent zero-bias convergence oscillations
CreateSiliconDriftDiffusion(device, "Silicon", **opts)
devsim.node_model(device=device, region="Silicon", name="LogElectrons", equation="log(Electrons + 1e-10) / log(10.0)")
devsim.node_model(device=device, region="Silicon", name="LogHoles", equation="log(Holes + 1e-10) / log(10.0)")
# Re-setup Silicon contacts for Drift-Diffusion
for c in silicon_contacts:
CreateSiliconDriftDiffusionContact(device, "Silicon", c, opts['Jn'], opts['Jp'])
# Skip solving initial zero-bias Drift-Diffusion
print("Skipping initial Drift-Diffusion solve, preparing to restore from checkpoint...")
# Switch continuity and potential equations for the bias sweep
print("Configuring continuity and potential equations for the bias sweep (min_error=1e5, positive update)...")
# Instantiate Avalanche and BTBT generation models
CreateAvalancheGeneration(device, "Silicon", opts['Jn'], opts['Jp'])
CreateBTBTGeneration(device, "Silicon")
av_model_n = "AvalancheGeneration" if is_avalanche_enabled else ""
av_model_p = "AvalancheGeneration_p" if is_avalanche_enabled else ""
btbt_model_n = "BTBTGeneration" if is_btbt_enabled else ""
btbt_model_p = "BTBTGeneration_p" if is_btbt_enabled else ""
if av_model_n and btbt_model_n:
main_edge_volume_model_n = f"{av_model_n} + {btbt_model_n}"
main_edge_volume_model_p = f"{av_model_p} + {btbt_model_p}"
elif av_model_n:
main_edge_volume_model_n = av_model_n
main_edge_volume_model_p = av_model_p
elif btbt_model_n:
main_edge_volume_model_n = btbt_model_n
main_edge_volume_model_p = btbt_model_p
else:
main_edge_volume_model_n = ""
main_edge_volume_model_p = ""
devsim.equation(device=device, region="Silicon", name="ElectronContinuityEquation", variable_name="Electrons",
time_node_model="NCharge", edge_model=opts['Jn'], edge_volume_model=main_edge_volume_model_n,
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=main_edge_volume_model_p,
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)
# Setup E_mag_log element models
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="E_mag_log", equation="asinh(Emag / 2.0) / log(10.0)")
# Save zero-bias tecplot and VTK
devsim.write_devices(file=f"{OUT_DIR}sweep_preview_0V.tec", type="tecplot")
# devsim.write_devices(file="sweep_preview_0V", type="vtk")
# 5. Define Sweep Parameters
v_target = 1000.0
sweep_target_env = os.environ.get("SWEEP_TARGET")
if sweep_target_env:
if sweep_target_env.endswith(('V', 'v')):
sweep_target_env = sweep_target_env[:-1]
try:
v_target = float(sweep_target_env)
print(f"Target voltage set from environment: {v_target} V")
except ValueError:
pass
if len(sys.argv) > 2:
try:
v_target = float(sys.argv[2])
print(f"Custom target voltage specified: {v_target} V")
except ValueError:
print(f"Warning: Invalid target voltage '{sys.argv[2]}', using default {v_target} V")
max_step = 50.0 # Maximum step size (V)
min_step = 1e-4 # Minimum step size (V)
compliance_current = 100e-3 # 100 mA compliance current
# Helper functions to save/restore state in case of convergence failure
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"])
# Variables were already loaded early from checkpoint
print(f"Resuming from V = {v_current:.4f} V, step = {step_size:.4f} V")
# 為了讓 log 連續，若 checkpoint 來自其他目錄，將舊日誌複製至當前的 output_this_run/ 中
cp_dir = os.path.dirname(latest_cp)
if cp_dir and cp_dir != OUT_DIR.rstrip("/"):
for log_name in ["simulation_time.log", "recon_avalanche.log", "recon_btbt.log", "recon_av_btbt.log"]:
old_log = os.path.join(cp_dir, log_name)
new_log = os.path.join(OUT_DIR, log_name)
if os.path.exists(old_log) and not os.path.exists(new_log):
try:
shutil.copy2(old_log, new_log)
print(f"Copied historical log {old_log} to {new_log} for continuity.")
except Exception as e:
print(f"Warning: Failed to copy historical log for continuity: {e}")
# Set initial contact biases to the resumed voltage v_current before restoring node values
# to avoid transient coordinate shock at the boundary during state loading
for c in ["MT1_Si", "MT1_P12_Si", "MT1_Ox", "MT1_Mold"]:
devsim.set_parameter(device=device, name=f"{c}_bias", value=v_current)
# Restore node values from checkpoint
restore_state(device, state)
# File logging setup (append mode for resume)
for log_name in ["recon_avalanche.log", "recon_btbt.log", "recon_av_btbt.log"]:
with open(f"{OUT_DIR}{log_name}", "a") as f:
pass
start_sweep_time = time.time()
print("Resuming adaptive bias sweep...")
iter_history = []
sn = 0
consecutive_fails = 0
# --- Adaptive step control thresholds ---
step_ctl_reduce = 12
step_ctl_enlarge = 8
just_refined = is_refined
if latest_cp and "seed_500V.pkl" in latest_cp:
step_size = 2.0528
just_refined = True
print(f"Resuming from seed_500V.pkl: forced initial step_size = {step_size:.4f} V, just_refined = {just_refined}")
# --- Output Milestones (Voltage & Current) ---
TEC_VOLTAGE_TARGETS = [5.0, 50.0, 250.0, 500.0]
TEC_CURRENT_TARGETS = [1e-8, 1e-7, 1e-6]
# Scan historical list to prevent duplication
saved_voltage_targets = {t for t in TEC_VOLTAGE_TARGETS if any(v >= t for v in voltage_list)}
saved_current_targets = {t for t in TEC_CURRENT_TARGETS if any(abs(i) >= t for i in current_list)}
# Pre-refinement seed save targets
seed_save_targets = [5.0, 25.0, 45.0, 95.0, 195.0, 395.0, 595.0, 795.0, 995.0, 1195.0]
saved_seeds = {t for t in seed_save_targets if v_current >= t}
# Open time log
time_log = open(f"{OUT_DIR}simulation_time.log", "a", buffering=1)
while v_current < v_target:
v_next = min(v_current + step_size, v_target)
# Apply new bias values to MT1 contacts
for c in ["MT1_Si", "MT1_P12_Si", "MT1_Ox", "MT1_Mold"]:
devsim.set_parameter(device=device, name=f"{c}_bias", value=v_next)
step_start_time = time.time()
try:
use_precondition = True
if use_precondition:
iters1 = 15
try:
# Always use log_damp preconditioning for Electron/Hole continuity equations in Stage 1
devsim.equation(device=device, region="Silicon", name="ElectronContinuityEquation", variable_name="Electrons",
time_node_model="NCharge", edge_model=opts['Jn'], edge_volume_model=main_edge_volume_model_n,
variable_update="log_damp", 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=main_edge_volume_model_p,
variable_update="log_damp", node_model="HoleGeneration", min_error=1e5)
res1 = devsim.solve(type="dc", absolute_error=1e10, relative_error=1e-3, charge_error=1e12, maximum_iterations=5, rollback=False, info=True)
iters1 = len(res1.get("iterations", []))
except devsim.error:
pass # Ignore non-convergence in pre-conditioning
finally:
# Always revert to positive variable update for Stage 2 precision Newton
devsim.equation(device=device, region="Silicon", name="ElectronContinuityEquation", variable_name="Electrons",
time_node_model="NCharge", edge_model=opts['Jn'], edge_volume_model=main_edge_volume_model_n,
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=main_edge_volume_model_p,
variable_update="positive", node_model="HoleGeneration", min_error=1e5)
import psutil
mem_stage1 = psutil.Process(os.getpid()).memory_info().rss / (1024**3)
print(f"Stage 1 (Precondition) Memory Usage: {mem_stage1:.1f} GB")
else:
iters1 = 0
mem_stage1 = 0.0
# Stage 2: Precision Newton (Strict Tolerance)
res = devsim.solve(type="dc", absolute_error=1e10, relative_error=1e-3, charge_error=1e12, maximum_iterations=10, info=True)
iters2 = len(res.get("iterations", []))
import psutil
mem_stage2 = psutil.Process(os.getpid()).memory_info().rss / (1024**3)
print(f"Stage 2 (Newton) Memory Usage: {mem_stage2:.1f} GB")
iters = f"{iters1}+{iters2}"
total_iters = iters1 + iters2
if not res.get("converged", False):
raise devsim.error("Convergence failure")
step_end_time = time.time()
time_taken = step_end_time - step_start_time
# Convergence succeeded! Compute current at MT1 terminal
# MT1 terminal current is the sum of currents on MT1_Si and MT1_P12_Si
i_n_si = devsim.get_contact_current(device=device, contact="MT1_Si", equation="ElectronContinuityEquation")
i_p_si = devsim.get_contact_current(device=device, contact="MT1_Si", equation="HoleContinuityEquation")
i_n_p12 = devsim.get_contact_current(device=device, contact="MT1_P12_Si", equation="ElectronContinuityEquation")
i_p_p12 = devsim.get_contact_current(device=device, contact="MT1_P12_Si", equation="HoleContinuityEquation")
total_curr = i_n_si + i_p_si + i_n_p12 + i_p_p12
# Update simulation status
v_current = v_next
state = save_state(device)
just_refined = False
voltage_list.append(v_current)
current_list.append(total_curr)
print(f"Step {step_count}: Converged at V = {v_current:.4f} V, I = {total_curr:.4e} A. Step size: {step_size:.4f} V. Iterations: {iters}. Time: {time_taken:.2f} s")
# Log to file
time_log.write(f"{time.strftime('%X')}\t{v_current:.4f}\t{step_size:.4f}\t{total_curr:.4e}\t{iters}\t{time_taken:.2f}\t{mem_stage1:.1f}GB\t{mem_stage2:.1f}GB\n")
# Check and save voltage milestones
for target in TEC_VOLTAGE_TARGETS:
if v_current >= target and target not in saved_voltage_targets:
filename = f"sweep_preview_{int(target)}V.tec"
print(f"Saving voltage milestone checkpoint at V = {v_current:.2f} V to {filename}...")
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="E_mag_log", equation="asinh(Emag / 2.0) / log(10.0)")
devsim.write_devices(file=f"{OUT_DIR}{filename}", type="tecplot")
saved_voltage_targets.add(target)
# Check and save current milestones
for target in TEC_CURRENT_TARGETS:
if abs(total_curr) >= target and target not in saved_current_targets:
filename = f"sweep_preview_current_{target:.1e}.tec"
print(f"Saving current milestone checkpoint at I = {total_curr:.4e} A to {filename}...")
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="E_mag_log", equation="asinh(Emag / 2.0) / log(10.0)")
devsim.write_devices(file=f"{OUT_DIR}{filename}", type="tecplot")
saved_current_targets.add(target)
# Save pre-refinement seeds when crossing target voltages
for target in seed_save_targets:
if v_current >= target and target not in saved_seeds:
seed_filename = f"{OUT_DIR}seed_{int(target)}V.pkl"
state_to_save = save_state(device)
seed_data = {"voltage": v_current, "step_size": step_size, "state": state_to_save}
with open(seed_filename, "wb") as f:
pickle.dump(seed_data, f)
print(f"Saved pre-refinement seed checkpoint at V = {v_current:.2f} V to {seed_filename}")
saved_seeds.add(target)
# Compliance check
if abs(total_curr) >= compliance_current:
print(f"Compliance current of {compliance_current:.1e} A reached at V = {v_current:.4f} V. Stopping sweep.")
time_log.write(f"Compliance current reached at V = {v_current:.4f} V.\n")
break
# --- RECONNAISSANCE PROBE & DYNAMIC REFINE ---
if is_refine_enabled and (v_current >= next_recon_v):
import dynamic_refine
try:
# 1. 執行網格自適應重劃與狀態插值
refined_device, refined_opts = dynamic_refine.refine_and_interpolate(device, v_current, is_avalanche_enabled=is_avalanche_enabled, is_btbt_enabled=is_btbt_enabled, time_log=time_log, out_dir=OUT_DIR)
device = refined_device
opts = refined_opts
just_refined = True
step_size = max(step_size / 5.0, 0.1) # 網格重建後，將步長縮小至 1/5 以防首步發散
# 2. 儲存優化後的狀態為 Seed
state = save_state(device)
seed_data = {"voltage": v_current, "step_size": step_size, "state": state}
seed_filename = f"{OUT_DIR}seed_{int(next_recon_v)}V.pkl"
with open(seed_filename, "wb") as f:
pickle.dump(seed_data, f)
print(f"\n--- RECON PROBE at {v_current:.2f} V ---")
print(f"Saved refined seed to {seed_filename}")
# Helper to run a recon probe
def run_recon_probe(probe_name, model_n, model_p, log_file):
print(f" Running Recon Probe ({probe_name})...")
# Temporarily set equations
devsim.equation(device=device, region="Silicon", name="ElectronContinuityEquation", variable_name="Electrons",
time_node_model="NCharge", edge_model=opts['Jn'], edge_volume_model=model_n,
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=model_p,
variable_update="positive", node_model="HoleGeneration", min_error=1e5)
converged = False
total_curr = 0.0
try:
# Stage 1 pre-conditioning
try:
devsim.solve(type="dc", absolute_error=1e10, relative_error=1e-3, charge_error=1e12, maximum_iterations=5, rollback=False, info=True)
except devsim.error:
pass
import psutil
mem_recon_stage1 = psutil.Process(os.getpid()).memory_info().rss / (1024**3)
print(f"Recon Stage 1 Memory Usage: {mem_recon_stage1:.1f} GB")
# Switch to positive update for Stage 2 precision Newton solve
devsim.equation(device=device, region="Silicon", name="ElectronContinuityEquation", variable_name="Electrons",
time_node_model="NCharge", edge_model=opts['Jn'], edge_volume_model=model_n,
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=model_p,
variable_update="positive", node_model="HoleGeneration", min_error=1e5)
# Stage 2 solve (Strict Tolerance)
res_av = devsim.solve(type="dc", absolute_error=1e10, relative_error=1e-3, charge_error=1e12, maximum_iterations=10, info=True)
mem_recon_stage2 = psutil.Process(os.getpid()).memory_info().rss / (1024**3)
print(f"Recon Stage 2 Memory Usage: {mem_recon_stage2:.1f} GB")
if res_av.get("converged", False):
ia_n_si = devsim.get_contact_current(device=device, contact="MT1_Si", equation="ElectronContinuityEquation")
ia_p_si = devsim.get_contact_current(device=device, contact="MT1_Si", equation="HoleContinuityEquation")
ia_n_p12 = devsim.get_contact_current(device=device, contact="MT1_P12_Si", equation="ElectronContinuityEquation")
ia_p_p12 = devsim.get_contact_current(device=device, contact="MT1_P12_Si", equation="HoleContinuityEquation")
total_curr = ia_n_si + ia_p_si + ia_n_p12 + ia_p_p12
converged = True
except devsim.error:
pass
# Restore state immediately after solve attempts
restore_state(device, state)
with open(log_file, "a") as f:
if converged:
print(f" {probe_name} Current at {v_current:.2f} V: {total_curr:.4e} A")
f.write(f"{v_current:.2f}\t{total_curr:.4e}\n")
else:
print(f" {probe_name} failed to converge.")
f.write(f"{v_current:.2f}\tFAILED\n")
# 1. Avalanche Only
run_recon_probe("Avalanche", "AvalancheGeneration", "AvalancheGeneration_p", f"{OUT_DIR}recon_avalanche.log")
# 2. BTBT Only
run_recon_probe("BTBT", "BTBTGeneration", "BTBTGeneration_p", f"{OUT_DIR}recon_btbt.log")
# 3. Avalanche + BTBT
run_recon_probe("Av+BTBT", "AvalancheGeneration + BTBTGeneration", "AvalancheGeneration_p + BTBTGeneration_p", f"{OUT_DIR}recon_av_btbt.log")
# Revert equations to main sweep configuration
devsim.equation(device=device, region="Silicon", name="ElectronContinuityEquation", variable_name="Electrons",
time_node_model="NCharge", edge_model=opts['Jn'], edge_volume_model=main_edge_volume_model_n,
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=main_edge_volume_model_p,
variable_update="positive", node_model="HoleGeneration", min_error=1e5)
else:
print("Avalanche probe skipped (AVALANCHE option is disabled).")
print("--- END RECON PROBE ---\n")
next_recon_v += refine_v_step
except devsim.error as ref_err:
print(f"\n!!! DYNAMIC REFINE FAILED at {v_current:.2f} V: {ref_err} !!!")
print("Restoring coarse device state and continuing simulation on the COARSE mesh.")
restore_state(device, state)
next_recon_v += refine_v_step
# Grow step size for next step adaptively based on Newton iterations (Rolling average)
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.6, min_step)
sn = 0
elif eff_iters < step_ctl_enlarge:
step_size = min(step_size * 1.2, max_step)
sn = 0
else:
sn += 1.25
step_count += 1
# Ping-pong rolling checkpoint
if step_count > 0 and step_count % 10 == 0:
checkpoint_file = f"{OUT_DIR}wsl_recovery_checkpoint_A.pkl" if (step_count // 10) % 2 != 0 else f"{OUT_DIR}wsl_recovery_checkpoint_B.pkl"
checkpoint_data = {
"v_current": v_current,
"step_size": step_size,
"step_count": step_count,
"voltage_list": voltage_list,
"current_list": current_list,
"next_recon_v": next_recon_v,
"state": state
}
with open(checkpoint_file, "wb") as f:
pickle.dump(checkpoint_data, f)
print(f"Rolling checkpoint saved to {checkpoint_file}")
time_log.write(f"Rolling checkpoint saved to {checkpoint_file}\n")
# Active garbage collection to prevent memory accumulation
gc.collect()
# Force Intel MKL to release thread buffers
try:
import ctypes
import glob
import sys
import os
mkl_libs = glob.glob(os.path.join(os.path.dirname(sys.executable), "../lib/libmkl_rt.so.*"))
if mkl_libs:
ctypes.CDLL(mkl_libs[0]).MKL_Free_Buffers()
# Force glibc to return freed memory to the OS (mitigate malloc fragmentation)
libc = ctypes.CDLL("libc.so.6")
libc.malloc_trim(0)
except Exception as e:
pass
except devsim.error as e:
# Convergence failure: restore last state and cut step size
step_end_time = time.time()
time_taken = step_end_time - step_start_time
consecutive_fails += 1
shrink_factor = 0.3 if consecutive_fails == 1 else 0.1
print(f"Convergence failure at V = {v_next:.4f} V. Restoring state and scaling step size by {shrink_factor} from {step_size:.4f} V (consecutive fails: {consecutive_fails}).")
time_log.write(f"{time.strftime('%X')}\t{v_next:.4f}\t{step_size:.4f}\tFAILED\t-\t{time_taken:.2f}\n")
restore_state(device, state)
step_size *= shrink_factor
iter_history = [step_ctl_reduce, step_ctl_reduce, step_ctl_reduce] # Option A virtual history to stay conservative
sn = 0 # Reset counter on failure
if step_size < min_step:
print("Step size has fallen below minimum limit. Aborting simulation.")
time_log.write(f"Aborted: step size fell below {min_step:.1e} V\n")
break
total_sweep_time = time.time() - start_sweep_time
print(f"Sweep completed in {total_sweep_time:.2f} s.")
time_log.write(f"Total Sweep Time: {total_sweep_time:.2f} s\n")
time_log.close()
# 6. Save final results and generate plots
# Save final tecplot and VTK at highest voltage
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="E_mag_log", equation="asinh(Emag / 2.0) / log(10.0)")
devsim.write_devices(file=f"{OUT_DIR}sweep_preview_final.tec", type="tecplot")
# devsim.write_devices(file="sweep_preview_final", type="vtk")
# Save I-V data to CSV
np.savetxt(f"{OUT_DIR}sweep_iv_2d.csv", np.column_stack((voltage_list, current_list)),
header="Voltage(V),Current(A/cm)", delimiter=",")
# Helper to read recon logs
def load_recon_log(filepath):
voltages = []
currents = []
if os.path.exists(filepath):
with open(filepath, "r") as f:
for line in f:
if "Voltage" in line or not line.strip():
continue
parts = line.strip().split()
if len(parts) == 2 and parts[1] != "FAILED":
try:
v = float(parts[0])
i = float(parts[1])
voltages.append(v)
currents.append(abs(i))
except ValueError:
pass
return voltages, currents
# Plot and save I-V curve
plt.figure(figsize=(8, 6))
plt.plot(voltage_list, np.abs(current_list), 'o-', color='#1f77b4', markersize=2, label='Main Sweep (Drift-Diffusion)')
v_av, i_av = load_recon_log(f"{OUT_DIR}recon_avalanche.log")
if v_av:
plt.plot(v_av, i_av, 's--', color='#ff7f0e', markersize=2, label='Avalanche Only (Recon)')
v_bt, i_bt = load_recon_log(f"{OUT_DIR}recon_btbt.log")
if v_bt:
plt.plot(v_bt, i_bt, '^--', color='#2ca02c', markersize=2, label='BTBT Only (Recon)')
v_av_bt, i_av_bt = load_recon_log(f"{OUT_DIR}recon_av_btbt.log")
if v_av_bt:
plt.plot(v_av_bt, i_av_bt, 'd--', color='#d62728', markersize=2, label='Avalanche + BTBT (Recon)')
plt.yscale('log')
plt.grid(True, which="both", ls="--")
plt.xlabel("Bias Voltage (V)")
plt.ylabel("Terminal Current (A/cm, Log Scale)")
plt.title(f"{SIM_NAME} (T={temp_val}K)")
plt.legend()
plt.tight_layout()
plt.savefig(f"{OUT_DIR}sweep_iv_2d.png", dpi=300)
plt.close()
# Generate potential & electric field plots at final converged bias
# Extract final node values
x_si = np.array(devsim.get_node_model_values(device=device, region="Silicon", name="x")) / um
y_si = np.array(devsim.get_node_model_values(device=device, region="Silicon", name="y")) / um
pot_si = np.array(devsim.get_node_model_values(device=device, region="Silicon", name="Potential"))
tri_si = np.array(devsim.get_element_node_list(device=device, region="Silicon"))
# Compute final Emag
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)")
emag_si = np.array(devsim.get_element_model_values(device=device, region="Silicon", name="Emag"))[::3]
x_ox = np.array(devsim.get_node_model_values(device=device, region="Oxide", name="x")) / um
y_ox = np.array(devsim.get_node_model_values(device=device, region="Oxide", name="y")) / um
pot_ox = np.array(devsim.get_node_model_values(device=device, region="Oxide", name="Potential"))
tri_ox = np.array(devsim.get_element_node_list(device=device, region="Oxide"))
emag_ox = np.array(devsim.get_element_model_values(device=device, region="Oxide", name="Emag"))[::3]
x_mold = np.array(devsim.get_node_model_values(device=device, region="Molding", name="x")) / um
y_mold = np.array(devsim.get_node_model_values(device=device, region="Molding", name="y")) / um
pot_mold = np.array(devsim.get_node_model_values(device=device, region="Molding", name="Potential"))
tri_mold = np.array(devsim.get_element_node_list(device=device, region="Molding"))
emag_mold = np.array(devsim.get_element_model_values(device=device, region="Molding", name="Emag"))[::3]
def draw_device_boundaries(ax):
ax.plot([-W_DEVICE/um, W_DEVICE/um], [-T_OX/um, -T_OX/um], color='black', linestyle='--', linewidth=0.8)
ax.plot([-W_DEVICE/um, W_DEVICE/um], [0, 0], color='black', linestyle='-', linewidth=0.8)
ax.plot([-W_DEVICE/um, -W_DEVICE/um], [0, H_SI/um], color='black', linestyle='-', linewidth=0.8)
ax.plot([W_DEVICE/um, W_DEVICE/um], [0, H_SI/um], color='black', linestyle='-', linewidth=0.8)
ax.plot([-W_SIM/um, W_SIM/um], [H_SI/um, H_SI/um], color='black', linestyle='-', linewidth=1.2)
fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(12, 14))
# Plot Potential
tcf1_si = ax1.tripcolor(x_si, y_si, tri_si, pot_si, cmap='RdYlBu_r', shading='gouraud')
tcf1_ox = ax1.tripcolor(x_ox, y_ox, tri_ox, pot_ox, cmap='RdYlBu_r', shading='gouraud')
tcf1_mold = ax1.tripcolor(x_mold, y_mold, tri_mold, pot_mold, cmap='RdYlBu_r', shading='gouraud')
fig.colorbar(tcf1_si, ax=ax1, label='Electrostatic Potential (V)')
draw_device_boundaries(ax1)
ax1.set_xlabel('X (μm)')
ax1.set_ylabel('Y (μm)')
ax1.set_title(f'2D Electrostatic Potential at V = {v_current:.2f} V')
ax1.set_xlim(-W_SIM / um, W_SIM / um)
ax1.set_ylim(H_SI/um + 15.0, -110.0)
# Plot EField Magnitude
tcf2_si = ax2.tripcolor(x_si, y_si, tri_si, facecolors=emag_si, cmap='inferno', shading='flat')
tcf2_ox = ax2.tripcolor(x_ox, y_ox, tri_ox, facecolors=emag_ox, cmap='inferno', shading='flat')
tcf2_mold = ax2.tripcolor(x_mold, y_mold, tri_mold, facecolors=emag_mold, cmap='inferno', shading='flat')
fig.colorbar(tcf2_si, ax=ax2, label='Electric Field Magnitude (V/cm)')
draw_device_boundaries(ax2)
ax2.set_xlabel('X (μm)')
ax2.set_ylabel('Y (μm)')
ax2.set_title(f'2D Electric Field Magnitude at V = {v_current:.2f} V')
ax2.set_xlim(-W_SIM / um, W_SIM / um)
ax2.set_ylim(H_SI/um + 15.0, -110.0)
plt.tight_layout()
plt.savefig(f"{OUT_DIR}sweep_potential_2d.png", dpi=300)
plt.close()
print(f"Sweep visualization plots saved: sweep_iv_2d.png and sweep_potential_2d.png.")
# Clean up rolling checkpoints upon successful completion
if v_current >= v_target:
print("Sweep completed successfully. Cleaning up rolling checkpoints...")
for cp in [f"{OUT_DIR}wsl_recovery_checkpoint_A.pkl", f"{OUT_DIR}wsl_recovery_checkpoint_B.pkl"]:
if os.path.exists(cp):
os.remove(cp)
else:
print(f"Sweep did not reach target voltage ({v_target} V). Keeping rolling checkpoints for recovery (current V = {v_current:.2f} V).")
import os; os.system(f'{sys.executable} plot_speed.py {OUT_DIR}simulation_time.log {OUT_DIR}simulation_speed.png')
Reference in New Issue View Git Blame Copy Permalink