主题020：静磁场仿真软件工程实践

场景描述

将"能运行的代码"升级为"可维护、可扩展、高质量的工程软件"，掌握科学计算软件的开发规范与最佳实践。

在前19个主题中，我们学习了静磁场仿真的各种数值方法和物理模型。然而，科研代码与工业级软件之间存在巨大鸿沟。本主题将系统介绍如何将零散的研究代码重构为模块化、可测试、易维护的专业仿真软件，涵盖软件架构设计、单元测试、性能优化、配置管理等工程实践。

---

1. 引言：从科研代码到工程软件

1.1 科研代码的常见问题

在科研和学术环境中，仿真代码往往具有以下特点：

1. 脚本化严重：所有功能堆砌在一个文件中，缺乏模块化
2. 硬编码参数：配置与逻辑混杂，修改参数需要改动源代码
3. 缺乏测试：没有单元测试，代码正确性难以保证
4. 文档缺失：仅有少量注释，缺乏系统文档
5. 性能未优化：算法效率低下，未考虑计算复杂度
6. 错误处理不足：异常情况处理不完善，容易崩溃

这些问题在原型开发阶段尚可接受，但当代码需要长期使用、多人协作或对外发布时，就会成为严重的技术债务。

1.2 工程软件的核心特征

与科研代码相比，工程级仿真软件应具备以下特征：

特征	说明	实现方式
模块化	功能分离，高内聚低耦合	类与模块设计
可配置	参数外部化，无需改代码	配置文件管理
可测试	自动化测试覆盖	单元测试框架
可观测	运行状态可监控	日志系统
高性能	算法优化，复杂度可控	性能分析与优化
可维护	代码清晰，文档完善	编码规范与注释

1.3 本主题学习目标

完成本主题学习后，您将能够：

1. 设计合理的软件架构，实现关注点分离
2. 使用面向对象方法组织仿真代码
3. 建立完善的配置管理系统
4. 编写有效的单元测试
5. 进行性能分析与优化
6. 应用Python工程化开发最佳实践

---

2. 软件架构设计

2.1 分层架构模式

科学计算软件通常采用分层架构，将系统划分为若干层次，每层只与相邻层交互：

┌─────────────────────────────────────────┐
│           应用层 (Application)           │
│         主程序、脚本、工作流              │
├─────────────────────────────────────────┤
│           业务逻辑层 (Business)          │
│      求解器、物理模型、算法实现           │
├─────────────────────────────────────────┤
│           数据访问层 (Data)              │
│      文件I/O、配置管理、数据持久化        │
├─────────────────────────────────────────┤
│           基础设施层 (Infrastructure)    │
│      数学库、可视化、日志、工具类         │
└─────────────────────────────────────────┘

这种架构的优势在于：

• 可替换性：每层可以独立替换实现
• 可测试性：每层可以独立测试
• 可维护性：修改一层不影响其他层

2.2 面向对象设计原则

在仿真软件开发中，应遵循SOLID原则：

2.2.1 单一职责原则 (SRP)

每个类应该只有一个引起变化的原因。例如：

# 不好的设计：一个类负责太多事情
class Simulation:
    def solve_electromagnetic(self): ...
    def solve_thermal(self): ...
    def plot_results(self): ...
    def save_to_file(self): ...

# 好的设计：职责分离
class ElectromagneticSolver:
    def solve(self): ...

class ThermalSolver:
    def solve(self): ...

class VisualizationManager:
    def plot(self): ...

class FileManager:
    def save(self): ...

2.2.2 开闭原则 (OCP)

软件实体应该对扩展开放，对修改关闭。通过抽象基类实现：

from abc import ABC, abstractmethod

class FieldSolver(ABC):
    @abstractmethod
    def solve(self) -> np.ndarray:
        pass
    
    @abstractmethod
    def get_field(self, field_type: str) -> np.ndarray:
        pass

class ElectromagneticSolver(FieldSolver):
    def solve(self) -> np.ndarray:
        # 电磁场求解实现
        ...
    
    def get_field(self, field_type: str) -> np.ndarray:
        # 获取电磁场量
        ...

class ThermalSolver(FieldSolver):
    def solve(self) -> np.ndarray:
        # 热场求解实现
        ...
    
    def get_field(self, field_type: str) -> np.ndarray:
        # 获取热场量
        ...

这样，新增求解器类型时无需修改现有代码。

2.2.3 依赖倒置原则 (DIP)

高层模块不应该依赖低层模块，两者都应该依赖抽象。

# 不好的设计：高层依赖具体实现
class SimulationWorkflow:
    def __init__(self):
        self.em_solver = ElectromagneticSolver()  # 依赖具体类
        self.thermal_solver = ThermalSolver()

# 好的设计：依赖抽象
class SimulationWorkflow:
    def __init__(self, em_solver: FieldSolver, thermal_solver: FieldSolver):
        self.em_solver = em_solver  # 依赖抽象接口
        self.thermal_solver = thermal_solver

2.3 设计模式应用

2.3.1 单例模式 - 日志管理器

日志管理器通常只需要一个实例：

class Logger:
    _instance = None
    
    def __new__(cls):
        if cls._instance is None:
            cls._instance = super().__new__(cls)
            cls._instance._setup_logger()
        return cls._instance

2.3.2 工厂模式 - 求解器创建

class SolverFactory:
    @staticmethod
    def create_solver(solver_type: str, config: SimulationConfig) -> FieldSolver:
        if solver_type == 'electromagnetic':
            return ElectromagneticSolver(config)
        elif solver_type == 'thermal':
            return ThermalSolver(config)
        else:
            raise ValueError(f"未知求解器类型: {solver_type}")

2.3.3 装饰器模式 - 性能分析

def timer_decorator(func):
    def wrapper(*args, **kwargs):
        start_time = time.time()
        result = func(*args, **kwargs)
        elapsed = time.time() - start_time
        logger.debug(f"{func.__name__} 执行时间: {elapsed:.4f} 秒")
        return result
    return wrapper

class ElectromagneticSolver:
    @timer_decorator
    def solve(self) -> np.ndarray:
        # 求解逻辑
        ...

---

3. 配置管理系统

3.1 配置外部化的重要性

将配置参数从代码中分离出来有以下好处：

1. 环境适应性：不同环境（开发、测试、生产）使用不同配置
2. 非侵入式修改：修改参数无需改动源代码
3. 版本控制：配置可以单独版本管理
4. 安全性：敏感信息（如API密钥）可以加密存储

3.2 使用dataclass管理配置

Python 3.7+ 引入的dataclass是管理配置的理想选择：

from dataclasses import dataclass, asdict

@dataclass
class SimulationConfig:
    # 网格参数
    nx: int = 101
    ny: int = 101
    Lx: float = 0.1
    Ly: float = 0.1
    
    # 求解器参数
    max_iter: int = 5000
    tol: float = 1e-6
    omega: float = 1.8
    
    # 物理参数
    current_density: float = 1e6
    
    # 计算属性
    @property
    def dx(self) -> float:
        return self.Lx / (self.nx - 1)

3.3 配置序列化

支持多种格式的配置持久化：

import json

class SimulationConfig:
    def to_dict(self) -> Dict:
        return asdict(self)
    
    def to_json(self, filepath: str):
        with open(filepath, 'w') as f:
            json.dump(self.to_dict(), f, indent=2)
    
    @classmethod
    def from_json(cls, filepath: str) -> 'SimulationConfig':
        with open(filepath, 'r') as f:
            return cls(**json.load(f))

3.4 配置验证

使用类型注解和验证器确保配置有效性：

@dataclass
class SimulationConfig:
    nx: int = 101
    
    def __post_init__(self):
        if self.nx < 3:
            raise ValueError("nx必须大于等于3")
        if self.tol <= 0:
            raise ValueError("tol必须为正数")

