55class AMSSolver :
66
77 def __init__ (self ):
8- """AMS预条件器求解类初始化。"""
9- # 初始化内部数据结构
10- self .A = None # 边空间刚度矩阵 (稀疏矩阵)
11- self .G = None # 离散梯度矩阵 G
8+
9+ self .A = None
10+ self .G = None
1211 self .Pi_x = None
1312 self .Pi_y = None
1413 self .Pi_z = None
15- self .dim = None # 空间维度 (2 或 3)
14+ self .dim = None
1615 self .n_nodes = 0
1716 self .n_edges = 0
18- self .vertices = None # 顶点坐标数组
19- self .edges = None # 每条边的顶点索引对 (取向已定)
20- # 平滑器与循环参数
17+ self .vertices = None
18+ self .edges = None
19+
2120 self .smoother = 'GaussSeidel'
2221 self .smooth_iters = 1
23- self .jacobi_weight = 1.0 # Jacobi 平滑的松弛因子
24- self .cycle_type = 'V' # 多重网格循环类型 ('V' 或 'W')
25- # 子空间矩阵
26- self .A_G = None # 梯度子空间矩阵 A_G = G^T * A * G
27- self .A_Pi = None # 插值子空间矩阵 A_Pi = Pi^T * A * Pi
22+ self .jacobi_weight = 1.0
23+ self .cycle_type = 'V'
24+
25+ self .A_G = None
26+ self .A_Pi = None
2827
2928 def AMSFEISetup (self , A , vertices , edges , dim = None ):
30- """
31- 初始化有限元结构,包括网格拓扑和自由度编号。
32- 参数:
33- - A: H(curl)空间的刚度矩阵 (scipy.sparse 矩阵).
34- - vertices: 顶点坐标数组,形状 (num_vertices, dim).
35- - edges: 边的顶点索引对列表 (长度 num_edges,每个元素为 (v1, v2)).
36- - dim: 空间维度 (2 或 3),可选参数,不提供则根据 vertices 列数推断。
37- """
38- # 设置维度
29+
3930 if dim is None :
4031 dim = vertices .shape [1 ]
4132 self .dim = dim
42- # 存储矩阵和几何信息
33+
4334 self .A = A
4435 self .vertices = vertices
4536 self .n_nodes = self .vertices .shape [0 ]
@@ -55,8 +46,8 @@ def AMSComputeGPi(self):
5546 num_nodes = self .n_nodes
5647
5748 edges = self .edges
58- ei = edges [:, 0 ].astype (int ) # 所有边的起点索引
59- ej = edges [:, 1 ].astype (int ) # 所有边的终点索引
49+ ei = edges [:, 0 ].astype (int )
50+ ej = edges [:, 1 ].astype (int )
6051 row = bm .repeat (bm .arange (num_edges , dtype = int ), 2 )
6152 col = bm .empty (2 * num_edges , dtype = int )
6253 col [0 ::2 ] = ei
@@ -69,9 +60,7 @@ def AMSComputeGPi(self):
6960 self .A_G = (self .G .T )@(self .A @(self .G ))
7061
7162 def AMSSetAGradient (self , A_grad ):
72- """
73- 直接设置梯度子空间矩阵 A_G。
74- """
63+
7564 self .A_G = A_grad
7665
7766 def AMSComputePi (self ):
@@ -117,7 +106,7 @@ def AMSSetCycle(self, cycle_type='V'):
117106
118107 def AMSSolve (self , b , tol = 1e-6 , maxiter = 100 , x0 = None ):
119108
120- # 初始解
109+
121110 if x0 is None :
122111 x = bm .zeros (self .n_edges ,dtype = self .A .dtype )
123112 else :
@@ -126,20 +115,19 @@ def AMSSolve(self, b, tol=1e-6, maxiter=100, x0=None):
126115 b = bm .array (b )
127116
128117 for it in range (maxiter ):
129- # 计算残差 r = b - A x
118+
130119 r = b - self .A @x
131120 res_norm = bm .linalg .norm (r )
132121 if res_norm < tol :
133- return x , it # 提前收敛
122+ return x , it
134123 el = spsolve (self .A .tril (), r ,'scipy' )
135124
136125 for _ in range (self .smooth_iters ):
137126 el += spsolve (self .A .tril (), r - self .A @ el ,'scipy' )
138127 x = x + el
139- r = b - self .A @x # 更新残差
140- # 2. 梯度子空间校正:求解 A_G e_g = G^T r,然后 x <- x + G e_g
128+ r = b - self .A @x
141129
142- r_g = self .G .T @r # 将残差限制到节点梯度空间
130+ r_g = self .G .T @r
143131
144132 e_g = spsolve (self .A_G , r_g ,'scipy' )
145133
@@ -149,7 +137,7 @@ def AMSSolve(self, b, tol=1e-6, maxiter=100, x0=None):
149137 for _ in range (self .smooth_iters ):
150138 el += spsolve (self .A .tril (), r - self .A @ el ,'scipy' )
151139 x = x + el
152- r = b - self .A @x # 更新残差
140+ r = b - self .A @x
153141 r_pi_combined = self .Pi .T @r
154142 e_pi_combined = bm .zeros (self .Pi .shape [0 ])
155143 ml = GAMGSolver (isolver = 'MG' , ptype = 'V' , sstep = 3 , theta = 0.25 )
@@ -162,11 +150,11 @@ def AMSSolve(self, b, tol=1e-6, maxiter=100, x0=None):
162150 for _ in range (self .smooth_iters ):
163151 el += spsolve (self .A .triu (), r - self .A @ el ,'scipy' )
164152 x = x + el
165- # 检查收敛
153+
166154 r = b - self .A @x
167155 res_norm = bm .linalg .norm (r )
168156 print ("iter:" , it , "res_norm:" , res_norm )
169157 if res_norm < tol :
170158 return x , it + 1
171- # 未在 maxiter 内收敛,返回当前解
159+
172160 return x , maxiter
0 commit comments