diff --git a/app/FuelRodSim/HeatEquationData.py b/app/FuelRodSim/HeatEquationData.py index 1cfcb2ea5..f3db92119 100644 --- a/app/FuelRodSim/HeatEquationData.py +++ b/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]) \ No newline at end of file + return bm.array([500]) \ No newline at end of file diff --git a/app/FuelRodSim/box_2dexample.py b/app/FuelRodSim/box_2dexample.py index 9f6977db5..5d2d26649 100644 --- a/app/FuelRodSim/box_2dexample.py +++ b/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() \ No newline at end of file diff --git a/app/FuelRodSim/box_3dexample.py b/app/FuelRodSim/box_3dexample.py index f346904a5..7b3a22e12 100644 --- a/app/FuelRodSim/box_3dexample.py +++ b/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() \ No newline at end of file diff --git a/app/FuelRodSim/heat_equation_solver.py b/app/FuelRodSim/heat_equation_solver.py index 9da58e11a..639e4d1ec 100644 --- a/app/FuelRodSim/heat_equation_solver.py +++ b/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,12 +94,12 @@ 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 @@ -110,27 +107,29 @@ def solve(self): 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: diff --git a/app/FuelRodSim/true_solution_2dexample.py b/app/FuelRodSim/true_solution_2dexample.py index 3e961f81d..0f964667c 100644 --- a/app/FuelRodSim/true_solution_2dexample.py +++ b/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() diff --git a/app/FuelRodSim/true_solution_3dexample.py b/app/FuelRodSim/true_solution_3dexample.py index 525351f65..96fde7c74 100644 --- a/app/FuelRodSim/true_solution_3dexample.py +++ b/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()