import devsim import os import shutil import pickle import numpy as np import sys sys.path.append("/home/pchan/devsim2026") from device_config import * from physics.model_create import * from physics.new_physics import * import generate_mesh_2d OUT_DIR = "output_this_run/" def setup_physics_for_device(device, is_avalanche_enabled=False, is_btbt_enabled=False): 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") # 1. Add GMSH region, contacts, and interfaces try: 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) except devsim.error as e: if "must not be finalized" in str(e) or "already exists" in str(e): pass else: raise e try: devsim.create_device(mesh=device, device=device) except devsim.error as e: if "already exists" in str(e): pass else: raise e # 2. Set up doping in Silicon region 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. Create solutions CreateSolution(device, "Silicon", "Potential") CreateSolution(device, "Silicon", "Electrons") CreateSolution(device, "Silicon", "Holes") sim_temp = os.environ.get("TEMP", "300") devsim.set_parameter(device=device, name="T", value=sim_temp) 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) 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) 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) # Redefine IntrinsicElectrons, IntrinsicHoles, and related models to avoid potential 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") # 全耦合 Drift-Diffusion 系統（Arora mobility + 連續方程式） opts = CreateAroraMobilityLF(device, "Silicon") CreateSiliconDriftDiffusion(device, "Silicon", **opts) # Drift diffusion contacts for c in silicon_contacts: CreateSiliconDriftDiffusionContact(device, "Silicon", c, opts['Jn'], opts['Jp']) # Avalanche generation model (enabled/disabled by refine_and_interpolate caller) if is_avalanche_enabled: CreateAvalancheGeneration(device, "Silicon", opts['Jn'], opts['Jp']) # BTBT generation model if is_btbt_enabled: 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: from physics.model_create import CreateEdgeModel, CreateEdgeModelDerivatives CreateEdgeModel(device, "Silicon", "CombinedGeneration", "AvalancheGeneration + BTBTGeneration") CreateEdgeModel(device, "Silicon", "CombinedGeneration_p", "AvalancheGeneration_p + BTBTGeneration_p") for i in ("Potential", "Electrons", "Holes"): CreateEdgeModelDerivatives(device, "Silicon", "CombinedGeneration", "AvalancheGeneration + BTBTGeneration", i) CreateEdgeModelDerivatives(device, "Silicon", "CombinedGeneration_p", "AvalancheGeneration_p + BTBTGeneration_p", i) gen_model_n = "CombinedGeneration" gen_model_p = "CombinedGeneration_p" elif av_model_n: gen_model_n = av_model_n gen_model_p = av_model_p elif btbt_model_n: gen_model_n = btbt_model_n gen_model_p = btbt_model_p else: gen_model_n = "" gen_model_p = "" # 預設以 positive (full Newton) 方式註冊連續方程式 # refine_and_interpolate 在插值後會臨時切換至 log_damp 做 Stage 1 預處理 devsim.equation(device=device, region="Silicon", name="ElectronContinuityEquation", variable_name="Electrons", time_node_model="NCharge", edge_model=opts['Jn'], edge_volume_model=gen_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=gen_model_p, variable_update="positive", node_model="HoleGeneration", min_error=1e5) 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)") 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)") return opts def solve_decoupled_carriers(device_name, opts): # We temporarily delete PotentialEquation to solve carrier equations individually for reg in ["Silicon", "Oxide", "Molding"]: try: devsim.delete_equation(device=device_name, region=reg, name="PotentialEquation") except Exception: pass for c in ["MT1_Si", "MT2_Si", "MT1_P12_Si", "MT2_P12_Si", "MT1_Ox", "MT2_Ox", "MT1_Mold", "MT2_Mold"]: try: devsim.delete_contact_equation(device=device_name, contact=c, name="PotentialEquation") except Exception: pass for i in ["Si_Ox", "Ox_Mold", "Si_Mold"]: try: devsim.delete_interface_equation(device=device_name, interface=i, name="PotentialEquation") except Exception: pass # Delete HoleContinuityEquation before solving ElectronContinuityEquation try: devsim.delete_equation(device=device_name, region="Silicon", name="HoleContinuityEquation") except Exception: pass for c in ["MT1_Si", "MT2_Si", "MT1_P12_Si", "MT2_P12_Si"]: try: devsim.delete_contact_equation(device=device_name, contact=c, name="HoleContinuityEquation") except Exception: pass # Solve ElectronContinuityEquation alone (linear in Electrons) print(" [Gummel] Solving Electron Continuity Equation alone...") devsim.equation(device=device_name, region="Silicon", name="ElectronContinuityEquation", variable_name="Electrons", time_node_model="NCharge", edge_model=opts['Jn'], edge_volume_model="", variable_update="default", node_model="ElectronGeneration", min_error=1e5) for c in ["MT1_Si", "MT2_Si", "MT1_P12_Si", "MT2_P12_Si"]: contact_electrons_name = f"{c}nodeelectrons" devsim.contact_equation(device=device_name, contact=c, name="ElectronContinuityEquation", node_model=contact_electrons_name, edge_current_model=opts['Jn']) try: devsim.solve(type="dc", absolute_error=1e10, relative_error=1e-5, charge_error=1e12, maximum_iterations=15, info=True) except Exception as e: print(f" [Gummel] Electron solve warning: {e}") devsim.delete_equation(device=device_name, region="Silicon", name="ElectronContinuityEquation") for c in ["MT1_Si", "MT2_Si", "MT1_P12_Si", "MT2_P12_Si"]: devsim.delete_contact_equation(device=device_name, contact=c, name="ElectronContinuityEquation") # Solve HoleContinuityEquation alone (linear in Holes) print(" [Gummel] Solving Hole Continuity Equation alone...") devsim.equation(device=device_name, region="Silicon", name="HoleContinuityEquation", variable_name="Holes", time_node_model="PCharge", edge_model=opts['Jp'], edge_volume_model="", variable_update="default", node_model="HoleGeneration", min_error=1e5) for c in ["MT1_Si", "MT2_Si", "MT1_P12_Si", "MT2_P12_Si"]: contact_holes_name = f"{c}nodeholes" devsim.contact_equation(device=device_name, contact=c, name="HoleContinuityEquation", node_model=contact_holes_name, edge_current_model=opts['Jp']) try: devsim.solve(type="dc", absolute_error=1e10, relative_error=1e-5, charge_error=1e12, maximum_iterations=15, info=True) except Exception as e: print(f" [Gummel] Hole solve warning: {e}") devsim.delete_equation(device=device_name, region="Silicon", name="HoleContinuityEquation") for c in ["MT1_Si", "MT2_Si", "MT1_P12_Si", "MT2_P12_Si"]: devsim.delete_contact_equation(device=device_name, contact=c, name="HoleContinuityEquation") def enforce_contact_boundary_conditions(device_name): import math print("[Boundary] Enforcing exact contact boundary conditions to prevent numerical shock...") # Silicon Contacts silicon_contacts = ["MT1_Si", "MT2_Si", "MT1_P12_Si", "MT2_P12_Si"] try: net_doping = np.array(devsim.get_node_model_values(device=device_name, region="Silicon", name="NetDoping")) vt = np.array(devsim.get_node_model_values(device=device_name, region="Silicon", name="V_t")) nie = np.array(devsim.get_node_model_values(device=device_name, region="Silicon", name="NIE")) pot = np.array(devsim.get_node_model_values(device=device_name, region="Silicon", name="Potential")) elec = np.array(devsim.get_node_model_values(device=device_name, region="Silicon", name="Electrons")) holes = np.array(devsim.get_node_model_values(device=device_name, region="Silicon", name="Holes")) changed_count = 0 for c in silicon_contacts: bias_name = f"{c}_bias" try: bias_val = devsim.get_parameter(device=device_name, name=bias_name) except Exception: bias_val = 0.0 try: nodes_list = devsim.get_element_node_list(device=device_name, region="Silicon", contact=c) contact_nodes = sorted(list(set(n for elem in nodes_list for n in elem))) except Exception as e_nodes: print(f" Warning: could not get nodes for contact {c}: {e_nodes}") continue for node in contact_nodes: doping = net_doping[node] ni = nie[node] v_t = vt[node] # Calculate celec and chole celec = 1e-10 + 0.5 * abs(doping + math.sqrt(doping**2 + 4.0 * ni**2)) chole = 1e-10 + 0.5 * abs(-doping + math.sqrt(doping**2 + 4.0 * ni**2)) # Calculate equilibrium values if doping > 0: elec_eq = celec holes_eq = (ni**2) / chole pot_eq = bias_val + v_t * math.log(celec / ni) else: elec_eq = (ni**2) / chole holes_eq = chole pot_eq = bias_val - v_t * math.log(chole / ni) pot[node] = pot_eq elec[node] = elec_eq holes[node] = holes_eq changed_count += 1 if changed_count > 0: devsim.set_node_values(device=device_name, region="Silicon", name="Potential", values=list(pot)) devsim.set_node_values(device=device_name, region="Silicon", name="Electrons", values=list(elec)) devsim.set_node_values(device=device_name, region="Silicon", name="Holes", values=list(holes)) print(f" Enforced boundary values on {changed_count} Silicon contact nodes.") except Exception as ex: print(f" Error enforcing Silicon contact boundary conditions: {ex}") # Oxide and Molding Contacts for reg, contacts in [("Oxide", ["MT1_Ox", "MT2_Ox"]), ("Molding", ["MT1_Mold", "MT2_Mold"])]: try: pot = np.array(devsim.get_node_model_values(device=device_name, region=reg, name="Potential")) changed_count = 0 for c in contacts: bias_name = f"{c}_bias" try: bias_val = devsim.get_parameter(device=device_name, name=bias_name) except Exception: bias_val = 0.0 try: nodes_list = devsim.get_element_node_list(device=device_name, region=reg, contact=c) contact_nodes = sorted(list(set(n for elem in nodes_list for n in elem))) except Exception: continue for node in contact_nodes: pot[node] = bias_val changed_count += 1 if changed_count > 0: devsim.set_node_values(device=device_name, region=reg, name="Potential", values=list(pot)) print(f" Enforced boundary values on {changed_count} {reg} contact nodes.") except Exception as ex: print(f" Error enforcing {reg} contact boundary conditions: {ex}") def refine_and_interpolate(device_old, v_bias, is_avalanche_enabled=False, is_btbt_enabled=False, time_log=None, out_dir=None): if out_dir is None: out_dir = OUT_DIR import time refine_start_time = time.time() device_new_name = f"device_refined_{int(v_bias)}V" print(f"\n--- DYNAMIC REFINE & INTERPOLATE at {v_bias:.2f} V ---") # 1. 確保舊元件中已計算出電場 for reg in ["Silicon", "Oxide", "Molding"]: devsim.element_from_edge_model(edge_model="EField", device=device_old, region=reg) devsim.element_model(device=device_old, region=reg, name="Emag", equation="(EField_x^2 + EField_y^2)^(0.5)") devsim.element_model(device=device_old, region=reg, name="E_mag_log", equation="asinh(Emag / 2.0) / log(10.0)") # 取得 Silicon 中各節點之座標與電場 x_si = np.array(devsim.get_node_model_values(device=device_old, region="Silicon", name="x")) y_si = np.array(devsim.get_node_model_values(device=device_old, region="Silicon", name="y")) # 由於 Emag 是一個 element model，我們手動將其平均到各節點上 emag_elements = np.array(devsim.get_element_model_values(device=device_old, region="Silicon", name="Emag")) triangles = devsim.get_element_node_list(device=device_old, region="Silicon") emag_si = np.zeros(len(x_si)) emag_count = np.zeros(len(x_si)) for i, tri in enumerate(triangles): n0, n1, n2 = tri[0], tri[1], tri[2] emag_si[n0] += emag_elements[3*i] emag_si[n1] += emag_elements[3*i+1] emag_si[n2] += emag_elements[3*i+2] emag_count[n0] += 1 emag_count[n1] += 1 emag_count[n2] += 1 emag_si = np.where(emag_count > 0, emag_si / emag_count, 0.0) # 2. 計算空乏區最深邊界 (E > 2e4 V/cm) dep_nodes = [y_val for y_val, e_val in zip(y_si, emag_si) if e_val > 2.0e4] y_dep = max(dep_nodes) if dep_nodes else 15.0 * um # 預測加密外推 10 µm (Guard Band) y_box_max = min(y_dep + 10.0 * um, H_SI) y_medium_max = min(y_dep + 30.0 * um, H_SI) print(f"Detected depletion max depth: {y_dep / um:.2f} um. Dynamic Box encryption boundary set to: {y_box_max / um:.2f} um, medium boundary set to: {y_medium_max / um:.2f} um.") # 3. 動態寫出 bgmesh.pos LcMin = 0.15 * um LcMax = 20.0 * um alpha = 1.0e-3 bgmesh_pos_path = f"{out_dir}device_bgmesh_refined.pos" with open(bgmesh_pos_path, "w") as f: f.write('View "background mesh" {\n') for reg in ["Silicon", "Oxide"]: x = np.array(devsim.get_node_model_values(device=device_old, region=reg, name="x")) y = np.array(devsim.get_node_model_values(device=device_old, region=reg, name="y")) triangles = np.array(devsim.get_element_node_list(device=device_old, region=reg)) emag = np.array(devsim.get_element_model_values(device=device_old, region=reg, name="Emag"))[::3] for i, tri in enumerate(triangles): n0, n1, n2 = tri[0], tri[1], tri[2] x0, y0 = x[n0], y[n0] x1, y1 = x[n1], y[n1] x2, y2 = x[n2], y[n2] e_val = emag[i] # 計算 lc_val lc_val = LcMax / (1.0 + alpha * e_val) if lc_val < LcMin: lc_val = LcMin f.write(f"ST({x0:.8e},{y0:.8e},0,{x1:.8e},{y1:.8e},0,{x2:.8e},{y2:.8e},0){{{lc_val:.8e},{lc_val:.8e},{lc_val:.8e}}};\n") f.write("};\n") # 4. 動態調用 Gmsh 重建網格 mesh_out_path = f"{out_dir}device_2d_refined.msh" generate_mesh_2d.create_mesh(y_box_max=y_box_max, y_medium_max=y_medium_max, mesh_out=mesh_out_path, bgmesh_pos=bgmesh_pos_path) # 5. 載入新網格並設定物理與 solutions devsim.create_gmsh_mesh(mesh=device_new_name, file=mesh_out_path) opts = setup_physics_for_device(device_new_name, is_avalanche_enabled=is_avalanche_enabled, is_btbt_enabled=is_btbt_enabled) # 6. Apply bias to contacts of the new device for c in ["MT1_Si", "MT1_P12_Si", "MT1_Ox", "MT1_Mold"]: devsim.set_parameter(device=device_new_name, name=f"{c}_bias", value=v_bias) for c in ["MT2_Si", "MT2_P12_Si", "MT2_Ox", "MT2_Mold"]: devsim.set_parameter(device=device_new_name, name=f"{c}_bias", value=0.0) # 7. Interpolate solutions from old device using scipy.interpolate.griddata print("Interpolating solutions to the refined mesh using scipy...") from scipy.interpolate import griddata # 7.1 Gather global old coordinates and potential to ensure 100% continuous interface potential interpolation x_old_si = np.array(devsim.get_node_model_values(device=device_old, region="Silicon", name="x")) y_old_si = np.array(devsim.get_node_model_values(device=device_old, region="Silicon", name="y")) pot_old_si = np.array(devsim.get_node_model_values(device=device_old, region="Silicon", name="Potential")) x_old_ox = np.array(devsim.get_node_model_values(device=device_old, region="Oxide", name="x")) y_old_ox = np.array(devsim.get_node_model_values(device=device_old, region="Oxide", name="y")) pot_old_ox = np.array(devsim.get_node_model_values(device=device_old, region="Oxide", name="Potential")) x_old_mold = np.array(devsim.get_node_model_values(device=device_old, region="Molding", name="x")) y_old_mold = np.array(devsim.get_node_model_values(device=device_old, region="Molding", name="y")) pot_old_mold = np.array(devsim.get_node_model_values(device=device_old, region="Molding", name="Potential")) points_old_global = np.column_stack(( np.concatenate((x_old_si, x_old_ox, x_old_mold)), np.concatenate((y_old_si, y_old_ox, y_old_mold)) )) pot_old_global = np.concatenate((pot_old_si, pot_old_ox, pot_old_mold)) # 7.2 Interpolate Potential to all regions of the new refined mesh using the global dataset for reg in ["Silicon", "Oxide", "Molding"]: x_new = np.array(devsim.get_node_model_values(device=device_new_name, region=reg, name="x")) y_new = np.array(devsim.get_node_model_values(device=device_new_name, region=reg, name="y")) points_new = np.column_stack((x_new, y_new)) pot_new = griddata(points_old_global, pot_old_global, points_new, method='linear') if np.isnan(pot_new).any(): pot_new_nearest = griddata(points_old_global, pot_old_global, points_new, method='nearest') pot_new = np.where(np.isnan(pot_new), pot_new_nearest, pot_new) devsim.set_node_values(device=device_new_name, region=reg, name="Potential", values=list(pot_new)) # 7.3 Volume-weighted conservative interpolation of Silicon carriers # Strategy: Instead of log-space re-sampling (Boltzmann-consistent but not charge-conserving), # use volume-weighted barycentric interpolation to conserve total charge from M0 to M1. # Formula: n_j = Σ_i (λ_i * n_i * V_i^M0) / Σ_i (λ_i * V_i^M0) # where λ_i are barycentric coords of M1 node j within M0 triangle, # and V_i^M0 are the Voronoi volumes of M0 nodes. print("Interpolating Electrons and Holes (log-space + global charge scaling)...") points_old_si = np.column_stack((x_old_si, y_old_si)) x_new_si = np.array(devsim.get_node_model_values(device=device_new_name, region="Silicon", name="x")) y_new_si = np.array(devsim.get_node_model_values(device=device_new_name, region="Silicon", name="y")) points_new_si = np.column_stack((x_new_si, y_new_si)) # Get Voronoi volumes from M0 (BEFORE deleting old device) V_old_si = np.array(devsim.get_node_model_values(device=device_old, region="Silicon", name="NodeVolume")) # Build Delaunay triangulation of M0 Silicon nodes from scipy.spatial import Delaunay tri_si = Delaunay(points_old_si) def log_space_interp(val_old, tri, points_new): """Log-space barycentric interpolation for Boltzmann consistency.""" log_val_old = np.log(np.maximum(val_old, 1e-40)) # Find containing simplex for each M1 node simplex_idx = tri.find_simplex(points_new) # Compute barycentric coordinates T = tri.transform # shape (N_simplices, ndim, ndim) r = points_new - T[simplex_idx, 2] # translate to simplex origin bary_partial = np.einsum('nij,nj->ni', T[simplex_idx, :2], r) # (N_new, 2) bary = np.concatenate([bary_partial, 1.0 - bary_partial.sum(axis=1, keepdims=True)], axis=1) # (N_new, 3) bary = np.clip(bary, 0.0, 1.0) bary /= bary.sum(axis=1, keepdims=True) # Vertex indices of each simplex verts = tri.simplices[simplex_idx] # (N_new, 3) log_val_verts = log_val_old[verts] # (N_new, 3) result_log = (bary * log_val_verts).sum(axis=1) # Handle M1 nodes outside M0 convex hull (simplex_idx == -1) outside = simplex_idx < 0 if outside.any(): from scipy.interpolate import griddata result_nearest = griddata(points_old_si, log_val_old, points_new[outside], method='nearest') result_log[outside] = result_nearest return np.exp(result_log) # Fetch old carrier values and new/old NodeVolumes to perform and verify interpolation elec_old_si = np.array(devsim.get_node_model_values(device=device_old, region="Silicon", name="Electrons")) holes_old_si = np.array(devsim.get_node_model_values(device=device_old, region="Silicon", name="Holes")) V_new_si = np.array(devsim.get_node_model_values(device=device_new_name, region="Silicon", name="NodeVolume")) # Perform log-space interpolation electrons_interp = log_space_interp(elec_old_si, tri_si, points_new_si) holes_interp = log_space_interp(holes_old_si, tri_si, points_new_si) # Scale interpolated concentrations globally to enforce exact charge conservation Q_m0_n = float(np.dot(elec_old_si, V_old_si)) Q_m0_p = float(np.dot(holes_old_si, V_old_si)) Q_before_n = float(np.dot(electrons_interp, V_new_si)) Q_before_p = float(np.dot(holes_interp, V_new_si)) if Q_before_n > 0: electrons_interp *= (Q_m0_n / Q_before_n) if Q_before_p > 0: holes_interp *= (Q_m0_p / Q_before_p) print(f" Interpolation (before scaling) — Electrons ratio: {Q_before_n/Q_m0_n:.6f} (M0: {Q_m0_n:.4e}, M1: {Q_before_n:.4e})") print(f" Interpolation (before scaling) — Holes ratio: {Q_before_p/Q_m0_p:.6f} (M0: {Q_m0_p:.4e}, M1: {Q_before_p:.4e})") print(f" After scaling charge conservation — Electrons ratio: {float(np.dot(electrons_interp, V_new_si))/Q_m0_n:.6f}") print(f" After scaling charge conservation — Holes ratio: {float(np.dot(holes_interp, V_new_si))/Q_m0_p:.6f}") def slope_preserving_conservative_smooth(n, V, pts, tri_m0, N_passes=5, alpha=0.3): """ Slope-preserving, charge-conserving smoothing (user's 11.3-11.4 algorithm). For each edge (i,j) in the Delaunay triangulation: 1. Estimate per-node local gradient [gx, gy] via inverse-distance-weighted LSQ 2. dn_expected = gradient_midpoint · (r_j - r_i) (what the slope predicts) 3. residual = (n_j - n_i) - dn_expected (deviation from local trend) 4. flux = alpha * residual * V_i*V_j/(V_i+V_j) (antisymmetric → conserved) Properties: - Exactly charge-conserving: Σ flux = 0 for each edge - Slope-preserving: if n follows local gradient, residual=0, no flux - Reduces 2nd-order deviations (peaks/valleys above/below local trend) - At depletion boundary: n_C=5e16 is ABOVE linear extrapolation → reduced toward A """ n = n.copy().astype(np.float64) N = len(n) # Build unique edge list from M0 Delaunay (used for M1 gradient; M1 tri not yet built) # For M1 smoothing we need M1 edges. Build M1 Delaunay from pts. from scipy.spatial import Delaunay as Delaunay_local tri_m1 = Delaunay_local(pts) edge_set = set() for s in tri_m1.simplices: a, b, c = s[0], s[1], s[2] edge_set.add((min(a,b), max(a,b))) edge_set.add((min(b,c), max(b,c))) edge_set.add((min(a,c), max(a,c))) ei_arr = np.array([e[0] for e in edge_set], dtype=np.int64) ej_arr = np.array([e[1] for e in edge_set], dtype=np.int64) dx = pts[ej_arr, 0] - pts[ei_arr, 0] dy = pts[ej_arr, 1] - pts[ei_arr, 1] d = np.sqrt(dx**2 + dy**2) inv_d = 1.0 / np.maximum(d, 1e-30) print(f" [smooth] M1 edges: {len(ei_arr)}, nodes: {N}, passes: {N_passes}, alpha: {alpha}") for pass_idx in range(N_passes): dn = n[ej_arr] - n[ei_arr] # actual concentration difference along each edge # --- Gradient estimation at each node via weighted LSQ --- # From edge (i,j): at node i, neighbor j is at (dx, dy) with value difference dn # gx_i * dx + gy_i * dy ≈ dn # From edge (i,j): at node j, neighbor i is at (-dx,-dy) with value difference -dn # gx_j * (-dx) + gy_j * (-dy) ≈ -dn → same normal equations (symmetric) w = inv_d A11 = np.bincount(ei_arr, w*dx**2, N) + np.bincount(ej_arr, w*dx**2, N) A12 = np.bincount(ei_arr, w*dx*dy, N) + np.bincount(ej_arr, w*dx*dy, N) A22 = np.bincount(ei_arr, w*dy**2, N) + np.bincount(ej_arr, w*dy**2, N) B1 = np.bincount(ei_arr, w*dx*dn, N) + np.bincount(ej_arr, w*dx*dn, N) B2 = np.bincount(ei_arr, w*dy*dn, N) + np.bincount(ej_arr, w*dy*dn, N) det = A11*A22 - A12**2 safe = np.abs(det) > 1e-60 det_s = np.where(safe, det, 1.0) gx = np.where(safe, ( A22*B1 - A12*B2) / det_s, 0.0) gy = np.where(safe, (-A12*B1 + A11*B2) / det_s, 0.0) # --- Edge-level gradient-corrected flux --- gx_mid = 0.5 * (gx[ei_arr] + gx[ej_arr]) gy_mid = 0.5 * (gy[ei_arr] + gy[ej_arr]) dn_expected = gx_mid * dx + gy_mid * dy # residual > 0 → n_j is HIGHER than gradient predicts → j should give charge to i residual = dn - dn_expected # Harmonic-volume weighted flux (antisymmetric → charge conserved) V_harm = V[ei_arr] * V[ej_arr] / (V[ei_arr] + V[ej_arr]) flux = alpha * residual * V_harm # charge [cm^-3 * cm^3] # i gains flux, j loses flux (antisymmetric) delta_charge = (np.bincount(ei_arr, flux, N) + np.bincount(ej_arr, -flux, N)) n = np.clip(n + delta_charge / V, 1e-40, None) return n # Apply slope-preserving conservative smooth on M1 grid # Bypassed because log-space interpolation is Boltzmann-consistent and inherently smooth. # Smoothing in linear space destroys the exponential carrier profiles across PN junctions. print(" Bypassing slope-preserving conservative smooth (log-space interpolation is Boltzmann-consistent)...") # Set on new mesh devsim.set_node_values(device=device_new_name, region="Silicon", name="Electrons", values=list(electrons_interp)) devsim.set_node_values(device=device_new_name, region="Silicon", name="Holes", values=list(holes_interp)) print(" Log-space interpolated and scaled Electrons/Holes set on new mesh.") enforce_contact_boundary_conditions(device_new_name) # 8. Destroy old device to release memory BEFORE solving print("Deleting old device to avoid simultaneous solving in DEVSIM...") devsim.delete_device(device=device_old) print("Old device deleted.") av_model = "AvalancheGeneration" if is_avalanche_enabled else "" # ========================================== # Stage 0.5: Nonlinear Poisson pre-solve with frozen carriers & Debye screening # ========================================== print("[Stage 0.5] Running nonlinear Poisson pre-solve with Debye screening...") # 1. Store reference solutions for Debye screening pot_ref_si = np.array(devsim.get_node_model_values(device=device_new_name, region="Silicon", name="Potential")) devsim.node_solution(device=device_new_name, region="Silicon", name="Electrons_ref") devsim.node_solution(device=device_new_name, region="Silicon", name="Holes_ref") devsim.node_solution(device=device_new_name, region="Silicon", name="Potential_ref") devsim.set_node_values(device=device_new_name, region="Silicon", name="Electrons_ref", values=list(electrons_interp)) devsim.set_node_values(device=device_new_name, region="Silicon", name="Holes_ref", values=list(holes_interp)) devsim.set_node_values(device=device_new_name, region="Silicon", name="Potential_ref", values=list(pot_ref_si)) # 2. Temporarily replace PotentialNodeCharge with nonlinear Gummel-style Poisson model devsim.node_model(device=device_new_name, region="Silicon", name="dV_gummel", equation="(Potential - Potential_ref) / V_t") devsim.node_model(device=device_new_name, region="Silicon", name="dV_gummel_clipped", equation="ifelse(dV_gummel > 50, 50, ifelse(dV_gummel < -50, -50, dV_gummel))") devsim.node_model(device=device_new_name, region="Silicon", name="PotentialNodeCharge", equation="-q * (Holes_ref * exp(-dV_gummel_clipped) - Electrons_ref * exp(dV_gummel_clipped) + NetDoping)") devsim.node_model(device=device_new_name, region="Silicon", name="dV_gummel_clipped:Potential", equation="ifelse(dV_gummel > 50, 0, ifelse(dV_gummel < -50, 0, 1 / V_t))") devsim.node_model(device=device_new_name, region="Silicon", name="PotentialNodeCharge:Potential", equation="-q * (-Holes_ref * exp(-dV_gummel_clipped) * dV_gummel_clipped:Potential - Electrons_ref * exp(dV_gummel_clipped) * dV_gummel_clipped:Potential)") # Since carrier equations are deleted, the derivatives w.r.t Electrons and Holes are set to 0 devsim.node_model(device=device_new_name, region="Silicon", name="PotentialNodeCharge:Electrons", equation="0") devsim.node_model(device=device_new_name, region="Silicon", name="PotentialNodeCharge:Holes", equation="0") # 3. Delete continuity equations from Silicon region and contacts devsim.delete_equation(device=device_new_name, region="Silicon", name="ElectronContinuityEquation") devsim.delete_equation(device=device_new_name, region="Silicon", name="HoleContinuityEquation") for c in ["MT1_Si", "MT2_Si", "MT1_P12_Si", "MT2_P12_Si"]: devsim.delete_contact_equation(device=device_new_name, contact=c, name="ElectronContinuityEquation") devsim.delete_contact_equation(device=device_new_name, contact=c, name="HoleContinuityEquation") # Temporarily set PotentialEquation variable update to log_damp devsim.equation(device=device_new_name, region="Silicon", name="PotentialEquation", variable_name="Potential", node_model="PotentialNodeCharge", edge_model="DField", variable_update="log_damp", min_error=1e-3) iters_s0 = 0 iters_s1 = 0 iters_s2 = 0 try: # Print Potential statistics before Poisson solve for reg in ["Silicon", "Oxide", "Molding"]: pot = np.array(devsim.get_node_model_values(device=device_new_name, region=reg, name="Potential")) print(f" [pre-solve] Potential in {reg} before solve: min={pot.min():.2e}, max={pot.max():.2e}, mean={pot.mean():.2e}") # 4. Solve nonlinear Poisson equation. We use a relaxed tolerance (2.0) and 6 iterations to ensure successful convergence and update. print(" Solving nonlinear Poisson with Debye screening...") res_s0 = devsim.solve(type="dc", absolute_error=1e10, relative_error=2.0, charge_error=1e12, maximum_iterations=6, info=True) iters_s0 = len(res_s0.get("iterations", [])) if (isinstance(res_s0, dict) and "iterations" in res_s0) else 6 print(" Nonlinear Poisson solve converged successfully.") # Print Potential statistics after Poisson solve for reg in ["Silicon", "Oxide", "Molding"]: pot = np.array(devsim.get_node_model_values(device=device_new_name, region=reg, name="Potential")) print(f" [pre-solve] Potential in {reg} after solve: min={pot.min():.2e}, max={pot.max():.2e}, mean={pot.mean():.2e}") except devsim.error as e: print(f" Nonlinear Poisson solve failed: {e}. Proceeding anyway...") # Bypass decoupled carrier solve and Stage 0.5.2 to preserve the high-quality interpolated and smoothed carriers. # Decoupled Gummel carriers solve is unstable at high reverse bias (500V+) and destroys the carriers by setting them to zero. elec_new = electrons_interp holes_new = holes_interp print("[Stage 0.5.2] Bypassed decoupled carrier initialization and Stage 0.5.2 to preserve charge-conserved carriers.") # 5. Restore original PotentialNodeCharge and its derivatives for drift-diffusion print(" Restoring original PotentialNodeCharge and derivatives...") devsim.node_model(device=device_new_name, region="Silicon", name="PotentialNodeCharge", equation="-q * kahan3(Holes, -Electrons, NetDoping)") devsim.node_model(device=device_new_name, region="Silicon", name="PotentialNodeCharge:Potential", equation="0") devsim.node_model(device=device_new_name, region="Silicon", name="PotentialNodeCharge:Electrons", equation="q") devsim.node_model(device=device_new_name, region="Silicon", name="PotentialNodeCharge:Holes", equation="-q") # Restore PotentialEquation variable update to log_damp devsim.equation(device=device_new_name, region="Silicon", name="PotentialEquation", variable_name="Potential", node_model="PotentialNodeCharge", edge_model="DField", variable_update="log_damp", min_error=1e-3) print("Solving refined mesh: Stage 1 fully-coupled log_damp (Poisson pre-solved), Stage 2 positive...") # ========================================== # Stage 1: Fully-coupled log_damp # ========================================== 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: from physics.model_create import CreateEdgeModel, CreateEdgeModelDerivatives CreateEdgeModel(device_new_name, "Silicon", "CombinedGeneration", "AvalancheGeneration + BTBTGeneration") CreateEdgeModel(device_new_name, "Silicon", "CombinedGeneration_p", "AvalancheGeneration_p + BTBTGeneration_p") for i in ("Potential", "Electrons", "Holes"): CreateEdgeModelDerivatives(device_new_name, "Silicon", "CombinedGeneration", "AvalancheGeneration + BTBTGeneration", i) CreateEdgeModelDerivatives(device_new_name, "Silicon", "CombinedGeneration_p", "AvalancheGeneration_p + BTBTGeneration_p", i) gen_model_n = "CombinedGeneration" gen_model_p = "CombinedGeneration_p" elif av_model_n: gen_model_n = av_model_n gen_model_p = av_model_p elif btbt_model_n: gen_model_n = btbt_model_n gen_model_p = btbt_model_p else: gen_model_n = "" gen_model_p = "" devsim.equation(device=device_new_name, region="Silicon", name="ElectronContinuityEquation", variable_name="Electrons", time_node_model="NCharge", edge_model=opts['Jn'], edge_volume_model=gen_model_n, variable_update="log_damp", node_model="ElectronGeneration", min_error=1e5) devsim.equation(device=device_new_name, region="Silicon", name="HoleContinuityEquation", variable_name="Holes", time_node_model="PCharge", edge_model=opts['Jp'], edge_volume_model=gen_model_p, variable_update="log_damp", node_model="HoleGeneration", min_error=1e5) for c in ["MT1_Si", "MT2_Si", "MT1_P12_Si", "MT2_P12_Si"]: contact_electrons_name = f"{c}nodeelectrons" contact_holes_name = f"{c}nodeholes" devsim.contact_equation(device=device_new_name, contact=c, name="ElectronContinuityEquation", node_model=contact_electrons_name, edge_current_model=opts['Jn']) devsim.contact_equation(device=device_new_name, contact=c, name="HoleContinuityEquation", node_model=contact_holes_name, edge_current_model=opts['Jp']) # Set carrier values again to ensure they are not reset/cleared by equation re-registration. # We update Electrons and Holes to their Gummel-adjusted Boltzmann values based on the solved Potential. try: pot = np.array(devsim.get_node_model_values(device=device_new_name, region="Silicon", name="Potential")) pot_ref = np.array(devsim.get_node_model_values(device=device_new_name, region="Silicon", name="Potential_ref")) vt = np.array(devsim.get_node_model_values(device=device_new_name, region="Silicon", name="V_t")) dv_clipped = np.clip((pot - pot_ref) / vt, -50.0, 50.0) electrons_updated = np.clip(elec_new * np.exp(dv_clipped), 1e-40, None) holes_updated = np.clip(holes_new * np.exp(-dv_clipped), 1e-40, None) print(f" Boltzmann carrier update applied. dV_gummel_clipped: min={dv_clipped.min():.4f}, max={dv_clipped.max():.4f}, mean={dv_clipped.mean():.4f}") print(f" Electrons sum: before={elec_new.sum():.4e}, after={electrons_updated.sum():.4e}") print(f" Holes sum: before={holes_new.sum():.4e}, after={holes_updated.sum():.4e}") except Exception as e_boltz: print(f" Failed to calculate Boltzmann updated carriers: {e_boltz}. Falling back to original carrier values.") electrons_updated = elec_new holes_updated = holes_new devsim.set_node_values(device=device_new_name, region="Silicon", name="Electrons", values=list(electrons_updated)) devsim.set_node_values(device=device_new_name, region="Silicon", name="Holes", values=list(holes_updated)) print(" Electrons/Holes values restored and Boltzmann-adjusted after equation registration.") enforce_contact_boundary_conditions(device_new_name) print("[Stage 1] Running fully-coupled log_damp solve...") # Change PotentialEquation variable_update to default for Stage 1 pre-conditioning stability devsim.equation(device=device_new_name, region="Silicon", name="PotentialEquation", variable_name="Potential", node_model="PotentialNodeCharge", edge_model="DField", variable_update="default", min_error=1e-3) for reg in ["Oxide", "Molding"]: devsim.equation(device=device_new_name, region=reg, name="PotentialEquation", variable_name="Potential", edge_model="PotentialEdgeFlux", variable_update="default", min_error=1e-3) try: res_s1 = devsim.solve(type="dc", absolute_error=1e10, relative_error=1e-3, charge_error=1e12, maximum_iterations=12, rollback=False, info=True) iters_s1 = len(res_s1.get("iterations", [])) if (isinstance(res_s1, dict) and "iterations" in res_s1) else 12 print(" Stage 1 Coupled log_damp solve finished.") except devsim.error as e: print(f" Stage 1 Coupled log_damp solve errored: {e}. Proceeding to Stage 2 anyway.") res_s1 = {"converged": False} iters_s1 = 20 # Clip any negative Electrons/Holes that log_damp Stage 1 may have introduced # (log_damp can sometimes produce negative values in extreme depletion regions due to # large Potential updates interacting with carrier boundary conditions) print("[Stage 1→2] Clipping negative Electrons/Holes to floor (1e-40) before Stage 2...") try: n_vals = np.array(devsim.get_node_model_values(device=device_new_name, region="Silicon", name="Electrons")) p_vals = np.array(devsim.get_node_model_values(device=device_new_name, region="Silicon", name="Holes")) n_clipped = np.clip(n_vals, 1e-40, None) p_clipped = np.clip(p_vals, 1e-40, None) if (n_clipped != n_vals).any() or (p_clipped != p_vals).any(): print(f" Clipped {(n_clipped != n_vals).sum()} Electrons nodes, {(p_clipped != p_vals).sum()} Holes nodes") devsim.set_node_values(device=device_new_name, region="Silicon", name="Electrons", values=list(n_clipped)) devsim.set_node_values(device=device_new_name, region="Silicon", name="Holes", values=list(p_clipped)) except Exception as e_clip: print(f" Clip step failed: {e_clip} (proceeding anyway)") # ========================================== # Stage 2: log_damp coupled precision solve # ========================================== 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: from physics.model_create import CreateEdgeModel, CreateEdgeModelDerivatives CreateEdgeModel(device_new_name, "Silicon", "CombinedGeneration", "AvalancheGeneration + BTBTGeneration") CreateEdgeModel(device_new_name, "Silicon", "CombinedGeneration_p", "AvalancheGeneration_p + BTBTGeneration_p") for i in ("Potential", "Electrons", "Holes"): CreateEdgeModelDerivatives(device_new_name, "Silicon", "CombinedGeneration", "AvalancheGeneration + BTBTGeneration", i) CreateEdgeModelDerivatives(device_new_name, "Silicon", "CombinedGeneration_p", "AvalancheGeneration_p + BTBTGeneration_p", i) gen_model_n = "CombinedGeneration" gen_model_p = "CombinedGeneration_p" elif av_model_n: gen_model_n = av_model_n gen_model_p = av_model_p elif btbt_model_n: gen_model_n = btbt_model_n gen_model_p = btbt_model_p else: gen_model_n = "" gen_model_p = "" devsim.equation(device=device_new_name, region="Silicon", name="ElectronContinuityEquation", variable_name="Electrons", time_node_model="NCharge", edge_model=opts['Jn'], edge_volume_model=gen_model_n, variable_update="positive", node_model="ElectronGeneration", min_error=1e5) devsim.equation(device=device_new_name, region="Silicon", name="HoleContinuityEquation", variable_name="Holes", time_node_model="PCharge", edge_model=opts['Jp'], edge_volume_model=gen_model_p, variable_update="positive", node_model="HoleGeneration", min_error=1e5) # Restore PotentialEquation variable update to default with min_error=1e-3 devsim.equation(device=device_new_name, region="Silicon", name="PotentialEquation", variable_name="Potential", node_model="PotentialNodeCharge", edge_model="DField", variable_update="default", min_error=1e-3) for reg in ["Oxide", "Molding"]: devsim.equation(device=device_new_name, region=reg, name="PotentialEquation", variable_name="Potential", edge_model="PotentialEdgeFlux", variable_update="default", min_error=1e-3) print("[Stage 2] Running log_damp coupled precision solve (max 30 iters)...") try: res_stage2 = devsim.solve(type="dc", absolute_error=1e10, relative_error=1e-3, charge_error=1e12, maximum_iterations=30, info=True) iters_s2 = len(res_stage2.get("iterations", [])) if (isinstance(res_stage2, dict) and "iterations" in res_stage2) else 30 # We check both absolute error and charge error convergence through devsim solve. # If solve completed without raising exception, we assume it achieved correct physical state. if res_stage2 is None: res_stage2 = {"converged": True} elif not res_stage2.get("converged", False): # Fallback check if info=True dict indicates failure raise devsim.error("Stage 2 solve did not achieve convergence status.") # Log successful refinement convergence to simulation_time.log try: import time time_taken = time.time() - refine_start_time if time_log is not None: time_log.write(f"{time.strftime('%X')}\t{v_bias:.4f}\t[Refine]\t-\t{iters_s0}+{iters_s1}+{iters_s2}\t{time_taken:.2f}s\tSuccess\n") time_log.flush() else: with open(f"{out_dir}simulation_time.log", "a") as local_log: local_log.write(f"{time.strftime('%X')}\t{v_bias:.4f}\t[Refine]\t-\t{iters_s0}+{iters_s1}+{iters_s2}\t{time_taken:.2f}s\tSuccess\n") except Exception as e_log: print(f"Warning: Failed to log refinement progress to simulation_time.log: {e_log}") except devsim.error as e: print(f"Refinement Precision Newton solve failed/errored: {e}! Recreating and restoring old device '{device_old}'...") # Log failure to simulation_time.log try: import time time_taken = time.time() - refine_start_time if time_log is not None: time_log.write(f"{time.strftime('%X')}\t{v_bias:.4f}\t[Refine]\t-\tStage 2 Failed\t{time_taken:.2f}s\tReverted to Coarse\n") time_log.flush() else: with open(f"{out_dir}simulation_time.log", "a") as local_log: local_log.write(f"{time.strftime('%X')}\t{v_bias:.4f}\t[Refine]\t-\tStage 2 Failed\t{time_taken:.2f}s\tReverted to Coarse\n") except Exception as e_log: pass try: if devsim.get_device_list() and device_new_name in devsim.get_device_list(): devsim.delete_device(device=device_new_name) except Exception: pass # Recreate old device to allow restore_state to work try: devsim.create_gmsh_mesh(mesh=device_old, file="device_2d.msh") except Exception: pass setup_physics_for_device(device_old, is_avalanche_enabled=is_avalanche_enabled, is_btbt_enabled=is_btbt_enabled) raise devsim.error(f"Precision Newton solve on refined mesh did not converge: {e}") print("Convergence on refined mesh achieved successfully!") return device_new_name, opts