---

4. 日志系统

4.1 日志的重要性

日志是软件可观测性的基础：

1. 调试：追踪程序执行流程
2. 监控：了解系统运行状态
3. 审计：记录关键操作
4. 故障排查：定位问题原因

4.2 Python logging模块

Python标准库logging提供了强大的日志功能：

import logging

# 创建logger
logger = logging.getLogger('MagneticSimulation')
logger.setLevel(logging.DEBUG)

# 控制台处理器
console_handler = logging.StreamHandler()
console_handler.setLevel(logging.INFO)

# 文件处理器
file_handler = logging.FileHandler('simulation.log')
file_handler.setLevel(logging.DEBUG)

# 格式化
formatter = logging.Formatter(
    '%(asctime)s - %(name)s - %(levelname)s - %(message)s'
)
console_handler.setFormatter(formatter)
file_handler.setFormatter(formatter)

logger.addHandler(console_handler)
logger.addHandler(file_handler)

4.3 日志级别使用规范

级别	使用场景	示例
DEBUG	详细调试信息	迭代误差、中间计算值
INFO	正常操作信息	求解开始/完成、配置加载
WARNING	潜在问题	求解未收敛、参数接近边界
ERROR	操作失败	文件读取失败、数值异常
CRITICAL	严重错误	程序无法继续运行

4.4 结构化日志

对于复杂应用，使用结构化日志便于后续分析：

import json

class StructuredLogFormatter(logging.Formatter):
    def format(self, record):
        log_data = {
            'timestamp': self.formatTime(record),
            'level': record.levelname,
            'message': record.getMessage(),
            'module': record.module,
            'function': record.funcName,
            'line': record.lineno
        }
        return json.dumps(log_data, ensure_ascii=False)

---

5. 单元测试

5.1 测试驱动开发 (TDD)

TDD的基本流程：

1. 编写测试：先写测试用例，定义期望行为
2. 运行测试：确认测试失败
3. 编写代码：实现功能使测试通过
4. 重构：优化代码，保持测试通过

5.2 单元测试框架

虽然Python有unittest和pytest等专业框架，但理解测试原理更重要。我们可以实现简单的测试套件：

class TestSuite:
    def __init__(self):
        self.passed = 0
        self.failed = 0
    
    def assert_equal(self, actual, expected, msg=""):
        if actual == expected:
            self.passed += 1
            return True
        else:
            self.failed += 1
            logger.error(f"测试失败: {msg}")
            return False
    
    def assert_almost_equal(self, actual, expected, tol=1e-6, msg=""):
        if abs(actual - expected) < tol:
            self.passed += 1
            return True
        else:
            self.failed += 1
            logger.error(f"测试失败: {msg}")
            return False

5.3 仿真代码的测试策略

科学计算代码的测试有其特殊性：

5.3.1 物理常数测试

def test_physics_constants(self):
    # 测试真空磁导率
    self.assert_almost_equal(
        PhysicsConstants.MU_0, 
        4*np.pi*1e-7, 
        msg="MU_0"
    )

5.3.2 数值算法测试

def test_solver_convergence(self):
    config = SimulationConfig(nx=21, ny=21, max_iter=1000)
    solver = ElectromagneticSolver(config)
    solver.initialize()
    solver.set_coil(0.05, 0.05, 0.02, 0.02, 1e6)
    solver.solve()
    
    # 验证收敛
    self.assert_true(solver.is_converged, msg="求解器应收敛")
    self.assert_true(solver.iterations < 1000, msg="迭代次数应小于最大值")

5.3.3 边界条件测试

def test_boundary_conditions(self):
    solver = ThermalSolver(config)
    solver.initialize()
    T = solver.solve(T_boundary=300)
    
    # 验证边界值
    self.assert_almost_equal(np.max(T[0, :]), 300, msg="上边界温度")
    self.assert_almost_equal(np.max(T[-1, :]), 300, msg="下边界温度")

5.4 测试覆盖率

理想情况下，测试应覆盖：

• 正常路径：典型输入的预期输出
• 边界条件：最小值、最大值、空值
• 异常情况：错误输入的处理
• 性能边界：大规模数据的处理能力

---

6. 性能分析与优化

6.1 性能分析工具

6.1.1 时间分析装饰器

def timer_decorator(func):
    def wrapper(*args, **kwargs):
        start_time = time.time()
        result = func(*args, **kwargs)
        elapsed = time.time() - start_time
        logger.debug(f"{func.__name__} 执行时间: {elapsed:.4f} 秒")
        return result
    return wrapper

6.1.2 性能分析器类

class PerformanceProfiler:
    def __init__(self):
        self.timings = {}
    
    def record_time(self, name: str, elapsed: float):
        if name not in self.timings:
            self.timings[name] = []
        self.timings[name].append(elapsed)
    
    def get_summary(self) -> Dict:
        summary = {}
        for name, times in self.timings.items():
            summary[name] = {
                'calls': len(times),
                'total': sum(times),
                'mean': np.mean(times),
                'min': min(times),
                'max': max(times)
            }
        return summary

6.2 算法复杂度分析

有限差分法求解泊松方程的时间复杂度：

• 空间复杂度：O(N)，存储网格数据
• 时间复杂度：
  • 单次迭代：O(N)，遍历所有网格点
  • 总求解时间：O(K×N)，K为迭代次数
  • 对于SOR方法，K与网格尺寸相关

6.3 优化策略

6.3.1 向量化计算

使用NumPy数组运算代替Python循环：

# 低效：Python循环
for i in range(ny):
    for j in range(nx):
        B[i, j] = np.sqrt(Bx[i, j]**2 + By[i, j]**2)

# 高效：向量化
B = np.sqrt(Bx**2 + By**2)

6.3.2 稀疏矩阵技术

对于大规模问题，使用稀疏矩阵存储：

from scipy.sparse import csr_matrix
from scipy.sparse.linalg import spsolve

# 构建稀疏矩阵
A_sparse = csr_matrix((data, (row_ind, col_ind)), shape=(N, N))

# 求解
x = spsolve(A_sparse, b)

6.3.3 并行计算

使用多线程或多进程加速：

from multiprocessing import Pool

def solve_subdomain(args):
    subdomain_id, config = args
    # 求解子域
    ...

with Pool(processes=4) as pool:
    results = pool.map(solve_subdomain, subdomain_tasks)

6.4 性能基准测试

建立性能基准，监控性能退化：

def benchmark_solver(grid_sizes=[21, 41, 61, 81]):
    results = []
    for n in grid_sizes:
        config = SimulationConfig(nx=n, ny=n)
        workflow = SimulationWorkflow(config)
        
        start = time.time()
        workflow.run_full_analysis()
        elapsed = time.time() - start
        
        results.append({
            'grid_size': n,
            'cells': n * n,
            'time': elapsed
        })
    
    return results

---

7. 完整代码实现

以下是主题020的完整代码实现，展示了软件工程最佳实践：

"""
主题020: 静磁场仿真软件工程实践
=====================================
模块化设计、性能优化与工程化开发演示

本代码演示：
1. 模块化架构设计
2. 面向对象编程最佳实践
3. 性能分析与优化
4. 单元测试框架
5. 配置驱动开发
6. 日志与错误处理
"""

import numpy as np
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
from matplotlib.patches import Rectangle, FancyBboxPatch
import time
import json
import logging
from abc import ABC, abstractmethod
from dataclasses import dataclass, asdict
from typing import Dict, List, Tuple, Optional, Union
from enum import Enum
import warnings
warnings.filterwarnings('ignore')

# 设置中文字体
plt.rcParams['font.sans-serif'] = ['SimHei', 'DejaVu Sans']
plt.rcParams['axes.unicode_minus'] = False

