Coverage for /usr/share/miniconda3/envs/dolfin/lib/python3.8/site-packages/block/algebraic/petsc/precond.py: 68%
176 statements
« prev ^ index » next coverage.py v7.2.1, created at 2023-03-20 13:03 +0000
1from __future__ import division
3from builtins import str
4from block.block_base import block_base
5from block.helpers import supports_mpi
6from block.object_pool import vec_pool
7from petsc4py import PETSc
8import dolfin as df
9import numpy as np
10import hashlib
11from time import time
12from dolfin import info
15class precond(block_base):
16 def __init__(self, A, prectype, pdes, nullspace, prefix, options, V=None, create_cb=None, defaults={}):
17 Ad = mat(A)
19 if nullspace:
20 from block.block_util import isscalar
21 ns = PETSc.NullSpace()
22 if isscalar(nullspace): 22 ↛ 23line 22 didn't jump to line 23, because the condition on line 22 was never true
23 ns.create(constant=True)
24 else:
25 ns.create(constant=False, vectors=[v.down_cast().vec() for v in nullspace])
26 try:
27 Ad.setNearNullSpace(ns)
28 except:
29 info('failed to set near null space (not supported in petsc4py version)')
31 # The PETSc options database is global, so we create a unique prefix
32 # that precisely matches the supplied options. The user can override
33 # this prefix if necessary.
34 if prefix is None: 34 ↛ 43line 34 didn't jump to line 43, because the condition on line 34 was never false
35 # Use class name instead of hash of defaults (we assume defaults
36 # are constant for a given class)
37 prefix = self.__class__.__name__
38 if options: 38 ↛ 39line 38 didn't jump to line 39, because the condition on line 38 was never true
39 sha = hashlib.sha256(str(tuple(sorted(options.items()))).encode())
40 prefix += sha.hexdigest()[:8]
41 prefix += ':'
43 self.optsDB = PETSc.Options(prefix)
44 params = defaults.copy()
45 if options: 45 ↛ 46line 45 didn't jump to line 46, because the condition on line 45 was never true
46 params.update(options)
47 for key, val in params.items():
48 # Allow overriding defaults by setting key='del'. Command-line
49 # options (already in optsDB) are never overwritten.
50 if val!='del' and not self.optsDB.hasName(key):
51 self.optsDB.setValue(key, val)
53 if prectype == PETSc.PC.Type.HYPRE and nullspace:
54 if not all(self.optsDB.hasName(n) for n in ['pc_hypre_boomeramg_nodal_coarsen', 'pc_hypre_boomeramg_vec_interp_variant']): 54 ↛ exit, 54 ↛ 572 missed branches: 1) line 54 didn't finish the generator expression on line 54, 2) line 54 didn't jump to line 57, because the condition on line 54 was never false
55 print('Near null space is ignored by hypre UNLESS you also supply pc_hypre_boomeramg_nodal_coarsen and pc_hypre_boomeramg_vec_interp_variant')
57 self.A = A
58 self.petsc_prec = PETSc.PC().create(V.mesh().mpi_comm() if V else None)
59 self.petsc_prec.setOptionsPrefix(prefix)
60 try:
61 self.petsc_prec.setType(prectype)
62 except PETSc.Error as e:
63 if e.ierr == 86:
64 raise ValueError(f'Unknown PETSc type "{prectype}"')
65 else:
66 raise
67 self.petsc_prec.setOperators(Ad)
68 self.petsc_prec.setFromOptions()
70 self.setup_time = -1.
72 def matvec(self, b):
73 # NOTE: we want to be lazy only setup when the action is needed
74 if self.setup_time < 0:
75 T = time()
76 # Create preconditioner based on the options database
77 self.petsc_prec.setUp()
79 setup_time = time()-T
80 comm = self.petsc_prec.getComm().tompi4py()
82 setup_time = comm.allreduce(setup_time)/comm.size
83 info('constructed %s preconditioner in %.2f s'%(self.__class__.__name__, setup_time))
85 self.setup_time = setup_time
87 from dolfin import GenericVector
88 if not isinstance(b, GenericVector): 88 ↛ 89line 88 didn't jump to line 89, because the condition on line 88 was never true
89 return NotImplemented
90 x = self.A.create_vec(dim=1)
91 if len(x) != len(b): 91 ↛ 92line 91 didn't jump to line 92, because the condition on line 91 was never true
92 raise RuntimeError(
93 'incompatible dimensions for PETSc matvec, %d != %d'%(len(x),len(b)))
95 self.petsc_prec.apply(b.down_cast().vec(), x.down_cast().vec())
96 return x
98 def down_cast(self):
99 return self.petsc_prec
101 def __str__(self):
102 return '<%s prec of %s>'%(self.__class__.__name__, str(self.A))
104 def create_vec(self, dim=1):
105 return self.A.create_vec(dim)
108class ML(precond):
109 def __init__(self, A, parameters=None, pdes=1, nullspace=None, prefix=None):
110 super().__init__(A, PETSc.PC.Type.ML,
111 pdes=pdes, nullspace=nullspace, options=parameters, prefix=prefix,
112 defaults={
113 # Symmetry- and PD-preserving smoother
114 'mg_levels_ksp_type': 'chebyshev',
115 'mg_levels_pc_type': 'jacobi',
116 # Fixed number of iterations to preserve linearity
117 'mg_levels_ksp_max_it': 2,
118 'mg_levels_ksp_check_norm_iteration': 9999,
119 # Exact inverse on coarse grid
120 'mg_coarse_ksp_type': 'preonly',
121 'mg_coarse_pc_type': 'lu',
122 })
124class ILU(precond):
125 def __init__(self, A, parameters=None, pdes=1, nullspace=None, prefix=None):
126 supports_mpi(False, 'PETSc ILU does not work in parallel', mat(A).comm.size)
127 super().__init__(A, PETSc.PC.Type.ILU,
128 pdes=pdes, nullspace=nullspace, options=parameters, prefix=prefix)
130class Cholesky(precond):
131 def __init__(self, A, parameters=None, prefix=None):
132 super().__init__(A, PETSc.PC.Type.CHOLESKY,
133 pdes=1, nullspace=None, options=parameters, prefix=prefix)
135class LU(precond):
136 def __init__(self, A, parameters=None, prefix=None, defaults={}):
137 super().__init__(A, PETSc.PC.Type.LU,
138 pdes=1, nullspace=None, options=parameters, prefix=prefix,
139 defaults=defaults)
142class MUMPS_LU(LU):
143 def __init__(self, A, parameters=None, prefix=None):
144 super().__init__(A,
145 parameters=parameters, prefix=prefix,
146 defaults={
147 'pc_factor_mat_solver_package': 'mumps',
148 })
151class SUPERLU_LU(LU):
152 def __init__(self, A, parameters=None, prefix=None):
153 super().__init__(A,
154 parameters=parameters, prefix=prefix,
155 defaults={
156 'pc_factor_mat_solver_package': 'superlu',
157 })
158 self.petsc_prec.setFactorPivot(1E-16)
160class AMG(precond):
161 """
162 BoomerAMG preconditioner from the Hypre Library
163 """
164 def __init__(self, A, parameters=None, pdes=1, nullspace=None, prefix=None):
165 super().__init__(A, PETSc.PC.Type.HYPRE,
166 pdes=pdes, nullspace=nullspace, options=parameters, prefix=prefix,
167 defaults={
168 "pc_hypre_type": "boomeramg",
169 #"pc_hypre_boomeramg_cycle_type": "V", # (V,W)
170 #"pc_hypre_boomeramg_max_levels": 25,
171 #"pc_hypre_boomeramg_max_iter": 1,
172 #"pc_hypre_boomeramg_tol": 0,
173 #"pc_hypre_boomeramg_truncfactor" : 0, # Truncation factor for interpolation
174 #"pc_hypre_boomeramg_P_max": 0, # Max elements per row for interpolation
175 #"pc_hypre_boomeramg_agg_nl": 0, # Number of levels of aggressive coarsening
176 #"pc_hypre_boomeramg_agg_num_paths": 1, # Number of paths for aggressive coarsening
177 #"pc_hypre_boomeramg_strong_threshold": .25,# Threshold for being strongly connected
178 #"pc_hypre_boomeramg_max_row_sum": 0.9,
179 #"pc_hypre_boomeramg_grid_sweeps_all": 1, # Number of sweeps for the up and down grid levels
180 #"pc_hypre_boomeramg_grid_sweeps_down": 1,
181 #"pc_hypre_boomeramg_grid_sweeps_up":1,
182 #"pc_hypre_boomeramg_grid_sweeps_coarse": 1,# Number of sweeps for the coarse level (None)
183 #"pc_hypre_boomeramg_relax_type_all": "symmetric-SOR/Jacobi", # (Jacobi, sequential-Gauss-Seidel, seqboundary-Gauss-Seidel,
184 # SOR/Jacobi, backward-SOR/Jacobi, symmetric-SOR/Jacobi,
185 # l1scaled-SOR/Jacobi Gaussian-elimination, CG, Chebyshev,
186 # FCF-Jacobi, l1scaled-Jacobi)
187 #"pc_hypre_boomeramg_relax_type_down": "symmetric-SOR/Jacobi",
188 #"pc_hypre_boomeramg_relax_type_up": "symmetric-SOR/Jacobi",
189 #"pc_hypre_boomeramg_relax_type_coarse": "Gaussian-elimination",
190 #"pc_hypre_boomeramg_relax_weight_all": 1, # Relaxation weight for all levels (0 = hypre estimates, -k = determined with k CG steps)
191 #"pc_hypre_boomeramg_relax_weight_level": (1,1), # Set the relaxation weight for a particular level
192 #"pc_hypre_boomeramg_outer_relax_weight_all": 1,
193 #"pc_hypre_boomeramg_outer_relax_weight_level": (1,1),
194 #"pc_hypre_boomeramg_no_CF": "", # Do not use CF-relaxation
195 #"pc_hypre_boomeramg_measure_type": "local", # (local global)
196 #"pc_hypre_boomeramg_coarsen_type": "Falgout", # (Ruge-Stueben, modifiedRuge-Stueben, Falgout, PMIS, HMIS)
197 #"pc_hypre_boomeramg_interp_type": "classical",# (classical, direct, multipass, multipass-wts, ext+i, ext+i-cc, standard, standard-wts, FF, FF1)
198 #"pc_hypre_boomeramg_print_statistics": "",
199 #"pc_hypre_boomeramg_print_debug": "",
200 #"pc_hypre_boomeramg_nodal_coarsen": "",
201 #"pc_hypre_boomeramg_nodal_relaxation": "",
202 })
205class SOR(precond):
206 def __init__(self, A, parameters=None, pdes=1, nullspace=None, prefix=None):
207 options = {
208 }
209 super().__init__(A, PETSc.PC.Type.SOR,
210 pdes=pdes, nullspace=nullspace, options=parameters, prefix=prefix,
211 defaults={
212 "pc_sor_omega": 1, # relaxation factor (0 < omega < 2, 1 is Gauss-Seidel)
213 "pc_sor_its": 1, # number of inner SOR iterations
214 "pc_sor_lits": 1, # number of local inner SOR iterations
215 "pc_sor_symmetric": "", # for SSOR
216 #"pc_sor_backward": "",
217 #"pc_sor_forward": "",
218 #"tmp_pc_sor_local_symmetric": "", # use SSOR separately on each processor
219 #"tmp_pc_sor_local_backward": "",
220 #"tmp_pc_sor_local_forward": "",
221 })
222 supports_mpi(not self.optsDB.hasName('pc_sor_symmetric'), 'PETSc symmetric SOR not supported in parallel', mat(A).comm.size)
225class Elasticity(precond):
226 def __init__(self, A, parameters=None, pdes=1, nullspace=None, prefix=None):
227 super().__init__(A, PETSc.PC.Type.GAMG,
228 pdes=pdes, nullspace=nullspace, options=parameters, prefix=prefix,
229 defaults={
230 'pc_mg_cycle_type': 'w',
231 'pc_mg_multiplicative_cycles': 2
232 })
234class ASM(precond):
235 """
236 Additive Scwharz Method.
237 Defaults (or should default, not tested) to one block per process.
238 """
239 def __init__(self, A, parameters=None, pdes=1, nullspace=None, prefix=None):
240 super().__init__(A, PETSc.PC.Type.ASM,
241 pdes=pdes, nullspace=nullspace, options=parameters, prefix=prefix,
242 defaults={
243 #"pc_asm_blocks": 1, # Number of subdomains
244 "pc_asm_overlap": 1, # Number of grid points overlap
245 "pc_asm_type": "RESTRICT", # (NONE, RESTRICT, INTERPOLATE, BASIC)
246 "sup_ksp_type": "preonly", # KSP solver for the subproblems
247 "sub_pc_type": "ilu" # Preconditioner for the subproblems
248 })
251class Jacobi(precond):
252 """
253 Actually this is only a diagonal scaling preconditioner; no support for relaxation or multiple iterations.
254 """
255 def __init__(self, A, parameters=None, pdes=1, nullspace=None, prefix=None):
256 super().__init__(A, PETSc.PC.Type.JACOBI,
257 pdes=pdes, nullspace=nullspace, options=parameters, prefix=prefix,
258 defaults={
259 #"pc_jacobi_rowmax": "", # Use row maximums for diagonal
260 #"pc_jacobi_rowsum": "", # Use row sums for diagonal
261 #"pc_jacobi_abs":, "", # Use absolute values of diagaonal entries
262 })
265def vec(x):
266 return df.as_backend_type(x).vec()
269def mat(A):
270 return df.as_backend_type(A).mat()
273class HypreAMS(precond):
274 '''
275 AMG auxiliary space preconditioner for H(curl) problems in 2d/3d and
276 H(div) problems in 2d. The bilinear form behind the matrix whose
277 approximate inverse is constructed is: Find u in V
279 alpha*inner(curl(u), curl(v)) + beta*inner(u, v)*dx for all v in V
281 where alpha > 0 and beta >= 0.
282 '''
283 def __init__(self, A, V, ams_zero_beta_poisson=False, prefix=None):
284 mesh = V.mesh()
285 if mesh.geometry().dim() == 2: 285 ↛ 288line 285 didn't jump to line 288, because the condition on line 285 was never false
286 assert V.ufl_element().family() in ('Raviart-Thomas', 'Nedelec 1st kind H(curl)')
287 else:
288 assert V.ufl_element().family() == 'Nedelec 1st kind H(curl)'
289 # FIXME: only tested this for lower orders
290 assert V.ufl_element().degree() == 1
292 # AMS setup requires auxiliary operators
293 Q = df.FunctionSpace(mesh, 'CG', 1)
294 G = df.DiscreteOperators.build_gradient(V, Q)
296 super().__init__(A, PETSc.PC.Type.HYPRE, V=V,
297 pdes=None, nullspace=None, options=None, prefix=prefix,
298 defaults={
299 'pc_hypre_type': 'ams',
300 })
301 self.petsc_prec.setHYPREDiscreteGradient(mat(G))
303 # Constants
304 constants = [vec(df.interpolate(df.Constant(c), V).vector())
305 for c in np.eye(mesh.geometry().dim())]
306 self.petsc_prec.setHYPRESetEdgeConstantVectors(*constants)
308 # NOTE: signal that we don't have lower order term
309 ams_zero_beta_poisson and self.petsc_prec.setHYPRESetBetaPoissonMatrix(None)
312class HypreADS(precond):
313 '''
314 AMG auxiliary space preconditioner for H(div) problems in 3d. The bilinear
315 form behind the matrix whose approximate inverse is constructed is: Find u in V
317 alpha*inner(div(u), div(v)) + beta*inner(u, v)*dx for all v in V
319 where alpha > 0 and beta >= 0.
320 '''
321 @staticmethod
322 def coordinates(mesh):
323 '''As P1 functions'''
324 V = df.FunctionSpace(mesh, 'CG', 1)
326 return V.tabulate_dof_coordinates()
328 @staticmethod
329 def discrete_gradient(mesh):
330 '''P1 to Ned1 map'''
331 Ned = df.FunctionSpace(mesh, 'Nedelec 1st kind H(curl)', 1)
332 P1 = df.FunctionSpace(mesh, 'CG', 1)
334 G = df.DiscreteOperators.build_gradient(Ned, P1)
336 return mat(G)
338 @staticmethod
339 def discrete_curl(mesh):
340 '''Ned1 to RT1 map'''
341 assert mesh.geometry().dim() == 3
342 assert mesh.topology().dim() == 3
344 Ned = df.FunctionSpace(mesh, 'Nedelec 1st kind H(curl)', 1)
345 RT = df.FunctionSpace(mesh, 'RT', 1)
347 RT_f2dof = np.array(RT.dofmap().entity_dofs(mesh, 2))
348 Ned_e2dof = np.array(Ned.dofmap().entity_dofs(mesh, 1))
350 # Alloc sparsity
351 a = df.Constant(0)*df.inner(df.TestFunction(RT), df.TrialFunction(Ned))*df.dx
352 C = df.PETScMatrix()
353 df.assemble(a, tensor=C)
355 RT_f2dof = np.array(RT.dofmap().entity_dofs(mesh, 2))
356 Ned_e2dof = np.array(Ned.dofmap().entity_dofs(mesh, 1))
358 # Facets in terms of edges
359 mesh.init(2, 1)
360 f2e = mesh.topology()(2, 1)
362 row_cols, row_values = np.zeros(3, dtype='int32'), np.array([-2, 2, -2])
363 Cmat = C.mat()
365 for facet, row in enumerate(RT_f2dof):
366 row_cols[:] = Ned_e2dof[f2e(facet)]
367 Cmat.setValues([row], row_cols, row_values, PETSc.InsertMode.INSERT_VALUES)
368 Cmat.assemble()
370 return Cmat
372 def __init__(self, A, V, prefix=None):
373 assert V.mesh().geometry().dim() == 3
374 assert V.ufl_element().family() == 'Raviart-Thomas'
375 assert V.ufl_element().degree() == 1
377 supports_mpi(False, 'ADS curl matrix needs to be fixed for parallel', mat(A).comm.size) # FIXME
379 mesh = V.mesh()
380 C = HypreADS.discrete_curl(mesh)
381 G = HypreADS.discrete_gradient(mesh)
382 xyz = HypreADS.coordinates(mesh)
384 super().__init__(A, PETSc.PC.Type.HYPRE, V=V,
385 pdes=None, nullspace=None, options=None, prefix=prefix,
386 defaults={
387 'pc_hypre_type': 'ads',
388 })
389 # Attach auxiliary operators
390 self.petsc_prec.setHYPREDiscreteGradient(G)
391 self.petsc_prec.setHYPREDiscreteCurl(G)
392 self.petsc_prec.setCoordinates(xyz)