更新新版本的热传导求解器

app/FuelRodSim/HeatEquationData.py

@@ -1,8 +1,11 @@
-import numpy as np
+
 from sympy import *
+from fealpy.backend import backend_manager as bm
+
 
 class Parabolic2dData: 
     def __init__(self,u, x, y, t, D=[0, 1, 0, 1], T=[0, 1]):
+
         self._domain = D 
         self._duration = T 
         self.u = lambdify([x,y,t], sympify(u))
@@ -19,15 +22,17 @@ def solution(self, p, t):
         y = p[..., 1]
         return self.u(x,y,t) 
 
+
     def init_solution(self, p):
         x = p[..., 0]
         y = p[..., 1]
         return self.u(x,y,self._duration[0])
-        
+
     def source(self, p, t):
         x = p[..., 0]
         y = p[..., 1]
-        return self.f(x,y,t)
+        val = self.f(x,y,t)
+        return val
 
     def gradient(self, p, t):
         x = p[..., 0]
@@ -38,20 +43,20 @@ def gradient(self, p, t):
     def dirichlet(self, p, t):
         return self.solution(p, t)
 
-
-    def source(self, p, t):
-        """
-        @brief 方程右端项 
 
-        @param[in] p numpy.ndarray, 空间点
-        @param[in] t float, 时间点 
+    # def source(self, p, t):
+    #     """
+    #     @brief 方程右端项 
 
-        @return 方程右端函数值
-        """
-        pi = np.pi
-        x = p[..., 0]
-        y = p[..., 1]
-        return np.zeros(x.shape)
+    #     @param[in] p numpy.ndarray, 空间点
+    #     @param[in] t float, 时间点 
+
+    #     @return 方程右端函数值
+    #     """
+    #     pi = np.pi
+    #     x = p[..., 0]
+    #     y = p[..., 1]
+    #     return np.zeros(x.shape)
 
 
     def dirichlet(self, p,t):
@@ -75,7 +80,9 @@ def source(self,p,t):
         x = p[..., 0]
         y = p[..., 1]
         z = p[..., 2]
-        return self.f(x,y,z,t)
+        val = self.f(x,y,z,t)
+        print(val.shape)
+        return val
 
     def solution(self, p, t):
         x = p[..., 0]
@@ -104,7 +111,7 @@ def source(self,p,t):
         return 0
 
     def dirichlet(self, p, t):
-        return np.array([500])
+        return bm.array([500])
 
 class FuelRod2dData:
     def domain(self):
@@ -117,4 +124,4 @@ def source(self, p, t):
         return 0
 
     def dirichlet(self, p, t):
-        return np.array([500])
+        return bm.array([500])

app/FuelRodSim/box_2dexample.py

@@ -8,7 +8,7 @@
 from fealpy.mesh import TriangleMesh
 mesh = TriangleMesh.from_box([0, 1, 0, 1],nx,ny)
 node = mesh.node
-isBdNode = mesh.ds.boundary_node_flag()
+isBdNode = mesh.boundary_node_flag()
 pde = FuelRod2dData()
 Boxslover = HeatEquationSolver(mesh,pde,120,isBdNode,300,alpha_caldding=0.3,layered=False,output='./rusult_boxtest_2d')
 Boxslover.solve()

app/FuelRodSim/box_3dexample.py

@@ -9,7 +9,7 @@
 from fealpy.mesh import TetrahedronMesh
 mesh = TetrahedronMesh.from_box([0, 1, 0, 1, 0, 1],nx,ny,nz)
 node = mesh.node
-isBdNode = mesh.ds.boundary_node_flag()
+isBdNode = mesh.boundary_node_flag()
 pde = FuelRod3dData()
 Boxslover = HeatEquationSolver(mesh,pde,120,isBdNode,300,alpha_caldding=0.08,layered=False,output='./rusult_boxtest_3d')
 Boxslover.solve()

app/FuelRodSim/heat_equation_solver.py

@@ -1,11 +1,11 @@
-import numpy as np
+from fealpy.backend import backend_manager as bm
 import os