# ============================================================
# 第一部分：配置管理模块
# ============================================================

class PhysicsConstants:
    """物理常数类 - 集中管理所有物理常量"""
    MU_0 = 4 * np.pi * 1e-7  # 真空磁导率 [H/m]
    EPSILON_0 = 8.854e-12    # 真空介电常数 [F/m]
    C = 2.998e8              # 光速 [m/s]
    
    # 材料属性库
    MATERIALS = {
        'copper': {
            'sigma': 5.8e7,      # 电导率 [S/m]
            'mu_r': 0.999991,    # 相对磁导率
            'rho': 8960,         # 密度 [kg/m³]
            'k': 400,            # 热导率 [W/(m·K)]
            'cp': 385,           # 比热容 [J/(kg·K)]
            'E': 110e9,          # 杨氏模量 [Pa]
            'nu': 0.34,          # 泊松比
            'alpha': 17e-6       # 热膨胀系数 [/K]
        },
        'iron': {
            'sigma': 1.0e7,
            'mu_r': 1000,
            'rho': 7874,
            'k': 80,
            'cp': 450,
            'E': 200e9,
            'nu': 0.30,
            'alpha': 12e-6
        },
        'air': {
            'sigma': 0,
            'mu_r': 1.0,
            'rho': 1.225,
            'k': 0.026,
            'cp': 1005,
            'E': 0,
            'nu': 0,
            'alpha': 0
        }
    }
    
    @classmethod
    def get_material(cls, name: str) -> Dict:
        """获取材料属性"""
        if name not in cls.MATERIALS:
            raise ValueError(f"未知材料: {name}。可用材料: {list(cls.MATERIALS.keys())}")
        return cls.MATERIALS[name].copy()


@dataclass
class SimulationConfig:
    """仿真配置数据类"""
    # 网格参数
    nx: int = 101
    ny: int = 101
    Lx: float = 0.1  # [m]
    Ly: float = 0.1  # [m]
    
    # 求解器参数
    max_iter: int = 5000
    tol: float = 1e-6
    omega: float = 1.8  # SOR松弛因子
    
    # 线圈参数
    coil_center_x: float = 0.05
    coil_center_y: float = 0.05
    coil_width: float = 0.02
    coil_height: float = 0.03
    current_density: float = 1e6  # [A/m²]
    
    # 铁芯参数
    core_center_x: float = 0.065
    core_center_y: float = 0.05
    core_width: float = 0.015
    core_height: float = 0.04
    core_mu_r: float = 1000
    
    # 输出参数
    output_dir: str = "."
    dpi: int = 150
    
    def to_dict(self) -> Dict:
        """转换为字典"""
        return asdict(self)
    
    @classmethod
    def from_dict(cls, config_dict: Dict) -> 'SimulationConfig':
        """从字典创建配置"""
        return cls(**config_dict)
    
    @classmethod
    def from_json(cls, filepath: str) -> 'SimulationConfig':
        """从JSON文件加载配置"""
        with open(filepath, 'r', encoding='utf-8') as f:
            return cls.from_dict(json.load(f))
    
    def to_json(self, filepath: str):
        """保存配置到JSON文件"""
        with open(filepath, 'w', encoding='utf-8') as f:
            json.dump(self.to_dict(), f, indent=2, ensure_ascii=False)
    
    @property
    def dx(self) -> float:
        """x方向网格步长"""
        return self.Lx / (self.nx - 1)
    
    @property
    def dy(self) -> float:
        """y方向网格步长"""
        return self.Ly / (self.ny - 1)


# ============================================================
# 第二部分：日志系统
# ============================================================

class Logger:
    """日志管理器"""
    
    _instance = None
    
    def __new__(cls):
        if cls._instance is None:
            cls._instance = super().__new__(cls)
            cls._instance._setup_logger()
        return cls._instance
    
    def _setup_logger(self):
        """设置日志配置"""
        self.logger = logging.getLogger('MagneticSimulation')
        self.logger.setLevel(logging.DEBUG)
        
        # 控制台处理器
        console_handler = logging.StreamHandler()
        console_handler.setLevel(logging.INFO)
        console_format = logging.Formatter(
            '%(asctime)s - %(levelname)s - %(message)s',
            datefmt='%H:%M:%S'
        )
        console_handler.setFormatter(console_format)
        
        # 文件处理器
        file_handler = logging.FileHandler('simulation.log', encoding='utf-8')
        file_handler.setLevel(logging.DEBUG)
        file_format = logging.Formatter(
            '%(asctime)s - %(name)s - %(levelname)s - %(funcName)s:%(lineno)d - %(message)s'
        )
        file_handler.setFormatter(file_format)
        
        self.logger.addHandler(console_handler)
        self.logger.addHandler(file_handler)
    
    def debug(self, msg: str):
        self.logger.debug(msg)
    
    def info(self, msg: str):
        self.logger.info(msg)
    
    def warning(self, msg: str):
        self.logger.warning(msg)
    
    def error(self, msg: str):
        self.logger.error(msg)
    
    def critical(self, msg: str):
        self.logger.critical(msg)


# 全局日志实例
logger = Logger()


# ============================================================
# 第三部分：性能分析装饰器
# ============================================================

def timer_decorator(func):
    """函数执行时间计时装饰器"""
    def wrapper(*args, **kwargs):
        start_time = time.time()
        result = func(*args, **kwargs)
        end_time = time.time()
        elapsed = end_time - start_time
        logger.debug(f"{func.__name__} 执行时间: {elapsed:.4f} 秒")
        return result
    return wrapper


def profile_decorator(func):
    """性能分析装饰器 - 记录调用次数和总时间"""
    call_count = 0
    total_time = 0.0
    
    def wrapper(*args, **kwargs):
        nonlocal call_count, total_time
        start_time = time.time()
        result = func(*args, **kwargs)
        elapsed = time.time() - start_time
        call_count += 1
        total_time += elapsed
        logger.debug(f"{func.__name__} 调用 #{call_count}: {elapsed:.4f}s, 累计: {total_time:.4f}s")
        return result
    
    wrapper.get_stats = lambda: {'calls': call_count, 'total_time': total_time}
    return wrapper


# ============================================================
# 第四部分：抽象基类与接口定义
# ============================================================

class FieldSolver(ABC):
    """场求解器抽象基类"""
    
    def __init__(self, config: SimulationConfig):
        self.config = config
        self.dx = config.dx
        self.dy = config.dy
        self.nx = config.nx
        self.ny = config.ny
        self._solution = None
        self._converged = False
        self._iterations = 0
    
    @abstractmethod
    def initialize(self):
        """初始化求解器"""
        pass
    
    @abstractmethod
    def solve(self) -> np.ndarray:
        """执行求解"""
        pass
    
    @abstractmethod
    def get_field(self, field_type: str) -> np.ndarray:
        """获取场量"""
        pass
    
    @property
    def is_converged(self) -> bool:
        return self._converged
    
    @property
    def iterations(self) -> int:
        return self._iterations


class PostProcessor(ABC):
    """后处理器抽象基类"""
    
    @abstractmethod
    def process(self, data: Dict) -> Dict:
        """处理数据"""
        pass


# ============================================================
# 第五部分：具体求解器实现
# ============================================================

