PINNs案例——二维磁场计算

news2025/7/22 11:37:14

基于物理信息的神经网络是一种解决偏微分方程计算问题的全新方法…

有关PINN基础详见：PINNs案例——中心热源温度场预测问题的torch代码

今日分享代码案例：二维带电流源磁场计算

该案例参考学习论文：[1]张宇娇，孙宏达，赵志涛，徐斌，黄雄峰．基于物理信息神经网络的电磁场计算方法[J/OL]．电工技术学报.

结果前瞻

磁场矢量

问题场景

  • ** 求解域** 案例选取边长为 10 的正方形区域作为研究对象，即定义域为
    Ω = [ 0 , 10 ] × [ 0 , 10 ] \Omega = [0, 10] \times [0, 10] Ω=[0,10]×[0,10]。以正方形左下角为原点，正方形区域内包含一个带电流密度为1的导体，导体区域为 Ω = [ 4 , 2 ] × [ 6 , 8 ] \Omega = [4, 2] \times [6, 8] Ω=[4,2]×[6,8]。
  • **参数设定：**导体处 J = 1 J = 1 J=1，导体材料为铜（copper），其余区域为空气（air），磁导率 γ \gamma γ 均为1
  • **边界条件：**方形四周边界的矢量磁位为0（理想导体边界）
  • **求解量：**区域内二维稳态磁场分布(矢量)

区域采样点

PINNs理论求解

在本案例中，我们引入矢量磁位 $ \mathbf{A} $ 来辅助磁场的计算。矢量磁位A设定为：
B = ∇ × A \mathbf{B} = \nabla \times \mathbf{A} B=∇×A
其中， B \mathbf{B} B 表示磁场强度，该关系表明磁场是矢量磁位的旋度。

在该正方形区域内，矢量磁位 $ A_z(x, y) $ 满足麦克斯韦方程，引入矢量磁位A后在二维稳态下，其数学表达式为：

( ∂ 2 A z ∂ x 2 + ∂ 2 A z ∂ y 2 ) = − μ J z \left( \frac{\partial^2 A_z}{\partial x^2} + \frac{\partial^2 A_z}{\partial y^2} \right) =-\mu J_z (∂x2∂2Az​​+∂y2∂2Az​​)=−μJz​

其中：
μ \mu μ 为材料的磁导率，表示材料对磁场的响应能力；
J z J_z Jz​ 为电流密度在 z z z 方向的分量，表示单位面积内电流的大小；
A z ( x , y ) A_z(x, y) Az​(x,y) 为矢量磁位在 z z z 方向的分量，通常用于描述二维平面内的磁场分布。

这种方程形式常见于静电学或静磁学中的泊松方程，用于描述在存在源（如电荷或电流）的情况下，势函数（如矢量磁位）的空间分布。矢量磁位的引入使得磁场的计算更为简便，尤其是在处理对称性问题或边界条件时。

本次案例中，由于存在铜导体区域（有源）和空气域（无源）,故采用PINNs案例——多介质分区热传导中相同的分区思想。此处不过多介绍

损失函数可以表示为：

  1. PDE 残差项： 控制方程如下：
    Loss PDE = 1 N c ∑ i = 1 N c ∣ f ( x c i , y c i ) ∣ 2 \text{Loss}_{\text{PDE}} = \frac{1}{N_c} \sum_{i=1}^{N_c} \left| f\left( x_c^i, y_c^i \right) \right|^2 LossPDE​=Nc​1​i=1∑Nc​​ ​f(xci​,yci​) ​2
    f ( x c i , y c i ) = ( ∂ 2 A z ∂ x 2 + ∂ 2 A z ∂ y 2 ) + μ J z f\left( x_c^i, y_c^i \right) = \left( \frac{\partial^2 A_z}{\partial x^2} + \frac{\partial^2 A_z}{\partial y^2} \right) +\mu J_z f(xci​,yci​)=(∂x2∂2Az​​+∂y2∂2Az​​)+μJz​

