Secure your code as it's written. Use Snyk Code to scan source code in minutes - no build needed - and fix issues immediately.
self._da = PETSc.DMDA().create(
[vs.nx, vs.ny],
stencil_width=1,
stencil_type='star',
comm=rs.mpi_comm,
proc_sizes=rs.num_proc,
boundary_type=boundary_type,
ownership_ranges=[
(vs.nx // rs.num_proc[0],) * rs.num_proc[0],
(vs.ny // rs.num_proc[1],) * rs.num_proc[1],
]
)
self._matrix, self._boundary_fac = self._assemble_poisson_matrix(vs)
petsc_options = PETSc.Options()
# setup krylov method
self._ksp = PETSc.KSP()
self._ksp.create(rs.mpi_comm)
self._ksp.setOperators(self._matrix)
self._ksp.setType('bcgs')
self._ksp.setTolerances(atol=0, rtol=vs.congr_epsilon, max_it=vs.congr_max_iterations)
# preconditioner
self._ksp.getPC().setType('hypre')
petsc_options['pc_hypre_type'] = 'boomeramg'
self._ksp.getPC().setFromOptions()
self._rhs_petsc = self._da.createGlobalVec()
self._sol_petsc = self._da.createGlobalVec()
def _setPETSc(self,petsc_file):
"""The function takes a file with petsc options and sets the options globally.
petsc_file : str
string with the location of the file
"""
# First, clear any existing PETSc options settings
for key in petsc4py.PETSc.Options().getAll():
petsc4py.PETSc.Options().delValue(key)
# Second, collect and add new PETSc options
petsc_options = []
with open(petsc_file) as petsc_file:
data = petsc_file.readlines()
def strip_comments(line):
if '#' in line:
line = line[:line.index('#')]
return line
stripped_data = [strip_comments(line) for line in data]
petsc_options = ''.join(stripped_data).split('\n')
new_petsc = []
for item in petsc_options:
if item != '':
new = item.split()
new[0] = new[0][1:]
t4 = t.time()
umfcon = A_lu.func_closure[1].cell_contents
err = scipy.linalg.norm(A.matvec(N.array(x_umf)) - RHS)
print 'LU decomp time: %f s' % (t2-t1)
print 'Solve time 1 : %f s' % (t3-t2)
print 'Solve time 2 : %f s' % (t4-t3)
print 'Err : %f' % err
print 'dofs: %d nnz: %d' % (A.shape[0], A.nnz)
else:
guess = int(sys.argv[3])
A = A.tocsr()
A.ensure_sorted_indices(inplace=True)
import petsc4py
petsc4py.init(sys.argv)
from petsc4py import PETSc
OptDB = PETSc.Options()
for k,v in precon_opts[precond_type].items():
OptDB[k]=v
A_petsc = PETSc.Mat().createAIJ(A.shape,
csr=(A.indptr,
A.indices,
A.data))
# obtain vectors for storing
# the solution and the rhs
x_petsc, b_petsc = A_petsc.getVecs()
# fill the rhs PETSc vector
# from the rhs numpy array
b_petsc[...] = RHS # note the syntax sugar
ksp = PETSc.KSP().create()
ksp.setFromOptions() # configure from OptDB
ksp.setInitialGuessNonzero(guess)
dF = block_jacobian(F, up)
# Setup H1 x H-0.5 preconditioner only once
B0 = ii_derivative(inner(grad(u), grad(v))*dx + inner(u, v)*dx, u)
# The Q norm via spectral
B1 = inverse(HsNorm(Q, s=-0.5)) # The norm is inverted exactly
# Preconditioner
B = block_diag_mat([AMG(ii_assemble(B0)), B1])
# Newton
eps = 1.0
tol = 1.0E-10
niter = 0.
maxiter = 25
OptDB = PETSc.Options()
# Since later b gets very small making relative too much work
OptDB.setValue('ksp_atol', 1E-6)
OptDB.setValue('ksp_monitor_true_residual', None)
dup = ii_Function(W)
x_vec = as_backend_type(dup.vector()).vec()
n_kspiters = []
while eps > tol and niter < maxiter:
niter += 1.
A, b = (ii_assemble(x) for x in (dF, F))
ksp = PETSc.KSP().create()
ksp.setType('minres')
ksp.setNormType(PETSc.KSP.NormType.NORM_PRECONDITIONED)
tao = PETSc.TAO().create(tao_comm)
tao.setType(user_specs['localopt_method'])
if user_specs['localopt_method'] == 'pounders':
f = PETSc.Vec().create(tao_comm)
f.setSizes(m)
f.setFromOptions()
if hasattr(tao, 'setResidual'):
tao.setResidual(lambda tao, x, f: tao_callback_fun_pounders(tao, x, f, comm_queue, child_can_read,
parent_can_read, user_specs), f)
else:
tao.setSeparableObjective(lambda tao, x, f: tao_callback_fun_pounders(tao, x, f, comm_queue, child_can_read,
parent_can_read, user_specs), f)
delta_0 = user_specs['dist_to_bound_multiple']*np.min([np.min(ub.array-x.array), np.min(x.array-lb.array)])
PETSc.Options().setValue('-tao_pounders_delta', str(delta_0))
elif user_specs['localopt_method'] == 'nm':
tao.setObjective(lambda tao, x: tao_callback_fun_nm(tao, x, comm_queue, child_can_read,
parent_can_read, user_specs))
elif user_specs['localopt_method'] == 'blmvm':
g = PETSc.Vec().create(tao_comm)
g.setSizes(n)
g.setFromOptions()
tao.setObjectiveGradient(lambda tao, x, g: tao_callback_fun_grad(tao, x, g, comm_queue, child_can_read,
parent_can_read, user_specs))
# Set everything for tao before solving
PETSc.Options().setValue('-tao_max_funcs', str(user_specs.get('run_max_eval', 1000*n)))
tao.setFromOptions()
tao.setVariableBounds((lb, ub))
f.setSizes(m)
f.setFromOptions()
delta_0 = gen_specs['dist_to_bound_multiple']*np.min([np.min(ub.array-x.array), np.min(x.array-lb.array)])
PETSc.Options().setValue('-tao_pounders_delta',str(delta_0))
# PETSc.Options().setValue('-pounders_subsolver_tao_type','bqpip')
tao.setSeparableObjective(lambda tao, x, f: pounders_obj_func(tao, x, f, Run_H), f)
# elif gen_specs['localopt_method'] == 'blmvm':
# g = PETSc.Vec().create(tao_comm)
# g.setSizes(n)
# g.setFromOptions()
# tao.setObjectiveGradient(lambda tao, x, g: blmvm_obj_func(tao, x, g, Run_H))
# Set everything for tao before solving
PETSc.Options().setValue('-tao_max_funcs',str(total_pts_in_run+1))
tao.setFromOptions()
tao.setVariableBounds((lb,ub))
# tao.setObjectiveTolerances(fatol=gen_specs['fatol'], frtol=gen_specs['frtol'])
# tao.setGradientTolerances(grtol=gen_specs['grtol'], gatol=gen_specs['gatol'])
tao.setTolerances(grtol=gen_specs['grtol'], gatol=gen_specs['gatol'])
tao.setInitial(x)
tao.solve(x)
x_opt = tao.getSolution().getArray()
exit_code = tao.getConvergedReason()
# print(exit_code)
# print(tao.view())
# print(x_opt)
# if gen_specs['localopt_method'] == 'pounders':
B : petsc4py matrix object
The discrete divergence operator.
Bt : petsc4py matrix object
The discrete gradient operator.
F : petsc4py matrix object
The A-block of the linear system.
"""
# TODO - Find a good way to assert that Qv is diagonal
self.Qv = Qv
self.B = B
self.Bt = Bt
self.F = F
self._constructBQinvBt()
self._options = p4pyPETSc.Options()
if self._options.hasName('innerLSCsolver_BTinvBt_ksp_constant_null_space'):
self._create_constant_nullspace()
self.BQinvBt.setNullSpace(self.const_null_space)
self.kspBQinvBt = p4pyPETSc.KSP().create()
self.kspBQinvBt.setOperators(self.BQinvBt,self.BQinvBt)
self.kspBQinvBt.setOptionsPrefix('innerLSCsolver_BTinvBt_')
self.kspBQinvBt.pc.setUp()
self.kspBQinvBt.setFromOptions()
self.kspBQinvBt.setUp()
# initialize solver for Qv
self.kspQv = p4pyPETSc.KSP().create()
self.kspQv.setOperators(self.Qv,self.Qv)
self.kspQv.setOptionsPrefix('innerLSCsolver_T_')
def krylov_solve(wh, AA, bb, BB, petsc_args):
'''The problem as hand is symmetric and by num evidence a positive definite'''
## AA and BB as block_mat
ksp = PETSc.KSP().create()
opts = PETSc.Options()
for key, value in zip(petsc_args[::2], petsc_args[1::2]):
opts.setValue(key, None if value == 'none' else value)
# FIXME: we do things differently for block mat and monolithic mat
if isinstance(AA, block_mat):
ksp.setOperators(ii_PETScOperator(AA))
ksp.setPC(ii_PETScPreconditioner(BB, ksp))
# b and x in A*x = b
b = as_petsc_nest(bb)
x = wh.petsc_vec()
else:
AA = as_backend_type(AA).mat()
# Monolithic
ksp.setOperators(AA)
assert BB is None # Use AA
def RunTest():
from petsc4py import PETSc
import example100
OptDB = PETSc.Options()
N = OptDB.getInt('N', 100)
draw = OptDB.getBool('draw', False)
A = PETSc.Mat()
A.create(comm=PETSc.COMM_WORLD)
A.setSizes([N,N])
A.setType(PETSc.Mat.Type.PYTHON)
A.setPythonContext(example100.Laplace1D())
A.setUp()
x, b = A.getVecs()
b.set(1)
ksp = PETSc.KSP()
ksp.create(comm=PETSc.COMM_WORLD)
ksp.setType(PETSc.KSP.Type.PYTHON)