class ElectromagneticSolver(FieldSolver):
    """电磁场求解器 - 实现FieldSolver接口"""
    
    def __init__(self, config: SimulationConfig):
        super().__init__(config)
        self.mu0 = PhysicsConstants.MU_0
        self.Az = np.zeros((self.ny, self.nx))
        self.Jz = np.zeros((self.ny, self.nx))
        self.mu = np.ones((self.ny, self.nx)) * self.mu0
        self._Bx = None
        self._By = None
        self._B_mag = None
        
    def initialize(self):
        """初始化场量"""
        self.Az.fill(0)
        self.Jz.fill(0)
        self.mu.fill(self.mu0)
        logger.info("电磁场求解器初始化完成")
    
    def set_coil(self, center_x: float, center_y: float, 
                 width: float, height: float, current_density: float):
        """设置线圈"""
        ix_start = int((center_x - width/2) / self.dx)
        ix_end = int((center_x + width/2) / self.dx)
        iy_start = int((center_y - height/2) / self.dy)
        iy_end = int((center_y + height/2) / self.dy)
        
        # 边界检查
        ix_start = max(1, min(ix_start, self.nx-2))
        ix_end = max(1, min(ix_end, self.nx-2))
        iy_start = max(1, min(iy_start, self.ny-2))
        iy_end = max(1, min(iy_end, self.ny-2))
        
        self.Jz[iy_start:iy_end, ix_start:ix_end] = current_density
        logger.info(f"线圈设置完成: 中心=({center_x*1000:.1f}, {center_y*1000:.1f})mm, "
                   f"尺寸={width*1000:.1f}x{height*1000:.1f}mm, J={current_density/1e6:.2f}MA/m²")
    
    def set_core(self, center_x: float, center_y: float,
                 width: float, height: float, mu_r: float):
        """设置铁芯"""
        ix_start = int((center_x - width/2) / self.dx)
        ix_end = int((center_x + width/2) / self.dx)
        iy_start = int((center_y - height/2) / self.dy)
        iy_end = int((center_y + height/2) / self.dy)
        
        ix_start = max(1, min(ix_start, self.nx-2))
        ix_end = max(1, min(ix_end, self.nx-2))
        iy_start = max(1, min(iy_start, self.ny-2))
        iy_end = max(1, min(iy_end, self.ny-2))
        
        self.mu[iy_start:iy_end, ix_start:ix_end] = self.mu0 * mu_r
        logger.info(f"铁芯设置完成: 中心=({center_x*1000:.1f}, {center_y*1000:.1f})mm, "
                   f"μr={mu_r}")
    
    @timer_decorator
    def solve(self) -> np.ndarray:
        """求解泊松方程"""
        omega = self.config.omega
        tol = self.config.tol
        max_iter = self.config.max_iter
        
        logger.info(f"开始求解电磁场: max_iter={max_iter}, tol={tol}")
        
        for iteration in range(max_iter):
            Az_old = self.Az.copy()
            
            # SOR迭代
            for i in range(1, self.ny-1):
                for j in range(1, self.nx-1):
                    mu_avg = self.mu[i, j]
                    rhs = -mu_avg * self.Jz[i, j]
                    
                    self.Az[i, j] = (1 - omega) * self.Az[i, j] + \
                                   omega / (2/self.dx**2 + 2/self.dy**2) * \
                                   ((self.Az[i+1, j] + self.Az[i-1, j])/self.dx**2 + \
                                    (self.Az[i, j+1] + self.Az[i, j-1])/self.dy**2 - rhs)
            
            # 收敛检查
            if iteration % 100 == 0:
                error = np.max(np.abs(self.Az - Az_old))
                if error < tol:
                    self._converged = True
                    self._iterations = iteration
                    logger.info(f"电磁场求解收敛于第 {iteration} 次迭代, 误差={error:.2e}")
                    break
        
        if not self._converged:
            logger.warning(f"电磁场求解未收敛，达到最大迭代次数 {max_iter}")
        
        self._solution = self.Az.copy()
        self._compute_magnetic_field()
        
        return self.Az
    
    def _compute_magnetic_field(self):
        """计算磁感应强度"""
        self._Bx = np.zeros_like(self.Az)
        self._By = np.zeros_like(self.Az)
        
        self._Bx[1:-1, :] = (self.Az[2:, :] - self.Az[:-2, :]) / (2 * self.dy)
        self._By[:, 1:-1] = -(self.Az[:, 2:] - self.Az[:, :-2]) / (2 * self.dx)
        
        self._B_mag = np.sqrt(self._Bx**2 + self._By**2)
    
    def get_field(self, field_type: str) -> np.ndarray:
        """获取场量"""
        if field_type == 'Az':
            return self.Az
        elif field_type == 'Bx':
            return self._Bx if self._Bx is not None else np.zeros_like(self.Az)
        elif field_type == 'By':
            return self._By if self._By is not None else np.zeros_like(self.Az)
        elif field_type == 'B_mag':
            return self._B_mag if self._B_mag is not None else np.zeros_like(self.Az)
        elif field_type == 'Jz':
            return self.Jz
        else:
            raise ValueError(f"未知场类型: {field_type}")
    
    def compute_joule_heating(self, material_name: str = 'copper') -> np.ndarray:
        """计算焦耳热"""
        material = PhysicsConstants.get_material(material_name)
        sigma = material['sigma']
        Q = self.Jz**2 / sigma
        return Q
    
    def compute_lorentz_force(self) -> Tuple[np.ndarray, np.ndarray]:
        """计算洛伦兹力"""
        if self._Bx is None or self._By is None:
            self._compute_magnetic_field()
        fx = self.Jz * self._By
        fy = -self.Jz * self._Bx
        return fx, fy


class ThermalSolver(FieldSolver):
    """热场求解器"""
    
    def __init__(self, config: SimulationConfig):
        super().__init__(config)
        self.T = np.ones((self.ny, self.nx)) * 300
        self.Q = np.zeros((self.ny, self.nx))
        self.k = np.ones((self.ny, self.nx)) * 400
        self.rho = np.ones((self.ny, self.nx)) * 8960
        self.cp = np.ones((self.ny, self.nx)) * 385
    
    def initialize(self):
        """初始化"""
        self.T.fill(300)
        self.Q.fill(0)
        logger.info("热场求解器初始化完成")
    
    def set_heat_source(self, Q: np.ndarray):
        """设置热源"""
        self.Q = Q.copy()
    
    def set_material_region(self, x_range: Tuple[float, float], 
                           y_range: Tuple[float, float],
                           material_name: str):
        """设置区域材料"""
        material = PhysicsConstants.get_material(material_name)
        
        ix_start = int(x_range[0] / self.dx)
        ix_end = int(x_range[1] / self.dx)
        iy_start = int(y_range[0] / self.dy)
        iy_end = int(y_range[1] / self.dy)
        
        ix_start = max(0, min(ix_start, self.nx-1))
        ix_end = max(0, min(ix_end, self.nx-1))
        iy_start = max(0, min(iy_start, self.ny-1))
        iy_end = max(0, min(iy_end, self.ny-1))
        
        self.k[iy_start:iy_end, ix_start:ix_end] = material['k']
        self.rho[iy_start:iy_end, ix_start:ix_end] = material['rho']
        self.cp[iy_start:iy_end, ix_start:ix_end] = material['cp']
    
    @timer_decorator
    def solve(self, T_boundary: float = 300) -> np.ndarray:
        """求解稳态热传导"""
        omega = 1.5
        tol = self.config.tol
        max_iter = self.config.max_iter
        
        logger.info("开始求解温度场")
        
        # 边界条件
        self.T[0, :] = T_boundary
        self.T[-1, :] = T_boundary
        self.T[:, 0] = T_boundary
        self.T[:, -1] = T_boundary
        
        for iteration in range(max_iter):
            T_old = self.T.copy()
            
            for i in range(1, self.ny-1):
                for j in range(1, self.nx-1):
                    k_p = self.k[i, j]
                    
                    k_e = (self.k[i, j+1] + k_p) / 2
                    k_w = (self.k[i, j-1] + k_p) / 2
                    k_n = (self.k[i+1, j] + k_p) / 2
                    k_s = (self.k[i-1, j] + k_p) / 2
                    
                    a_e = k_e / self.dx**2
                    a_w = k_w / self.dx**2
                    a_n = k_n / self.dy**2
                    a_s = k_s / self.dy**2
                    a_p = a_e + a_w + a_n + a_s
                    
                    rhs = self.Q[i, j]
                    
                    self.T[i, j] = (1 - omega) * self.T[i, j] + \
                                  omega / a_p * (a_e * self.T[i, j+1] + \
                                                a_w * self.T[i, j-1] + \
                                                a_n * self.T[i+1, j] + \
                                                a_s * self.T[i-1, j] + rhs)
            
            if iteration % 100 == 0:
                error = np.max(np.abs(self.T - T_old))
                if error < tol:
                    self._converged = True
                    self._iterations = iteration
                    logger.info(f"温度场求解收敛于第 {iteration} 次迭代")
                    break
        
        self._solution = self.T.copy()
        return self.T
    
    def get_field(self, field_type: str) -> np.ndarray:
        if field_type == 'T':
            return self.T
        elif field_type == 'Q':
            return self.Q
        else:
            raise ValueError(f"未知场类型: {field_type}")