注意： 由于PDE损失中需加入 μ 0 \mu_{0} μ0​，而 μ 0 = 4 π ∗ 10 − 7 \mu_{0} = 4\pi * 10^{-7} μ0​=4π∗10−7 数量级过小，会对神经网络的训练带来较大困难。故需对PDE方程进行简单的无量纲化操作，令 A ′ = A / μ 0 A' = A/\mu_0 A′=A/μ0​, 于是可以得到 f ( x c i , y c i ) = ( ∂ 2 A z ′ ∂ x 2 + ∂ 2 A z ′ ∂ y 2 ) + μ r J z f\left( x_c^i, y_c^i \right) = \left( \frac{\partial^2 A_z'}{\partial x^2} + \frac{\partial^2 A_z'}{\partial y^2} \right) +\mu_r J_z f(xci​,yci​)=(∂x2∂2Az′​​+∂y2∂2Az′​​)+μr​Jz​.预测后再将A乘 μ 0 \mu_0 μ0​以还原A

  1. 边界条件项：在多介质问题中，边界条件除矢量磁位A=0的方形边界条件，还要添加不同介质交界面的连续性边界条件。

  • 固定温度边界：
    Loss BC = 1 N b ∑ i = 1 N b ∣ A ( x b i , y b i ) ∣ 2 \text{Loss}_{\text{BC}} = \frac{1}{N_b} \sum_{i=1}^{N_b} \left|A\left( x_b^i, y_b^i \right) \right|^2 LossBC​=Nb​1​i=1∑Nb​​ ​A(xbi​,ybi​) ​2

  • 介质交界面：
    此案例中以A的第一类和第二类边界条件作为介质交界面的连续性条件

第一类边界条件：
Loss I1 = 1 N I ∑ j = 1 N I ∣ A p ( x I j ) − A p + 1 ( x I j ) ∣ 2 \text{Loss}_{\text{I1}} = \frac{1}{N_{I}}\sum_{j=1}^{N_{I}}\left| A_{p}\left(x_{I}^{j}\right)-A_{p+1}\left(x_{I}^{j}\right)\right|^{2} LossI1​=NI​1​j=1∑NI​​ ​Ap​(xIj​)−Ap+1​(xIj​) ​2
第二类边界条件（不同介质在边界面上矢量磁位法向导数与磁导率的乘积连续）：
Loss I2 = 1 N I ∑ j = 1 N I ∣ γ p ∂ A p ( x I j ) ∂ n − γ p + 1 ∂ A p + 1 ( x I j ) ∂ n ∣ 2 \text{Loss}_{\text{I2}} = \frac{1}{N_{I}}\sum_{j=1}^{N_{I}}\left|\gamma_{p}\frac{\partial A_{p}\left(x_{I}^{j}\right)}{\partial n}-\gamma_{p+1}\frac{\partial A_{p+1}\left(x_{I}^{j}\right)}{\partial n}\right|^{2} LossI2​=NI​1​j=1∑NI​​ ​γp​∂n∂Ap​(xIj​)​−γp+1​∂n∂Ap+1​(xIj​)​ ​2
其中：

  • N I N_I NI​ 表示介质边界点的点数
  • A p A_p Ap​ 表示第P个介质的矢量磁位
  • γ p \gamma_ {p} γp​ 表示第P个介质的热导率

训练

采用10000次adam优化与lbfgs优化器进行优化。各项损失收敛过程如下，收敛效果良好：

损失曲线

复现结果

矢量磁位

磁场标量

磁场矢量

comsol仿真结果

矢量磁位

磁场

可见结果吻合良好

pytorch代码

为避免繁杂的点生成，作者编写了一组点生成函数，使用时单独建立points_generator.py，并在训练代码中引入。

points_generator.py

import numpy as np
from pyDOE import lhs

def generate_boundary_points(x_min, y_min, x_max, y_max, n):
    # 计算每条边上的点数
    num_points_x = int((x_max - x_min) * n) + 1
    num_points_y = int((y_max - y_min) * n) + 1

    # 生成每条边上的点
    bottom_edge = np.array([np.linspace(x_min, x_max, num_points_x), np.full(num_points_x, y_min)]).T
    top_edge = np.array([np.linspace(x_min, x_max, num_points_x), np.full(num_points_x, y_max)]).T
    left_edge = np.array([np.full(num_points_y, x_min), np.linspace(y_min, y_max, num_points_y)]).T
    right_edge = np.array([np.full(num_points_y, x_max), np.linspace(y_min, y_max, num_points_y)]).T

    # 合并所有边上的点，去掉重复的角点
    boundary_points = np.concatenate([bottom_edge, top_edge[1:], left_edge[1:-1], right_edge[1:-1]])
    boundary_vertical = np.concatenate([left_edge[1:-1], right_edge[1:-1]])
    boundary_horizontal = np.concatenate([bottom_edge, top_edge[1:]])

    return boundary_points,boundary_vertical,boundary_horizontal

def generate_internal_points(x_min, y_min, x_max, y_max, n):
    # 计算内部点的网格
    x = np.linspace(x_min, x_max, int((x_max - x_min) * n) + 1)
    y = np.linspace(y_min, y_max, int((y_max - y_min) * n) + 1)
    x_grid, y_grid = np.meshgrid(x, y)

    # 去掉边界点
    x_grid = x_grid[1:-1, 1:-1]
    y_grid = y_grid[1:-1, 1:-1]

    # 将网格点转换为一维数组
    internal_points = np.array([x_grid.flatten(), y_grid.flatten()]).T

    return internal_points

def generate_mask(points, x_min, y_min, x_max, y_max):
    # 检查每个点是否在矩形区域内
    mask = (points[:, 0] >= x_min) & (points[:, 0] <= x_max) & (points[:, 1] >= y_min) & (points[:, 1] <= y_max)

    return mask


def lhs_internal_points(x_min, y_min, x_max, y_max, n):  
    '''
    使用超拉丁立方采样矩形区域内的点
    '''
    ub = np.array([x_max, y_max])
    lb = np.array([x_min, y_min])
    internal_points = lb + (ub - lb) * lhs(2, n)

    return internal_points

训练代码

import numpy as np
import matplotlib.pyplot as plt
import torch
import torch.nn as nn
from torch.autograd import grad
import time
import datetime
import points_generator as pg

device = torch.device("cuda") if torch.cuda.is_available() else torch.device("cpu")


# 设置空间尺寸
length = 10.0
width = 10.0

#真空磁导率
mu0 = 4*np.pi*1e-7
mu_air = 1
mu_copper = 1

# Generate boundary points using NumPy
# boundary_bottom = np.stack([np.linspace(0, length, N_each_boundary), np.zeros(N_each_boundary)], axis=1)  # 下边界
# boundary_top = np.stack([np.linspace(0, length, N_each_boundary), np.ones(N_each_boundary) * width], axis=1)  # 上边界
# boundary_left = np.stack([np.zeros(N_each_boundary), np.linspace(0, width, N_each_boundary)], axis=1)  # 左边界
# boundary_right = np.stack([np.ones(N_each_boundary) * length, np.linspace(0, width, N_each_boundary)], axis=1)  # 右边界

boundary,_,_ = pg.generate_boundary_points(0,0,width,length,8)

# Concatenate boundary points
boundary = torch.tensor(boundary,dtype=torch.float32).to(device)

# 生成公共界面上的点
interface_points,interface_vertical,interface_horizontal = pg.generate_boundary_points(4,2,6,8,8)

# Convert NumPy arrays to Torch tensors
interface_points = torch.tensor(interface_points,dtype=torch.float32).to(device)
interface_vertical = torch.tensor(interface_vertical,dtype=torch.float32).to(device)
interface_horizontal = torch.tensor(interface_horizontal,dtype=torch.float32).to(device)


# 生成板内部的点
air_points = pg.generate_internal_points(0,0,width,length,8)
mask = pg.generate_mask(air_points,4,2,6,8)
air_points = air_points[~mask]
copper_points = pg.generate_internal_points(4,2,6,8,8)


# Convert NumPy arrays to Torch tensors
air_points = torch.tensor(air_points,dtype=torch.float32).to(device)
copper_points = torch.tensor(copper_points,dtype=torch.float32).to(device)

# Concatenate all points
internal_points = torch.cat([copper_points,air_points],dim=0)
all_points = torch.cat([boundary, interface_points, copper_points,air_points], dim=0)


def pointsplot():
    plt.figure(figsize=(8, 8))
    plt.scatter(air_points[:, 0].numpy(), air_points[:, 1].numpy(), s=10, c='blue', alpha=0.5,label='air_points')
    plt.scatter(copper_points[:, 0].numpy(), copper_points[:, 1].numpy(), s=10, c='red', alpha=0.5,label='copper_points')
    plt.scatter(boundary[:, 0].numpy(), boundary[:, 1].numpy(), s=10, c='black', alpha=0.5,label='boundary')
    plt.scatter(interface_points[:, 0].numpy(), interface_points[:, 1].numpy(), s=10, c='orange', alpha=0.5,
                label='interface_points')
    plt.xlabel('X ')
    plt.ylabel('Y')
    plt.title('points visualize')
    plt.legend()
    plt.grid(True)
    plt.show()

# 打印各部分点数和总点数
print(f"边界点数: {boundary.shape[0]}")
print(f"界面点数: {interface_points.shape[0]}")
print(f"内部点数: {internal_points.shape[0]}")
print(f"总点数: {all_points.shape[0]}")
pointsplot()


class layer(nn.Module):

    def __init__(self, n_in, n_out, activation):
        super().__init__()
        self.layer = nn.Linear(n_in, n_out)
        self.activation = activation

    def forward(self, x):
        x = self.layer(x)
        if self.activation:
            x = self.activation(x)
        return x


def weights_init(m):
    if isinstance(m, nn.Linear):
        torch.nn.init.xavier_uniform_(m.weight.data)
        torch.nn.init.zeros_(m.bias.data)
class DNN(nn.Module):
    # dim_in 和 dim_out：分别为网络的输入和输出维度，n_layer 和 n_node：n_layer 表示隐藏层的数量， activation为激活函数，ub 和 lb：分别是区域的上界和下界，用于输入的归一化处理。
    def __init__(self, dim_in, dim_out, n_layer, n_node, activation = nn.Tanh()):
        super().__init__()
        self.net = nn.ModuleList()
        self.net.append(layer(dim_in, n_node, activation))  #第一个输入层
        for _ in range(n_layer):
            self.net.append(layer(n_node, n_node, activation))
        self.net.append(layer(n_node, dim_out, activation=None))   #最后一个输出层

        self.net.apply(weights_init)  # weights_init 函数对网络的权重和偏置进行初始化（例如 Xavier 初始化），从而提供较好的初始收敛性。

    def forward(self, x):
        out = x
        for layer in self.net:
            out = layer(out)
        return out


class PINN:

    def __init__(self) -> None:
        # 定义每个介质的网络
        self.nets = [
            DNN(dim_in=2, dim_out=1, n_layer=6, n_node=80).to(device),
            DNN(dim_in=2, dim_out=1, n_layer=6, n_node=80).to(device),
            # DNN(dim_in=2, dim_out=1, n_layer=6, n_node=50).to(device)
        ]

        self.lbfgs = torch.optim.LBFGS(
            params=[p for net in self.nets for p in net.parameters()],
            lr=0.1,
            max_iter=20000,
            max_eval=20000,
            tolerance_grad=1e-5,
            tolerance_change= 1 * np.finfo(float).eps,
            history_size=50,
            line_search_fn="strong_wolfe",
        )

        #分别定义了 LBFGS 和 Adam 两种优化器。LBFGS 是一种基于二阶信息的优化算法，适用于物理信息神经网络的优化；Adam 是一种常用的一阶优化方法，用于更新网络权重。

        # self.adam = [
        #     torch.optim.Adam(net.parameters(), lr=1e-4) for net in self.nets
        # ]
        self.adam = torch.optim.Adam(
            params=[p for net in self.nets for p in net.parameters()], lr=1e-5
        )
        self.losses = {"bc": [], "pde": [], "interface": [], "loss": []}
        self.iter = 0

    def predict(self, xy):
        output = torch.zeros(xy.shape[0], 1) # Ensure output has the right shape

        mask = pg.generate_mask(xy,4,2,6,8)

        output[mask] = self.nets[0](xy[mask])

        output[~mask] = self.nets[1](xy[~mask])

        return output


    def bc_loss(self):

        boundary.requires_grad = True

        A = self.predict(boundary)
        # 边界条件处理
        mse_bc = torch.mean(A**2)
        #边界损失
        l_bc = mse_bc
        return l_bc

    def interface_bc_loss(self):
        """
            计算公共边界条件损失函数
        """

        interface_vertical.requires_grad = True
        interface_horizontal.requires_grad = True



        A_interface_1_vertical = self.nets[0](interface_vertical)
        A_interface_2_vertical = self.nets[1](interface_vertical)
        A_interface_1_horizontal = self.nets[0](interface_horizontal)
        A_interface_2_horizontal = self.nets[1](interface_horizontal)
        # 第一边界条件
        l_1 = torch.mean(torch.square(A_interface_1_vertical - A_interface_2_vertical)) + torch.mean(torch.square(A_interface_1_horizontal - A_interface_2_horizontal))


        # 第二边界条件
        grad_1_vertical = grad(A_interface_1_vertical.sum(), interface_vertical, create_graph=True)[0][:, 0:1]
        grad_2_vertical = grad(A_interface_2_vertical.sum(), interface_vertical, create_graph=True)[0][:, 0:1]

        grad_1_horizontal = grad(A_interface_1_horizontal.sum(), interface_horizontal, create_graph=True)[0][:, 1:2]
        grad_2_horizontal = grad(A_interface_2_horizontal.sum(), interface_horizontal, create_graph=True)[0][:, 1:2]
        grad_1 = torch.cat([grad_1_vertical, grad_1_horizontal])
        grad_2 = torch.cat([grad_2_vertical, grad_2_horizontal])
        l_2 = torch.mean(torch.square(mu_copper*grad_1-mu_air*grad_2))



        return l_1+l_2



    def pde_loss(self, xy):
        xy = xy.clone().detach().requires_grad_(True)

        mask = pg.generate_mask(xy, 4, 2, 6, 8)

        # 铜导体区域
        xy_copper = xy[mask]
        A_copper = self.nets[0](xy_copper)
        grad_A_copper = grad(A_copper.sum(), xy_copper, create_graph=True)[0]
        dA_dx_copper = grad_A_copper[:, 0:1]
        dA_dy_copper = grad_A_copper[:, 1:2]

        d2A_dx2_copper = grad(dA_dx_copper.sum(), xy_copper, create_graph=True)[0][:, 0:1]
        d2A_dy2_copper = grad(dA_dy_copper.sum(), xy_copper, create_graph=True)[0][:, 1:2]

        nu_copper = 1 / (mu_copper)
        pde_copper = d2A_dx2_copper + d2A_dy2_copper + 1
        loss_copper = torch.mean((nu_copper * pde_copper) ** 2)

        # 空气区域
        xy_air = xy[~mask]
        A_air = self.nets[1](xy_air)
        grad_A_air = grad(A_air.sum(), xy_air, create_graph=True)[0]
        dA_dx_air = grad_A_air[:, 0:1]
        dA_dy_air = grad_A_air[:, 1:2]

        d2A_dx2_air = grad(dA_dx_air.sum(), xy_air, create_graph=True)[0][:, 0:1]
        d2A_dy2_air = grad(dA_dy_air.sum(), xy_air, create_graph=True)[0][:, 1:2]

        nu_air = 1 / (mu_air)
        pde_air = d2A_dx2_air + d2A_dy2_air
        loss_air = torch.mean((nu_air * pde_air) ** 2)

        return loss_copper + loss_air

    def closure(self):

        self.lbfgs.zero_grad()
        self.adam.zero_grad()

        mse_bc = self.bc_loss()
        mse_pde = self.pde_loss(internal_points)

        mse_interface = self.interface_bc_loss()

        loss =  mse_pde +mse_bc+mse_interface

        loss.backward()

        self.losses["bc"].append(mse_bc.detach().cpu().item())
        self.losses["pde"].append(mse_pde.detach().cpu().item())
        self.losses["interface"].append(mse_interface.detach().cpu().item())
        self.losses["loss"].append(loss.cpu().item())

        self.iter += 1

        print(
            f"\r It: {self.iter} Loss: {loss.item():.5e} BC: {mse_bc.item():.3e} PDE: {mse_pde.item():.3e} Interface: {mse_interface.item():.3e}",
            end="",
        )

        if self.iter % 100 == 0:
            print("")

        return loss

    def save_model(self, path):
        """
        保存模型和训练状态
        :param path: 保存文件路径
        """
        save_dict = {
            "nets": [net.state_dict() for net in self.nets],  # 保存所有子网络的参数
            "optimizers": [self.adam.state_dict(), self.lbfgs.state_dict()],  # 保存优化器状态
            "losses": self.losses,  # 保存损失记录
            "iteration": self.iter,  # 保存当前迭代次数
        }
        torch.save(save_dict, path)
        print(f"模型已保存到 {path}")

    def load_model(self, path):
        """
        加载模型和训练状态
        :param path: 加载文件路径
        """
        checkpoint = torch.load(path, map_location=device)

        # 加载每个网络的参数
        for net, state_dict in zip(self.nets, checkpoint["nets"]):
            net.load_state_dict(state_dict)

        # 加载每个优化器的状态
        for opt, state_dict in zip([self.adam, self.lbfgs], checkpoint["optimizers"]):
            opt.load_state_dict(state_dict)

        # 加载其他训练状态
        self.losses = checkpoint["losses"]
        self.iter = checkpoint["iteration"]
        print(f"模型已从 {path} 加载")


def plotLoss(losses_dict, path, info=["BC", "PDE","Interface","Loss"]):
    fig, axes = plt.subplots(1, 4, sharex=True, sharey=True, figsize=(12, 6))
    axes[0].set_yscale("log")
    for i, j in zip(range(4), info):
        axes[i].plot(losses_dict[j.lower()])
        axes[i].set_title(j)
    plt.show()
    fig.savefig(path)



if __name__ == "__main__":

    pinn = PINN()
    start_time = time.time()

    for i in range(10000):
        pinn.closure()
        pinn.adam.step()
    pinn.lbfgs.step(pinn.closure)

    print("--- %s seconds ---" % (time.time() - start_time))
    print(f'-- {(time.time() - start_time)/60} mins --')
    now = datetime.datetime.now()
    # 格式化时间为字符串，例如：2024-11-28_15-30-00
    time_str = now.strftime("%m%d-%H%M")
    pinn.save_model( f"Param{time_str}.pt")
    plotLoss(pinn.losses, f"LossCurve{time_str}.png", )

预测代码

import numpy as np
import matplotlib.pyplot as plt
import torch
from main import PINN, mu0
from scipy.interpolate import griddata

device = torch.device("cuda") if torch.cuda.is_available() else torch.device("cpu")



pinn = PINN()
# 定义网格的分辨率
grid_size_x = 200
grid_size_y = 200
x_min = 0
x_max = 10
y_min =0
y_max = 10


# 生成x和y坐标的等距网格点
x = np.linspace(x_min, x_max, grid_size_x)
y = np.linspace(y_min, y_max, grid_size_y)
X, Y = np.meshgrid(x, y)

# 将网格点转化为输入格式（grid_size * grid_size, 2）
xy_grid = np.hstack((X.reshape(-1, 1), Y.reshape(-1, 1)))
xy_grid = torch.tensor(xy_grid, dtype=torch.float32).to(device)
time_str = '0531-1718'
pinn.load_model(f"Param{time_str}.pt")

# 预测温度
with torch.no_grad():
    A_pred = pinn.predict(xy_grid)*mu0

# 将预测的温度值转化为numpy格式，并恢复为网格形状
A_pred = A_pred.detach().cpu().numpy().reshape(Y.shape)



plt.figure(figsize=(8, 6))
plt.contourf(x,y,A_pred,levels=50, cmap='jet')
plt.colorbar(label='A')
plt.xlabel('X')
plt.ylabel('Y')
plt.title('Predicted A Distribution')
plt.savefig(f"Sol{time_str}.jpg")
plt.show()

# 计算磁场分量 Bx 和 By
Ay, Ax = np.gradient(A_pred, y, x)
Bx = Ay
By = -Ax
B = np.sqrt(Bx**2+By**2)

skip = 10  # 每隔 skip 个点绘制一个箭头

# 绘制磁场矢量图
plt.figure(figsize=(8, 6))
# 使用 B_magnitude 作为颜色的依据
quiver = plt.quiver(X[::skip, ::skip], Y[::skip, ::skip], Bx[::skip, ::skip], By[::skip, ::skip],
                    B[::skip, ::skip], cmap='rainbow')
#设置图的背景色
# plt.gca().set_facecolor('black')
plt.colorbar(quiver, label='Magnetic Field Magnitude')
plt.title('Magnetic Field Vector Plot')
plt.xlabel('x')
plt.ylabel('y')
plt.grid()
plt.show()

如果代码运行有问题或对内容有疑问，都可以后台私信哦！

期待下一期PINNs案例代码分享吧！

基于物理信息的神经网络是一种解决偏微分方程计算问题的全新方法…

有关PINN基础详见：PINNs案例——中心热源温度场预测问题的torch代码

今日分享代码案例：二维带电流源磁场计算

该案例参考学习论文：[1]张宇娇，孙宏达，赵志涛，徐斌，黄雄峰．基于物理信息神经网络的电磁场计算方法[J/OL]．电工技术学报.

结果前瞻

磁场矢量

问题场景

  • **求解域：**案例选取边长为 10 的正方形区域作为研究对象，即定义域为
    $ \Omega = [0, 10] \times [0, 10] 。以正方形左下角为原点，正方形区域内包含一个带电流密度为 1 的导体，导体区域为 。以正方形左下角为原点，正方形区域内包含一个带电流密度为1的导体，导体区域为 。以正方形左下角为原点，正方形区域内包含一个带电流密度为1的导体，导体区域为 \Omega = [4, 2] \times [6, 8] $。
  • **参数设定：**导体处 J = 1 J = 1 J=1，导体材料为铜（copper），其余区域为空气（air），磁导率 γ \gamma γ 均为1
  • **边界条件：**方形四周边界的矢量磁位为0（理想导体边界）
  • **求解量：**区域内二维稳态磁场分布(矢量)

区域采样点

PINNs理论求解

在本案例中，我们引入矢量磁位 $ \mathbf{A} $ 来辅助磁场的计算。矢量磁位A设定为：
B = ∇ × A \mathbf{B} = \nabla \times \mathbf{A} B=∇×A
其中， B \mathbf{B} B 表示磁场强度，该关系表明磁场是矢量磁位的旋度。

在该正方形区域内，矢量磁位 $ A_z(x, y) $ 满足麦克斯韦方程，引入矢量磁位A后在二维稳态下，其数学表达式为：

( ∂ 2 A z ∂ x 2 + ∂ 2 A z ∂ y 2 ) = − μ J z \left( \frac{\partial^2 A_z}{\partial x^2} + \frac{\partial^2 A_z}{\partial y^2} \right) =-\mu J_z (∂x2∂2Az​​+∂y2∂2Az​​)=−μJz​

其中：

  • $ \mu $ 为材料的磁导率，表示材料对磁场的响应能力；
  • $ J_z $ 为电流密度在 $ z $ 方向的分量，表示单位面积内电流的大小；
  • $ A_z(x, y) $ 为矢量磁位在 $ z $ 方向的分量，通常用于描述二维平面内的磁场分布。

这种方程形式常见于静电学或静磁学中的泊松方程，用于描述在存在源（如电荷或电流）的情况下，势函数（如矢量磁位）的空间分布。矢量磁位的引入使得磁场的计算更为简便，尤其是在处理对称性问题或边界条件时。

本次案例中，由于存在铜导体区域（有源）和空气域（无源）,故采用PINNs案例——多介质分区热传导中相同的分区思想。此处不过多介绍

损失函数可以表示为：

  1. **PDE 残差项：**控制方程如下：
    Loss PDE = 1 N c ∑ i = 1 N c ∣ f ( x c i , y c i ) ∣ 2 \text{Loss}_{\text{PDE}} = \frac{1}{N_c} \sum_{i=1}^{N_c} \left| f\left( x_c^i, y_c^i \right) \right|^2 LossPDE​=Nc​1​i=1∑Nc​​ ​f(xci​,yci​) ​2
    f ( x c i , y c i ) = ( ∂ 2 A z ∂ x 2 + ∂ 2 A z ∂ y 2 ) + μ J z f\left( x_c^i, y_c^i \right) = \left( \frac{\partial^2 A_z}{\partial x^2} + \frac{\partial^2 A_z}{\partial y^2} \right) +\mu J_z f(xci​,yci​)=(∂x2∂2Az​​+∂y2∂2Az​​)+μJz​

注意： 由于PDE损失中需加入 μ 0 \mu_{0} μ0​，而 μ 0 = 4 π ∗ 10 − 7 \mu_{0} = 4\pi * 10^{-7} μ0​=4π∗10−7 数量级过小，会对神经网络的训练带来较大困难。故需对PDE方程进行简单的无量纲化操作，令 A ′ = A / μ 0 A' = A/\mu_0 A′=A/μ0​, 于是可以得到 f ( x c i , y c i ) = ( ∂ 2 A z ′ ∂ x 2 + ∂ 2 A z ′ ∂ y 2 ) + μ r J z f\left( x_c^i, y_c^i \right) = \left( \frac{\partial^2 A_z'}{\partial x^2} + \frac{\partial^2 A_z'}{\partial y^2} \right) +\mu_r J_z f(xci​,yci​)=(∂x2∂2Az′​​+∂y2∂2Az′​​)+μr​Jz​.预测后再将A乘 μ 0 \mu_0 μ0​以还原A

  1. 边界条件项：在多介质问题中，边界条件除矢量磁位A=0的方形边界条件，还要添加不同介质交界面的连续性边界条件。

  • 固定温度边界：
    Loss BC = 1 N b ∑ i = 1 N b ∣ A ( x b i , y b i ) ∣ 2 \text{Loss}_{\text{BC}} = \frac{1}{N_b} \sum_{i=1}^{N_b} \left|A\left( x_b^i, y_b^i \right) \right|^2 LossBC​=Nb​1​i=1∑Nb​​ ​A(xbi​,ybi​) ​2

  • 介质交界面：
    此案例中以A的第一类和第二类边界条件作为介质交界面的连续性条件

第一类边界条件：
Loss I1 = 1 N I ∑ j = 1 N I ∣ A p ( x I j ) − A p + 1 ( x I j ) ∣ 2 \text{Loss}_{\text{I1}} = \frac{1}{N_{I}}\sum_{j=1}^{N_{I}}\left| A_{p}\left(x_{I}^{j}\right)-A_{p+1}\left(x_{I}^{j}\right)\right|^{2} LossI1​=NI​1​j=1∑NI​​ ​Ap​(xIj​)−Ap+1​(xIj​) ​2
第二类边界条件（不同介质在边界面上矢量磁位法向导数与磁导率的乘积连续）：
Loss I2 = 1 N I ∑ j = 1 N I ∣ γ p ∂ A p ( x I j ) ∂ n − γ p + 1 ∂ A p + 1 ( x I j ) ∂ n ∣ 2 \text{Loss}_{\text{I2}} = \frac{1}{N_{I}}\sum_{j=1}^{N_{I}}\left|\gamma_{p}\frac{\partial A_{p}\left(x_{I}^{j}\right)}{\partial n}-\gamma_{p+1}\frac{\partial A_{p+1}\left(x_{I}^{j}\right)}{\partial n}\right|^{2} LossI2​=NI​1​j=1∑NI​​ ​γp​∂n∂Ap​(xIj​)​−γp+1​∂n∂Ap+1​(xIj​)​ ​2
其中：

  • $ N_I $ 表示介质边界点的点数
  • $ A_p $ 表示第P个介质的矢量磁位
  • $ \gamma_ {p} $ 表示第P个介质的热导率

训练

采用10000次adam优化与lbfgs优化器进行优化。各项损失收敛过程如下，收敛效果良好：

损失曲线

复现结果

矢量磁位

磁场标量

磁场矢量

comsol仿真结果

矢量磁位

磁场

可见结果吻合良好

pytorch代码

为避免繁杂的点生成，作者编写了一组点生成函数，使用时单独建立points_generator.py，并在训练代码中引入。

points_generator.py

import numpy as np
from pyDOE import lhs

def generate_boundary_points(x_min, y_min, x_max, y_max, n):
    # 计算每条边上的点数
    num_points_x = int((x_max - x_min) * n) + 1
    num_points_y = int((y_max - y_min) * n) + 1

    # 生成每条边上的点
    bottom_edge = np.array([np.linspace(x_min, x_max, num_points_x), np.full(num_points_x, y_min)]).T
    top_edge = np.array([np.linspace(x_min, x_max, num_points_x), np.full(num_points_x, y_max)]).T
    left_edge = np.array([np.full(num_points_y, x_min), np.linspace(y_min, y_max, num_points_y)]).T
    right_edge = np.array([np.full(num_points_y, x_max), np.linspace(y_min, y_max, num_points_y)]).T

    # 合并所有边上的点，去掉重复的角点
    boundary_points = np.concatenate([bottom_edge, top_edge[1:], left_edge[1:-1], right_edge[1:-1]])
    boundary_vertical = np.concatenate([left_edge[1:-1], right_edge[1:-1]])
    boundary_horizontal = np.concatenate([bottom_edge, top_edge[1:]])

    return boundary_points,boundary_vertical,boundary_horizontal

def generate_internal_points(x_min, y_min, x_max, y_max, n):
    # 计算内部点的网格
    x = np.linspace(x_min, x_max, int((x_max - x_min) * n) + 1)
    y = np.linspace(y_min, y_max, int((y_max - y_min) * n) + 1)
    x_grid, y_grid = np.meshgrid(x, y)

    # 去掉边界点
    x_grid = x_grid[1:-1, 1:-1]
    y_grid = y_grid[1:-1, 1:-1]

    # 将网格点转换为一维数组
    internal_points = np.array([x_grid.flatten(), y_grid.flatten()]).T

    return internal_points

def generate_mask(points, x_min, y_min, x_max, y_max):
    # 检查每个点是否在矩形区域内
    mask = (points[:, 0] >= x_min) & (points[:, 0] <= x_max) & (points[:, 1] >= y_min) & (points[:, 1] <= y_max)

    return mask


def lhs_internal_points(x_min, y_min, x_max, y_max, n):  
    '''
    使用超拉丁立方采样矩形区域内的点
    '''
    ub = np.array([x_max, y_max])
    lb = np.array([x_min, y_min])
    internal_points = lb + (ub - lb) * lhs(2, n)

    return internal_points

训练代码与预测代码请关注：

本文来自互联网用户投稿，该文观点仅代表作者本人，不代表本站立场。本站仅提供信息存储空间服务，不拥有所有权，不承担相关法律责任。如若转载，请注明出处：http://www.coloradmin.cn/o/2395522.html

如若内容造成侵权/违法违规/事实不符，请联系多彩编程网进行投诉反馈，一经查实，立即删除！

相关文章

算法打开13天

算法打开13天

41.前 K 个高频元素 &#xff08;力扣347题&#xff09; 给你一个整数数组 nums 和一个整数 k &#xff0c;请你返回其中出现频率前 k 高的元素。你可以按 任意顺序 返回答案。 示例 1: 输入: nums [1,1,1,2,2,3], k 2 输出: [1,2]示例 2: 输入: nums [1], k 1 输出: …
阅读更多...
广告拦截器：全方位拦截，畅享无广告体验

广告拦截器：全方位拦截，畅享无广告体验

在数字时代&#xff0c;广告无处不在。无论是浏览网页、使用社交媒体&#xff0c;还是观看视频&#xff0c;广告的频繁弹出常常打断我们的体验&#xff0c;让人不胜其烦。更令人担忧的是&#xff0c;一些广告可能包含恶意软件&#xff0c;威胁我们的设备安全和个人隐私。AdGuar…
阅读更多...
主数据编码体系全景解析：从基础到高级的编码策略全指南

主数据编码体系全景解析：从基础到高级的编码策略全指南

在数字化转型的浪潮中&#xff0c;主数据管理&#xff08;MDM&#xff09;已成为企业数字化转型的基石。而主数据编码作为MDM的核心环节&#xff0c;其设计质量直接关系到数据管理的效率、系统的可扩展性以及业务决策的准确性。本文将系统性地探讨主数据编码的七大核心策略&…
阅读更多...
Selenium操作指南（全）

Selenium操作指南（全）

&#x1f345; 点击文末小卡片&#xff0c;免费获取软件测试全套资料&#xff0c;资料在手&#xff0c;涨薪更快 大家好&#xff0c;今天带大家一起系统的学习下模拟浏览器运行库Selenium&#xff0c;它是一个用于Web自动化测试及爬虫应用的重要工具。 Selenium测试直接运行在…
阅读更多...
智绅科技——科技赋能健康养老，构建智慧晚年新生态

智绅科技——科技赋能健康养老，构建智慧晚年新生态

当老龄化浪潮与数字技术深度碰撞&#xff0c;智绅科技以 “科技赋能健康&#xff0c;智慧守护晚年” 为核心理念&#xff0c;锚定数字健康与养老服务赛道&#xff0c;通过人工智能、物联网、大数据等技术集成&#xff0c;为亚健康群体与中老年人群构建 “监测 - 预防 - 辅助 - …
阅读更多...
STM32通过KEIL pack包轻松移植LVGL，并学会使用GUI guider

STM32通过KEIL pack包轻松移植LVGL，并学会使用GUI guider

先展示最终实现的功能效果如下&#xff1a; 1.目的与意义 之前在学习STM32移植LVGL图形库的时候&#xff0c;搜到的很多教程都是在官网下载LVGL的文件包&#xff0c;然后一个个文件包含进去&#xff0c;还要添加路径&#xff0c;还要给文件改名字&#xff0c;最后才能修改程序…
阅读更多...
SQL的查询优化

SQL的查询优化

1. 查询优化器 1.1. SQL语句执行需要经历的环节 解析阶段&#xff1a;语法分析和语义检查&#xff0c;确保语句正确&#xff1b;优化阶段&#xff1a;通过优化器生成查询计划&#xff1b;执行阶段&#xff1a;由执行器根据查询计划实际执行操作。 1.2. 查询优化器 查询优化器…
阅读更多...
MCU如何从向量表到中断服务

MCU如何从向量表到中断服务

目录 1、中断向量表 2、编写中断服务例程 中断处理的核心是中断向量表&#xff08;IVT&#xff09;&#xff0c;它是一个存储中断服务例程&#xff08;ISR&#xff09;地址的内存结构。当中断发生时&#xff0c;MCU通过IVT找到对应的ISR地址并跳转执行。本文将深入探讨MCU&am…
阅读更多...
Linux线程同步实战：多线程程序的同步与调度

Linux线程同步实战：多线程程序的同步与调度

个人主页&#xff1a;chian-ocean 文章专栏-Linux Linux线程同步实战&#xff1a;多线程程序的同步与调度 个人主页&#xff1a;chian-ocean文章专栏-Linux 前言&#xff1a;为什么要实现线程同步线程饥饿&#xff08;Thread Starvation&#xff09;示例&#xff1a;抢票问题 …
阅读更多...
【MySQL】事务及隔离性

【MySQL】事务及隔离性

目录 一、什么是事务 &#xff08;一&#xff09;概念 &#xff08;二&#xff09;事务的四大属性 &#xff08;三&#xff09;事务的作用 &#xff08;四&#xff09;事务的提交方式 二、事务的启动、回滚与提交 &#xff08;一&#xff09;事务的启动、回滚与提交 &am…
阅读更多...
yolo目标检测助手：具有模型预测、图像标注功能

yolo目标检测助手：具有模型预测、图像标注功能

在人工智能浪潮席卷各行各业的今天&#xff0c;计算机视觉模型&#xff08;如 YOLO&#xff09;已成为目标检测领域的标杆。然而&#xff0c;模型的强大能力需要直观的界面和便捷的工具才能充分发挥其演示、验证与迭代优化的价值。为此&#xff0c;我开发了一款基于 WPF 的桌面…
阅读更多...
2022 RoboCom 世界机器人开发者大赛(睿抗 caip) -高职组（国赛）解题报告 | 科学家

2022 RoboCom 世界机器人开发者大赛(睿抗 caip) -高职组（国赛）解题报告 | 科学家

前言 题解 2022 RoboCom 世界机器人开发者大赛(睿抗 caip) -高职组&#xff08;国赛&#xff09;。 最后一题还考验能力&#xff0c;需要找到合适的剪枝。 RC-v1 智能管家 分值: 20分 签到题&#xff0c;map的简单实用 #include <bits/stdc.h>using namespace std;int…
阅读更多...
基于物联网（IoT）的电动汽车（EVs）智能诊断

基于物联网（IoT）的电动汽车（EVs）智能诊断

我是穿拖鞋的汉子&#xff0c;魔都中坚持长期主义的汽车电子工程师。 老规矩&#xff0c;分享一段喜欢的文字&#xff0c;避免自己成为高知识低文化的工程师&#xff1a; 做到欲望极简&#xff0c;了解自己的真实欲望&#xff0c;不受外在潮流的影响&#xff0c;不盲从&#x…
阅读更多...
JDBC+HTML+AJAX实现登陆和单表的CRUD

JDBC+HTML+AJAX实现登陆和单表的CRUD

JDBCHTMLAJAX实现登陆和单表的CRUD 导入maven依赖 <?xml version"1.0" encoding"UTF-8"?><project xmlns"http://maven.apache.org/POM/4.0.0" xmlns:xsi"http://www.w3.org/2001/XMLSchema-instance" xsi:schemaLocatio…
阅读更多...
【C++】位图详解（一文彻底搞懂位图的使用方法与底层原理）

【C++】位图详解（一文彻底搞懂位图的使用方法与底层原理）

&#x1f308; 个人主页&#xff1a;谁在夜里看海. &#x1f525; 个人专栏&#xff1a;《C系列》《Linux系列》 ⛰️ 天高地阔&#xff0c;欲往观之。 目录 1.位图的概念 2.位图的使用方法 定义与创建 设置和清除 位访问和检查 转换为其他格式 3.位图的使用场景 1.快速…
阅读更多...
【笔记】开源通用人工智能代理 Suna 部署全流程准备清单（Windows 系统）

【笔记】开源通用人工智能代理 Suna 部署全流程准备清单（Windows 系统）

#工作记录 一、基础工具与环境 开发工具 Git 或 GitHub Desktop&#xff08;代码管理&#xff09;Docker Desktop&#xff08;需启用 WSL2&#xff0c;容器化部署&#xff09;Python 3.11&#xff08;推荐版本&#xff0c;需添加到系统环境变量&#xff09;Node.js LTS&#xf…
阅读更多...
海康工业相机SDK二次开发（VS+QT+海康SDK+C++）

海康工业相机SDK二次开发（VS+QT+海康SDK+C++）

前言 工业相机在现代制造和工业自动化中扮演了至关重要的角色&#xff0c;尤其是在高精度、高速度检测中。海康威视工业相机以其性能稳定、图像质量高、兼容性强而受到广泛青睐。特别是搞机器视觉的小伙伴们跟海康打交道肯定不在少数&#xff0c;笔者在平常项目中跟海康相关人…
阅读更多...
深度学习|pytorch基本运算-乘除法和幂运算

深度学习|pytorch基本运算-乘除法和幂运算

【1】引言 前序学习进程中&#xff0c;已经对pytorch张量数据的生成和广播做了详细探究&#xff0c;文章链接为&#xff1a; 深度学习|pytorch基本运算-CSDN博客 深度学习|pytorch基本运算-广播失效-CSDN博客 上述探索的内容还止步于张量的加减法&#xff0c;在此基础上&am…
阅读更多...
4.2.4 Spark SQL 数据写入模式

4.2.4 Spark SQL 数据写入模式

在本节实战中&#xff0c;我们详细探讨了Spark SQL中数据写入的四种模式&#xff1a;ErrorIfExists、Append、Overwrite和Ignore。通过具体案例&#xff0c;我们演示了如何使用mode()方法结合SaveMode枚举类来控制数据写入行为。我们首先读取了一个JSON文件生成DataFrame&#…
阅读更多...
论文笔记： Urban Region Embedding via Multi-View Contrastive Prediction

论文笔记： Urban Region Embedding via Multi-View Contrastive Prediction

AAAI 2024 1 INTRO 之前基于多视图的region embedding工作大多遵循相同的模式 单独的单视图表示多视图融合 但这种方法存在明显的局限性&#xff1a;忽略了不同视图之间的信息一致性 一个区域的多个视图所携带的信息是高度相关的&#xff0c;因此它们的表示应该是一致的如果能…
阅读更多...
推荐文章
最新文章

Copyright @ 2022~2023