Secure your code as it's written. Use Snyk Code to scan source code in minutes - no build needed - and fix issues immediately.
def petsc_serial_matrix(test_space, trial_space, nnz=None):
'''
PETsc.Mat from trial_space to test_space to be filled in the
with block. The spaces can be represented by intergers meaning
generic R^n.
'''
# Decide local to global map
# For our custom case everything is serial
if is_number(test_space) and is_number(trial_space):
comm = mpi_comm_world().tompi4py()
# Local same as global
sizes = [[test_space, test_space], [trial_space, trial_space]]
row_map = PETSc.IS().createStride(test_space, 0, 1, comm)
col_map = PETSc.IS().createStride(trial_space, 0, 1, comm)
# With function space this can be extracted
else:
mesh = test_space.mesh()
comm = mesh.mpi_comm().tompi4py()
row_map = test_space.dofmap()
col_map = trial_space.dofmap()
sizes = [[row_map.index_map().size(IndexMap.MapSize_OWNED),
row_map.index_map().size(IndexMap.MapSize_GLOBAL)],
[col_map.index_map().size(IndexMap.MapSize_OWNED),
col_map.index_map().size(IndexMap.MapSize_GLOBAL)]]
row_map = map(int, row_map.tabulate_local_to_global_dofs())
col_map = map(int, col_map.tabulate_local_to_global_dofs())
self.pressureDOF = numpy.arange(start=L_range[0],
stop=L_range[0]+neqns,
step=nc,
dtype="i")
velocityDOF = []
for start in range(1,nc):
velocityDOF.append(numpy.arange(start=L_range[0]+start,
stop=L_range[0]+neqns,
step=nc,
dtype="i"))
self.velocityDOF = numpy.vstack(velocityDOF).transpose().flatten()
self.pc = p4pyPETSc.PC().create()
self.pc.setType('fieldsplit')
self.isp = p4pyPETSc.IS()
self.isp.createGeneral(self.pressureDOF,comm=p4pyPETSc.COMM_WORLD)
self.isv = p4pyPETSc.IS()
self.isv.createGeneral(self.velocityDOF,comm=p4pyPETSc.COMM_WORLD)
self.pc.setFieldSplitIS(('velocity',self.isv),('pressure',self.isp))
rows.extend(shift + np.array(bc.keys(), dtype='int32'))
x_values.extend(bc.values())
rows = np.hstack(rows)
x_values = np.array(x_values)
x = bb.copy()
x.zeroEntries()
x.setValues(rows, x_values)
# Apply to monolithic
len(rows) and AA.zeroRowsColumns(rows, diag=diag_val, x=x, b=bb)
blocks = []
for first, last in zip(offsets[:-1], offsets[1:]):
blocks.append(PETSc.IS().createStride(last-first, first, 1))
# Reasamble
comm = mpi_comm_world()
b = [PETSc.Vec().createWithArray(bb.getValues(block), comm=comm) for block in blocks]
for bi in b:
bi.assemblyBegin()
bi.assemblyEnd()
b = block_vec(map(PETScVector, b))
# PETSc 3.7.x
try:
AA.getSubMatrix
A = [[PETScMatrix(AA.getSubMatrix(block_row, block_col)) for block_col in blocks]
for block_row in blocks]
# NOTE: 3.8+ does not ahve getSubMatrix, there is getLocalSubMatrix
# but cannot get that to work - petsc error 73. So for now everything
dofmap = bc.function_space().dofmap()
local_range = dofmap.ownership_range()
local_dofs.update(list(range(local_range[0], local_range[1])))
constrained_local_dofs.update([
dofmap.local_to_global_index(local_dof_index) for local_dof_index in bc.get_boundary_values().keys()
])
# List all unconstrained dofs
unconstrained_local_dofs = local_dofs.difference(constrained_local_dofs)
unconstrained_local_dofs = list(unconstrained_local_dofs)
# Generate IS accordingly
comm = bcs[0].function_space().mesh().mpi_comm()
for bc in bcs:
assert comm == bc.function_space().mesh().mpi_comm()
self._is = PETSc.IS().createGeneral(unconstrained_local_dofs, comm)
def field_ises(self):
"""A list of PETSc ISes defining the global indices for each set in
the DataSet.
Used when extracting blocks from matrices for solvers."""
ises = []
nlocal_rows = 0
for dset in self:
nlocal_rows += dset.size * dset.cdim
offset = self.comm.scan(nlocal_rows)
offset -= nlocal_rows
for dset in self:
nrows = dset.size * dset.cdim
iset = PETSc.IS().createStride(nrows, first=offset, step=1,
comm=self.comm)
iset.setBlockSize(dset.cdim)
ises.append(iset)
offset += nrows
return tuple(ises)
step=3,
dtype="i")
velocityDOF = []
for start in range(1,3):
velocityDOF.append(numpy.arange(start=L_range[0]+start,
stop=L_range[0]+neqns,
step=3,
dtype="i"))
self.velocityDOF = numpy.vstack(velocityDOF).transpose().flatten()
self.pc = p4pyPETSc.PC().create()
if prefix:
self.pc.setOptionsPrefix(prefix)
self.pc.setType('fieldsplit')
self.isp = p4pyPETSc.IS()
self.isp.createGeneral(self.pressureDOF,comm=p4pyPETSc.COMM_WORLD)
self.isv = p4pyPETSc.IS()
self.isv.createGeneral(self.velocityDOF,comm=p4pyPETSc.COMM_WORLD)
self.pc.setFieldSplitIS(('velocity',self.isv),('pressure',self.isp))
self.pc.setFromOptions()
def setUp(self, global_ksp=None):
if(i == i0n + In - 1 and j0 == j0n and k0 == k0n):
ighosts += list(range(gindex+2*In*Jn+2*In*Kn+Jn*Kn,
gindex+2*In*Jn+2*In*Kn+2*Jn*Kn))
# zmax boundary
i = i0+I
if(i > Nz-1): i = 0
for pd in pos_data:
gindex, i0n, j0n, k0n, In, Jn, Kn = pd
if(i == i0n and j0 == j0n and k0 == k0n):
ighosts += list(range(gindex+2*In*Jn+2*In*Kn,
gindex+2*In*Jn+2*In*Kn+Jn*Kn))
# Set the ghost indices = finish constructing ghost vector
self._ighosts = np.array(ighosts, dtype=np.int32)
ghosts = PETSc.IS().createGeneral(self._ighosts)
self._gvec.setMPIGhost(ghosts)