# ============================================================
# 第六部分：仿真工作流管理
# ============================================================

class SimulationWorkflow:
    """仿真工作流管理器"""
    
    def __init__(self, config: SimulationConfig):
        self.config = config
        self.em_solver = ElectromagneticSolver(config)
        self.thermal_solver = ThermalSolver(config)
        self.results = {}
        logger.info("仿真工作流初始化完成")
    
    def setup_geometry(self):
        """设置几何模型"""
        cfg = self.config
        
        # 设置线圈
        self.em_solver.set_coil(
            cfg.coil_center_x, cfg.coil_center_y,
            cfg.coil_width, cfg.coil_height,
            cfg.current_density
        )
        
        # 设置铁芯
        self.em_solver.set_core(
            cfg.core_center_x, cfg.core_center_y,
            cfg.core_width, cfg.core_height,
            cfg.core_mu_r
        )
        
        logger.info("几何模型设置完成")
    
    def run_electromagnetic_analysis(self) -> Dict:
        """执行电磁分析"""
        logger.info("=" * 50)
        logger.info("开始电磁场分析")
        
        self.em_solver.initialize()
        Az = self.em_solver.solve()
        
        B_mag = self.em_solver.get_field('B_mag')
        Q_joule = self.em_solver.compute_joule_heating()
        fx, fy = self.em_solver.compute_lorentz_force()
        
        results = {
            'Az': Az,
            'B_mag': B_mag,
            'B_max': np.max(B_mag),
            'Q_joule': Q_joule,
            'P_joule': np.sum(Q_joule) * self.config.dx * self.config.dy,
            'fx': fx,
            'fy': fy,
            'iterations': self.em_solver.iterations,
            'converged': self.em_solver.is_converged
        }
        
        logger.info(f"电磁场分析完成: B_max={results['B_max']:.4e} T")
        logger.info(f"焦耳热功率: {results['P_joule']:.4f} W")
        
        return results
    
    def run_thermal_analysis(self, Q_joule: np.ndarray) -> Dict:
        """执行热分析"""
        logger.info("=" * 50)
        logger.info("开始热场分析")
        
        self.thermal_solver.initialize()
        self.thermal_solver.set_heat_source(Q_joule)
        
        # 设置材料区域
        cfg = self.config
        self.thermal_solver.set_material_region(
            (cfg.coil_center_x - cfg.coil_width/2, cfg.coil_center_x + cfg.coil_width/2),
            (cfg.coil_center_y - cfg.coil_height/2, cfg.coil_center_y + cfg.coil_height/2),
            'copper'
        )
        
        T = self.thermal_solver.solve()
        
        results = {
            'T': T,
            'T_max': np.max(T),
            'T_avg': np.mean(T),
            'iterations': self.thermal_solver.iterations,
            'converged': self.thermal_solver.is_converged
        }
        
        logger.info(f"热场分析完成: T_max={results['T_max']:.2f} K")
        
        return results
    
    def run_full_analysis(self) -> Dict:
        """执行完整分析"""
        logger.info("=" * 60)
        logger.info("开始完整仿真分析")
        
        start_time = time.time()
        
        # 设置几何
        self.setup_geometry()
        
        # 电磁分析
        em_results = self.run_electromagnetic_analysis()
        
        # 热分析
        thermal_results = self.run_thermal_analysis(em_results['Q_joule'])
        
        # 汇总结果
        self.results = {
            'electromagnetic': em_results,
            'thermal': thermal_results,
            'config': self.config.to_dict(),
            'runtime': time.time() - start_time
        }
        
        logger.info("=" * 60)
        logger.info(f"完整分析完成，总耗时: {self.results['runtime']:.2f} 秒")
        
        return self.results


# ============================================================
# 第七部分：可视化模块
# ============================================================

class VisualizationManager:
    """可视化管理器"""
    
    def __init__(self, config: SimulationConfig):
        self.config = config
        self.output_dir = config.output_dir
        self.dpi = config.dpi
    
    def create_field_plot(self, field: np.ndarray, title: str, 
                         cmap: str = 'viridis', 
                         colorbar_label: str = '',
                         filename: str = None):
        """创建单场图"""
        fig, ax = plt.subplots(figsize=(10, 8))
        
        x = np.linspace(0, self.config.Lx, self.config.nx) * 1000
        y = np.linspace(0, self.config.Ly, self.config.ny) * 1000
        X, Y = np.meshgrid(x, y)
        
        im = ax.contourf(X, Y, field, levels=50, cmap=cmap)
        ax.set_xlabel('x (mm)')
        ax.set_ylabel('y (mm)')
        ax.set_title(title)
        
        cbar = plt.colorbar(im, ax=ax)
        if colorbar_label:
            cbar.set_label(colorbar_label)
        
        if filename:
            filepath = f"{self.output_dir}/{filename}"
            plt.savefig(filepath, dpi=self.dpi, bbox_inches='tight')
            logger.info(f"图像已保存: {filepath}")
        
        plt.close()
    
    def create_multifield_plot(self, results: Dict, filename: str = None):
        """创建多物理场对比图"""
        fig, axes = plt.subplots(2, 2, figsize=(14, 12))
        
        x = np.linspace(0, self.config.Lx, self.config.nx) * 1000
        y = np.linspace(0, self.config.Ly, self.config.ny) * 1000
        X, Y = np.meshgrid(x, y)
        
        em = results['electromagnetic']
        thermal = results['thermal']
        
        # 磁场
        ax = axes[0, 0]
        im = ax.contourf(X, Y, em['B_mag']*1000, levels=50, cmap='viridis')
        ax.set_title('磁感应强度 (mT)')
        ax.set_xlabel('x (mm)')
        ax.set_ylabel('y (mm)')
        plt.colorbar(im, ax=ax)
        
        # 焦耳热
        ax = axes[0, 1]
        im = ax.contourf(X, Y, em['Q_joule']/1e6, levels=30, cmap='hot')
        ax.set_title('焦耳热功率密度 (MW/m³)')
        ax.set_xlabel('x (mm)')
        plt.colorbar(im, ax=ax)
        
        # 温度
        ax = axes[1, 0]
        im = ax.contourf(X, Y, thermal['T']-273.15, levels=50, cmap='jet')
        ax.set_title('温度 (°C)')
        ax.set_xlabel('x (mm)')
        ax.set_ylabel('y (mm)')
        plt.colorbar(im, ax=ax)
        
        # 汇总信息
        ax = axes[1, 1]
        ax.axis('off')
        
        summary = f"""
仿真结果汇总:

电磁场:
  B_max = {em['B_max']*1000:.4f} mT
  P_joule = {em['P_joule']:.4f} W
  迭代次数 = {em['iterations']}

热场:
  T_max = {thermal['T_max']-273.15:.2f} °C
  T_avg = {thermal['T_avg']-273.15:.2f} °C
  迭代次数 = {thermal['iterations']}

运行时间: {results['runtime']:.2f} s
"""
        ax.text(0.1, 0.5, summary, transform=ax.transAxes, fontsize=11,
               verticalalignment='center', fontfamily='monospace',
               bbox=dict(boxstyle='round', facecolor='wheat', alpha=0.5))
        
        plt.tight_layout()
        
        if filename:
            filepath = f"{self.output_dir}/{filename}"
            plt.savefig(filepath, dpi=self.dpi, bbox_inches='tight')
            logger.info(f"多物理场图已保存: {filepath}")
        
        plt.close()


