From 983cc660785e52fb7cdfb49aba9a5bf502e33f5f Mon Sep 17 00:00:00 2001 From: Joachim B Haga Date: Thu, 9 Feb 2023 10:01:27 +0100 Subject: [PATCH 01/14] Initial attempt (not working) --- block/algebraic/petsc/__init__.py | 1 + block/algebraic/petsc/precond.py | 84 +++++++++++----------- block/algebraic/petsc/solver.py | 116 ++++++++++++++++++++++++++++++ demo/petsc.py | 49 +++++++++++++ 4 files changed, 210 insertions(+), 40 deletions(-) create mode 100644 block/algebraic/petsc/solver.py create mode 100644 demo/petsc.py diff --git a/block/algebraic/petsc/__init__.py b/block/algebraic/petsc/__init__.py index 98cbfce..3419ebb 100644 --- a/block/algebraic/petsc/__init__.py +++ b/block/algebraic/petsc/__init__.py @@ -18,3 +18,4 @@ def __call__(self): from .precond import * from .matrix import * +from .solver import * diff --git a/block/algebraic/petsc/precond.py b/block/algebraic/petsc/precond.py index 8a8eee3..d32cb4e 100644 --- a/block/algebraic/petsc/precond.py +++ b/block/algebraic/petsc/precond.py @@ -11,55 +11,63 @@ 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): + @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) - self.petsc_prec.setType(prectype) - 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) + self.petsc_op.setType(prectype) + self.petsc_op.setOperators(Ad) + self.petsc_op.setFromOptions() self.setup_time = -1. @@ -68,10 +76,10 @@ 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)) @@ -86,19 +94,15 @@ def matvec(self, 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, parameters, pdes, nullspace, options=parameters, prefix=prefix, @@ -142,7 +146,7 @@ def __init__(self, A, parameters=None): defaults={ 'pc_factor_mat_solver_package': 'superlu', }) - self.petsc_prec.setFactorPivot(1E-16) + self.petsc_op.setFactorPivot(1E-16) class AMG(precond): """ @@ -279,15 +283,15 @@ def __init__(self, A, V, ams_zero_beta_poisson=False): 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) class HypreADS(precond): @@ -367,6 +371,6 @@ def __init__(self, A, V): 'pc_hypre_type': 'ads', }) # 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) diff --git a/block/algebraic/petsc/solver.py b/block/algebraic/petsc/solver.py new file mode 100644 index 0000000..9a6663f --- /dev/null +++ b/block/algebraic/petsc/solver.py @@ -0,0 +1,116 @@ +from block import block_mat, block_vec +from block.block_base import block_base +from block.object_pool import vec_pool +from dolfin import GenericVector, Matrix, PETScVector +import numpy as np +from petsc4py import PETSc +from .precond import petsc_base, mat as single_mat, vec as single_vec + +def mat(A, _sizes=None): + if isinstance(A, block_mat): + # We need sizes upfront in case there are matrix-free blocks, creating + # vectors is wasteful but easy + row_sz = [(v.local_size(), v.size()) for v in A.create_vec(dim=1)] + col_sz = [(v.local_size(), v.size()) for v in A.create_vec(dim=0)] + m = PETSc.Mat().createNest([[mat(m, (row_sz[r],col_sz[c])) for (c,m) in enumerate(row)] for (r,row) in enumerate(A.blocks)]) + return m + elif isinstance(A, block_base) or np.isscalar(A): + if not _sizes: + raise ValueError('Cannot create unsized PETSc matrix') + Ad = PETSc.Mat().createPython(_sizes) + Ad.setPythonContext(wrap_mat(A)) + Ad.setUp() + return Ad + elif isinstance(A, Matrix): + return single_mat(A) + else: + raise TypeError(str(type(A))) + +def vec(v): + if isinstance(v, block_vec): + return PETSc.Vec().createNest([vec(x) for x in v]) + elif isinstance(v, GenericVector): + return single_vec(v) + else: + raise TypeError(str(type(v))) + +def Vec(v): + if isinstance(v, PETSc.Vec): + if v.type == PETSc.Vec.Type.NEST: + return block_vec([Vec(x) for x in v]) + else: + return PETScVector(v) + else: + raise TypeError(str(type(v))) + +class wrap_mat: + def __init__(self, A): + self.A = A + + def mult(self, mat, x, y): + new_y = vec(self.A * Vec(x)) + y.aypx(0, new_y) + +class wrap_precond: + def __init__(self, P): + self.P = P + + def setUp(self, pc): + #B, P = pc.getOperators() + #ctx = B.getPythonContext() + pass + + def apply(self, pc, x, y): + # x and y are petsc vectors + new_y = vec(self.P * Vec(x)) + y.aypx(0, new_y) + +class petsc_solver(petsc_base): + def __init__(self, A, precond, V, ksp_type, prefix, options, defaults={}): + self.A = A + self.Ad = mat(A) + print(f'# petsc matrix type: {self.Ad.getType()}') + if self.Ad.getType() == 'nest': + for i in range(self.A.blocks.shape[0]): + for j in range(self.A.blocks.shape[1]): + print(f'# {i,j} -> {self.Ad.getNestSubMatrix(0,0).getType()}') + + prefix, self.optsDB = self._merge_options(prefix=prefix, options=options, defaults=defaults) + + self.petsc_op = PETSc.KSP().create(V.mesh().mpi_comm() if V else None) + self.petsc_op.setOptionsPrefix(prefix) + if ksp_type is not None: + self.petsc_op.setType(ksp_type) + self.petsc_op.setOperators(self.Ad) + self.petsc_op.setFromOptions() + if precond is not None: + pc = self.petsc_op.getPC() + pc.setType(PETSc.PC.Type.PYTHON) + pc.setPythonContext(wrap_precond(precond)) + + def __str__(self): + return '<%s ksp of %s>'%(self.__class__.__name__, str(self.A)) + + def matvec(self, b): + self.petsc_op.setUp() + x = self.A.create_vec(dim=1) + print('# petsc sizes (A, b, x):', self.petsc_op.getOperators()[0].getSizes(), vec(b).getSizes(), vec(x).getSizes()) + self.petsc_op.solve(vec(b), vec(x)) + return x + + @vec_pool + def create_vec(self, dim=1): + from dolfin import PETScVector + if dim == 0: + m = self.Ad.createVecRight() + elif dim == 1: + m = self.Ad.createVecLeft() + else: + raise ValueError('dim must be 0 or 1') + return PETScVector(m) + +class KSP(petsc_solver): + def __init__(self, A, precond=None, prefix=None, **parameters): + super().__init__(A, precond=precond, V=None, ksp_type=None, prefix=prefix, options=parameters, + defaults={ + }) diff --git a/demo/petsc.py b/demo/petsc.py new file mode 100644 index 0000000..bf069ce --- /dev/null +++ b/demo/petsc.py @@ -0,0 +1,49 @@ +from block import * +from block.algebraic.hazmath import AMG +from block.algebraic.petsc import KSP +from block.dolfin_util import BoxBoundary +from block.iterative import Richardson +from dolfin import * + +######################### +# Basic Poisson (poisson.py) +######################### + +mesh = UnitSquareMesh(16,16) +V = FunctionSpace(mesh, "CG", 1) + +f = Expression("sin(pi*x[0])", degree=2) +u, v = TrialFunction(V), TestFunction(V) + +a = u*v*dx + dot(grad(u), grad(v))*dx +L = f*v*dx + +A = assemble(a) +b = assemble(L) + +u = Function(V) +solve(A, u.vector(), b) + +######################### +# Solve as single-matrix +######################### + +Ainv = KSP(A, precond=AMG(A), ksp_type='cg') +# Alternatively: +#Ainv = KSP(A, ksp_type='cg', pc_type='hypre', pc_hypre_type='boomeramg') +x = Ainv*b + +check_expected('x', x, expected=u.vector(), show=True) + +######################### +# Solve as (trivial) block-matrix +######################### + +AA = block_mat([[A]]) +bb = block_vec([b]) +AAinv = KSP(AA, ksp_type='cg', pc_type='none') + +xx = AAinv*bb # PETSC_ERR_ARG_WRONG 62 /* wrong argument (but object probably ok) */ + +check_expected('x', xx[0], expected=u2.vector(), show=True) + From ec9f266114d39fa06f8d22834c61a7440fab80c9 Mon Sep 17 00:00:00 2001 From: Joachim B Haga Date: Thu, 9 Feb 2023 14:23:54 +0100 Subject: [PATCH 02/14] Working, but ugly/inefficient --- block/algebraic/petsc/solver.py | 85 ++++++++++++++++++----------- data/regression/petsc.Sigma.pickle | Bin 0 -> 5201 bytes data/regression/petsc.U.pickle | Bin 0 -> 1816 bytes demo/petsc.py | 82 ++++++++++++++++++++++++---- 4 files changed, 125 insertions(+), 42 deletions(-) create mode 100644 data/regression/petsc.Sigma.pickle create mode 100644 data/regression/petsc.U.pickle diff --git a/block/algebraic/petsc/solver.py b/block/algebraic/petsc/solver.py index 9a6663f..7d4f5f4 100644 --- a/block/algebraic/petsc/solver.py +++ b/block/algebraic/petsc/solver.py @@ -1,24 +1,34 @@ from block import block_mat, block_vec -from block.block_base import block_base +from block.block_base import block_base, block_container from block.object_pool import vec_pool from dolfin import GenericVector, Matrix, PETScVector import numpy as np from petsc4py import PETSc from .precond import petsc_base, mat as single_mat, vec as single_vec +# VecNest is a perfect fit for block_vec, but I can't get it to work (gives +# error code 86), and I don't know how to get proper error messages or traces +# from PETSc. When set to False, it works but +# - is very inefficient (lots of copying of vectors) +# - doesn't work in parallel +USE_VECNEST=False + def mat(A, _sizes=None): if isinstance(A, block_mat): # We need sizes upfront in case there are matrix-free blocks, creating - # vectors is wasteful but easy + # a vector to find the size is wasteful but easy row_sz = [(v.local_size(), v.size()) for v in A.create_vec(dim=1)] col_sz = [(v.local_size(), v.size()) for v in A.create_vec(dim=0)] - m = PETSc.Mat().createNest([[mat(m, (row_sz[r],col_sz[c])) for (c,m) in enumerate(row)] for (r,row) in enumerate(A.blocks)]) + m = PETSc.Mat().createNest([[mat(m, (row_sz[r],col_sz[c])) for (c,m) in enumerate(row)] + for (r,row) in enumerate(A.blocks)]) + if USE_VECNEST: + m.setVecType(PETSc.Vec.Type.NEST) return m elif isinstance(A, block_base) or np.isscalar(A): if not _sizes: raise ValueError('Cannot create unsized PETSc matrix') Ad = PETSc.Mat().createPython(_sizes) - Ad.setPythonContext(wrap_mat(A)) + Ad.setPythonContext(petsc_py_wrapper(A)) Ad.setUp() return Ad elif isinstance(A, Matrix): @@ -28,52 +38,66 @@ def mat(A, _sizes=None): def vec(v): if isinstance(v, block_vec): - return PETSc.Vec().createNest([vec(x) for x in v]) + if USE_VECNEST: + return PETSc.Vec().createNest([single_vec(x) for x in v]) + else: + vecs = [vec(x) for x in v] + sizes = [x.getSizes() for x in vecs] + vs = PETSc.Vec().create(comm=vecs[0].comm) + vs.setType(PETSc.Vec.Type.STANDARD) + vs.setSizes(sum(s[1] for s in sizes)) + arr = vs.getArray() + i0 = 0 + for vv in vecs: + arr2 = vv.getArray() + arr[i0:i0+len(arr2)] = arr2 + i0 += len(arr2) + return vs elif isinstance(v, GenericVector): return single_vec(v) else: raise TypeError(str(type(v))) -def Vec(v): +def Vec(v, creator): if isinstance(v, PETSc.Vec): - if v.type == PETSc.Vec.Type.NEST: - return block_vec([Vec(x) for x in v]) + if USE_VECNEST: + if v.type == PETSc.Vec.Type.NEST: + return block_vec([Vec(x) for x in v]) + else: + return PETScVector(v) else: - return PETScVector(v) + if isinstance(creator, block_container): + ret = creator.create_vec(dim=1) + arr = v.getArray() + i0 = 0 + for vv in ret: + vv.set_local(arr[i0:i0+len(vv)]) + i0 += len(vv) + return ret + else: + return PETScVector(v) else: raise TypeError(str(type(v))) -class wrap_mat: +class petsc_py_wrapper: + # Python context for cbc.block actions (PC.Type.PYTHON, Mat.Type.PYTHON) def __init__(self, A): self.A = A def mult(self, mat, x, y): - new_y = vec(self.A * Vec(x)) + new_y = vec(self.A * Vec(x, self.A)) y.aypx(0, new_y) -class wrap_precond: - def __init__(self, P): - self.P = P + def apply(self, pc, x, y): + return self.mult(None, x, y) def setUp(self, pc): - #B, P = pc.getOperators() - #ctx = B.getPythonContext() pass - def apply(self, pc, x, y): - # x and y are petsc vectors - new_y = vec(self.P * Vec(x)) - y.aypx(0, new_y) - class petsc_solver(petsc_base): def __init__(self, A, precond, V, ksp_type, prefix, options, defaults={}): self.A = A self.Ad = mat(A) - print(f'# petsc matrix type: {self.Ad.getType()}') - if self.Ad.getType() == 'nest': - for i in range(self.A.blocks.shape[0]): - for j in range(self.A.blocks.shape[1]): - print(f'# {i,j} -> {self.Ad.getNestSubMatrix(0,0).getType()}') prefix, self.optsDB = self._merge_options(prefix=prefix, options=options, defaults=defaults) @@ -86,17 +110,16 @@ def __init__(self, A, precond, V, ksp_type, prefix, options, defaults={}): if precond is not None: pc = self.petsc_op.getPC() pc.setType(PETSc.PC.Type.PYTHON) - pc.setPythonContext(wrap_precond(precond)) + pc.setPythonContext(petsc_py_wrapper(precond)) def __str__(self): return '<%s ksp of %s>'%(self.__class__.__name__, str(self.A)) def matvec(self, b): self.petsc_op.setUp() - x = self.A.create_vec(dim=1) - print('# petsc sizes (A, b, x):', self.petsc_op.getOperators()[0].getSizes(), vec(b).getSizes(), vec(x).getSizes()) - self.petsc_op.solve(vec(b), vec(x)) - return x + x = self.Ad.createVecLeft() + self.petsc_op.solve(vec(b), x) + return Vec(x, self.A) @vec_pool def create_vec(self, dim=1): diff --git a/data/regression/petsc.Sigma.pickle b/data/regression/petsc.Sigma.pickle new file mode 100644 index 0000000000000000000000000000000000000000..9e62a78186eaa9197b17a80504e3bd45f9b5e501 GIT binary patch literal 5201 zcmX9?c|26@`##nZAzOt?C3})0p}B@6Az9*;kZf7P7-Oq282dKNShBoLB`GbI-qNO# zcu|T-DnyB-qLPs6H^2A${Bh3be4hJ!p6guCeP7r89J*k>xeW2&h+(T`3I>NzgV?+A zE%5a!2O2w5Ai~Fsohcg}9z==O_6ni;Xa|KKrTHJGQV&P7GsO;3eY`?~!)VlSFB+T0 z&J_ADke#_G*!zDWEVdPq#NNs0$xdP0v(>El_R{#1*nYy9Lf*7!iqHSQ!|!XrPD}mY zv_K}GG26z*#_|tu{CBY;Xl&0kQ;sx2ydQNv+&00%Gu`&F=QTyZ=j`sB$ZR~VY;(~w z&&LL(M$cMi0&1s}*edx+=qDu4PpTwhQRrJAzQ?%;O)Fh4j3f9k^+kMQDjuM=h;Ppi z#<8mf0=riyqTx5!G}SU+Hv{j@^p>OaDqrEO!Cf$cLVuVb$q@crf1GtMHv9wtO;?-V=o@e^=JB zy~E%cQQy7l3k8$5tDW}E2g5~UWssa@0_a6Cr%YW1p%neg%hgdC#2#(^?e-f%)Yp~X zNcloq7ynF-|3)Zmo(lDe$@*jTXtsU+o-j0Cw{6-py7%*!F=^+AxtKYn1IEEPqmw6epXrIS zopaB-$wwhL92k9X;-z>a8LHH$! zu>KWU<1CNjk17XqU3~CT$W7cT!W}C0Z=T+D2}81+4s-`x|_6a>Ywok5k=E}KJ$P?m>YPW&gQNASWj85Rsi%z=M>^RWI%aVm{Bsl z8H8)XnhA|@68Az~irZr;Si~91#TnaS@!s{GguER*UbQ-zH3VR$Az*5*F9Hi$iH);a z6tL|Qzlp97MY-?||INMbaPBDicp=Rbop)Bh*fZ%4f>68K(?{~77-Ux}%DXF!J zohK1ot?0Q!k|2JMUhMKv0a-8R;(YN+5czqx@2Sm`xXT8%L1H^4{b`yGS(X1o54PsyJmhZgg0H1No8`MfaFn0vpH+VwsxE5NP%TaHRTReXe*EPub{?5D*0X0LShtd zEK%*b_Wm%6R2to;#%&;Z;Kn0MQG0aCx~eo+>;uvH`oM?DTO^{j-Do86IS<#1rA_wz zy!%}xvT6s3-n8jOjkOHC_ZnO_IkppneE-Wxc(^7`Ze@NV5$SuW zHLBl9#4A^3#SHgxcnggd_PsC0x+-g7D2%i!TO z$(KEqh`u#ZU+D$`pdXEX@ZH7)Oumu0VS8gtG3-5@224Qs-BGEj{Fg*jKYIEZJe>07 z8(V~Vbuf7HOs<3oxZd;;@6&5RIkjx~T%!TRWJ)@QUwLA${l$j&rAhczz1K3MG!cx@ zP3ynyKLX)fx;F)uc!GRo;{q*uKjb$PPja z|IPS{AqHIAwz)9(Cqd`4tnO$e1J4Dzo8L(Jqf^nHyCZKqJT#rx1>|f3;dm}b_oD&N zS7b+imht#K!-+kh^n*m`nmlJ*7X)!*Qd7^Kr^lY*dc%vQ9w1Mamw? z_A$73pOOH*G`7O)z-%zMv;?YSBC=)1Ja=^jK+ga9!unoskR796^cAlFIkobaefJU& zTeSHK>Un(Rc9q%*FN5nfX)(Xj6DT;glnf|MH*YFMyl(&4gydp0keO@)w9feez!K57{87Td6`vWbFi4e2(!=^+M zyuTPPSz{al<&V|e5shFdt=AMtfeg5}jyb1L7J>M{zbwM-3(2E`y24C-9YSe@Vnvy;+CGLONtr6knGdI0%a+hmKz^eBXQRq?>0+MDyuf|6@x*mpeXA5nc}RJO2+q_}79Q!*EpW-h!dEo#e&U zrkMNmGvuO!Hx|x!C_eco9sxfr44PYV@M)8ix;!fv;)mnL@jexdpmyi=+<1IwdU4~^ z3mT>*h(j$dav;`zh_X)jNg|^9T|SESlE|@@!gaSd;9A&^h(Eo@vF+o1FG1-fX!r!T zxppPMLYdg=U!0EBfDwc~aF9!Q4rN^yIY)Oc|I z+MSpumX6*PZ#{FEiHPR}TZ>;}fU#zDi%?rAW=k180}fso5-k(Ht78Fj=&?V^Q)VDM zMn`Z-^UlE?cvW3x0l2>*JuX;8xrf{!ZCMd~NxRBND&-R)Ei zZ+}*kR*`@ONxc)I;^|n*)cU2lED=vH>=C!-Pr`me?DUWa4Lxse54D~PhOyco*Oi~0 zU^dYiQs?6UB54)Zv{Dp9rddK?!<0bx(N6o>whr_s0~2;#t9Uv>Up73j21@d#ipJ|U z;E|AUK@LM78CNVgdS>A` zE~@+A(m1ZU*z?qW;d4hi$O~GhqHpu+-S)ulY|(#6Tuy;^F=r{Hnvz_9JhF$hj8L+) zNCX-QL$7*H4hGWqv~~?MP`X#v?PhBl1S}5)_nb^bPYK~2=9LC{{!&$u+jOX_Z)eAq z#luaYT0`HJ0())gs?6+Q{5Jia=5f{ou}g8p5Y^+dHGE46BY;;<>k@*Z216l zj-;&k;UeqoR$TnqIPI>R?vLiF z_9f!j+y+bWp)@GHz83#mCIRP9>(1-SCLtz<^`vk-5)7G+rU!S9LA>JO27NZKpK;ea z(D2>^bWTD&HC6{h`_#849@PO6Air|{pgM><+iPX3d3kQI+NkE%BHlb-NmPjzMzswL ze3kx!VW#X>hnlWK1)^&FFVIhKtUHXIiSg?A%)p6);9vc5#SfwGzLMb!M zu5b|xe*CPJXGaU5+iaqs+miu~;fG4_1GxMoq^4-P?{DugJj2+M2A9;Y1Wj7n4IyreSB7 zfP-*DJU0Gx-@W%_5;z9Y7X&UJLl3=fg>z6mF8wf;~r*<&0rzMKDq6oaUr_o zU-sErvB8w59!kyH2IBh3()bw(Ft(eBKAX3N&2A_5Cx-+u`eV+UaWgPj^>0Dn8YXJH zr>(C}vyduxE|Hv{hG4Z;lg|E3(3ae9kzz3*neJ#M6rGOBoMX))riln?kb8fqDiu>E z3)_+x=s17I^x&7aWXMJzeP=fshZb&Lblq4yQblUg4QH&;88Kj|b8;uhW#XR8w;O`n zAvMH5xfXo(rtROxx{-bSUQApQGh7@xiCZH?SldmeY$_J+my> z95rKZ^b2r?D#towx9Uy;-06W z)LHe;zJUagpBSxA`ZpbQhVymYe;5!-%AIokmSuMVL5|R}w-z!1 zdDpFbq2@dv>N#z*QFN9xHEg!Z&w{s)xAF(A{0@Lw9{u`(e-i2+YL052%7ylR3#}zx zWuP3A|1BxYLBid^YUZ}n$oQ?jvFLjtUZ$^Y|Fn{g-9N3Qf*33e^^1qCJH|$)#EJdu zlTJeNj665Ws|ZB9RAYBP3FI|1dm6s0!2QgQ<37T?ep1l)lSXq64r7&hpJzVwwc_1% z+cWXOH&gMXQxPtyXr0gM$wG3K*ZHJvOvpNCCA3K;!{XVqm+#{;ppg8cc1)K6wW*e! z+iKDwW2do@yO;sm#<3NnH#0CbrR-QeoP_KD1o}!T=fcnZ6fA#a;YX~_+nBLP7)YGz z;y>?-p?AyPSFPCvkFmLw`fT30)NZ%8)yOh%bsO`u%QO(WZ>vv<;$gVl^M9ba^8}9S zs%1;R&4=2VIe+%K5_ITW;hMo&oU49%eYe6{TzfKR`*4~Ay@M%{ams}VZ5X|{`4Q{y%j_bau|Yg5lf3_y=i6(IdVfkALGV@Iz3*F2 z;PPsc)xi1!=HW%*!NVofmi;Lm*$b`wdtw=gC=u&gQucbyWl zBM--8s~vnQczkCN7QH=8=-hc_SGXn^IsyVIJ+Agp(l@yA^DOUtu_UzWzBJEg14!Ch z8&*Ksemc^!))8{sTi@zSh9Q7Y(8}g63mP+LC#DCAK{i|dXk$S+oSodRJYFb=eTnWv zOTG$(Qinr|e{tZ-S@`_Xq7Z^sbSbALY-l&jq-NEyu)}*$agl#MRwQ4#BjHgB=fw(& z7M8^zOx;yyh5wK`3$IQH%o5Ud>R4zrpOHZv6UyiKv%* P3i=bg!RGy`3)lWX5%LmX literal 0 HcmV?d00001 diff --git a/data/regression/petsc.U.pickle b/data/regression/petsc.U.pickle new file mode 100644 index 0000000000000000000000000000000000000000..994128d7db96389935dff014cbe5ddf0079908c7 GIT binary patch literal 1816 zcmX9;& zS&T%Ym|*LCClYMk*w8;I#e^&61mUO=LSzwM1kF{$Re6Glpt5yB6O-e^{x8uyw3*00 z^=DU$t+AEx@bGY%p27FxN=YO_vK=C*6w1|yA2aURUWv$KwYy&4lfrwK#cQXOE=Z2t^vS)h1jlab zc!IkSzN3$eWI83dx|+THv}rsJZO9s+r?eoA+<)Um=RE`!SwG$Ar9#(l1e!Jl z+q@R@AWmt_wF{4czu@c01~nIz#d}=m=1V}a&}F+fEAaAcWGHJJi3iK5d(t~v&^tJ< z;i8&@9i7M5PrS)Pab)L5k)#!`mp&iP2qE#jt!lffUXHp2wmGr8?}9V(HqLt&4@U1k z>n^a!L7_-P@?KdBZqGHnpzd}=&+$xL;VegoZfVt~UMcKPa!Q720)#dwc2N5(kvtjv zV#uHzNtXU2_ZPLnEu;3%&zoiVD$!R)Qb_cf2ucQ5iLq%cfBG+(JWv-qmd{#I2=m*M zj(&y>_-#d7uhoP08U0k{L?#L5@#T7d=twd4j~5~CQO$VbTW)A*Qw4E<#!E4^8qX~g zS@TZ{(H#{q(HJU&rF+ruf~PVpUE}uLr&R=%6GLA)i-&@(KCz|pVpL4|=G~2LgmIJa z>`0jey)2f~f1Knf-KxJ_(?E(O(-6f|;x_6cMUyvfR3P9X-E3s%b*OfmyT$5@q1W|U z@6bCE!$PY|wqg=0^QHy2VnnDI+iMwMSA#K=B}yfE4rf1|&oTK)fckgCS8+o!!O;MrNciQf@ zM1oC%LW{A`Ti}a@j&qNeLg7Z|zxL zz5_o-@;`Hasz%b6RkW19a?!dmdqm4t3~I*xN7Nn#dn5KZO=!m+hn#DV`;-WJqfT!uP=H&(;?nX&*la~#*nIUmMzohJ zai9_t0~Xrqi)~nwT%vXkk>m1FW-(nu1<7y2fp2I0U3f6mZ2HhGIHViZwHKG-?_nnA ziBvvpD<~JrY@`?r%Q5+6phQ6LLep1fB;t)}BM%m~q9~xb)yc043$g~kTu^ag5IA7y zIVnQ!{E1K9xpMI4808kemm?{&YhHtm7_~Y_vUB&>qW5Bwu}4H1g1wi8`=tw@nHI6o zIJ{#f$K%3E{h3_ztLHg-O7L;bBl}}heE4Mt_WW!8IjGNlx!uLD!{IX@PPrSlgXVY6 z(?UxQOVP8iCF~CD7@5#~e54gCy^aHKE`xPg*L1|N8kFDH?mm?%#I`*_pTmP?__i(m zZ;PMFaJj~$??<5sq7QMLvN=2qXqs;fC?s&H|B%#Wx(OM+PfYdvq;Pv^u|MP(iSuSe zvLLM;vV51keSAKM?&_t@W#>`yYM7B!S_5j!D%&ynT`&qZXeJes7!^qN-kp@ffl<|1 z#A|_F^;NpT&dYFHlYPv&)DxOOG7J%YGHXDAj*` zKfMOyp6vF~!)y$Be>)JU(S-ZXv}w6S z^B{1*X>QIYz3cE(^xn*3OR&=Xb826c0ws&yPmIpY+0s=V#gVn3G;zo=ua$$PoR4v< z4;DckJg~w4j0lO|T2t3ulvt~#J7oka;Wgisthc?3Uk>KJ7;3x0_KGX zl?#avH!CukOYXyQc;AUnWG7DX*QNf~Qi!208;0rHEYuD-KPnjB2V&;qu&n Date: Thu, 9 Feb 2023 14:56:21 +0100 Subject: [PATCH 03/14] Mark as MPI-not-supported --- block/algebraic/petsc/solver.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/block/algebraic/petsc/solver.py b/block/algebraic/petsc/solver.py index 7d4f5f4..ebfa75c 100644 --- a/block/algebraic/petsc/solver.py +++ b/block/algebraic/petsc/solver.py @@ -1,4 +1,4 @@ -from block import block_mat, block_vec +from block import block_mat, block_vec, supports_mpi from block.block_base import block_base, block_container from block.object_pool import vec_pool from dolfin import GenericVector, Matrix, PETScVector @@ -96,6 +96,8 @@ def setUp(self, pc): class petsc_solver(petsc_base): def __init__(self, A, precond, V, ksp_type, prefix, options, defaults={}): + supports_mpi(False, 'PETSc solver interface currently does not work in parallel') + self.A = A self.Ad = mat(A) From 88e1f26998b17475ca91bf32eb60c71439662929 Mon Sep 17 00:00:00 2001 From: Joachim B Haga Date: Thu, 9 Feb 2023 10:01:27 +0100 Subject: [PATCH 04/14] Initial attempt (not working) --- block/algebraic/petsc/__init__.py | 1 + block/algebraic/petsc/precond.py | 84 +++++++++++----------- block/algebraic/petsc/solver.py | 116 ++++++++++++++++++++++++++++++ demo/petsc.py | 49 +++++++++++++ 4 files changed, 210 insertions(+), 40 deletions(-) create mode 100644 block/algebraic/petsc/solver.py create mode 100644 demo/petsc.py diff --git a/block/algebraic/petsc/__init__.py b/block/algebraic/petsc/__init__.py index 98cbfce..3419ebb 100644 --- a/block/algebraic/petsc/__init__.py +++ b/block/algebraic/petsc/__init__.py @@ -18,3 +18,4 @@ def __call__(self): from .precond import * from .matrix import * +from .solver import * diff --git a/block/algebraic/petsc/precond.py b/block/algebraic/petsc/precond.py index 873cbc7..d2232e7 100644 --- a/block/algebraic/petsc/precond.py +++ b/block/algebraic/petsc/precond.py @@ -11,61 +11,69 @@ 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): + @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) + self.petsc_op = PETSc.PC().create(V.mesh().mpi_comm() if V else None) + self.petsc_op.setOptionsPrefix(prefix) try: - self.petsc_prec.setType(prectype) + 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_prec.setOperators(Ad) - self.petsc_prec.setFromOptions() + self.petsc_op.setOperators(Ad) + self.petsc_op.setFromOptions() self.setup_time = -1. @@ -74,10 +82,10 @@ 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)) @@ -92,19 +100,15 @@ def matvec(self, 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, @@ -155,7 +159,7 @@ def __init__(self, A, parameters=None, prefix=None): defaults={ 'pc_factor_mat_solver_package': 'superlu', }) - self.petsc_prec.setFactorPivot(1E-16) + self.petsc_op.setFactorPivot(1E-16) class AMG(precond): """ @@ -298,15 +302,15 @@ 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) class HypreADS(precond): @@ -387,6 +391,6 @@ def __init__(self, A, V, prefix=None): 'pc_hypre_type': 'ads', }) # 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) diff --git a/block/algebraic/petsc/solver.py b/block/algebraic/petsc/solver.py new file mode 100644 index 0000000..9a6663f --- /dev/null +++ b/block/algebraic/petsc/solver.py @@ -0,0 +1,116 @@ +from block import block_mat, block_vec +from block.block_base import block_base +from block.object_pool import vec_pool +from dolfin import GenericVector, Matrix, PETScVector +import numpy as np +from petsc4py import PETSc +from .precond import petsc_base, mat as single_mat, vec as single_vec + +def mat(A, _sizes=None): + if isinstance(A, block_mat): + # We need sizes upfront in case there are matrix-free blocks, creating + # vectors is wasteful but easy + row_sz = [(v.local_size(), v.size()) for v in A.create_vec(dim=1)] + col_sz = [(v.local_size(), v.size()) for v in A.create_vec(dim=0)] + m = PETSc.Mat().createNest([[mat(m, (row_sz[r],col_sz[c])) for (c,m) in enumerate(row)] for (r,row) in enumerate(A.blocks)]) + return m + elif isinstance(A, block_base) or np.isscalar(A): + if not _sizes: + raise ValueError('Cannot create unsized PETSc matrix') + Ad = PETSc.Mat().createPython(_sizes) + Ad.setPythonContext(wrap_mat(A)) + Ad.setUp() + return Ad + elif isinstance(A, Matrix): + return single_mat(A) + else: + raise TypeError(str(type(A))) + +def vec(v): + if isinstance(v, block_vec): + return PETSc.Vec().createNest([vec(x) for x in v]) + elif isinstance(v, GenericVector): + return single_vec(v) + else: + raise TypeError(str(type(v))) + +def Vec(v): + if isinstance(v, PETSc.Vec): + if v.type == PETSc.Vec.Type.NEST: + return block_vec([Vec(x) for x in v]) + else: + return PETScVector(v) + else: + raise TypeError(str(type(v))) + +class wrap_mat: + def __init__(self, A): + self.A = A + + def mult(self, mat, x, y): + new_y = vec(self.A * Vec(x)) + y.aypx(0, new_y) + +class wrap_precond: + def __init__(self, P): + self.P = P + + def setUp(self, pc): + #B, P = pc.getOperators() + #ctx = B.getPythonContext() + pass + + def apply(self, pc, x, y): + # x and y are petsc vectors + new_y = vec(self.P * Vec(x)) + y.aypx(0, new_y) + +class petsc_solver(petsc_base): + def __init__(self, A, precond, V, ksp_type, prefix, options, defaults={}): + self.A = A + self.Ad = mat(A) + print(f'# petsc matrix type: {self.Ad.getType()}') + if self.Ad.getType() == 'nest': + for i in range(self.A.blocks.shape[0]): + for j in range(self.A.blocks.shape[1]): + print(f'# {i,j} -> {self.Ad.getNestSubMatrix(0,0).getType()}') + + prefix, self.optsDB = self._merge_options(prefix=prefix, options=options, defaults=defaults) + + self.petsc_op = PETSc.KSP().create(V.mesh().mpi_comm() if V else None) + self.petsc_op.setOptionsPrefix(prefix) + if ksp_type is not None: + self.petsc_op.setType(ksp_type) + self.petsc_op.setOperators(self.Ad) + self.petsc_op.setFromOptions() + if precond is not None: + pc = self.petsc_op.getPC() + pc.setType(PETSc.PC.Type.PYTHON) + pc.setPythonContext(wrap_precond(precond)) + + def __str__(self): + return '<%s ksp of %s>'%(self.__class__.__name__, str(self.A)) + + def matvec(self, b): + self.petsc_op.setUp() + x = self.A.create_vec(dim=1) + print('# petsc sizes (A, b, x):', self.petsc_op.getOperators()[0].getSizes(), vec(b).getSizes(), vec(x).getSizes()) + self.petsc_op.solve(vec(b), vec(x)) + return x + + @vec_pool + def create_vec(self, dim=1): + from dolfin import PETScVector + if dim == 0: + m = self.Ad.createVecRight() + elif dim == 1: + m = self.Ad.createVecLeft() + else: + raise ValueError('dim must be 0 or 1') + return PETScVector(m) + +class KSP(petsc_solver): + def __init__(self, A, precond=None, prefix=None, **parameters): + super().__init__(A, precond=precond, V=None, ksp_type=None, prefix=prefix, options=parameters, + defaults={ + }) diff --git a/demo/petsc.py b/demo/petsc.py new file mode 100644 index 0000000..bf069ce --- /dev/null +++ b/demo/petsc.py @@ -0,0 +1,49 @@ +from block import * +from block.algebraic.hazmath import AMG +from block.algebraic.petsc import KSP +from block.dolfin_util import BoxBoundary +from block.iterative import Richardson +from dolfin import * + +######################### +# Basic Poisson (poisson.py) +######################### + +mesh = UnitSquareMesh(16,16) +V = FunctionSpace(mesh, "CG", 1) + +f = Expression("sin(pi*x[0])", degree=2) +u, v = TrialFunction(V), TestFunction(V) + +a = u*v*dx + dot(grad(u), grad(v))*dx +L = f*v*dx + +A = assemble(a) +b = assemble(L) + +u = Function(V) +solve(A, u.vector(), b) + +######################### +# Solve as single-matrix +######################### + +Ainv = KSP(A, precond=AMG(A), ksp_type='cg') +# Alternatively: +#Ainv = KSP(A, ksp_type='cg', pc_type='hypre', pc_hypre_type='boomeramg') +x = Ainv*b + +check_expected('x', x, expected=u.vector(), show=True) + +######################### +# Solve as (trivial) block-matrix +######################### + +AA = block_mat([[A]]) +bb = block_vec([b]) +AAinv = KSP(AA, ksp_type='cg', pc_type='none') + +xx = AAinv*bb # PETSC_ERR_ARG_WRONG 62 /* wrong argument (but object probably ok) */ + +check_expected('x', xx[0], expected=u2.vector(), show=True) + From 607ba70158cb14bce0a5e5fc4c14143e2df83c9f Mon Sep 17 00:00:00 2001 From: Joachim B Haga Date: Thu, 9 Feb 2023 14:23:54 +0100 Subject: [PATCH 05/14] Working, but ugly/inefficient --- block/algebraic/petsc/solver.py | 85 ++++++++++++++++++----------- data/regression/petsc.Sigma.pickle | Bin 0 -> 5201 bytes data/regression/petsc.U.pickle | Bin 0 -> 1816 bytes demo/petsc.py | 82 ++++++++++++++++++++++++---- 4 files changed, 125 insertions(+), 42 deletions(-) create mode 100644 data/regression/petsc.Sigma.pickle create mode 100644 data/regression/petsc.U.pickle diff --git a/block/algebraic/petsc/solver.py b/block/algebraic/petsc/solver.py index 9a6663f..7d4f5f4 100644 --- a/block/algebraic/petsc/solver.py +++ b/block/algebraic/petsc/solver.py @@ -1,24 +1,34 @@ from block import block_mat, block_vec -from block.block_base import block_base +from block.block_base import block_base, block_container from block.object_pool import vec_pool from dolfin import GenericVector, Matrix, PETScVector import numpy as np from petsc4py import PETSc from .precond import petsc_base, mat as single_mat, vec as single_vec +# VecNest is a perfect fit for block_vec, but I can't get it to work (gives +# error code 86), and I don't know how to get proper error messages or traces +# from PETSc. When set to False, it works but +# - is very inefficient (lots of copying of vectors) +# - doesn't work in parallel +USE_VECNEST=False + def mat(A, _sizes=None): if isinstance(A, block_mat): # We need sizes upfront in case there are matrix-free blocks, creating - # vectors is wasteful but easy + # a vector to find the size is wasteful but easy row_sz = [(v.local_size(), v.size()) for v in A.create_vec(dim=1)] col_sz = [(v.local_size(), v.size()) for v in A.create_vec(dim=0)] - m = PETSc.Mat().createNest([[mat(m, (row_sz[r],col_sz[c])) for (c,m) in enumerate(row)] for (r,row) in enumerate(A.blocks)]) + m = PETSc.Mat().createNest([[mat(m, (row_sz[r],col_sz[c])) for (c,m) in enumerate(row)] + for (r,row) in enumerate(A.blocks)]) + if USE_VECNEST: + m.setVecType(PETSc.Vec.Type.NEST) return m elif isinstance(A, block_base) or np.isscalar(A): if not _sizes: raise ValueError('Cannot create unsized PETSc matrix') Ad = PETSc.Mat().createPython(_sizes) - Ad.setPythonContext(wrap_mat(A)) + Ad.setPythonContext(petsc_py_wrapper(A)) Ad.setUp() return Ad elif isinstance(A, Matrix): @@ -28,52 +38,66 @@ def mat(A, _sizes=None): def vec(v): if isinstance(v, block_vec): - return PETSc.Vec().createNest([vec(x) for x in v]) + if USE_VECNEST: + return PETSc.Vec().createNest([single_vec(x) for x in v]) + else: + vecs = [vec(x) for x in v] + sizes = [x.getSizes() for x in vecs] + vs = PETSc.Vec().create(comm=vecs[0].comm) + vs.setType(PETSc.Vec.Type.STANDARD) + vs.setSizes(sum(s[1] for s in sizes)) + arr = vs.getArray() + i0 = 0 + for vv in vecs: + arr2 = vv.getArray() + arr[i0:i0+len(arr2)] = arr2 + i0 += len(arr2) + return vs elif isinstance(v, GenericVector): return single_vec(v) else: raise TypeError(str(type(v))) -def Vec(v): +def Vec(v, creator): if isinstance(v, PETSc.Vec): - if v.type == PETSc.Vec.Type.NEST: - return block_vec([Vec(x) for x in v]) + if USE_VECNEST: + if v.type == PETSc.Vec.Type.NEST: + return block_vec([Vec(x) for x in v]) + else: + return PETScVector(v) else: - return PETScVector(v) + if isinstance(creator, block_container): + ret = creator.create_vec(dim=1) + arr = v.getArray() + i0 = 0 + for vv in ret: + vv.set_local(arr[i0:i0+len(vv)]) + i0 += len(vv) + return ret + else: + return PETScVector(v) else: raise TypeError(str(type(v))) -class wrap_mat: +class petsc_py_wrapper: + # Python context for cbc.block actions (PC.Type.PYTHON, Mat.Type.PYTHON) def __init__(self, A): self.A = A def mult(self, mat, x, y): - new_y = vec(self.A * Vec(x)) + new_y = vec(self.A * Vec(x, self.A)) y.aypx(0, new_y) -class wrap_precond: - def __init__(self, P): - self.P = P + def apply(self, pc, x, y): + return self.mult(None, x, y) def setUp(self, pc): - #B, P = pc.getOperators() - #ctx = B.getPythonContext() pass - def apply(self, pc, x, y): - # x and y are petsc vectors - new_y = vec(self.P * Vec(x)) - y.aypx(0, new_y) - class petsc_solver(petsc_base): def __init__(self, A, precond, V, ksp_type, prefix, options, defaults={}): self.A = A self.Ad = mat(A) - print(f'# petsc matrix type: {self.Ad.getType()}') - if self.Ad.getType() == 'nest': - for i in range(self.A.blocks.shape[0]): - for j in range(self.A.blocks.shape[1]): - print(f'# {i,j} -> {self.Ad.getNestSubMatrix(0,0).getType()}') prefix, self.optsDB = self._merge_options(prefix=prefix, options=options, defaults=defaults) @@ -86,17 +110,16 @@ def __init__(self, A, precond, V, ksp_type, prefix, options, defaults={}): if precond is not None: pc = self.petsc_op.getPC() pc.setType(PETSc.PC.Type.PYTHON) - pc.setPythonContext(wrap_precond(precond)) + pc.setPythonContext(petsc_py_wrapper(precond)) def __str__(self): return '<%s ksp of %s>'%(self.__class__.__name__, str(self.A)) def matvec(self, b): self.petsc_op.setUp() - x = self.A.create_vec(dim=1) - print('# petsc sizes (A, b, x):', self.petsc_op.getOperators()[0].getSizes(), vec(b).getSizes(), vec(x).getSizes()) - self.petsc_op.solve(vec(b), vec(x)) - return x + x = self.Ad.createVecLeft() + self.petsc_op.solve(vec(b), x) + return Vec(x, self.A) @vec_pool def create_vec(self, dim=1): diff --git a/data/regression/petsc.Sigma.pickle b/data/regression/petsc.Sigma.pickle new file mode 100644 index 0000000000000000000000000000000000000000..9e62a78186eaa9197b17a80504e3bd45f9b5e501 GIT binary patch literal 5201 zcmX9?c|26@`##nZAzOt?C3})0p}B@6Az9*;kZf7P7-Oq282dKNShBoLB`GbI-qNO# zcu|T-DnyB-qLPs6H^2A${Bh3be4hJ!p6guCeP7r89J*k>xeW2&h+(T`3I>NzgV?+A zE%5a!2O2w5Ai~Fsohcg}9z==O_6ni;Xa|KKrTHJGQV&P7GsO;3eY`?~!)VlSFB+T0 z&J_ADke#_G*!zDWEVdPq#NNs0$xdP0v(>El_R{#1*nYy9Lf*7!iqHSQ!|!XrPD}mY zv_K}GG26z*#_|tu{CBY;Xl&0kQ;sx2ydQNv+&00%Gu`&F=QTyZ=j`sB$ZR~VY;(~w z&&LL(M$cMi0&1s}*edx+=qDu4PpTwhQRrJAzQ?%;O)Fh4j3f9k^+kMQDjuM=h;Ppi z#<8mf0=riyqTx5!G}SU+Hv{j@^p>OaDqrEO!Cf$cLVuVb$q@crf1GtMHv9wtO;?-V=o@e^=JB zy~E%cQQy7l3k8$5tDW}E2g5~UWssa@0_a6Cr%YW1p%neg%hgdC#2#(^?e-f%)Yp~X zNcloq7ynF-|3)Zmo(lDe$@*jTXtsU+o-j0Cw{6-py7%*!F=^+AxtKYn1IEEPqmw6epXrIS zopaB-$wwhL92k9X;-z>a8LHH$! zu>KWU<1CNjk17XqU3~CT$W7cT!W}C0Z=T+D2}81+4s-`x|_6a>Ywok5k=E}KJ$P?m>YPW&gQNASWj85Rsi%z=M>^RWI%aVm{Bsl z8H8)XnhA|@68Az~irZr;Si~91#TnaS@!s{GguER*UbQ-zH3VR$Az*5*F9Hi$iH);a z6tL|Qzlp97MY-?||INMbaPBDicp=Rbop)Bh*fZ%4f>68K(?{~77-Ux}%DXF!J zohK1ot?0Q!k|2JMUhMKv0a-8R;(YN+5czqx@2Sm`xXT8%L1H^4{b`yGS(X1o54PsyJmhZgg0H1No8`MfaFn0vpH+VwsxE5NP%TaHRTReXe*EPub{?5D*0X0LShtd zEK%*b_Wm%6R2to;#%&;Z;Kn0MQG0aCx~eo+>;uvH`oM?DTO^{j-Do86IS<#1rA_wz zy!%}xvT6s3-n8jOjkOHC_ZnO_IkppneE-Wxc(^7`Ze@NV5$SuW zHLBl9#4A^3#SHgxcnggd_PsC0x+-g7D2%i!TO z$(KEqh`u#ZU+D$`pdXEX@ZH7)Oumu0VS8gtG3-5@224Qs-BGEj{Fg*jKYIEZJe>07 z8(V~Vbuf7HOs<3oxZd;;@6&5RIkjx~T%!TRWJ)@QUwLA${l$j&rAhczz1K3MG!cx@ zP3ynyKLX)fx;F)uc!GRo;{q*uKjb$PPja z|IPS{AqHIAwz)9(Cqd`4tnO$e1J4Dzo8L(Jqf^nHyCZKqJT#rx1>|f3;dm}b_oD&N zS7b+imht#K!-+kh^n*m`nmlJ*7X)!*Qd7^Kr^lY*dc%vQ9w1Mamw? z_A$73pOOH*G`7O)z-%zMv;?YSBC=)1Ja=^jK+ga9!unoskR796^cAlFIkobaefJU& zTeSHK>Un(Rc9q%*FN5nfX)(Xj6DT;glnf|MH*YFMyl(&4gydp0keO@)w9feez!K57{87Td6`vWbFi4e2(!=^+M zyuTPPSz{al<&V|e5shFdt=AMtfeg5}jyb1L7J>M{zbwM-3(2E`y24C-9YSe@Vnvy;+CGLONtr6knGdI0%a+hmKz^eBXQRq?>0+MDyuf|6@x*mpeXA5nc}RJO2+q_}79Q!*EpW-h!dEo#e&U zrkMNmGvuO!Hx|x!C_eco9sxfr44PYV@M)8ix;!fv;)mnL@jexdpmyi=+<1IwdU4~^ z3mT>*h(j$dav;`zh_X)jNg|^9T|SESlE|@@!gaSd;9A&^h(Eo@vF+o1FG1-fX!r!T zxppPMLYdg=U!0EBfDwc~aF9!Q4rN^yIY)Oc|I z+MSpumX6*PZ#{FEiHPR}TZ>;}fU#zDi%?rAW=k180}fso5-k(Ht78Fj=&?V^Q)VDM zMn`Z-^UlE?cvW3x0l2>*JuX;8xrf{!ZCMd~NxRBND&-R)Ei zZ+}*kR*`@ONxc)I;^|n*)cU2lED=vH>=C!-Pr`me?DUWa4Lxse54D~PhOyco*Oi~0 zU^dYiQs?6UB54)Zv{Dp9rddK?!<0bx(N6o>whr_s0~2;#t9Uv>Up73j21@d#ipJ|U z;E|AUK@LM78CNVgdS>A` zE~@+A(m1ZU*z?qW;d4hi$O~GhqHpu+-S)ulY|(#6Tuy;^F=r{Hnvz_9JhF$hj8L+) zNCX-QL$7*H4hGWqv~~?MP`X#v?PhBl1S}5)_nb^bPYK~2=9LC{{!&$u+jOX_Z)eAq z#luaYT0`HJ0())gs?6+Q{5Jia=5f{ou}g8p5Y^+dHGE46BY;;<>k@*Z216l zj-;&k;UeqoR$TnqIPI>R?vLiF z_9f!j+y+bWp)@GHz83#mCIRP9>(1-SCLtz<^`vk-5)7G+rU!S9LA>JO27NZKpK;ea z(D2>^bWTD&HC6{h`_#849@PO6Air|{pgM><+iPX3d3kQI+NkE%BHlb-NmPjzMzswL ze3kx!VW#X>hnlWK1)^&FFVIhKtUHXIiSg?A%)p6);9vc5#SfwGzLMb!M zu5b|xe*CPJXGaU5+iaqs+miu~;fG4_1GxMoq^4-P?{DugJj2+M2A9;Y1Wj7n4IyreSB7 zfP-*DJU0Gx-@W%_5;z9Y7X&UJLl3=fg>z6mF8wf;~r*<&0rzMKDq6oaUr_o zU-sErvB8w59!kyH2IBh3()bw(Ft(eBKAX3N&2A_5Cx-+u`eV+UaWgPj^>0Dn8YXJH zr>(C}vyduxE|Hv{hG4Z;lg|E3(3ae9kzz3*neJ#M6rGOBoMX))riln?kb8fqDiu>E z3)_+x=s17I^x&7aWXMJzeP=fshZb&Lblq4yQblUg4QH&;88Kj|b8;uhW#XR8w;O`n zAvMH5xfXo(rtROxx{-bSUQApQGh7@xiCZH?SldmeY$_J+my> z95rKZ^b2r?D#towx9Uy;-06W z)LHe;zJUagpBSxA`ZpbQhVymYe;5!-%AIokmSuMVL5|R}w-z!1 zdDpFbq2@dv>N#z*QFN9xHEg!Z&w{s)xAF(A{0@Lw9{u`(e-i2+YL052%7ylR3#}zx zWuP3A|1BxYLBid^YUZ}n$oQ?jvFLjtUZ$^Y|Fn{g-9N3Qf*33e^^1qCJH|$)#EJdu zlTJeNj665Ws|ZB9RAYBP3FI|1dm6s0!2QgQ<37T?ep1l)lSXq64r7&hpJzVwwc_1% z+cWXOH&gMXQxPtyXr0gM$wG3K*ZHJvOvpNCCA3K;!{XVqm+#{;ppg8cc1)K6wW*e! z+iKDwW2do@yO;sm#<3NnH#0CbrR-QeoP_KD1o}!T=fcnZ6fA#a;YX~_+nBLP7)YGz z;y>?-p?AyPSFPCvkFmLw`fT30)NZ%8)yOh%bsO`u%QO(WZ>vv<;$gVl^M9ba^8}9S zs%1;R&4=2VIe+%K5_ITW;hMo&oU49%eYe6{TzfKR`*4~Ay@M%{ams}VZ5X|{`4Q{y%j_bau|Yg5lf3_y=i6(IdVfkALGV@Iz3*F2 z;PPsc)xi1!=HW%*!NVofmi;Lm*$b`wdtw=gC=u&gQucbyWl zBM--8s~vnQczkCN7QH=8=-hc_SGXn^IsyVIJ+Agp(l@yA^DOUtu_UzWzBJEg14!Ch z8&*Ksemc^!))8{sTi@zSh9Q7Y(8}g63mP+LC#DCAK{i|dXk$S+oSodRJYFb=eTnWv zOTG$(Qinr|e{tZ-S@`_Xq7Z^sbSbALY-l&jq-NEyu)}*$agl#MRwQ4#BjHgB=fw(& z7M8^zOx;yyh5wK`3$IQH%o5Ud>R4zrpOHZv6UyiKv%* P3i=bg!RGy`3)lWX5%LmX literal 0 HcmV?d00001 diff --git a/data/regression/petsc.U.pickle b/data/regression/petsc.U.pickle new file mode 100644 index 0000000000000000000000000000000000000000..994128d7db96389935dff014cbe5ddf0079908c7 GIT binary patch literal 1816 zcmX9;& zS&T%Ym|*LCClYMk*w8;I#e^&61mUO=LSzwM1kF{$Re6Glpt5yB6O-e^{x8uyw3*00 z^=DU$t+AEx@bGY%p27FxN=YO_vK=C*6w1|yA2aURUWv$KwYy&4lfrwK#cQXOE=Z2t^vS)h1jlab zc!IkSzN3$eWI83dx|+THv}rsJZO9s+r?eoA+<)Um=RE`!SwG$Ar9#(l1e!Jl z+q@R@AWmt_wF{4czu@c01~nIz#d}=m=1V}a&}F+fEAaAcWGHJJi3iK5d(t~v&^tJ< z;i8&@9i7M5PrS)Pab)L5k)#!`mp&iP2qE#jt!lffUXHp2wmGr8?}9V(HqLt&4@U1k z>n^a!L7_-P@?KdBZqGHnpzd}=&+$xL;VegoZfVt~UMcKPa!Q720)#dwc2N5(kvtjv zV#uHzNtXU2_ZPLnEu;3%&zoiVD$!R)Qb_cf2ucQ5iLq%cfBG+(JWv-qmd{#I2=m*M zj(&y>_-#d7uhoP08U0k{L?#L5@#T7d=twd4j~5~CQO$VbTW)A*Qw4E<#!E4^8qX~g zS@TZ{(H#{q(HJU&rF+ruf~PVpUE}uLr&R=%6GLA)i-&@(KCz|pVpL4|=G~2LgmIJa z>`0jey)2f~f1Knf-KxJ_(?E(O(-6f|;x_6cMUyvfR3P9X-E3s%b*OfmyT$5@q1W|U z@6bCE!$PY|wqg=0^QHy2VnnDI+iMwMSA#K=B}yfE4rf1|&oTK)fckgCS8+o!!O;MrNciQf@ zM1oC%LW{A`Ti}a@j&qNeLg7Z|zxL zz5_o-@;`Hasz%b6RkW19a?!dmdqm4t3~I*xN7Nn#dn5KZO=!m+hn#DV`;-WJqfT!uP=H&(;?nX&*la~#*nIUmMzohJ zai9_t0~Xrqi)~nwT%vXkk>m1FW-(nu1<7y2fp2I0U3f6mZ2HhGIHViZwHKG-?_nnA ziBvvpD<~JrY@`?r%Q5+6phQ6LLep1fB;t)}BM%m~q9~xb)yc043$g~kTu^ag5IA7y zIVnQ!{E1K9xpMI4808kemm?{&YhHtm7_~Y_vUB&>qW5Bwu}4H1g1wi8`=tw@nHI6o zIJ{#f$K%3E{h3_ztLHg-O7L;bBl}}heE4Mt_WW!8IjGNlx!uLD!{IX@PPrSlgXVY6 z(?UxQOVP8iCF~CD7@5#~e54gCy^aHKE`xPg*L1|N8kFDH?mm?%#I`*_pTmP?__i(m zZ;PMFaJj~$??<5sq7QMLvN=2qXqs;fC?s&H|B%#Wx(OM+PfYdvq;Pv^u|MP(iSuSe zvLLM;vV51keSAKM?&_t@W#>`yYM7B!S_5j!D%&ynT`&qZXeJes7!^qN-kp@ffl<|1 z#A|_F^;NpT&dYFHlYPv&)DxOOG7J%YGHXDAj*` zKfMOyp6vF~!)y$Be>)JU(S-ZXv}w6S z^B{1*X>QIYz3cE(^xn*3OR&=Xb826c0ws&yPmIpY+0s=V#gVn3G;zo=ua$$PoR4v< z4;DckJg~w4j0lO|T2t3ulvt~#J7oka;Wgisthc?3Uk>KJ7;3x0_KGX zl?#avH!CukOYXyQc;AUnWG7DX*QNf~Qi!208;0rHEYuD-KPnjB2V&;qu&n Date: Thu, 9 Feb 2023 14:56:21 +0100 Subject: [PATCH 06/14] Mark as MPI-not-supported --- block/algebraic/petsc/solver.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/block/algebraic/petsc/solver.py b/block/algebraic/petsc/solver.py index 7d4f5f4..ebfa75c 100644 --- a/block/algebraic/petsc/solver.py +++ b/block/algebraic/petsc/solver.py @@ -1,4 +1,4 @@ -from block import block_mat, block_vec +from block import block_mat, block_vec, supports_mpi from block.block_base import block_base, block_container from block.object_pool import vec_pool from dolfin import GenericVector, Matrix, PETScVector @@ -96,6 +96,8 @@ def setUp(self, pc): class petsc_solver(petsc_base): def __init__(self, A, precond, V, ksp_type, prefix, options, defaults={}): + supports_mpi(False, 'PETSc solver interface currently does not work in parallel') + self.A = A self.Ad = mat(A) From dd426e608e271c54049a8a39fddd260dbbf32055 Mon Sep 17 00:00:00 2001 From: mirok Date: Thu, 23 Mar 2023 18:39:38 +0100 Subject: [PATCH 07/14] bring API closer to that of block.iterative --- block/algebraic/petsc/precond.py | 10 ++++++---- block/algebraic/petsc/solver.py | 17 ++++++++++++++--- 2 files changed, 20 insertions(+), 7 deletions(-) diff --git a/block/algebraic/petsc/precond.py b/block/algebraic/petsc/precond.py index d32cb4e..034ebbb 100644 --- a/block/algebraic/petsc/precond.py +++ b/block/algebraic/petsc/precond.py @@ -124,16 +124,16 @@ def __init__(self, A, parameters=None, pdes=1, nullspace=None, prefix=None): super().__init__(A, PETSc.PC.Type.ILU, pdes, nullspace, options=parameters, prefix=prefix) class Cholesky(precond): - def __init__(self, A, parameters=None): + 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): + def __init__(self, A, parameters=None, prefix=None): super().__init__(A, PETSc.PC.Type.LU, 1, None, options=parameters, prefix=prefix) class MUMPS_LU(LU): - def __init__(self, A, parameters=None): + def __init__(self, A, parameters=None, prefix=None): super().__init__(A, PETSc.PC.Type.LU, 1, None, options=parameters, prefix=prefix, defaults={ 'pc_factor_mat_solver_package': 'mumps', @@ -141,7 +141,7 @@ def __init__(self, A, parameters=None): class SUPERLU_LU(LU): - def __init__(self, A, parameters=None): + def __init__(self, A, parameters=None, prefix=None): super().__init__(A, PETSc.PC.Type.LU, 1, None, options=parameters, prefix=prefix, defaults={ 'pc_factor_mat_solver_package': 'superlu', @@ -293,7 +293,9 @@ def __init__(self, A, V, ams_zero_beta_poisson=False): # NOTE: signal that we don't have lower order term 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 diff --git a/block/algebraic/petsc/solver.py b/block/algebraic/petsc/solver.py index ebfa75c..93bf296 100644 --- a/block/algebraic/petsc/solver.py +++ b/block/algebraic/petsc/solver.py @@ -1,7 +1,7 @@ from block import block_mat, block_vec, supports_mpi from block.block_base import block_base, block_container from block.object_pool import vec_pool -from dolfin import GenericVector, Matrix, PETScVector +from dolfin import GenericVector, Matrix, PETScVector, PETScMatrix import numpy as np from petsc4py import PETSc from .precond import petsc_base, mat as single_mat, vec as single_vec @@ -31,7 +31,7 @@ def mat(A, _sizes=None): Ad.setPythonContext(petsc_py_wrapper(A)) Ad.setUp() return Ad - elif isinstance(A, Matrix): + elif isinstance(A, (PETScMatrix, Matrix)): return single_mat(A) else: raise TypeError(str(type(A))) @@ -102,11 +102,14 @@ def __init__(self, A, precond, V, ksp_type, prefix, options, defaults={}): self.Ad = mat(A) prefix, self.optsDB = self._merge_options(prefix=prefix, options=options, defaults=defaults) - self.petsc_op = PETSc.KSP().create(V.mesh().mpi_comm() if V else None) self.petsc_op.setOptionsPrefix(prefix) if ksp_type is not None: self.petsc_op.setType(ksp_type) + + self.petsc_op.setConvergenceHistory() # Record residuals + self.residuals = [] + self.petsc_op.setOperators(self.Ad) self.petsc_op.setFromOptions() if precond is not None: @@ -121,6 +124,9 @@ def matvec(self, b): self.petsc_op.setUp() x = self.Ad.createVecLeft() self.petsc_op.solve(vec(b), x) + + self.residuals = self.petsc_op.getConvergenceHistory() + return Vec(x, self.A) @vec_pool @@ -134,6 +140,11 @@ def create_vec(self, dim=1): raise ValueError('dim must be 0 or 1') return PETScVector(m) + @property + def iterations(self): + return len(self.residuals)-1 + + class KSP(petsc_solver): def __init__(self, A, precond=None, prefix=None, **parameters): super().__init__(A, precond=precond, V=None, ksp_type=None, prefix=prefix, options=parameters, From 00672663d918f6e88805391e7efc31a4df731aa8 Mon Sep 17 00:00:00 2001 From: mirok Date: Sat, 25 Mar 2023 14:29:06 +0100 Subject: [PATCH 08/14] eigenvalue_estimates from petsc --- block/algebraic/petsc/solver.py | 11 +++++++++++ 1 file changed, 11 insertions(+) diff --git a/block/algebraic/petsc/solver.py b/block/algebraic/petsc/solver.py index 5ce9a64..856da78 100644 --- a/block/algebraic/petsc/solver.py +++ b/block/algebraic/petsc/solver.py @@ -111,6 +111,8 @@ def __init__(self, A, precond, V, ksp_type, prefix, options, defaults={}): self.petsc_op.setConvergenceHistory() # Record residuals self.residuals = [] + self.petsc_op.setComputeEigenvalues(True) + self.petsc_op.setOperators(self.Ad) self.petsc_op.setFromOptions() if precond is not None: @@ -145,6 +147,15 @@ def create_vec(self, dim=1): def iterations(self): return len(self.residuals)-1 + def eigenvalue_estimates(self): + if self.petsc_op.getComputeEigenvalues(): + return self.petsc_op.computeEigenvalues() + print(self.petsc_op.getComputeSingularValues()) + exit() + if self.petsc_op.getComputeSingularValues(): + return self.petsc_op.computeExtremeSingularValues() + + return np.array([]) class KSP(petsc_solver): def __init__(self, A, precond=None, prefix=None, **parameters): From 43a2efaeba85414662204e263090f3e5bae78b7c Mon Sep 17 00:00:00 2001 From: mirok Date: Fri, 7 Apr 2023 07:34:23 +0200 Subject: [PATCH 09/14] Handle scalar 0 blocks not on diagonal as None --- block/algebraic/petsc/solver.py | 10 +++++++--- 1 file changed, 7 insertions(+), 3 deletions(-) diff --git a/block/algebraic/petsc/solver.py b/block/algebraic/petsc/solver.py index 856da78..a81d859 100644 --- a/block/algebraic/petsc/solver.py +++ b/block/algebraic/petsc/solver.py @@ -19,7 +19,8 @@ def mat(A, _sizes=None): # a vector to find the size is wasteful but easy row_sz = [(v.local_size(), v.size()) for v in A.create_vec(dim=1)] col_sz = [(v.local_size(), v.size()) for v in A.create_vec(dim=0)] - m = PETSc.Mat().createNest([[mat(m, (row_sz[r],col_sz[c])) for (c,m) in enumerate(row)] + + m = PETSc.Mat().createNest([[mat(m, (row_sz[r], col_sz[c])) for (c,m) in enumerate(row)] for (r,row) in enumerate(A.blocks)]) if USE_VECNEST: m.setVecType(PETSc.Vec.Type.NEST) @@ -27,6 +28,10 @@ def mat(A, _sizes=None): elif isinstance(A, block_base) or np.isscalar(A): if not _sizes: raise ValueError('Cannot create unsized PETSc matrix') + + if A == 0: + return None + Ad = PETSc.Mat().createPython(_sizes) Ad.setPythonContext(petsc_py_wrapper(A)) Ad.setUp() @@ -150,8 +155,7 @@ def iterations(self): def eigenvalue_estimates(self): if self.petsc_op.getComputeEigenvalues(): return self.petsc_op.computeEigenvalues() - print(self.petsc_op.getComputeSingularValues()) - exit() + if self.petsc_op.getComputeSingularValues(): return self.petsc_op.computeExtremeSingularValues() From 897bb5fa626a64f172e964996e161b6b6766aa1c Mon Sep 17 00:00:00 2001 From: mirok Date: Fri, 4 Aug 2023 11:38:32 +0200 Subject: [PATCH 10/14] update --- block/algebraic/petsc/precond.py | 36 +++++++++++++++++--------------- block/algebraic/petsc/solver.py | 10 ++++++--- block/block_vec.py | 4 ++++ 3 files changed, 30 insertions(+), 20 deletions(-) diff --git a/block/algebraic/petsc/precond.py b/block/algebraic/petsc/precond.py index 411735d..1330af2 100644 --- a/block/algebraic/petsc/precond.py +++ b/block/algebraic/petsc/precond.py @@ -66,16 +66,17 @@ def __init__(self, A, prectype, pdes, nullspace, prefix, options, V=None, defaul self.petsc_op = PETSc.PC().create(V.mesh().mpi_comm() if V else None) self.petsc_op.setOptionsPrefix(prefix) - try: - self.petsc_op.setType(prectype) - except PETSc.Error as e: - if e.ierr == 86: - raise ValueError(f'Unknown PETSc type "{prectype}"') - else: - raise + 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.petsc_op.setOperators(Ad) + self.petsc_op.setFromOptions() self.setup_time = -1. @@ -94,8 +95,8 @@ def matvec(self, b): 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): @@ -384,7 +385,6 @@ 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() @@ -392,12 +392,14 @@ def __init__(self, A, V, prefix=None): 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_op.setHYPREDiscreteGradient(G) self.petsc_op.setHYPREDiscreteCurl(G) self.petsc_op.setCoordinates(xyz) + self.petsc_op.setUp() + diff --git a/block/algebraic/petsc/solver.py b/block/algebraic/petsc/solver.py index a81d859..2089b7d 100644 --- a/block/algebraic/petsc/solver.py +++ b/block/algebraic/petsc/solver.py @@ -58,7 +58,7 @@ def vec(v): arr[i0:i0+len(arr2)] = arr2 i0 += len(arr2) return vs - elif isinstance(v, GenericVector): + elif isinstance(v, (PETScVector, GenericVector)): return single_vec(v) else: raise TypeError(str(type(v))) @@ -76,7 +76,9 @@ def Vec(v, creator): arr = v.getArray() i0 = 0 for vv in ret: - vv.set_local(arr[i0:i0+len(vv)]) + # vv.set_local(arr[i0:i0+len(vv)]) + vv *= 0. + vv.axpy(1., PETScVector(PETSc.Vec().createWithArray(arr[i0:i0+len(vv)]))) i0 += len(vv) return ret else: @@ -91,7 +93,9 @@ def __init__(self, A): def mult(self, mat, x, y): new_y = vec(self.A * Vec(x, self.A)) - y.aypx(0, new_y) + # y.aypx(0, new_y) + y *= 0 + y += new_y def apply(self, pc, x, y): return self.mult(None, x, y) diff --git a/block/block_vec.py b/block/block_vec.py index 8b2df96..9b15148 100644 --- a/block/block_vec.py +++ b/block/block_vec.py @@ -43,8 +43,12 @@ def allocate(self, template, dim=None, alternative_templates=[]): except Exception: pass if not isinstance(self[i], GenericVector): + from IPython import embed + embed() + raise ValueError( f"Can't allocate vector - no usable template for block {i}.\n" + f'{self[i]}' "Consider calling something like bb.allocate([V, Q]) to initialise the block_vec." ) self[i].zero() From 6484a27e8b533212f1289b97b5a01a26d8c71fea Mon Sep 17 00:00:00 2001 From: Joachim B Haga Date: Wed, 17 Apr 2024 09:50:40 +0200 Subject: [PATCH 11/14] Work around a liveness issue in fenics --- demo/petsc.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/demo/petsc.py b/demo/petsc.py index 3e52266..f3ab95c 100644 --- a/demo/petsc.py +++ b/demo/petsc.py @@ -72,10 +72,10 @@ def eval_cell(self, values, x, ufc_cell): def value_shape(self): return (2,) -bcs = block_bc([ - [DirichletBC(BDM, BoundarySource(mesh), BoxBoundary(mesh).ns)], - None -], symmetric=True) +G = BoundarySource(mesh) +boundary = BoxBoundary(mesh).ns +bcs_BDM = [DirichletBC(BDM, G, boundary)] +bcs = block_bc([bcs_BDM, None], symmetric=True) f = Expression("10*exp(-(pow(x[0] - 0.5, 2) + pow(x[1] - 0.5, 2)) / 0.02)", degree=4) From 1ce12484607e5907dea06eddd6c29db10563e5f6 Mon Sep 17 00:00:00 2001 From: mirok Date: Wed, 8 May 2024 08:56:48 +0200 Subject: [PATCH 12/14] pass nullspace to KSP --- block/algebraic/petsc/precond.py | 14 ++++++++++++-- block/algebraic/petsc/solver.py | 10 ++++++++-- block/block_mat.py | 15 +++++++++++++-- 3 files changed, 33 insertions(+), 6 deletions(-) diff --git a/block/algebraic/petsc/precond.py b/block/algebraic/petsc/precond.py index 1330af2..ae1687c 100644 --- a/block/algebraic/petsc/precond.py +++ b/block/algebraic/petsc/precond.py @@ -149,12 +149,22 @@ def __init__(self, A, parameters=None, prefix=None, defaults={}): pdes=1, nullspace=None, options=parameters, prefix=prefix, 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', }) @@ -163,7 +173,7 @@ 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_op.setFactorPivot(1E-16) diff --git a/block/algebraic/petsc/solver.py b/block/algebraic/petsc/solver.py index 2089b7d..e4dcad4 100644 --- a/block/algebraic/petsc/solver.py +++ b/block/algebraic/petsc/solver.py @@ -104,7 +104,7 @@ def setUp(self, pc): pass class petsc_solver(petsc_base): - def __init__(self, A, precond, V, ksp_type, prefix, options, defaults={}): + def __init__(self, A, precond, V, ksp_type, prefix, options, defaults={}, nullspace=None): supports_mpi(False, 'PETSc solver interface currently does not work in parallel') self.A = A @@ -122,8 +122,13 @@ def __init__(self, A, precond, V, ksp_type, prefix, options, defaults={}): self.petsc_op.setComputeEigenvalues(True) + if nullspace is not None: + ns = PETSc.NullSpace().create(vectors=[vec(fi) for fi in nullspace]) + self.Ad.setNearNullSpace(ns) self.petsc_op.setOperators(self.Ad) + self.petsc_op.setFromOptions() + if precond is not None: pc = self.petsc_op.getPC() pc.setType(PETSc.PC.Type.PYTHON) @@ -166,7 +171,8 @@ def eigenvalue_estimates(self): return np.array([]) class KSP(petsc_solver): - def __init__(self, A, precond=None, prefix=None, **parameters): + def __init__(self, A, precond=None, prefix=None, nullspace=None, **parameters): super().__init__(A, precond=precond, V=None, ksp_type=None, prefix=prefix, options=parameters, + nullspace=nullspace, defaults={ }) diff --git a/block/block_mat.py b/block/block_mat.py index 1c55d1d..8dc6213 100644 --- a/block/block_mat.py +++ b/block/block_mat.py @@ -2,6 +2,8 @@ from builtins import range from .block_base import block_container from .block_vec import block_vec +import numpy as np + class block_mat(block_container): """Block of matrices or other operators. Empty blocks doesn't need to be set @@ -18,7 +20,7 @@ def __init__(self, m, n=None, blocks=0): block_container.__init__(self, (m,n), blocks) def matvec(self, x): - from dolfin import GenericVector, Matrix, PETScMatrix + from dolfin import GenericVector, Matrix, PETScMatrix, PETScVector from .block_util import isscalar import numpy m,n = self.blocks.shape @@ -45,7 +47,16 @@ def matvec(self, x): # Do the block multiply if isinstance(self[i,j], PETScMatrix): z = self[i,j].create_vec(dim=0) - self[i,j].mult(x[j], z) + + if isinstance(x[j], (GenericVector, PETScVector)): + self[i,j].mult(x[j], z) + else: + xx = self.create_vec() + offsets = np.cumsum([0] + [xi.size() for xi in xx]) + xj = xx[j] + + xj.set_local(x.get_local()[offsets[j]:offsets[j+1]]) + self[i, j].mult(xj, z) else: z = self[i,j] * x[j] if isinstance(z, type(NotImplemented)): return NotImplemented From 9e8d1c123f3e80ae45de0acd48114fcf9f3a35ad Mon Sep 17 00:00:00 2001 From: mirok Date: Tue, 25 Jun 2024 09:26:10 +0200 Subject: [PATCH 13/14] KSP with shells --- block/algebraic/petsc/solver.py | 19 ++++++++++++++++--- 1 file changed, 16 insertions(+), 3 deletions(-) diff --git a/block/algebraic/petsc/solver.py b/block/algebraic/petsc/solver.py index e4dcad4..cdd1020 100644 --- a/block/algebraic/petsc/solver.py +++ b/block/algebraic/petsc/solver.py @@ -24,7 +24,12 @@ def mat(A, _sizes=None): for (r,row) in enumerate(A.blocks)]) if USE_VECNEST: m.setVecType(PETSc.Vec.Type.NEST) + return m + + elif isinstance(A, (PETScMatrix, Matrix)): + return single_mat(A) + elif isinstance(A, block_base) or np.isscalar(A): if not _sizes: raise ValueError('Cannot create unsized PETSc matrix') @@ -36,10 +41,18 @@ def mat(A, _sizes=None): Ad.setPythonContext(petsc_py_wrapper(A)) Ad.setUp() return Ad - elif isinstance(A, (PETScMatrix, Matrix)): - return single_mat(A) else: - raise TypeError(str(type(A))) + try: + row_vec, col_vec = A.create_vec(dim=1), A.create_vec(dim=0) + _sizes = (row_vec.local_size(), col_vec.local_size()) + + Ad = PETSc.Mat().createPython(_sizes) + Ad.setPythonContext(petsc_py_wrapper(A)) + Ad.setUp() + + return Ad + except: + raise TypeError(str(type(A))) def vec(v): if isinstance(v, block_vec): From 2554c792cafc7f95377aaea3c8f474dfa8c0f57f Mon Sep 17 00:00:00 2001 From: mirok Date: Thu, 11 Jul 2024 10:52:44 +0200 Subject: [PATCH 14/14] Vec for block vectors --- block/algebraic/petsc/solver.py | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/block/algebraic/petsc/solver.py b/block/algebraic/petsc/solver.py index cdd1020..079e11a 100644 --- a/block/algebraic/petsc/solver.py +++ b/block/algebraic/petsc/solver.py @@ -84,8 +84,12 @@ def Vec(v, creator): else: return PETScVector(v) else: - if isinstance(creator, block_container): + try: ret = creator.create_vec(dim=1) + except: + pass + + if isinstance(ret, block_vec): arr = v.getArray() i0 = 0 for vv in ret: