feat: Stabilizing High-Voltage BJT Simulation

feat: Remove .msh file generation for Paraview
feat: Add PROJECT GUIDE.md.
2026-03-27 16:13:23 +08:00 · 2026-03-19 16:30:19 +08:00 · 2026-03-16 14:48:02 +08:00 · 2026-03-16 09:51:13 +08:00 · 2026-03-13 16:01:52 +08:00 · 2026-03-13 15:13:51 +08:00
38 changed files with 1044 additions and 78 deletions
@@ -0,0 +1,280 @@
+# DEVSIM TCAD 專案 — 完整安裝與執行流程
+
+## 專案概覽
+
+**DEVSIM** 是一套開源的 TCAD（Technology Computer-Aided Design）元件模擬器，採用有限體積法（Finite Volume Method）求解半導體方程式，透過 Python 腳本驅動模擬。
+
+本專案路徑為 `/home/wsx52114/devsim`，目錄結構如下：
+
+```
+devsim/                           # 專案根目錄
+├── denv/                         # Python 虛擬環境（已建立）
+├── devsim/                       # DEVSIM 原始碼 + 自訂模擬
+│   ├── wisetop_bjt/              # ★ 自訂 BJT 模擬專案（主要工作區）
+│   ├── wisetop_opto/             # 光電元件模擬
+│   ├── examples/                 # 官方範例（diode, capacitance 等）
+│   ├── testing/                  # 官方測試腳本
+│   ├── src/                      # C++ 原始碼
+│   ├── scripts/                  # 建置腳本
+│   └── INSTALL.md, BUILD.md ...  # 官方文件
+└── devsim_bjt_example-main/      # 官方 BJT 範例（論文用）
+```
+
+---
+
+## 一、環境安裝
+
+### 1.1 系統需求
+
+| 軟體 | 最低版本 | 用途 | 目前已安裝版本 |
+|------|---------|------|---------------|
+| **Python** | ≥ 3.9 | 模擬核心引擎 | 3.12.3 |
+| **Gmsh** | ≥ 4.0 | 有限元素網格生成 | 4.15.0 |
+| **ParaView** | 選用 | `.tec` 結果視覺化 | — |
+
+### 1.2 安裝步驟（從零開始）
+
+```bash
+# ── Step 1：安裝系統套件 ──
+sudo apt update
+sudo apt install python3 python3-venv gmsh
+
+# ── Step 2：Clone 專案 ──
+git clone https://gitlab.com/wisetop/devsim-tcad/devsim.git
+cd devsim
+
+# ── Step 3：建立 Python 虛擬環境 ──
+python3 -m venv denv
+
+# ── Step 4：啟用虛擬環境 ──
+source denv/bin/activate
+
+# ── Step 5：安裝 Python 套件 ──
+pip install devsim numpy
+# 選用：安裝 MKL（提升線性代數運算效能）
+pip install mkl
+```
+
+> [!IMPORTANT]
+> **每次開啟新終端都必須先啟用虛擬環境**：`source denv/bin/activate`
+
+### 1.3 目前已安裝的套件
+
+本專案虛擬環境 `denv/` 中已安裝：
+
+| 套件 | 版本 | 說明 |
+|------|------|------|
+| `devsim` | 2.10.0 | TCAD 模擬核心 |
+| `gmsh` | 4.15.0 | Gmsh Python API |
+| `numpy` | 2.4.0 | 數值計算 |
+| `mkl` | 2025.3.0 | Intel Math Kernel Library |
+| `intel_openmp` | 2025.3.1 | 多執行緒支援 |
+
+### 1.4 驗證安裝
+
+```bash
+source denv/bin/activate
+
+# 測試 devsim
+python3 -c "import devsim; print('devsim OK')"
+
+# 測試 gmsh
+gmsh --version
+
+# 跑官方範例
+cd devsim/testing
+python3 cap2.py
+```
+
+---
+
+## 二、模擬執行流程
+
+### 2.1 主要工作區：`wisetop_bjt/`
+
+該目錄包含自訂的平面 NPN BJT 模擬，已封裝為 **Makefile** 自動化流程。
+
+```bash
+# 進入工作目錄
+cd /home/wsx52114/devsim/devsim/wisetop_bjt
+
+# 啟用虛擬環境（若尚未啟用）
+source /home/wsx52114/devsim/denv/bin/activate
+```
+
+### 2.2 模擬流程圖
+
+```
+┌─────────┐    ┌─────────┐    ┌─────────┐    ┌──────────────────┐
+│  mesh   │───▶│ refine  │───▶│  init   │───▶│ gummel / output  │
+│ (Gmsh)  │    │ (迴圈)  │    │ (DC=0V) │    │   / cb / ac      │
+└─────────┘    └─────────┘    └─────────┘    └──────────────────┘
+  bjt.msh    bjt_refined.msh  bjt_dd_0.msh       *.csv
+```
+
+### 2.3 Make 指令一覽
+
+#### 完整流程指令
+
+| 指令 | 說明 | 包含步驟 |
+|------|------|---------|
+| `make all` | **完整模擬**（推薦） | refine → init → gummel → output → cb → ac |
+| `make quick` | 快速模擬（無網格細化） | mesh → init → gummel |
+| `make clean` | 清除所有生成檔案 | 刪除 `*.msh *.pos *.tec *.csv *.log` |
+
+#### 個別步驟指令
+
+| 指令 | 說明 | 輸入 | 輸出 |
+|------|------|------|------|
+| `make mesh` | 生成初始網格 | `bjt.geo` | `bjt.msh` |
+| `make refine` | 網格細化迴圈 | `bjt.geo` | `bjt_refined.msh` |
+| `make init` | 零偏壓 DD 初始化 | `bjt_refined.msh` | `bjt_dd_0.msh`, `bjt_dd_0.tec` |
+| `make gummel` | Gummel 曲線 (Ic,Ib vs Vb) | `bjt_dd_0.msh` | `gummel.csv` |
+| `make output` | 輸出特性 (Ic vs Vc) | `bjt_dd_0.msh` | `output_Vb0.70.csv` |
+| `make cb` | 共基極特性 (Ic vs Ve) | `bjt_dd_0.msh` | `cb_Vc0.50.csv` |
+| `make ac` | AC 小信號分析 | `bjt_dd_0.msh` | `ac_Vc0.5_Ve-0.7.csv` |
+| `make batch` | 批次多參數掃描 | `bjt_dd_0.msh` | `data/*.log` |
+
+#### 環境變數調整
+
+```bash
+# 設定網格細化迭代次數（預設 1）
+make all LOOPS=3
+
+# 指定 Gmsh 執行緒數
+make all NPROC=8
+
+# 修改偏壓參數
+make output VB_DEFAULT=0.8
+make cb VC_DEFAULT=1.0
+```
+
+---
+
+## 三、各步驟詳細說明
+
+### 3.1 網格生成 (`make mesh`)
+
+```bash
+gmsh -nt <NPROC> -2 -format msh2 bjt.geo -o bjt.msh
+```
+
+- 讀取 `bjt.geo` 幾何定義檔
+- 生成 2D 網格，格式為 MSH2
+- `-nt` 啟用多執行緒
+
+### 3.2 網格細化 (`make refine`)
+
+迭代流程（重複 LOOPS 次）：
+
+1. 用 `bjt_refine.py` 分析電場，計算各區域需要的網格密度
+2. 輸出背景場檔 `bjt_bg.pos`
+3. Gmsh 依據背景場重新生成網格
+4. 最終結果複製為 `bjt_refined.msh`
+
+### 3.3 零偏壓初始化 (`make init`)
+
+```bash
+python3 bjt_dd.py
+```
+
+- 載入細化後的網格
+- 設定 ERFC 摻雜模型
+- 求解 Poisson 方程（Potential Only → Drift-Diffusion）
+- 零偏壓 DC 求解
+- 輸出 `bjt_dd_0.msh`（可續算狀態）和 `bjt_dd_0.tec`（視覺化用）
+
+### 3.4 電氣特性分析
+
+| 分析 | 腳本 | 說明 |
+|------|------|------|
+| **Gummel** | `bjt_circuit3.py` | 固定 Vc=2V，掃描 Vb 0→0.8V（20mV 步進） |
+| **輸出特性** | `bjt_circuit2.py <Vb>` | 固定 Vb，掃描 Vc 0→2V（0.1V 步進） |
+| **共基極** | `bjt_circuit4.py <Vc>` | 固定 Vc，掃描 Ve 0→-1V（-50mV 步進） |
+| **AC 分析** | `bjt_circuit5.py <Vc> <fmin> <fmax> <ppd>` | AC 小信號，頻率 1kHz→100GHz |
+
+---
+
+## 四、結果檢視
+
+```bash
+# 網格視覺化
+gmsh bjt_refined.msh
+
+# 摻雜/電位分佈（需 ParaView）
+paraview bjt_dd_0.tec
+# 載入後選擇 Variable: LogNetDoping 或 Potential
+
+# CSV 數據
+# 可用 Python/matplotlib 繪圖，或用試算表軟體開啟
+```
+
+---
+
+## 五、官方範例
+
+### 5.1 內建範例
+
+位於 [devsim/examples/](file:///home/wsx52114/devsim/devsim/examples)：
+- `diode/` — 二極體模擬
+- `capacitance/` — 電容模擬
+- `mobility/` — 遷移率模型
+- `bioapp1/` — 生物應用
+- `plotting/` — 繪圖範例
+
+### 5.2 BJT 論文範例
+
+位於 [devsim_bjt_example-main/](file:///home/wsx52114/devsim/devsim_bjt_example-main)，對應論文 *Semiconductor Device Simulation Using DEVSIM*。
+
+---
+
+## 六、從原始碼建置（進階）
+
+若需修改 DEVSIM 核心，可從原始碼建置：
+
+```bash
+# Clone
+git clone https://github.com/devsim/devsim
+cd devsim
+git submodule init
+git submodule update
+
+# Linux 建置
+bash scripts/build_manylinux_2_28.sh <version>
+
+# 或使用 Docker
+bash scripts/build_docker_manylinux_2_28.sh <version>
+
+# 建置後的 .whl 檔在 dist/ 目錄
+pip install dist/devsim-<version>*.whl
+```
+
+---
+
+## 七、疑難排解
+
+| 問題 | 解決方案 |
+|------|---------|
+| `ModuleNotFoundError: devsim` | 確認虛擬環境已啟用：`source denv/bin/activate` |
+| `devsim` 未安裝 | `pip install devsim` |
+| `gmsh: command not found` | `sudo apt install gmsh` 或 `pip install gmsh` |
+| 模擬收斂失敗 | `make clean && make refine && make init` 重新開始 |
+| 摻雜分佈異常 | 用 ParaView 檢視 `bjt_dd_0.tec` 的 `LogNetDoping` |
+| 記憶體不足 | 減少網格細化次數：`make all LOOPS=1` |
+
+---
+
+## 八、快速開始指令摘要
+
+```bash
+# ── 完整流程（一行搞定）──
+cd /home/wsx52114/devsim/devsim/wisetop_bjt
+source /home/wsx52114/devsim/denv/bin/activate
+make all
+
+# ── 快速測試 ──
+make quick
+
+# ── 清除重來 ──
+make clean && make all
+```
@@ -1,21 +1,21 @@
 # =============================================================================
 # Makefile for DEVSIM BJT Simulation
 # =============================================================================
-# 本 Makefile 詳細列出每個模擬步驟所需的指令
+# -----------------------------------------------------------------------------
+# 快速操作:
+#   make all        - 執行完整模擬：網格細化 → 零偏壓初始化 → 各項特性分析
+#   make quick      - 快速模擬驗證：跳過網格細化，直接分析
 #
-# 使用方式:
-#   make all        - 完整模擬流程 (refine → init → 所有分析)
-#   make quick      - 快速模擬 (無網格細化)
-#
-# 細部步驟:
-#   make mesh       - 生成初始網格
-#   make refine     - 執行網格細化迴圈 (LOOPS 次)
-#   make init       - 零偏壓初始化 (drift-diffusion)
-#   make gummel     - Gummel 曲線掃描
-#   make output     - 輸出特性掃描
-#   make cb         - 共基極特性掃描
-#   make ac         - AC 小信號分析
-#   make batch      - 批次參數掃描 (多組參數)
+# 個別分析項目:
+#   make mesh       - 生成初始網格 (.msh)
+#   make refine     - 依據電場強度執行網格細化迴圈 (次數由 LOOPS 參數控制)
+#   make init       - 零偏壓初始化 (Drift-Diffusion 模型)
+#   make gummel     - 繪製 Gummel 曲線 (Ic, Ib 對 Vb)
+#   make output     - 輸出特性曲線 (Ic 對 Vc)
+#   make cb         - 共基極特性曲線 (Ic 對 Ve)
+#   make ac         - 小信號交流響應分析 (頻率掃描)
+#   make static     - 靜態偏壓分析（輸出場分佈供 ParaView 檢視）
+#   make batch      - 批次執行所有模擬條件
 # =============================================================================

 # -----------------------------------------------------------------------------
@@ -45,12 +45,18 @@ LOOPS ?= 1
 # -----------------------------------------------------------------------------
 VB_DEFAULT ?= 0.7
 VC_DEFAULT ?= 0.5
+VE_DEFAULT ?= 0.0
 VC_GUMMEL  ?= 2.0

+# static 命令預設環境變數
+VB ?= 0.7
+VC ?= 30.0
+VE ?= 0.0
+
 # -----------------------------------------------------------------------------
 # PHONY 目標宣告
 # -----------------------------------------------------------------------------
-.PHONY: help all quick mesh refine init gummel output cb ac batch clean
+.PHONY: help all quick mesh refine init gummel output cb ac static batch clean

 # =============================================================================
 # HELP - 顯示可用指令
@@ -74,16 +80,17 @@ help:
 	@echo "    make output           輸出特性 (Ic vs Vc)"
 	@echo "    make cb               共基極特性 (Ic vs Ve)"
 	@echo "    make ac               AC 小信號分析"
+	@echo "    make static           靜態偏壓分析 (場分佈輸出)"
 	@echo "    make batch            批次參數掃描"
 	@echo ""
 	@echo "  其他:"
 	@echo "    make clean            清除所有生成檔案"
 	@echo ""
 	@echo "  環境變數:"
-	@echo "    LOOPS=N               細化迭代次數 (預設: 3)"
-	@echo "    NPROC=N               設定 Gmsh 執行緒數 (預設: 自動偵測)"
-	@echo "    VB_DEFAULT=x          output 的 Vb 值 (預設: 0.7)"
-	@echo "    VC_DEFAULT=x          cb/ac 的 Vc 值 (預設: 0.5)"
+	@echo "    LOOPS=N               網格細化迭代次數 (預設: 1，數值越大網格越密)"
+	@echo "    NPROC=N               設定 Gmsh 處理器核心數 (預設: 自動偵測系統規格)"
+	@echo "    VB_DEFAULT=x          設定 output 分析的基極電壓 Vb (預設: 0.7V)"
+	@echo "    VC_DEFAULT=x          設定 cb/ac 分析的集極電壓 Vc (預設: 0.5V)"
 	@echo ""
 	@echo "  目前設定: NPROC=$(NPROC), LOOPS=$(LOOPS)"
 	@echo ""
@@ -93,8 +100,8 @@ help:
 # =============================================================================

 # all: 完整模擬流程
-# 流程: refine → init → gummel → output → cb → ac
-all: refine init gummel output cb ac
+# 流程: refine → init → gummel → output → cb → ac → static
+all: refine init gummel output cb ac static
 	@echo "=========================================="
 	@echo ">>> [ALL] 完整模擬流程完成！"
 	@echo "=========================================="
@@ -171,7 +178,7 @@ refine: mesh
 # 指令: python bjt_dd.py
 # 輸入: bjt_refined.msh (或 bjt.msh)
 # 輸出: bjt_dd_0.msh, bjt_dd_0.tec
-# 說明:
+# 說明: 
 #   1. 載入網格並設定摻雜
 #   2. 求解 Potential Only
 #   3. 設定 Drift-Diffusion 模型
@@ -245,10 +252,10 @@ output:
 	@echo "    參數:"
 	@echo "      - Vb = $(VB_DEFAULT)V (固定)"
 	@echo "      - Vc: 0 → 2V (0.1V 步進)"
-	@echo "    輸出: output_Vb$(VB_DEFAULT).csv"
+	@echo "    輸出: output_Vb$(VB_DEFAULT)0.csv"
 	$(PYTHON) bjt_circuit2.py $(VB_DEFAULT)
 	@echo ""
-	@echo ">>> [output] 完成: output_Vb$(VB_DEFAULT).csv"
+	@echo ">>> [output] 完成: output_Vb$(VB_DEFAULT)0.csv"
 	@echo ""

 # -----------------------------------------------------------------------------
@@ -272,10 +279,10 @@ cb:
 	@echo "    參數:"
 	@echo "      - Vc = $(VC_DEFAULT)V (固定)"
 	@echo "      - Ve: 0 → -1V (-50mV 步進)"
-	@echo "    輸出: cb_Vc$(VC_DEFAULT).csv"
+	@echo "    輸出: cb_Vc$(VC_DEFAULT)0.csv"
 	$(PYTHON) bjt_circuit4.py $(VC_DEFAULT)
 	@echo ""
-	@echo ">>> [cb] 完成: cb_Vc$(VC_DEFAULT).csv"
+	@echo ">>> [cb] 完成: cb_Vc$(VC_DEFAULT)0.csv"
 	@echo ""

 # -----------------------------------------------------------------------------
@@ -307,6 +314,29 @@ ac:
 	@echo ">>> [ac] 完成"
 	@echo ""

+# -----------------------------------------------------------------------------
+# static: 靜態偏壓分析 (輸出場分佈 .tec)
+# 指令: python bjt_circuit6.py <Vb> <Vc> [Ve]
+# 輸入: bjt_dd_0.msh
+# 輸出: static_Vb*_Vc*_Ve*.tec
+# 說明:
+#   - 給定特定偏壓，DC 求解後輸出完整場分佈
+#   - 用 ParaView 檢視 Potential, Electrons, Holes, EField 等
+# -----------------------------------------------------------------------------
+static:
+	@if [ ! -f $(INIT_MSH) ]; then \
+		echo ">>> [static] 未找到初始狀態，先執行 init..."; \
+		$(MAKE) init; \
+	fi
+	@echo ">>> [static] 執行靜態偏壓分析..."
+	@echo "    指令: $(PYTHON) bjt_circuit6.py $(VB) $(VC) $(VE)"
+	@echo "    輸入: $(INIT_MSH)"
+	@echo "    參數: Vb=$(VB)V, Vc=$(VC)V, Ve=$(VE)V"
+	$(PYTHON) bjt_circuit6.py $(VB) $(VC) $(VE)
+	@echo ""
+	@echo ">>> [static] 完成: 用 paraview static_*.tec 檢視結果"
+	@echo ""
+
 # -----------------------------------------------------------------------------
 # batch: 批次參數掃描
 # 指令: bash sims.sh
@@ -341,9 +371,14 @@ batch:
 # CLEAN - 清除所有生成檔案
 # =============================================================================
 clean:
-	@echo ">>> [clean] 清除所有生成檔案..."
+	@echo ">>> [clean] 清除所有生成檔案與暫存..."
 	@echo "    rm -f *.msh *.pos *.tec *.csv *.log"
 	@echo "    rm -rf data"
+	@echo "    rm -rf __pycache__ physics/__pycache__"
+	@echo "    rm -f *.pyc physics/*.pyc"
 	rm -f *.msh *.pos *.tec *.csv *.log
 	rm -rf data
+	rm -rf __pycache__ physics/__pycache__
+	find . -name '*.pyc' -delete 2>/dev/null || true
+	find . -name '__pycache__' -type d -exec rm -rf {} + 2>/dev/null || true
 	@echo ">>> [clean] 完成"
@@ -94,21 +94,37 @@ make all    # 執行完整模擬流程

 | 指令 | 說明 | 輸出檔案 |
 |------|------|----------|
-| `make all` | 執行完整模擬流程 | 所有 CSV 檔案 |
+| `make all` | 執行完整模擬流程 | 所有 CSV 與 TEC 檔案 |
 | `make quick` | 快速模擬 (無網格細化) | — |
-| `make clean` | 清除所有生成檔案 | — |
+| `make clean` | 清除所有生成檔案與暫存檔 (`__pycache__`, `*.pyc`) | — |

 ### 個別步驟

 | 指令 | 說明 | 輸出檔案 |
 |------|------|----------|
 | `make mesh` | 產生初始網格 | `bjt.msh` |
-| `make refine` | 網格細化 (2 次迭代) | `bjt_refined.msh` |
+| `make refine` | 網格細化 (LOOPS 次迭代) | `bjt_refined.msh` |
 | `make init` | 零偏壓初始化 | `bjt_dd_0.msh`, `bjt_dd_0.tec` |
 | `make gummel` | Gummel 曲線 | `gummel.csv` |
 | `make output` | 輸出特性 (Vb=0.7V) | `output_Vb0.70.csv` |
 | `make cb` | 共基極掃描 (Vc=0.5V) | `cb_Vc0.50.csv` |
 | `make ac` | AC 小信號分析 | `ac_Vc0.5_Ve-0.7.csv` |
+| `make static` | 靜態偏壓分析 (場分佈輸出) | `static_Vb*_Vc*_Ve*.tec`, `*.msh` |
+| `make batch` | 批次參數掃描 | `data/*.log` |
+
+### 靜態偏壓分析 (`bjt_circuit6.py`)
+
+給定 Vb、Vc、Ve，DC 求解後輸出場分佈 `.tec` 檔，用 ParaView 觀察電位、電場、載子濃度等（檔案內輸出與註解為全英文）。
+
+```bash
+make static VB=0.7 VC=2.0              # Vb=0.7V, Vc=2.0V, Ve=0V
+make static VB=0.7 VC=2.0 VE=-0.5      # 三端都指定
+make static                            # 使用預設值 (全部 0V)
+```
+
+輸出 `static_Vb*_Vc*_Ve*.tec`，用 `paraview static_*.tec` 開啟。
+
+可觀察變數：`Potential`、`Electrons`/`Holes`、`LogElectrons`/`LogHoles`、`NetDoping`、`LogNetDoping`、`EField`

 ---

@@ -0,0 +1,134 @@
+#!/usr/bin/env python3
+"""
+bjt_circuit6.py - Static DC Bias Analysis (export .tec for ParaView)
+"""
+import sys
+import os
+import math
+from devsim import (
+    load_devices, set_parameter, solve, write_devices,
+    node_model, edge_from_node_model, edge_model
+)
+
+from bjt_common import DOPING_PARAMS, get_contact_currents, restore_doping_params
+from physics.new_physics_c6 import GetContactBiasName, SetSiliconParameters, CreateNodeModel, CreateNodeModelDerivative
+
+
+def _ramp_contact(device, contact, target, step_v=0.1, max_retries=3):
+    """Ramp a single contact from 0 to target voltage in fixed step_v increments."""
+    if target == 0.0:
+        return
+    n_steps = max(1, math.ceil(abs(target) / step_v))
+    print(f"  Ramping {contact}: 0 -> {target}V ({n_steps} steps, ~{abs(target)/n_steps:.3f}V/step)")
+    for i in range(1, n_steps + 1):
+        v = target * i / n_steps
+        set_parameter(device=device, name=GetContactBiasName(contact), value=v)
+        try:
+            solve(type="dc", absolute_error=1e6, relative_error=1.0, maximum_iterations=100)
+        except Exception:
+            # Retry with finer sub-steps
+            ok = target * (i - 1) / n_steps
+            success = False
+            for retry in range(1, max_retries + 1):
+                n_sub = 2 ** retry
+                try:
+                    for j in range(1, n_sub + 1):
+                        sv = ok + (v - ok) * j / n_sub
+                        set_parameter(device=device, name=GetContactBiasName(contact), value=sv)
+                        solve(type="dc", absolute_error=1e6, relative_error=1.0, maximum_iterations=100)
+                    success = True
+                    break
+                except Exception:
+                    set_parameter(device=device, name=GetContactBiasName(contact), value=ok)
+                    try:
+                        solve(type="dc", absolute_error=1e6, relative_error=1.0, maximum_iterations=100)
+                    except Exception:
+                        pass
+            if not success:
+                raise RuntimeError(f"Failed to ramp {contact} to {v:.3f}V")
+    print(f"  {contact} -> {target}V OK")
+
+
+def run_static(Vb=0.7, Vc=2.0, Ve=0.0):
+    """Apply static DC bias and export field distributions."""
+    device = "bjt"
+    region = "Silicon"
+    init_file = "bjt_dd_0.msh"
+
+    if not os.path.exists(init_file):
+        sys.exit("Error: Run 'make init' first")
+
+    print(f"Loading: {init_file}")
+    load_devices(file=init_file)
+
+    # Override IntrinsicElectrons to prevent overflow at high voltage
+    # Clamp exp argument to ±80 (corresponds to 10^35 upper limit)
+    ie_exp = "NIE*exp(ifelse(Potential/V_t > 80, 80, ifelse(Potential/V_t < -80, -80, Potential/V_t)))"
+    CreateNodeModel(device, region, "IntrinsicElectrons", ie_exp)
+    CreateNodeModelDerivative(device, region, "IntrinsicElectrons", ie_exp, 'Potential')
+    CreateNodeModel(device, region, "IntrinsicHoles", "NIE^2/IntrinsicElectrons")
+    CreateNodeModelDerivative(device, region, "IntrinsicHoles", "NIE^2/IntrinsicElectrons", 'Potential')
+
+    restore_doping_params(device, region)
+    SetSiliconParameters(device, region)
+
+    # Initialize contact bias
+    for contact in ["collector", "emitter", "base"]:
+        set_parameter(device=device, name=GetContactBiasName(contact), value=0.0)
+
+    # Ramp bias sequentially: Vc -> Vb -> Ve
+    print(f"Ramping bias: Vc={Vc}V, Vb={Vb}V, Ve={Ve}V")
+    try:
+        _ramp_contact(device, "collector", Vc)
+        _ramp_contact(device, "base", Vb)
+        _ramp_contact(device, "emitter", Ve)
+    except RuntimeError as e:
+        print(f"FAILED: {e}")
+        sys.exit(1)
+    print("DC Solve OK")
+
+    # Setup additional models for visualization
+    node_model(device=device, region=region, name="BuiltinPotential", equation="V_t * asinh(NetDoping / (2 * NIE))")
+    node_model(device=device, region=region, name="AppliedPotential", equation="Potential - BuiltinPotential")
+    edge_from_node_model(device=device, region=region, node_model="Potential")
+    edge_model(device=device, region=region, name="EField", equation="(Potential@n0 - Potential@n1) * EdgeInverseLength")
+    node_model(device=device, region=region, name="LogElectrons", equation="log(Electrons)/log(10)")
+    node_model(device=device, region=region, name="LogHoles", equation="log(Holes)/log(10)")
+    node_model(device=device, region=region, name="NetCarrier", equation="Electrons - Holes")
+
+    try:
+        node_model(device=device, region=region, name="LogNetDoping", equation="asinh(NetDoping/2)/log(10)")
+    except Exception:
+        pass
+
+    # Print terminal currents
+    currents = get_contact_currents(device)
+    print("\n[Terminal Currents]")
+    print(f"Ic = {currents['collector']:+.6e} A/cm")
+    print(f"Ib = {currents['base']:+.6e} A/cm")
+    print(f"Ie = {currents['emitter']:+.6e} A/cm")
+    print(f"KCL: Ic+Ib+Ie = {sum(currents.values()):+.4e} A/cm")
+
+    if abs(currents['collector']) > 1e-30:
+        beta = abs(currents['collector'] / currents['base']) if abs(currents['base']) > 1e-30 else float('inf')
+        print(f"β ≈ {beta:.1f}")
+
+    # Output files
+    tag = f"static_Vb{Vb:.2f}_Vc{Vc:.2f}_Ve{Ve:.2f}"
+    tec_file = f"{tag}.tec"
+
+    write_devices(file=tec_file, type="tecplot")
+
+    print(f"\nDone: {tec_file}")
+
+
+if __name__ == "__main__":
+    if len(sys.argv) < 3:
+        print("Usage: python3 bjt_circuit6.py <Vb> <Vc> [Ve]")
+        sys.exit(0)
+
+    Vb = float(sys.argv[1])
+    Vc = float(sys.argv[2])
+    Ve = float(sys.argv[3]) if len(sys.argv) > 3 else 0.0
+    run_static(Vb=Vb, Vc=Vc, Ve=Ve)
+
@@ -0,0 +1,484 @@
+# Copyright 2013 Devsim LLC
+#
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+#
+#     http://www.apache.org/licenses/LICENSE-2.0
+#
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+
+from devsim import *
+from .model_create import *
+
+def SetUniversalParameters(device, region):
+    universal = {
+        'q' :          1.6e-19,          #, 'coul'),
+      'k' :          1.3806503e-23,    #, 'J/K'),
+      'Permittivity_0' :  8.85e-14     #, 'F/cm^2')
+    }
+    for k, v in universal.items():
+        set_parameter(device=device, region=region, name=k, value=v)
+
+
+
+def SetSiliconParameters(device, region):
+    '''
+      Sets Silicon device parameters on the specified region.
+    '''
+
+    SetUniversalParameters(device, region)
+
+##D. B. M. Klaassen, J. W. Slotboom, and H. C. de Graaff, "Unified apparent bandgap narrowing in n- and p-type Silicon," Solid-State Electronics, vol. 35, no. 2, pp. 125-29, 1992.
+    par = {
+        'Permittivity'     : 11.1*get_parameter(device=device, region=region, name='Permittivity_0'),
+      'NC300'       : 2.8e19,  # '1/cm^3'
+      'NV300'       : 3.1e19, # '1/cm^3'
+      'EG300'       : 1.12,    # 'eV'
+      'EGALPH'      : 2.73e-4, # 'eV/K'
+      'EGBETA'      : 0      , # 'K'
+      'Affinity'    : 4.05      , # 'K'
+      # Canali model
+      'BETAN0'      : 2.57e-2, # '1'
+      'BETANE'      : 0.66,    # '1'
+      'BETAP0'      : 0.46,    # '1'
+      'BETAPE'      : 0.17,    # '1'
+      'VSATN0'      : 1.43e9,
+      'VSATNE'      : -0.87,
+      'VSATP0'      : 1.62e8,
+      'VSATPE'      : -0.52,
+      # Arora model
+      'MUMN'  : 88,
+      'MUMEN' : -0.57,
+      'MU0N'  : 7.4e8,
+      'MU0EN' : -2.33,
+      'NREFN' : 1.26e17,
+      'NREFNE' : 2.4,
+      'ALPHA0N' : 0.88,
+      'ALPHAEN' : -0.146,
+      'MUMP'  : 54.3,
+      'MUMEP' : -0.57,
+      'MU0P'  : 1.36e8,
+      'MU0EP' : -2.23,
+      'NREFP' : 2.35e17,
+      'NREFPE' : 2.4,
+      'ALPHA0P' : 0.88,
+      'ALPHAEP' : -0.146,
+      # SRH
+      "taun" : 1e-5,
+      "taup" : 1e-5,
+      "n1" : 1e10,
+      "p1" : 1e10,
+      # TEMP
+      "T" : 300
+    }
+
+    for k, v in par.items():
+        set_parameter(device=device, region=region, name=k, value=v)
+
+def CreateQuasiFermiLevels(device, region, electron_model, hole_model, variables):
+    '''
+    Creates the models for the quasi-Fermi levels.  Assuming Boltzmann statistics.
+    '''
+    eq = (
+        ('EFN', 'EC + V_t * log(%s/NC)' % electron_model, ('Potential', 'Electrons')),
+      ('EFP', 'EV - V_t * log(%s/NV)' % hole_model, ('Potential', 'Holes')),
+    )
+    for (model, equation, variable_list) in eq:
+        #print "MODEL: " + model + " equation " + equation
+        CreateNodeModel(device, region, model, equation)
+        vset = set(variable_list)
+        for v in variables:
+            if v in vset:
+                CreateNodeModelDerivative(device, region, model, equation, v)
+
+def CreateDensityOfStates(device, region, variables):
+    '''
+      Set up models for density of states.
+      Neglects Bandgap narrowing.
+    '''
+    eq = (
+        ('NC', 'NC300 * (T/300)^1.5', ('T',)),
+      ('NV', 'NV300 * (T/300)^1.5', ('T',)),
+      ('NTOT', 'Donors + Acceptors', ()),
+      # Band Gap Narrowing
+      ('DEG', '0', ()),
+      #('DEG', 'V0.BGN * (log(NTOT/N0.BGN) + ((log(NTOT/N0.BGN)^2 + CON.BGN)^(0.5)))', ()),
+      ('EG', 'EG300 + EGALPH*((300^2)/(300+EGBETA) - (T^2)/(T+EGBETA)) - DEG', ('T')),
+      ('NIE', '((NC * NV)^0.5) * exp(-EG/(2*V_t))*exp(DEG)', ('T')),
+      ('EC', '-Potential - Affinity - DEG/2', ('Potential',)),
+      ('EV', 'EC - EG + DEG/2', ('Potential', 'T')),
+      ('EI', '0.5 * (EC + EV + V_t*log(NC/NV))', ('Potential', 'T')),
+    )
+
+    for (model, equation, variable_list) in eq:
+        #print "MODEL: " + model + " equation " + equation
+        CreateNodeModel(device, region, model, equation)
+        vset = set(variable_list)
+        for v in variables:
+            if v in vset:
+                CreateNodeModelDerivative(device, region, model, equation, v)
+
+
+def GetContactBiasName(contact):
+    return "{0}_bias".format(contact)
+
+def GetContactNodeModelName(contact):
+    return "{0}nodemodel".format(contact)
+
+
+def CreateVT(device, region, variables):
+    '''
+      Calculates the thermal voltage, based on the temperature.
+      V_t : node model
+      V_t_edge : edge model from arithmetic mean
+    '''
+    CreateNodeModel(device, region, 'V_t', "k*T/q")
+    CreateArithmeticMean(device, region, 'V_t', 'V_t_edge')
+    if 'T' in variables:
+        CreateArithmeticMeanDerivative(device, region, 'V_t', 'V_t_edge', 'T')
+
+
+def CreateEField(device, region):
+    '''
+      Creates the EField and DField.
+    '''
+    edge_average_model(device=device, region=region, node_model="Potential",
+                       edge_model="EField", average_type="negative_gradient")
+    edge_average_model(device=device, region=region, node_model="Potential",
+                       edge_model="EField", average_type="negative_gradient", derivative="Potential")
+
+def CreateDField(device, region):
+    CreateEdgeModel(device, region, "DField", "Permittivity * EField")
+    CreateEdgeModel(device, region, "DField:Potential@n0", "Permittivity * EField:Potential@n0")
+    CreateEdgeModel(device, region, "DField:Potential@n1", "Permittivity * EField:Potential@n1")
+
+def CreateSiliconPotentialOnly(device, region):
+    '''
+      Creates the physical models for a Silicon region for equilibrium simulation.
+    '''
+
+    variables = ("Potential",)
+    CreateVT(device, region, variables)
+    CreateDensityOfStates(device, region, variables)
+
+    SetSiliconParameters(device, region)
+
+    # require NetDoping
+    for i in (
+        ("IntrinsicElectrons",       "NIE*exp(ifelse(Potential/V_t > 80, 80, ifelse(Potential/V_t < -80, -80, Potential/V_t)))"),
+         ("IntrinsicHoles",           "NIE^2/IntrinsicElectrons"),
+         ("IntrinsicCharge",          "kahan3(IntrinsicHoles, -IntrinsicElectrons, NetDoping)"),
+         ("PotentialIntrinsicCharge", "-q * IntrinsicCharge")
+    ):
+        n = i[0]
+        e = i[1]
+        CreateNodeModel(device, region, n, e)
+        CreateNodeModelDerivative(device, region, n, e, 'Potential')
+
+    CreateQuasiFermiLevels(device, region, 'IntrinsicElectrons', 'IntrinsicHoles', variables)
+
+    CreateEField(device, region)
+    CreateDField(device, region)
+
+    equation(device=device, region=region, name="PotentialEquation", variable_name="Potential",
+             node_model="PotentialIntrinsicCharge", edge_model="DField", variable_update="log_damp")
+
+def CreateSiliconPotentialOnlyContact(device, region, contact, is_circuit=False):
+    '''
+      Creates the potential equation at the contact
+      if is_circuit is true, than use node given by GetContactBiasName
+    '''
+    if not InNodeModelList(device, region, "contactcharge_node"):
+        CreateNodeModel(device, region, "contactcharge_node", "q*IntrinsicCharge")
+
+    celec_model = "(1e-10 + 0.5*abs(NetDoping+(NetDoping^2 + 4 * NIE^2)^(0.5)))"
+    chole_model = "(1e-10 + 0.5*abs(-NetDoping+(NetDoping^2 + 4 * NIE^2)^(0.5)))"
+    contact_model = "Potential -{0} + ifelse(NetDoping > 0, \
+    -V_t*log({1}/NIE), \
+    V_t*log({2}/NIE))".format(GetContactBiasName(contact), celec_model, chole_model)
+
+    contact_model_name = GetContactNodeModelName(contact)
+    CreateContactNodeModel(device, contact, contact_model_name, contact_model)
+    CreateContactNodeModel(device, contact, "{0}:{1}".format(contact_model_name,"Potential"), "1")
+    if is_circuit:
+        CreateContactNodeModel(device, contact, "{0}:{1}".format(contact_model_name,GetContactBiasName(contact)), "-1")
+
+    if is_circuit:
+        contact_equation(device=device, contact=contact, name="PotentialEquation",
+                         node_model=contact_model_name, edge_model="",
+                         node_charge_model="contactcharge_node", edge_charge_model="DField",
+                         node_current_model="", edge_current_model="", circuit_node=GetContactBiasName(contact))
+    else:
+        contact_equation(device=device, contact=contact, name="PotentialEquation",
+                         node_model=contact_model_name, edge_model="",
+                         node_charge_model="contactcharge_node", edge_charge_model="DField",
+                         node_current_model="", edge_current_model="")
+
+
+def CreateSRH(device, region, variables):
+    '''
+      Shockley Read hall recombination model in terms of generation.
+    '''
+    USRH="(Electrons*Holes - NIE^2)/(taup*(Electrons + n1) + taun*(Holes + p1))"
+    Gn = "-q * USRH"
+    Gp = "+q * USRH"
+    CreateNodeModel(device, region, "USRH", USRH)
+    CreateNodeModel(device, region, "ElectronGeneration", Gn)
+    CreateNodeModel(device, region, "HoleGeneration", Gp)
+    for i in ("Electrons", "Holes", "T"):
+        if i in variables:
+            CreateNodeModelDerivative(device, region, "USRH", USRH, i)
+            CreateNodeModelDerivative(device, region, "ElectronGeneration", Gn, i)
+            CreateNodeModelDerivative(device, region, "HoleGeneration", Gp, i)
+
+def CreateECE(device, region, Jn):
+    '''
+      Electron Continuity Equation using specified equation for Jn
+    '''
+    NCharge = "q * Electrons"
+    CreateNodeModel(device, region, "NCharge", NCharge)
+    CreateNodeModelDerivative(device, region, "NCharge", NCharge, "Electrons")
+
+    equation(device=device, region=region, name="ElectronContinuityEquation", variable_name="Electrons",
+             time_node_model = "NCharge",
+             edge_model=Jn, variable_update="positive", node_model="ElectronGeneration")
+
+def CreateHCE(device, region, Jp):
+    '''
+      Hole Continuity Equation using specified equation for Jp
+    '''
+    PCharge = "-q * Holes"
+    CreateNodeModel(device, region, "PCharge", PCharge)
+    CreateNodeModelDerivative(device, region, "PCharge", PCharge, "Holes")
+
+    equation(device=device, region=region, name="HoleContinuityEquation", variable_name="Holes",
+             time_node_model = "PCharge",
+             edge_model=Jp, variable_update="positive", node_model="HoleGeneration")
+
+def CreatePE(device, region):
+    '''
+      Create Poisson Equation assuming the Electrons and Holes as solution variables
+    '''
+    pne = "-q*kahan3(Holes, -Electrons, NetDoping)"
+    CreateNodeModel(device, region,           "PotentialNodeCharge", pne)
+    CreateNodeModelDerivative(device, region, "PotentialNodeCharge", pne, "Electrons")
+    CreateNodeModelDerivative(device, region, "PotentialNodeCharge", pne, "Holes")
+
+    equation(device=device, region=region, name="PotentialEquation", variable_name="Potential",
+             node_model="PotentialNodeCharge", edge_model="DField",
+             time_node_model="", variable_update="log_damp")
+
+
+def CreateSiliconDriftDiffusion(device, region, mu_n="mu_n", mu_p="mu_p", Jn='Jn', Jp='Jp'):
+    '''
+      Instantiate all equations for drift diffusion simulation
+    '''
+    CreateDensityOfStates(device, region, ("Potential",))
+    CreateQuasiFermiLevels(device, region, "Electrons", "Holes", ("Electrons", "Holes", "Potential"))
+    CreatePE(device, region)
+    CreateSRH(device, region, ("Electrons", "Holes", "Potential"))
+    CreateECE(device, region, Jn)
+    CreateHCE(device, region, Jp)
+
+
+def CreateSiliconDriftDiffusionContact(device, region, contact, Jn, Jp, is_circuit=False): 
+    '''
+      Restrict electrons and holes to their equilibrium values
+      Integrates current into circuit
+    '''
+    CreateSiliconPotentialOnlyContact(device, region, contact, is_circuit)
+
+    celec_model = "(1e-10 + 0.5*abs(NetDoping+(NetDoping^2 + 4 * NIE^2)^(0.5)))"
+    chole_model = "(1e-10 + 0.5*abs(-NetDoping+(NetDoping^2 + 4 * NIE^2)^(0.5)))"
+    contact_electrons_model = "Electrons - ifelse(NetDoping > 0, {0}, NIE^2/{1})".format(celec_model, chole_model)
+    contact_holes_model = "Holes - ifelse(NetDoping < 0, +{1}, +NIE^2/{0})".format(celec_model, chole_model)
+    contact_electrons_name = "{0}nodeelectrons".format(contact)
+    contact_holes_name = "{0}nodeholes".format(contact)
+
+    CreateContactNodeModel(device, contact, contact_electrons_name, contact_electrons_model)
+    CreateContactNodeModel(device, contact, "{0}:{1}".format(contact_electrons_name, "Electrons"), "1")
+
+    CreateContactNodeModel(device, contact, contact_holes_name, contact_holes_model)
+    CreateContactNodeModel(device, contact, "{0}:{1}".format(contact_holes_name, "Holes"), "1")
+
+    if is_circuit:
+        contact_equation(device=device, contact=contact, name="ElectronContinuityEquation",
+                         node_model=contact_electrons_name,
+                         edge_current_model=Jn, circuit_node=GetContactBiasName(contact))
+
+        contact_equation(device=device, contact=contact, name="HoleContinuityEquation",
+                         node_model=contact_holes_name,
+                         edge_current_model=Jp, circuit_node=GetContactBiasName(contact))
+
+    else:
+        contact_equation(device=device, contact=contact, name="ElectronContinuityEquation",
+                         node_model=contact_electrons_name,
+                         edge_current_model=Jn)
+
+        contact_equation(device=device, contact=contact, name="HoleContinuityEquation",
+                         node_model=contact_holes_name,
+                         edge_current_model=Jp)
+
+
+def CreateBernoulliString (Potential="Potential", scaling_variable="V_t", sign=-1):
+    '''
+    Creates the Bernoulli function for Scharfetter Gummel
+    sign -1 for potential
+    sign +1 for energy
+    scaling variable should be V_t
+    Potential should be scaled by V_t in V
+    Ec, Ev should scaled by V_t in eV
+
+    returns the Bernoulli expression and its argument
+    Caller should understand that B(-x) = B(x) + x
+    '''
+
+    tdict = {
+        "Potential" : Potential,
+      "V_t"       : scaling_variable
+    }
+    #### test for requisite models here
+    if sign == -1:
+        vdiff="(%(Potential)s@n0 - %(Potential)s@n1)/%(V_t)s" % tdict
+    elif sign == 1:
+        vdiff="(%(Potential)s@n1 - %(Potential)s@n0)/%(V_t)s" % tdict
+    else:
+        raise NameError("Invalid Sign %s" % sign)
+
+    Bern01 = "B(%s)" % vdiff
+    return (Bern01, vdiff)
+
+
+def CreateElectronCurrent(device, region, mu_n, Potential="Potential", sign=-1, ElectronCurrent="ElectronCurrent", V_t="V_t_edge"):
+    '''
+    Electron current
+    mu_n = mobility name
+    Potential is the driving potential
+    '''
+    EnsureEdgeFromNodeModelExists(device, region, "Potential")
+    EnsureEdgeFromNodeModelExists(device, region, "Electrons")
+    EnsureEdgeFromNodeModelExists(device, region, "Holes")
+    if Potential == "Potential":
+        (Bern01, vdiff) = CreateBernoulliString(scaling_variable=V_t, Potential=Potential, sign=sign)
+    else:
+        raise NameError("Implement proper call")
+
+    tdict = {
+        'Bern01' : Bern01,
+      'vdiff'  : vdiff,
+      'mu_n'   : mu_n,
+      'V_t'    : V_t
+    }
+
+    Jn = "q*%(mu_n)s*EdgeInverseLength*%(V_t)s*kahan3(Electrons@n1*%(Bern01)s,  Electrons@n1*%(vdiff)s,  -Electrons@n0*%(Bern01)s)" % tdict
+
+    CreateEdgeModel(device, region, ElectronCurrent, Jn)
+    for i in ("Electrons", "Potential", "Holes"):
+        CreateEdgeModelDerivatives(device, region, ElectronCurrent, Jn, i)
+
+def CreateHoleCurrent(device, region, mu_p, Potential="Potential", sign=-1, HoleCurrent="HoleCurrent", V_t="V_t_edge"):
+    '''
+    Hole current
+    '''
+    EnsureEdgeFromNodeModelExists(device, region, "Potential")
+    EnsureEdgeFromNodeModelExists(device, region, "Electrons")
+    EnsureEdgeFromNodeModelExists(device, region, "Holes")
+    # Make sure the bernoulli functions exist
+    if Potential == "Potential":
+        (Bern01, vdiff) = CreateBernoulliString(scaling_variable=V_t, Potential=Potential, sign=sign)
+    else:
+        raise NameError("Implement proper call for " + Potential)
+
+    tdict = {
+        'Bern01' : Bern01,
+      'vdiff'  : vdiff,
+      'mu_p'   : mu_p,
+      'V_t'    : V_t
+    }
+
+    Jp ="-q*%(mu_p)s*EdgeInverseLength*%(V_t)s*kahan3(Holes@n1*%(Bern01)s, -Holes@n0*%(Bern01)s, -Holes@n0*%(vdiff)s)" % tdict
+    CreateEdgeModel(device, region, HoleCurrent, Jp)
+    for i in ("Holes", "Potential", "Electrons"):
+        CreateEdgeModelDerivatives(device, region, HoleCurrent, Jp, i)
+
+def CreateAroraMobilityLF(device, region):
+    '''
+      Creates node mobility models and then averages them on edge
+      Uses model from Muller and Kamins
+      Add T derivative dependence later
+    '''
+    models = (
+        ('Tn', 'T/300'),
+      ('mu_arora_n_node',
+            'MUMN * pow(Tn, MUMEN) + (MU0N * pow(T, MU0EN))/(1 + pow((NTOT/(NREFN*pow(Tn, NREFNE))), ALPHA0N*pow(Tn, ALPHAEN)))'),
+      ('mu_arora_p_node',
+            'MUMP * pow(Tn, MUMEP) + (MU0P * pow(T, MU0EP))/(1 + pow((NTOT/(NREFP*pow(Tn, NREFPE))), ALPHA0P*pow(Tn, ALPHAEP)))')
+    )
+
+    for k, v in models:
+        CreateNodeModel(device, region, k, v)
+    CreateArithmeticMean(device, region, 'mu_arora_n_node', 'mu_arora_n_lf')
+    CreateArithmeticMean(device, region, 'mu_arora_p_node', 'mu_arora_p_lf')
+    CreateElectronCurrent(device, region, mu_n = 'mu_arora_n_lf', Potential="Potential", sign=-1, ElectronCurrent="Jn_arora_lf", V_t="V_t_edge")
+    CreateHoleCurrent(device, region, mu_p = 'mu_arora_p_lf', Potential="Potential", sign=-1, HoleCurrent="Jp_arora_lf", V_t="V_t_edge")
+    return {
+        'mu_n' : 'mu_arora_n_lf',
+      'mu_p' : 'mu_arora_p_lf',
+      'Jn'   : 'Jn_arora_lf',
+      'Jp'   : 'Jp_arora_lf',
+    }
+
+
+def CreateHFMobility(device, region, mu_n, mu_p, Jn, Jp):
+    '''
+      Add T derivatives when debugged
+      use parameters to set model flags
+      Caughey Thomas
+    '''
+
+    tdict = {
+        'Jn' : Jn,
+      'mu_n' : mu_n,
+      'Jp' : Jp,
+      'mu_p' : mu_p
+    }
+    tlist = (
+        ("vsat_n", "VSATN0 * pow(T, VSATNE)" % tdict, ('T')),
+      ("beta_n", "BETAN0 * pow(T, BETANE)" % tdict, ('T')),
+      ("Epar_n",
+            "ifelse((%(Jn)s * EField) > 0, abs(EField), 1e-15)" % tdict, ('Potential')),
+      ("mu_n", "%(mu_n)s * pow(1 + pow((%(mu_n)s*Epar_n/vsat_n), beta_n), -1/beta_n)"
+            % tdict, ('Electrons', 'Holes', 'Potential', 'T')),
+      ("vsat_p", "VSATP0 * pow(T, VSATPE)" % tdict, ('T')),
+      ("beta_p", "BETAP0 * pow(T, BETAPE)" % tdict, ('T')),
+      ("Epar_p",
+            "ifelse((%(Jp)s * EField) > 0, abs(EField), 1e-15)" % tdict, ('Potential')),
+      ("mu_p", "%(mu_p)s * pow(1 + pow(%(mu_p)s*Epar_p/vsat_p, beta_p), -1/beta_p)"
+            % tdict, ('Electrons', 'Holes', 'Potential', 'T')),
+    )
+
+    variable_list = ('Electrons', 'Holes', 'Potential')
+    for (model, equation, variables) in tlist:
+        CreateEdgeModel(device, region, model, equation)
+        for v in variable_list:
+            if v in variables:
+                CreateEdgeModelDerivatives(device, region, model, equation, v)
+
+    # This create derivatives automatically
+    CreateElectronCurrent(device, region, mu_n='mu_n', Potential="Potential", sign=-1, ElectronCurrent="Jn", V_t="V_t_edge")
+    CreateHoleCurrent(    device, region, mu_p='mu_p', Potential="Potential", sign=-1, HoleCurrent="Jp", V_t="V_t_edge")
+    return {
+        'mu_n' : 'mu_n',
+      'mu_p' : 'mu_p',
+      'Jn'   : 'Jn',
+      'Jp'   : 'Jp',
+    }
+
+
+
+
@@ -139,13 +139,13 @@ refine2d:
 # sim2d: 執行 2D 模擬 (使用細化後的網格)
 # 指令: python3 opto_simplified_run.py opto_refined.msh
 # 輸入: opto_refined.msh
-# 輸出: opto_simplified_result.msh, opto_simplified_result.tec
+# 輸出: opto_simplified_result.tec
 # -----------------------------------------------------------------------------
 sim2d: $(REFINED_MESH_2D)
 	@echo ">>> [sim2d] 執行 2D 模擬..."
 	@echo "    指令: $(PYTHON) opto_simplified_run.py $(REFINED_MESH_2D) $(RESULT_2D)"
 	$(PYTHON) opto_simplified_run.py $(REFINED_MESH_2D) $(RESULT_2D)
-	@echo ">>> [sim2d] 完成: $(RESULT_2D).msh, $(RESULT_2D).tec"
+	@echo ">>> [sim2d] 完成: $(RESULT_2D).tec"
 	@echo ""

 # -----------------------------------------------------------------------------
@@ -158,7 +158,7 @@ sim2d: $(REFINED_MESH_2D)
 #   - $(MESH_2D)         - 初始網格
 #   - $(BG_POS_2D)       - 背景場
 #   - $(REFINED_MESH_2D) - 細化後網格
-#   - $(RESULT_2D).msh/.tec - 模擬結果
+#   - $(RESULT_2D).tec      - 模擬結果
 # -----------------------------------------------------------------------------
 2d:
 	@echo "=========================================="
@@ -189,7 +189,7 @@ sim2d: $(REFINED_MESH_2D)
 	@echo "    初始網格:   $(MESH_2D)"
 	@echo "    細化網格:   $(REFINED_MESH_2D)"
 	@echo "    背景場:     $(BG_POS_2D)"
-	@echo "    結果檔案:   $(RESULT_2D).msh, $(RESULT_2D).tec"
+	@echo "    結果檔案:   $(RESULT_2D).tec"
 	@echo "=========================================="


@@ -254,13 +254,13 @@ refine3d:
 # sim3d: 執行 3D 模擬 (使用細化後的網格)
 # 指令: python3 opto_3d_run.py opto_3d_refined.msh
 # 輸入: opto_3d_refined.msh
-# 輸出: opto_3d_result.msh, opto_3d_result.tec
+# 輸出: opto_3d_result.tec
 # -----------------------------------------------------------------------------
 sim3d: $(REFINED_MESH_3D)
 	@echo ">>> [sim3d] 執行 3D 模擬..."
 	@echo "    指令: $(PYTHON) opto_3d_run.py $(REFINED_MESH_3D) $(RESULT_3D)"
 	$(PYTHON) opto_3d_run.py $(REFINED_MESH_3D) $(RESULT_3D)
-	@echo ">>> [sim3d] 完成: $(RESULT_3D).msh, $(RESULT_3D).tec"
+	@echo ">>> [sim3d] 完成: $(RESULT_3D).tec"
 	@echo ""

 # -----------------------------------------------------------------------------
@@ -273,7 +273,7 @@ sim3d: $(REFINED_MESH_3D)
 #   - $(MESH_3D)         - 初始網格
 #   - $(BG_POS_3D)       - 背景場
 #   - $(REFINED_MESH_3D) - 細化後網格
-#   - $(RESULT_3D).msh/.tec - 模擬結果
+#   - $(RESULT_3D).tec      - 模擬結果
 # -----------------------------------------------------------------------------
 3d:
 	@echo "=========================================="
@@ -304,7 +304,7 @@ sim3d: $(REFINED_MESH_3D)
 	@echo "    初始網格:   $(MESH_3D)"
 	@echo "    細化網格:   $(REFINED_MESH_3D)"
 	@echo "    背景場:     $(BG_POS_3D)"
-	@echo "    結果檔案:   $(RESULT_3D).msh, $(RESULT_3D).tec"
+	@echo "    結果檔案:   $(RESULT_3D).tec"
 	@echo "=========================================="


@@ -313,6 +313,8 @@ sim3d: $(REFINED_MESH_3D)
 # =============================================================================
 clean:
 	@echo ">>> [clean] 清除所有生成檔案..."
-	@echo "    rm -f *.msh *.tec *.pos *.vtu *.vtm *.vtk current_mesh*.msh"
-	rm -f *.msh *.tec *.pos *.vtu *.vtm *.vtk current_mesh*.msh
+	@echo "    rm -f *.msh *.tec *.pos *.vtu *.vtm *.vtk current_mesh*.msh *.log"
+	rm -f *.msh *.tec *.pos *.vtu *.vtm *.vtk current_mesh*.msh *.log
+	@echo "    rm -rf __pycache__ *.pyc"
+	rm -rf __pycache__ *.pyc
 	@echo ">>> [clean] 完成"
@@ -106,12 +106,19 @@ def generate_mesh(geo_file, mesh_file, dimension=2, force=False):
 # --- Refinement Parameters ---
 MESH_SIZE_MIN_2D = 2.0e-4   # 2 um
 MESH_SIZE_MAX_2D = 5.0e-3   # 50 um
-MESH_SIZE_MIN_3D = 2.0e-4   # 2 um
-MESH_SIZE_MAX_3D = 5.0e-3   # 50 um
+MESH_SIZE_MIN_3D = 10.0e-4   # 10 um
+MESH_SIZE_MAX_3D = 2.0e-3   # 200 um

-# E-field thresholds
-E_VERY_HIGH = 1.0e4   # 10 kV/cm  -> 1x
-E_HIGH = 5.0e3        # 5 kV/cm   -> 2x
-E_MEDIUM = 1.0e3      # 1 kV/cm   -> 4x
-E_LOW = 5.0e2         # 500 V/cm  -> 8x
-E_VERYLOW = 1.0e2     # 100 V/cm  -> 16x
+# E-field thresholds (2D)
+E_VERY_HIGH_2D = 1.0e4   # 10 kV/cm  -> 1x
+E_HIGH_2D = 5.0e3        # 5 kV/cm   -> 2x
+E_MEDIUM_2D = 1.0e3      # 1 kV/cm   -> 4x
+E_LOW_2D = 5.0e2         # 500 V/cm  -> 8x
+E_VERYLOW_2D = 1.0e2     # 100 V/cm  -> 16x
+
+# E-field thresholds (3D)
+E_VERY_HIGH_3D = 2.5e3   # 2.5 kV/cm -> 1x
+E_HIGH_3D = 2.0e3        # 2 kV/cm   -> 2x
+E_MEDIUM_3D = 1.5e3      # 1.5 kV/cm -> 4x
+E_LOW_3D = 1.0e3         # 1 kV/cm   -> 8x
+E_VERYLOW_3D = 5.0e2     # 500 V/cm  -> 16x
@@ -131,17 +131,28 @@ class OptoDevice:
                    pass
                    
    def create_element_models(self):
-        """Create element models for 3D visualization."""
-        if self.dimension == 3:
-            print("--- 3D Element Models ---")
-            for region in self.regions:
-                try:
-                    element_from_edge_model(device=self.name, region=region, edge_model="ElectricField")
-                except:
-                    pass
+        """Create element models for visualization (Emag, etc.)."""
+        from devsim import element_model
+        
+        print(f"--- Creating Element Models ({self.dimension}D) ---")
+        for region in self.regions:
+            try:
+                # Convert edge ElectricField to element model
+                element_from_edge_model(device=self.name, region=region, edge_model="ElectricField")
+                
+                # Calculate Emag (electric field magnitude)
+                if self.dimension == 2:
+                    element_model(device=self.name, region=region, name="Emag",
+                                  equation="(ElectricField_x^2 + ElectricField_y^2)^(0.5)")
+                else:  # 3D
+                    element_model(device=self.name, region=region, name="Emag",
+                                  equation="(ElectricField_x^2 + ElectricField_y^2 + ElectricField_z^2)^(0.5)")
+                
+                print(f"  {region}: Emag created")
+            except Exception as e:
+                print(f"  {region}: Failed - {e}")
                    
    def output(self, prefix):
        """Write output files."""
-        write_devices(file=f"{prefix}.msh", type="devsim")
        write_devices(file=f"{prefix}.tec", type="tecplot")
-        print(f"Done: {prefix}.msh, {prefix}.tec")
+        print(f"Done: {prefix}.tec")
@@ -20,11 +20,11 @@ def emag_refinement(device, region):
    element_model(device=device, region=region, name="Emag",
                  equation="(ElectricField_x^2 + ElectricField_y^2)^(0.5)")
    element_model(device=device, region=region, name="Enorm", equation=f'''
-        ifelse(Emag > {opto_common.E_VERY_HIGH}, 1.0,
-          ifelse(Emag > {opto_common.E_HIGH}, 2.0,
-            ifelse(Emag > {opto_common.E_MEDIUM}, 4.0,
-              ifelse(Emag > {opto_common.E_LOW}, 8.0,
-                if(Emag > {opto_common.E_VERYLOW}, 16)))))
+        ifelse(Emag > {opto_common.E_VERY_HIGH_2D}, 1.0,
+          ifelse(Emag > {opto_common.E_HIGH_2D}, 2.0,
+            ifelse(Emag > {opto_common.E_MEDIUM_2D}, 4.0,
+              ifelse(Emag > {opto_common.E_LOW_2D}, 8.0,
+                if(Emag > {opto_common.E_VERYLOW_2D}, 16)))))
    ''')
    return "Enorm"

@@ -20,11 +20,11 @@ def emag_refinement_3d(device, region):
    element_model(device=device, region=region, name="Emag",
                  equation="(ElectricField_x^2 + ElectricField_y^2 + ElectricField_z^2)^(0.5)")
    element_model(device=device, region=region, name="Enorm", equation=f'''
-        ifelse(Emag > {opto_common.E_VERY_HIGH}, 1.0,
-          ifelse(Emag > {opto_common.E_HIGH}, 2.0,
-            ifelse(Emag > {opto_common.E_MEDIUM}, 4.0,
-              ifelse(Emag > {opto_common.E_LOW}, 8.0,
-                if(Emag > {opto_common.E_VERYLOW}, 16)))))
+        ifelse(Emag > {opto_common.E_VERY_HIGH_3D}, 1.0,
+          ifelse(Emag > {opto_common.E_HIGH_3D}, 2.0,
+            ifelse(Emag > {opto_common.E_MEDIUM_3D}, 4.0,
+              ifelse(Emag > {opto_common.E_LOW_3D}, 8.0,
+                if(Emag > {opto_common.E_VERYLOW_3D}, 16)))))
    ''')
    return "Enorm"

@@ -10,8 +10,7 @@ Mesh.Algorithm = 6;
 /* Scale & Mesh */
 sf = 0.1;                         // 1 mm = 0.1 cm
 um = 1.0e-4;                      // 1 µm = 1e-4 cm
-lc_fine = 5.0 * um;
-lc_base = 50.0 * um;
+lc_coarse = 50.0 * um;

 /* Dimensions */
 W_cu = 3.5*sf;    H_cu = 0.2*sf;      cu_ext = 1.0*sf;
@@ -137,8 +136,7 @@ Physical Curve("interface_encap_led_l") = {ifc_enc_led_l()};
 Physical Curve("interface_encap_led_r") = {ifc_enc_led_r()};

 /* === Mesh Settings === */
-// Initial mesh uses lc_base, refinement can go down to lc_fine
-Mesh.CharacteristicLengthMin = lc_fine;
-Mesh.CharacteristicLengthMax = lc_base;
+Mesh.CharacteristicLengthExtendFromBoundary = 0;
+Mesh.CharacteristicLengthMax = lc_coarse;


@@ -17,8 +17,7 @@ Mesh.Algorithm3D = 10;
 /* Scale & Mesh */
 sf = 0.1;                         // 1 mm = 0.1 cm
 um = 1.0e-4;                      // 1 µm = 1e-4 cm
-lc_fine = 20.0 * um;
-lc_base = 50.0 * um;
+lc_coarse = 50.0 * um;

 /* XY Dimensions (SAME as 2D version) */
 W_cu = 3.5*sf;    H_cu = 0.2*sf;      cu_ext = 1.0*sf;
@@ -214,6 +213,5 @@ ifc_enc_led_z_back() = Surface In BoundingBox{x_led_l-d, y_led_bot-d, -z_led-d,
 Physical Surface("interface_encap_led_z") = {ifc_enc_led_z_front(), ifc_enc_led_z_back()};

 /* === Mesh Settings === */
-// Initial mesh uses lc_base, refinement can go down to lc_fine
-Mesh.CharacteristicLengthMin = lc_fine;
-Mesh.CharacteristicLengthMax = lc_base;
+Mesh.CharacteristicLengthExtendFromBoundary = 0;
+Mesh.CharacteristicLengthMax = lc_coarse;
@@ -77,6 +77,7 @@ def run_simulation():
    # Solve and output
    device.solve()
    device.print_results()
+    device.create_element_models()
    device.output(output_prefix)
Author	SHA1	Message	Date
Dean Huang	be5107731d	feat: Stabilizing High-Voltage BJT Simulation	2026-03-27 16:13:23 +08:00
Dean Huang	90c0da0194	feat: Remove .msh file generation for Paraview	2026-03-19 16:30:19 +08:00
Dean Huang	e0f941ad13	feat: Add PROJECT GUIDE.md.	2026-03-16 14:48:02 +08:00
Dean Huang	0d65a7a0d1	refactor: 更新程式碼以及Makefile	2026-03-16 09:51:13 +08:00
Dean Huang	00afb79e86	docs: 更新註解與README。	2026-03-13 16:01:52 +08:00
Dean Huang	5df4b3432b	feat: 新增靜態直流偏壓分析，更新 README 和 Makefile。	2026-03-13 15:13:51 +08:00
Dean Huang	dd44ab08c4	feat: 新增並更新了多個3D繪圖相關圖檔。	2025-12-29 16:33:20 +08:00
Dean Huang	fe66120592	feat: 新增 P 與 E 的 3D 繪圖	2025-12-29 16:25:41 +08:00
Dean Huang	be38fce34a	Update	2025-12-29 15:07:45 +08:00
Dean Huang	3be76542f6	Update	2025-12-29 12:31:47 +08:00
Dean Huang	d75b08dcc2	Update	2025-12-29 12:23:32 +08:00
Dean Huang	b92e0ae19d	Update	2025-12-29 12:09:32 +08:00
Dean Huang	a664fca506	Update	2025-12-29 11:29:22 +08:00
Dean Huang	eb849b3e39	Update	2025-12-29 10:58:29 +08:00
Dean Huang	b767114c1d	Update	2025-12-29 10:28:21 +08:00
Dean Huang	f1299dadeb	feat: 調整 3D 網格最大尺寸並區分 2D/3D 電場閾值。	2025-12-26 22:02:32 +08:00
Dean Huang	8f12dc4bc1	Merge branch 'refacter'	2025-12-26 16:14:16 +08:00