# ============================================================
# 第八部分：单元测试
# ============================================================

class TestSuite:
    """单元测试套件"""
    
    def __init__(self):
        self.passed = 0
        self.failed = 0
        self.tests = []
    
    def assert_equal(self, actual, expected, msg=""):
        """断言相等"""
        try:
            if actual == expected:
                self.passed += 1
                return True
            else:
                self.failed += 1
                logger.error(f"测试失败: {msg}, 期望={expected}, 实际={actual}")
                return False
        except Exception as e:
            self.failed += 1
            logger.error(f"测试异常: {msg}, 错误={e}")
            return False
    
    def assert_almost_equal(self, actual, expected, tol=1e-6, msg=""):
        """断言近似相等"""
        try:
            if abs(actual - expected) < tol:
                self.passed += 1
                return True
            else:
                self.failed += 1
                logger.error(f"测试失败: {msg}, 期望={expected}, 实际={actual}, 差值={abs(actual-expected)}")
                return False
        except Exception as e:
            self.failed += 1
            logger.error(f"测试异常: {msg}, 错误={e}")
            return False
    
    def assert_true(self, condition, msg=""):
        """断言为真"""
        if condition:
            self.passed += 1
            return True
        else:
            self.failed += 1
            logger.error(f"测试失败: {msg}")
            return False
    
    def run_tests(self):
        """运行所有测试"""
        logger.info("=" * 60)
        logger.info("开始单元测试")
        
        # 测试物理常数
        self.test_physics_constants()
        
        # 测试配置类
        self.test_simulation_config()
        
        # 测试求解器
        self.test_electromagnetic_solver()
        
        logger.info("=" * 60)
        logger.info(f"测试完成: 通过={self.passed}, 失败={self.failed}")
        
        return self.failed == 0
    
    def test_physics_constants(self):
        """测试物理常数"""
        logger.info("测试物理常数...")
        
        # 测试真空磁导率
        self.assert_almost_equal(PhysicsConstants.MU_0, 4*np.pi*1e-7, msg="MU_0")
        
        # 测试材料库
        copper = PhysicsConstants.get_material('copper')
        self.assert_true('sigma' in copper, msg="铜材料应有sigma属性")
        self.assert_true(copper['sigma'] > 0, msg="铜电导率应大于0")
    
    def test_simulation_config(self):
        """测试配置类"""
        logger.info("测试配置类...")
        
        config = SimulationConfig(nx=51, ny=51)
        
        # 测试网格步长计算
        self.assert_almost_equal(config.dx, config.Lx/(config.nx-1), msg="dx计算")
        self.assert_almost_equal(config.dy, config.Ly/(config.ny-1), msg="dy计算")
        
        # 测试字典转换
        config_dict = config.to_dict()
        self.assert_true('nx' in config_dict, msg="配置字典应有nx")
        
        # 测试从字典恢复
        config2 = SimulationConfig.from_dict(config_dict)
        self.assert_equal(config2.nx, config.nx, msg="从字典恢复nx")
    
    def test_electromagnetic_solver(self):
        """测试电磁场求解器"""
        logger.info("测试电磁场求解器...")
        
        config = SimulationConfig(nx=21, ny=21, max_iter=100)
        solver = ElectromagneticSolver(config)
        
        # 测试初始化
        solver.initialize()
        self.assert_equal(solver.Az.shape, (config.ny, config.nx), msg="Az形状")
        
        # 测试线圈设置
        solver.set_coil(0.05, 0.05, 0.02, 0.02, 1e6)
        J_max = np.max(solver.Jz)
        self.assert_true(J_max > 0, msg="线圈电流密度应大于0")


# ============================================================
# 第九部分：性能分析器
# ============================================================

class PerformanceProfiler:
    """性能分析器"""
    
    def __init__(self):
        self.timings = {}
        self.memory_usage = {}
    
    def record_time(self, name: str, elapsed: float):
        """记录时间"""
        if name not in self.timings:
            self.timings[name] = []
        self.timings[name].append(elapsed)
    
    def get_summary(self) -> Dict:
        """获取性能摘要"""
        summary = {}
        for name, times in self.timings.items():
            summary[name] = {
                'calls': len(times),
                'total': sum(times),
                'mean': np.mean(times),
                'min': min(times),
                'max': max(times)
            }
        return summary
    
    def print_report(self):
        """打印性能报告"""
        logger.info("=" * 60)
        logger.info("性能分析报告")
        logger.info("=" * 60)
        
        summary = self.get_summary()
        for name, stats in summary.items():
            logger.info(f"{name}:")
            logger.info(f"  调用次数: {stats['calls']}")
            logger.info(f"  总时间: {stats['total']:.4f} s")
            logger.info(f"  平均时间: {stats['mean']:.4f} s")
            logger.info(f"  最短时间: {stats['min']:.4f} s")
            logger.info(f"  最长时间: {stats['max']:.4f} s")


# ============================================================
# 第十部分：主程序
# ============================================================

