Skip to content

Jobh/petsc solver #10

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Draft
jobh wants to merge 16 commits into
base: master
Choose a base branch
from
Draft

Jobh/petsc solver #10

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
16 commits
Select commit Hold shift + click to select a range
983cc66
Initial attempt (not working)
jobh Feb 9, 2023
ec9f266
Working, but ugly/inefficient
jobh Feb 9, 2023
7c2e4c5
Mark as MPI-not-supported
jobh Feb 9, 2023
88e1f26
Initial attempt (not working)
jobh Feb 9, 2023
607ba70
Working, but ugly/inefficient
jobh Feb 9, 2023
396861d
Mark as MPI-not-supported
jobh Feb 9, 2023
dd426e6
bring API closer to that of block.iterative
MiroK Mar 23, 2023
77ca85c
bring API closer to that of block.iterative
MiroK Mar 23, 2023
0067266
eigenvalue_estimates from petsc
MiroK Mar 25, 2023
43a2efa
Handle scalar 0 blocks not on diagonal as None
MiroK Apr 7, 2023
897bb5f
update
MiroK Aug 4, 2023
6484a27
Work around a liveness issue in fenics
jobh Apr 17, 2024
1ce1248
pass nullspace to KSP
MiroK May 8, 2024
83a4ba3
Merge branch 'jobh/petsc-solver' of github.com:blocknics/cbc.block in…
MiroK May 8, 2024
9e8d1c1
KSP with shells
MiroK Jun 25, 2024
2554c79
Vec for block vectors
MiroK Jul 11, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
1 change: 1 addition & 0 deletions block/algebraic/petsc/__init__.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -18,3 +18,4 @@ def __call__(self):

from .precond import *
from .matrix import *
from .solver import *
135 changes: 79 additions & 56 deletions block/algebraic/petsc/precond.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -11,61 +11,72 @@
from time import time
from dolfin import info


class precond(block_base):
def __init__(self, A, prectype, pdes, nullspace, prefix, options, V=None, create_cb=None, defaults={}):
Ad = mat(A)

if nullspace:
from block.block_util import isscalar
ns = PETSc.NullSpace()
if isscalar(nullspace):
ns.create(constant=True)
else:
ns.create(constant=False, vectors=[v.down_cast().vec() for v in nullspace])
try:
Ad.setNearNullSpace(ns)
except:
info('failed to set near null space (not supported in petsc4py version)')

class petsc_base(block_base):
Copy link
Collaborator

@MiroK MiroK Feb 21, 2023

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Some of the methods, e.g. LU, Cholesky miss the prefix argument in their __init__ signature.

Copy link
Collaborator Author

@jobh jobh Feb 21, 2023

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think it must be a problem introduced in PR #9, not this one? I'll look into it.

Copy link
Collaborator Author

@jobh jobh Feb 26, 2023

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Rebased on top of bugfix

@classmethod
def _merge_options(cls, prefix, options, defaults):
# The PETSc options database is global, so we create a unique prefix
# that precisely matches the supplied options. The user can override
# this prefix if necessary.
if prefix is None:
# Use class name instead of hash of defaults (we assume defaults
# are constant for a given class)
prefix = self.__class__.__name__
prefix = cls.__name__
if options:
sha = hashlib.sha256(str(tuple(sorted(options.items()))).encode())
prefix += sha.hexdigest()[:8]
prefix += ':'

self.optsDB = PETSc.Options(prefix)
optsDB = PETSc.Options(prefix)
params = defaults.copy()
if options:
params.update(options)
for key, val in params.items():
# Allow overriding defaults by setting key='del'. Command-line
# options (already in optsDB) are never overwritten.
if val!='del' and not self.optsDB.hasName(key):
self.optsDB.setValue(key, val)
if val!='del' and not optsDB.hasName(key):
optsDB.setValue(key, val)
return prefix, optsDB

def down_cast(self):
return self.petsc_op

class precond(petsc_base):
def __init__(self, A, prectype, pdes, nullspace, prefix, options, V=None, defaults={}):
Ad = mat(A)

if nullspace:
from block.block_util import isscalar
ns = PETSc.NullSpace()
if isscalar(nullspace):
ns.create(constant=True)
else:
ns.create(constant=False, vectors=[vec(v) for v in nullspace])
try:
Ad.setNearNullSpace(ns)
except:
info('failed to set near null space (not supported in petsc4py version)')

prefix, self.optsDB = self._merge_options(prefix=prefix, options=options, defaults=defaults)

if prectype == PETSc.PC.Type.HYPRE and nullspace:
if not all(self.optsDB.hasName(n) for n in ['pc_hypre_boomeramg_nodal_coarsen', 'pc_hypre_boomeramg_vec_interp_variant']):
print('Near null space is ignored by hypre UNLESS you also supply pc_hypre_boomeramg_nodal_coarsen and pc_hypre_boomeramg_vec_interp_variant')

self.A = A
self.petsc_prec = PETSc.PC().create(V.mesh().mpi_comm() if V else None)
self.petsc_prec.setOptionsPrefix(prefix)
try:
self.petsc_prec.setType(prectype)
except PETSc.Error as e:
if e.ierr == 86:
raise ValueError(f'Unknown PETSc type "{prectype}"')
else:
raise
self.petsc_prec.setOperators(Ad)
self.petsc_prec.setFromOptions()
self.petsc_op = PETSc.PC().create(V.mesh().mpi_comm() if V else None)
self.petsc_op.setOptionsPrefix(prefix)

if prectype:
try:
self.petsc_op.setType(prectype)
except PETSc.Error as e:
if e.ierr == 86:
raise ValueError(f'Unknown PETSc type "{prectype}"')
else:
raise

self.petsc_op.setOperators(Ad)
self.petsc_op.setFromOptions()