-from scipy.sparse.linalg import spsolve
 from fealpy.mesh import TriangleMesh,TetrahedronMesh
 from fealpy.functionspace import LagrangeFESpace
-from fealpy.fem import DiffusionIntegrator, BilinearForm, ScalarMassIntegrator, LinearForm, ScalarSourceIntegrator
+from fealpy.fem import ScalarDiffusionIntegrator, BilinearForm, ScalarMassIntegrator, LinearForm, ScalarSourceIntegrator
 from fealpy.fem.dirichlet_bc import DirichletBC
 import matplotlib.pyplot as plt
+from fealpy.solver import cg
 
 class FuelRod3dData:
     def domain(self):
@@ -18,7 +18,7 @@ def source(self,p,t):
         return 0
 
     def dirichlet(self, p, t):
-        return np.array([500])
+        return bm.array([500])
 
 class HeatEquationSolver:
     def __init__(self, mesh: TriangleMesh, pde:FuelRod3dData, nt, bdnidx, p0, alpha_caldding=4e-4, alpha_inner=8e-4, layered=True, ficdx=None ,cacidx=None,output: str = './result', filename: str = 'temp'):
@@ -53,39 +53,36 @@ def __init__(self, mesh: TriangleMesh, pde:FuelRod3dData, nt, bdnidx, p0, alpha_
         self.threshold = self.create_threshold()
         self.errors = []  # 用于存储每个时间步的误差
 
-    def initialize_output_directory(self):
-        if not os.path.exists(self.output):
-            os.makedirs(self.output)
-
     def create_threshold(self):
         """
         可以分辨布尔数组和编号索引，如果是布尔数组直接传入，如果是索引编号转换为布尔数组
         """
-        if isinstance(self.bdnidx, np.ndarray) and self.bdnidx.dtype == bool:
+        if isinstance(self.bdnidx, bm.DATA_CLASS) and self.bdnidx.dtype == bool:
             return self.bdnidx
         else:
             NN = len(self.mesh.entity('node'))
-            isbdnidx = np.full(NN, False, dtype=bool)
+            isbdnidx = bm.full(NN, False, dtype=bool)
             isbdnidx[self.bdnidx] = True
             return isbdnidx
 
 
     def solve(self):
-        self.p = self.space.function()
+        d = self.space.function()
+        self.p = bm.zeros_like(d)
         self.p += self.p0  
         for n in range(self.nt):
             t = self.duration[0] + n * self.tau
             bform3 = LinearForm(self.space)
             # 这里应该传入后的得到的为qc
             #f=self.pde.source(node,t)
-            source=lambda p: self.pde.source(p, t)
-            print(source)
-            bform3.add_domain_integrator(ScalarSourceIntegrator(source, q=3))
+            source=lambda p, index: self.pde.source(p, t)
+            # source = pde.source
+            bform3.add_integrator(ScalarSourceIntegrator(source))
             self.F = bform3.assembly()
 
             # 组装刚度矩阵
             NC = self.mesh.number_of_cells()
-            alpha = np.zeros(NC)
+            alpha = bm.zeros(NC)
 
             if self.layered:
                 # 假设 ficdx 和 cacidx 是定义好的两个索引列表
@@ -97,40 +94,42 @@ def solve(self):
                 alpha += self.alpha_caldding
 
             bform = BilinearForm(self.space)
-            bform.add_domain_integrator(DiffusionIntegrator(alpha, q=3))
+            bform.add_integrator(ScalarDiffusionIntegrator(alpha, q=3))
             self.K = bform.assembly()
 
             # 组装质量矩阵
             bform2 = BilinearForm(self.space)
-            bform2.add_domain_integrator(ScalarMassIntegrator(q=3))
+            bform2.add_integrator(ScalarMassIntegrator(q=3))
             self.M = bform2.assembly()
             A = self.M + self.alpha_caldding * self.K * self.tau
             b = self.M @ self.p + self.tau * self.F
             if n == 0:
                 A = A
                 b = b
             else:
-                gD=lambda x : self.pde.dirichlet(x,t)
+                gd=lambda x : self.pde.dirichlet(x,t)
                 ipoints = self.space.interpolation_points()
-                gD=gD(ipoints[self.threshold])
-                bc = DirichletBC(space=self.space, gD=gD.reshape(-1,1), threshold=self.threshold)
+                gd=gd(ipoints[self.threshold])
+                bc = DirichletBC(space=self.space, gd=gd, threshold=self.threshold)
                 A, b = bc.apply(A, b)
-                self.p = spsolve(A, b)
+                self.p = cg(A, b, maxiter=5000, atol=1e-14, rtol=1e-14)
             print(self.p)
             # 计算并存储误差，如果 solution 方法存在
             if hasattr(self.pde, 'solution'):
                 exact_solution = self.pde.solution(self.mesh.node, t)
-                error = np.linalg.norm(exact_solution - self.p.flatten('F'))
+                error = bm.linalg.norm(exact_solution - self.p.flatten('F'))
                 self.errors.append(error)
             self.mesh.nodedata['temp'] = self.p.flatten('F')
             name = os.path.join(self.output, f'{self.filename}_{n:010}.vtu')
             self.mesh.to_vtk(fname=name)
-        print('self.p',self.p)
-        print(self.p.shape)
+
+    def initialize_output_directory(self):
+        if not os.path.exists(self.output):
+            os.makedirs(self.output)
 
     def plot_error_over_time(self):
         plt.figure(figsize=(10, 6))
-        plt.plot(np.linspace(self.duration[0], self.duration[1], self.nt), self.errors, marker='o', linestyle='-')
+        plt.plot(bm.linspace(self.duration[0], self.duration[1], self.nt), self.errors, marker='o', linestyle='-')
         plt.title('Error over Time')
         plt.xlabel('Time')
         plt.ylabel('L2 Norm of Error')
@@ -168,7 +167,7 @@ def plot_error(self):
         t = self.duration[1]
         exact_solution = self.pde.solution(self.mesh.node, t)
         numerical_solution = self.p.flatten('F')
-        error = np.abs(exact_solution - numerical_solution)
+        error = bm.abs(exact_solution - numerical_solution)
         print('error',error)
 
         if self.GD == 2:
@@ -218,7 +217,7 @@ def plot_error_heatmap(self):
         t = self.duration[1]
         exact_solution = self.pde.solution(self.mesh.node, t)
         numerical_solution = self.p.flatten('F')
-        error = np.abs(exact_solution - numerical_solution)
+        error = bm.abs(exact_solution - numerical_solution)
         print('error', error)
 
         if self.GD == 2:

app/FuelRodSim/true_solution_2dexample.py

@@ -8,7 +8,7 @@
 mesh = TriangleMesh.from_box([0, 1, 0, 1], nx,ny)
 node = mesh.node
 print(node.shape)
-isBdNode = mesh.ds.boundary_node_flag()
+isBdNode = mesh.boundary_node_flag()
 p0=pde.init_solution(node) #准备一个初值
 Box2dsolver = HeatEquationSolver(mesh,pde,160,isBdNode,p0=p0,alpha_caldding=1,layered=False,output='./rusult_box2dtesttest_to')
 Box2dsolver.solve()

app/FuelRodSim/true_solution_3dexample.py

@@ -1,4 +1,5 @@
 # 三维带真解的测试
+# TODO: 1. 三维真解的测试(其他五个已经通过)
 from fealpy.mesh import TetrahedronMesh
 from app.FuelRodSim.HeatEquationData import Parabolic3dData
 from app.FuelRodSim.heat_equation_solver import HeatEquationSolver
@@ -8,7 +9,7 @@
 nz = 5
 mesh = TetrahedronMesh.from_box([0, 1, 0, 1, 0, 1], nx, ny, nz)
 node = mesh.node
-isBdNode = mesh.ds.boundary_node_flag()
+isBdNode = mesh.boundary_node_flag()
 p0=pde.init_solution(node) #准备一个初值
 Box3DSolver = HeatEquationSolver(mesh, pde, 160, isBdNode, p0=p0, alpha_caldding=1, layered=False, output='./result_box3dtest')
 Box3DSolver.solve()