def main():
    """主函数"""
    print("=" * 70)
    print("主题020: 静磁场仿真软件工程实践")
    print("=" * 70)
    
    # 1. 运行单元测试
    print("\n[1/5] 运行单元测试...")
    test_suite = TestSuite()
    tests_passed = test_suite.run_tests()
    
    if not tests_passed:
        print("警告: 部分单元测试未通过，继续执行...")
    
    # 2. 创建配置
    print("\n[2/5] 创建仿真配置...")
    config = SimulationConfig(
        nx=101,
        ny=101,
        Lx=0.1,
        Ly=0.1,
        max_iter=3000,
        tol=1e-6,
        omega=1.8,
        current_density=1e6,
        output_dir=".",
        dpi=150
    )
    
    # 保存配置
    config.to_json('simulation_config.json')
    print("  配置已保存到 simulation_config.json")
    
    # 3. 执行仿真
    print("\n[3/5] 执行仿真分析...")
    workflow = SimulationWorkflow(config)
    results = workflow.run_full_analysis()
    
    # 4. 生成可视化
    print("\n[4/5] 生成可视化...")
    viz = VisualizationManager(config)
    
    # 单场图
    viz.create_field_plot(
        results['electromagnetic']['B_mag'] * 1000,
        '磁感应强度分布',
        cmap='viridis',
        colorbar_label='B (mT)',
        filename='fig1_magnetic_field.png'
    )
    
    viz.create_field_plot(
        results['electromagnetic']['Q_joule'] / 1e6,
        '焦耳热功率密度',
        cmap='hot',
        colorbar_label='Q (MW/m³)',
        filename='fig2_joule_heating.png'
    )
    
    viz.create_field_plot(
        results['thermal']['T'] - 273.15,
        '温度分布',
        cmap='jet',
        colorbar_label='T (°C)',
        filename='fig3_temperature.png'
    )
    
    # 多物理场对比图
    viz.create_multifield_plot(results, filename='fig4_multiphysics_summary.png')
    
    # 5. 性能分析
    print("\n[5/5] 性能分析...")
    
    # 对比不同网格大小的性能
    grid_sizes = [21, 41, 61, 81]
    performance_data = []
    
    for n in grid_sizes:
        test_config = SimulationConfig(nx=n, ny=n, max_iter=500, tol=1e-5)
        test_workflow = SimulationWorkflow(test_config)
        
        start = time.time()
        test_workflow.setup_geometry()
        test_workflow.em_solver.solve()
        elapsed = time.time() - start
        
        performance_data.append({
            'grid_size': n,
            'total_cells': n * n,
            'time': elapsed
        })
        print(f"  网格 {n}x{n}: {elapsed:.4f} 秒")
    
    # 绘制性能曲线
    fig, ax = plt.subplots(figsize=(10, 6))
    
    grid_sizes_arr = [d['grid_size'] for d in performance_data]
    times_arr = [d['time'] for d in performance_data]
    cells_arr = [d['total_cells'] for d in performance_data]
    
    ax.plot(cells_arr, times_arr, 'bo-', linewidth=2, markersize=8)
    ax.set_xlabel('网格单元数')
    ax.set_ylabel('求解时间 (s)')
    ax.set_title('求解器性能分析')
    ax.grid(True, alpha=0.3)
    
    # 添加复杂度参考线 O(N^2)
    if len(times_arr) > 1:
        ref_time = times_arr[-1]
        ref_cells = cells_arr[-1]
        ref_line = [ref_time * (c / ref_cells)**2 for c in cells_arr]
        ax.plot(cells_arr, ref_line, 'r--', label='O(N²) 参考', alpha=0.7)
        ax.legend()
    
    plt.tight_layout()
    plt.savefig('fig5_performance_analysis.png', dpi=150, bbox_inches='tight')
    print("  性能分析图已保存: fig5_performance_analysis.png")
    plt.close()
    
    # 软件架构图
    fig, ax = plt.subplots(figsize=(14, 10))
    ax.set_xlim(0, 10)
    ax.set_ylim(0, 10)
    ax.axis('off')
    
    ax.text(5, 9.5, '静磁场仿真软件架构', fontsize=16, ha='center', fontweight='bold')
    
    # 配置层
    rect1 = FancyBboxPatch((0.5, 8), 9, 0.8, boxstyle="round,pad=0.1", 
                           facecolor='lightblue', edgecolor='blue', linewidth=2)
    ax.add_patch(rect1)
    ax.text(5, 8.4, '配置管理层 (SimulationConfig)', ha='center', fontsize=11)
    
    # 核心求解器层
    rect2 = FancyBboxPatch((0.5, 6.5), 4, 1, boxstyle="round,pad=0.1",
                           facecolor='lightgreen', edgecolor='green', linewidth=2)
    ax.add_patch(rect2)
    ax.text(2.5, 7.2, '电磁场求解器', ha='center', fontsize=10)
    ax.text(2.5, 6.8, 'ElectromagneticSolver', ha='center', fontsize=9)
    
    rect3 = FancyBboxPatch((5.5, 6.5), 4, 1, boxstyle="round,pad=0.1",
                           facecolor='lightyellow', edgecolor='orange', linewidth=2)
    ax.add_patch(rect3)
    ax.text(7.5, 7.2, '热场求解器', ha='center', fontsize=10)
    ax.text(7.5, 6.8, 'ThermalSolver', ha='center', fontsize=9)
    
    # 工作流层
    rect4 = FancyBboxPatch((0.5, 5), 9, 1, boxstyle="round,pad=0.1",
                           facecolor='lightcoral', edgecolor='red', linewidth=2)
    ax.add_patch(rect4)
    ax.text(5, 5.7, '仿真工作流管理器 (SimulationWorkflow)', ha='center', fontsize=11)
    ax.text(5, 5.3, '协调各求解器执行', ha='center', fontsize=9)
    
    # 工具层
    rect5 = FancyBboxPatch((0.5, 3.5), 2.8, 1, boxstyle="round,pad=0.1",
                           facecolor='plum', edgecolor='purple', linewidth=2)
    ax.add_patch(rect5)
    ax.text(1.9, 4.2, '可视化', ha='center', fontsize=10)
    ax.text(1.9, 3.8, 'Visualization', ha='center', fontsize=9)
    
    rect6 = FancyBboxPatch((3.6, 3.5), 2.8, 1, boxstyle="round,pad=0.1",
                           facecolor='plum', edgecolor='purple', linewidth=2)
    ax.add_patch(rect6)
    ax.text(5, 4.2, '日志系统', ha='center', fontsize=10)
    ax.text(5, 3.8, 'Logger', ha='center', fontsize=9)
    
    rect7 = FancyBboxPatch((6.7, 3.5), 2.8, 1, boxstyle="round,pad=0.1",
                           facecolor='plum', edgecolor='purple', linewidth=2)
    ax.add_patch(rect7)
    ax.text(8.1, 4.2, '性能分析', ha='center', fontsize=10)
    ax.text(8.1, 3.8, 'Profiler', ha='center', fontsize=9)
    
    # 基础层
    rect8 = FancyBboxPatch((0.5, 2), 4, 1, boxstyle="round,pad=0.1",
                           facecolor='lightgray', edgecolor='gray', linewidth=2)
    ax.add_patch(rect8)
    ax.text(2.5, 2.7, '物理常数', ha='center', fontsize=10)
    ax.text(2.5, 2.3, 'PhysicsConstants', ha='center', fontsize=9)
    
    rect9 = FancyBboxPatch((5.5, 2), 4, 1, boxstyle="round,pad=0.1",
                           facecolor='lightgray', edgecolor='gray', linewidth=2)
    ax.add_patch(rect9)
    ax.text(7.5, 2.7, '单元测试', ha='center', fontsize=10)
    ax.text(7.5, 2.3, 'TestSuite', ha='center', fontsize=9)
    
    # 箭头
    ax.annotate('', xy=(5, 8), xytext=(5, 7.5),
                arrowprops=dict(arrowstyle='->', color='black', lw=2))
    ax.annotate('', xy=(5, 6.5), xytext=(5, 6),
                arrowprops=dict(arrowstyle='->', color='black', lw=2))
    ax.annotate('', xy=(5, 5), xytext=(5, 4.5),
                arrowprops=dict(arrowstyle='->', color='black', lw=2))
    
    # 设计原则
    principles = [
        '设计原则:',
        '• 单一职责原则 (SRP)',
        '• 开闭原则 (OCP)',
        '• 依赖倒置原则 (DIP)',
        '• 接口隔离原则 (ISP)'
    ]
    ax.text(0.5, 1, '\n'.join(principles), fontsize=10, 
           verticalalignment='top', fontfamily='monospace',
           bbox=dict(boxstyle='round', facecolor='wheat', alpha=0.5))
    
    plt.savefig('fig6_software_architecture.png', dpi=150, bbox_inches='tight')
    print("  软件架构图已保存: fig6_software_architecture.png")
    plt.close()
    
    # 最终总结
    print("\n" + "=" * 70)
    print("主题020仿真完成！")
    print("=" * 70)
    print("\n生成的文件:")
    print("  - simulation_config.json: 仿真配置文件")
    print("  - simulation.log: 运行日志")
    print("  - fig1_magnetic_field.png: 磁场分布")
    print("  - fig2_joule_heating.png: 焦耳热分布")
    print("  - fig3_temperature.png: 温度分布")
    print("  - fig4_multiphysics_summary.png: 多物理场汇总")
    print("  - fig5_performance_analysis.png: 性能分析")
    print("  - fig6_software_architecture.png: 软件架构")
    print("=" * 70)