self.setup_time = -1.

Expand All @@ -74,37 +85,33 @@ def matvec(self, b):
if self.setup_time < 0:
T = time()
# Create preconditioner based on the options database
self.petsc_prec.setUp()
self.petsc_op.setUp()

setup_time = time()-T
comm = self.petsc_prec.getComm().tompi4py()
comm = self.petsc_op.getComm().tompi4py()

setup_time = comm.allreduce(setup_time)/comm.size
info('constructed %s preconditioner in %.2f s'%(self.__class__.__name__, setup_time))

self.setup_time = setup_time

from dolfin import GenericVector
if not isinstance(b, GenericVector):
from dolfin import GenericVector, PETScVector
if not isinstance(b, (PETScVector, GenericVector)):
return NotImplemented
x = self.A.create_vec(dim=1)
if len(x) != len(b):
raise RuntimeError(
'incompatible dimensions for PETSc matvec, %d != %d'%(len(x),len(b)))

self.petsc_prec.apply(b.down_cast().vec(), x.down_cast().vec())
self.petsc_op.apply(vec(b), vec(x))
return x

def down_cast(self):
return self.petsc_prec

def __str__(self):
return '<%s prec of %s>'%(self.__class__.__name__, str(self.A))

def create_vec(self, dim=1):
return self.A.create_vec(dim)


class ML(precond):
def __init__(self, A, parameters=None, pdes=1, nullspace=None, prefix=None):
super().__init__(A, PETSc.PC.Type.ML,
Expand All @@ -128,6 +135,10 @@ def __init__(self, A, parameters=None, pdes=1, nullspace=None, prefix=None):
pdes=pdes, nullspace=nullspace, options=parameters, prefix=prefix)

class Cholesky(precond):
def __init__(self, A, parameters=None, prefix=None):
super().__init__(A, PETSc.PC.Type.CHOLESKY, 1, None, options=parameters, prefix=prefix)

class LU(precond):
def __init__(self, A, parameters=None, prefix=None):
super().__init__(A, PETSc.PC.Type.CHOLESKY,
pdes=1, nullspace=None, options=parameters, prefix=prefix)
Expand All @@ -139,12 +150,21 @@ def __init__(self, A, parameters=None, prefix=None, defaults={}):
defaults=defaults)


class UMFPACK_LU(LU):
def __init__(self, A, parameters=None, prefix=None):
super().__init__(A,
parameters=parameters, prefix=prefix,
defaults={
'pc_factor_mat_solver_type': 'umfpack',
})


class MUMPS_LU(LU):
def __init__(self, A, parameters=None, prefix=None):
super().__init__(A,
parameters=parameters, prefix=prefix,
defaults={
'pc_factor_mat_solver_package': 'mumps',
'pc_factor_mat_solver_type': 'mumps',
})


Expand All @@ -153,9 +173,9 @@ def __init__(self, A, parameters=None, prefix=None):
super().__init__(A,
parameters=parameters, prefix=prefix,
defaults={
'pc_factor_mat_solver_package': 'superlu',
'pc_factor_mat_solver_type': 'superlu',
})
self.petsc_prec.setFactorPivot(1E-16)
self.petsc_op.setFactorPivot(1E-16)

class AMG(precond):
"""
Expand Down Expand Up @@ -298,17 +318,19 @@ def __init__(self, A, V, ams_zero_beta_poisson=False, prefix=None):
defaults={
'pc_hypre_type': 'ams',
})
self.petsc_prec.setHYPREDiscreteGradient(mat(G))
self.petsc_op.setHYPREDiscreteGradient(mat(G))

# Constants
constants = [vec(df.interpolate(df.Constant(c), V).vector())
for c in np.eye(mesh.geometry().dim())]
self.petsc_prec.setHYPRESetEdgeConstantVectors(*constants)
self.petsc_op.setHYPRESetEdgeConstantVectors(*constants)

# NOTE: signal that we don't have lower order term
ams_zero_beta_poisson and self.petsc_prec.setHYPRESetBetaPoissonMatrix(None)
ams_zero_beta_poisson and self.petsc_op.setHYPRESetBetaPoissonMatrix(None)

self.setup_time = -1


class HypreADS(precond):
'''
AMG auxiliary space preconditioner for H(div) problems in 3d. The bilinear
Expand Down Expand Up @@ -373,20 +395,21 @@ def __init__(self, A, V, prefix=None):
assert V.mesh().geometry().dim() == 3
assert V.ufl_element().family() == 'Raviart-Thomas'
assert V.ufl_element().degree() == 1

supports_mpi(False, 'ADS curl matrix needs to be fixed for parallel', mat(A).comm.size) # FIXME

mesh = V.mesh()
C = HypreADS.discrete_curl(mesh)
G = HypreADS.discrete_gradient(mesh)
xyz = HypreADS.coordinates(mesh)

super().__init__(A, PETSc.PC.Type.HYPRE, V=V,
pdes=None, nullspace=None, options=None, prefix=prefix,
defaults={
'pc_hypre_type': 'ads',
})
super().__init__(A, '', V=V,
pdes=None, nullspace=None, options=None, prefix=prefix)
self.petsc_op.setType('hypre')
self.petsc_op.setHYPREType('ads')
self.petsc_op.setOperators(mat(A))
# Attach auxiliary operators
self.petsc_prec.setHYPREDiscreteGradient(G)
self.petsc_prec.setHYPREDiscreteCurl(G)
self.petsc_prec.setCoordinates(xyz)
self.petsc_op.setHYPREDiscreteGradient(G)
self.petsc_op.setHYPREDiscreteCurl(G)
self.petsc_op.setCoordinates(xyz)
self.petsc_op.setUp()

Loading