if __name__ == "__main__":
    main()

---

8. 代码深度解析

8.1 物理常数类设计

PhysicsConstants类采用类变量和类方法设计，集中管理所有物理常量和材料属性：

class PhysicsConstants:
    MU_0 = 4 * np.pi * 1e-7  # 类变量
    
    MATERIALS = {
        'copper': {'sigma': 5.8e7, ...},
        'iron': {'sigma': 1.0e7, ...}
    }
    
    @classmethod
    def get_material(cls, name: str) -> Dict:
        if name not in cls.MATERIALS:
            raise ValueError(f"未知材料: {name}")
        return cls.MATERIALS[name].copy()  # 返回副本防止修改

设计要点：

• 使用类变量确保全局唯一
• 返回副本防止外部修改
• 异常处理提供清晰的错误信息

8.2 配置类设计

SimulationConfig使用Python 3.7+的@dataclass装饰器：

@dataclass
class SimulationConfig:
    nx: int = 101
    ny: int = 101
    
    @property
    def dx(self) -> float:
        return self.Lx / (self.nx - 1)

优势：

• 自动生成__init__、__repr__等方法
• 类型注解提供IDE支持
• 计算属性动态派生相关参数

8.3 抽象基类设计

FieldSolver定义了求解器的通用接口：

class FieldSolver(ABC):
    @abstractmethod
    def solve(self) -> np.ndarray:
        pass
    
    @property
    def is_converged(self) -> bool:
        return self._converged

设计意图：

• 强制子类实现必要方法
• 统一接口便于多态使用
• 公共逻辑在基类中实现

8.4 日志单例模式

Logger类确保全局只有一个日志实例：

class Logger:
    _instance = None
    
    def __new__(cls):
        if cls._instance is None:
            cls._instance = super().__new__(cls)
            cls._instance._setup_logger()
        return cls._instance

实现细节：

• 重写__new__控制实例创建
• 延迟初始化减少启动开销
• 双处理器分别输出到控制台和文件

8.5 性能装饰器

timer_decorator为函数添加计时功能：

def timer_decorator(func):
    def wrapper(*args, **kwargs):
        start_time = time.time()
        result = func(*args, **kwargs)
        elapsed = time.time() - start_time
        logger.debug(f"{func.__name__} 执行时间: {elapsed:.4f} 秒")
        return result
    return wrapper

使用方式：

class ElectromagneticSolver:
    @timer_decorator
    def solve(self) -> np.ndarray:
        ...

---

9. 运行结果预期

执行代码后将生成以下输出和文件：

9.1 控制台输出

======================================================================
主题020: 静磁场仿真软件工程实践
======================================================================

[1/5] 运行单元测试...
08:36:54 - INFO - 开始单元测试
08:36:54 - INFO - 测试完成: 通过=9, 失败=0

[2/5] 创建仿真配置...
  配置已保存到 simulation_config.json

[3/5] 执行仿真分析...
08:36:54 - INFO - 开始完整仿真分析
08:36:54 - INFO - 电磁场求解收敛于第 0 次迭代, 误差=0.00e+00
08:36:54 - INFO - 完整分析完成，总耗时: 0.04 秒

[4/5] 生成可视化...
[5/5] 性能分析...
  网格 21x21: 0.0578 秒
  网格 41x41: 0.2641 秒
  网格 61x61: 1.3721 秒
  网格 81x81: 3.5902 秒

9.2 生成的文件

文件	说明
simulation_config.json	JSON格式的仿真配置
simulation.log	详细运行日志
fig1_magnetic_field.png	磁感应强度分布云图
fig2_joule_heating.png	焦耳热功率密度分布
fig3_temperature.png	温度场分布
fig4_multiphysics_summary.png	多物理场汇总对比图
fig5_performance_analysis.png	求解器性能曲线
fig6_software_architecture.png	软件架构示意图

9.3 性能分析结果

从性能测试可以看出：

• 21×21网格：约0.06秒
• 41×41网格：约0.26秒（4.5倍）
• 61×61网格：约1.37秒（23.7倍）
• 81×81网格：约3.59秒（62.1倍）

时间复杂度接近O(N²)，与理论预期一致。

---

10. 进阶挑战

挑战1：扩展材料库

添加更多材料到PhysicsConstants.MATERIALS：

'aluminum': {
    'sigma': 3.5e7,
    'mu_r': 1.00002,
    'rho': 2700,
    'k': 237,
    ...
}

挑战2：实现配置文件热加载

添加文件监控，当配置文件修改时自动重新加载：

import os
import time

class ConfigWatcher:
    def __init__(self, filepath):
        self.filepath = filepath
        self.last_mtime = os.path.getmtime(filepath)
    
    def check_reload(self):
        current_mtime = os.path.getmtime(self.filepath)
        if current_mtime != self.last_mtime:
            self.last_mtime = current_mtime
            return SimulationConfig.from_json(self.filepath)
        return None

挑战3：添加命令行参数支持

使用argparse支持命令行配置：

import argparse

def parse_args():
    parser = argparse.ArgumentParser(description='静磁场仿真')
    parser.add_argument('--nx', type=int, default=101, help='x方向网格数')
    parser.add_argument('--current', type=float, default=1e6, help='电流密度')
    parser.add_argument('--config', type=str, help='配置文件路径')
    return parser.parse_args()

挑战4：并行求解器

使用multiprocessing实现多网格并行求解：

from multiprocessing import Pool

def solve_grid_size(n):
    config = SimulationConfig(nx=n, ny=n)
    workflow = SimulationWorkflow(config)
    return workflow.run_full_analysis()

with Pool(4) as pool:
    results = pool.map(solve_grid_size, [41, 61, 81, 101])

---

11. 总结与展望

11.1 核心知识点回顾

本主题系统介绍了科学计算软件工程化的关键实践：

1. 架构设计：分层架构、SOLID原则、设计模式
2. 配置管理：dataclass、序列化、外部化配置
3. 日志系统：分级日志、多处理器、结构化日志
4. 单元测试：测试策略、断言方法、覆盖率
5. 性能优化：复杂度分析、向量化、并行计算

11.2 最佳实践清单

开发科学计算软件时，建议遵循以下清单：

• 使用类型注解增强代码可读性
• 配置与代码分离，支持外部配置
• 建立完善的日志系统
• 编写单元测试覆盖核心功能
• 使用抽象基类定义接口
• 性能关键路径添加计时装饰器
• 异常处理提供清晰的错误信息
• 代码注释说明"为什么"而非"做什么"

11.3 未来扩展方向

1. Web界面：使用Flask/Django开发Web前端
2. 数据库集成：使用SQLAlchemy持久化仿真结果
3. 分布式计算：使用MPI或Ray实现集群并行
4. 机器学习集成：使用PyTorch/TensorFlow加速求解
5. CI/CD：建立自动化测试和部署流程

---

附录：常见报错与解决方案

错误1：ModuleNotFoundError

ModuleNotFoundError: No module named 'matplotlib'

解决方案：

pip install matplotlib numpy

错误2：日志文件权限错误

PermissionError: [Errno 13] Permission denied: 'simulation.log'

解决方案：

• 检查文件是否被其他程序占用
• 修改输出目录到可写位置
• 使用管理员权限运行

错误3：内存不足

MemoryError: Unable to allocate array

解决方案：

• 减小网格尺寸
• 使用稀疏矩阵存储
• 分块处理大数据

错误4：收敛失败

WARNING: 电磁场求解未收敛

解决方案：

• 增加最大迭代次数
• 调整SOR松弛因子
• 检查边界条件设置

---

主题020完成！ 您已掌握了科学计算软件工程化的核心技能，可以开始构建专业的仿真软件系统了。