[petsc-users] petsc4py snes behavior
Nicola Creati
ncreati at inogs.it
Mon Feb 8 02:41:35 CST 2016
Hello,
I'm trying to convert some PETSc examples from C to Python. I noted a
strange behavior of the snes solver. It seems that the same example in
Python and C are characterized by a different number of function
evaluations. I converted ex25.c example in the "
snes/examples/tutorials" folder. The number of function evaluations is
always bigger in Python and increases with the number of grid nodes. The
increased number of function evaluations affects the performance and in
other examples can drastically increase the execution time (for example
in ex30.c " Steady-state 2D subduction...").
This is my Python version of the simple ex25.c:
## # SNES/Tutorials/ex25
import sys
import petsc4py
petsc4py.init(sys.argv)
from petsc4py import PETSc
import numpy as np
class UserContext(object):
def __init__(self, da):
self.da = da
self.localX = da.createLocalVec()
def formFunction(self, snes, X, F):
self.da.globalToLocal(X, self.localX)
x = self.da.getVecArray(self.localX)
f = self.da.getVecArray(F)
(x0, x1), (y0, y1)= self.da.getRanges()
hx = 1/float(x1-1)
hy = 1/float(y1-1)
for j in range(y0, y1):
for i in range(x0, x1):
if i == 0 or i == (x1-1) or j == 0 or j == (y1-1 ):
f[j,i] = x[j,i] - (1.0 - (2.0 * hx * i - 1.0) *
(2.0 * hx * i -1.0))
else:
gradup = (x[j+1, i] -x[j, i])/hy
graddown = (x[j, i] -x[j-1, i])/hy
gradright = (x[j, i+1] -x[j, i])/hx
gradleft = (x[j, i] -x[j, i-1])/hx
gradx = 0.5 * (x[j, i+1] - x[j, i-1])/hx
grady = 0.5 * (x[j+1, i] - x[j-1, i])/hy
coeffup = 1.0/np.sqrt(1.0 + gradup * gradup + gradx
* gradx)
coeffdown = 1.0/np.sqrt(1.0 + graddown * graddown +
gradx * gradx)
coeffleft = 1.0/np.sqrt(1.0 + gradleft * gradleft
+ grady * grady)
coeffright = 1.0/np.sqrt(1.0 + gradright *
gradright + grady * grady)
f[j, i] = (coeffup * gradup - coeffdown * graddown)
* hx + (coeffright * gradright - coeffleft * gradleft) * hy
snes = PETSc.SNES().create()
da = PETSc.DMDA().create(dim=(-5, -5),
stencil_type = PETSc.DMDA.StencilType.STAR,
stencil_width =1,
setup=False)
da.setFromOptions()
da.setUp()
snes.setDM(da)
ctx = UserContext(da)
da.setAppCtx(ctx)
F = da.createGlobalVec()
snes.setFunction(ctx.formFunction, F)
snes.setFromOptions()
x = da.createGlobalVector()
snes.solve(None, x)
its = snes.getIterationNumber()
lits = snes.getLinearSolveIterations()
print "Number of SNES iterations = %d" % its
print "Number of Linear iterations = %d" % lits
litspit = lits/float(its)
print "Average Linear its / SNES = %e" % float(litspit)
#####################################
I executed the code with the following options:
python ex25.py -snes_view -snes_monitor -da_grid_x 50 -da_grid_y 50
This is the output from Python:
Konsole output
0 SNES Function norm 7.229568285876e+00
1 SNES Function norm 1.970854029907e-02
2 SNES Function norm 5.721344031051e-03
3 SNES Function norm 1.179738644817e-03
4 SNES Function norm 1.028879149924e-04
5 SNES Function norm 7.472012969660e-07
6 SNES Function norm 4.744184684658e-11
SNES Object: 1 MPI processes
type: newtonls
maximum iterations=50, maximum function evaluations=10000
tolerances: relative=1e-08, absolute=1e-50, solution=1e-08
total number of linear solver iterations=223
total number of function evaluations=37
norm schedule ALWAYS
SNESLineSearch Object: 1 MPI processes
type: bt
interpolation: cubic
alpha=1.000000e-04
maxstep=1.000000e+08, minlambda=1.000000e-12
tolerances: relative=1.000000e-08, absolute=1.000000e-15,
lambda=1.000000e-08
maximum iterations=40
KSP Object: 1 MPI processes
type: gmres
GMRES: restart=30, using Classical (unmodified) Gram-Schmidt
Orthogonalization with no iterative refinement
GMRES: happy breakdown tolerance 1e-30
maximum iterations=10000, initial guess is zero
tolerances: relative=1e-05, absolute=1e-50, divergence=10000
left preconditioning
using PRECONDITIONED norm type for convergence test
PC Object: 1 MPI processes
type: ilu
ILU: out-of-place factorization
0 levels of fill
tolerance for zero pivot 2.22045e-14
matrix ordering: natural
factor fill ratio given 1, needed 1
Factored matrix follows:
Mat Object: 1 MPI processes
type: seqaij
rows=2500, cols=2500
package used to perform factorization: petsc
total: nonzeros=12300, allocated nonzeros=12300
total number of mallocs used during MatSetValues calls =0
not using I-node routines
linear system matrix = precond matrix:
Mat Object: 1 MPI processes
type: seqaij
rows=2500, cols=2500
total: nonzeros=12300, allocated nonzeros=12300
total number of mallocs used during MatSetValues calls =0
not using I-node routines
Number of SNES iterations = 6
Number of Linear iterations = 223
Average Linear its / SNES = 3.716667e+01
This the output from the original C code with the same options (./ex25
-snes_view -snes_monitor -da_grid_x 50 -da_grid_y 50):
Konsole output
0 SNES Function norm 7.229568285876e+00
1 SNES Function norm 1.970854029907e-02
2 SNES Function norm 5.721344035538e-03
3 SNES Function norm 1.179738636852e-03
4 SNES Function norm 1.028879138425e-04
5 SNES Function norm 7.472012249912e-07
6 SNES Function norm 4.743884624164e-11
SNES Object: 1 MPI processes
type: newtonls
maximum iterations=50, maximum function evaluations=10000
tolerances: relative=1e-08, absolute=1e-50, solution=1e-08
total number of linear solver iterations=223
total number of function evaluations=7
norm schedule ALWAYS
SNESLineSearch Object: 1 MPI processes
type: bt
interpolation: cubic
alpha=1.000000e-04
maxstep=1.000000e+08, minlambda=1.000000e-12
tolerances: relative=1.000000e-08, absolute=1.000000e-15,
lambda=1.000000e-08
maximum iterations=40
KSP Object: 1 MPI processes
type: gmres
GMRES: restart=30, using Classical (unmodified) Gram-Schmidt
Orthogonalization with no iterative refinement
GMRES: happy breakdown tolerance 1e-30
maximum iterations=10000, initial guess is zero
tolerances: relative=1e-05, absolute=1e-50, divergence=10000
left preconditioning
using PRECONDITIONED norm type for convergence test
PC Object: 1 MPI processes
type: ilu
ILU: out-of-place factorization
0 levels of fill
tolerance for zero pivot 2.22045e-14
matrix ordering: natural
factor fill ratio given 1, needed 1
Factored matrix follows:
Mat Object: 1 MPI processes
type: seqaij
rows=2500, cols=2500
package used to perform factorization: petsc
total: nonzeros=12300, allocated nonzeros=12300
total number of mallocs used during MatSetValues calls =0
not using I-node routines
linear system matrix = precond matrix:
Mat Object: 1 MPI processes
type: seqaij
rows=2500, cols=2500
total: nonzeros=12300, allocated nonzeros=12300
total number of mallocs used during MatSetValues calls =0
not using I-node routines
Number of SNES iterations = 6
Number of Linear iterations = 223
Average Linear its / SNES = 3.716667e+01
Why the number of function evaluations is so different?
I'm using last version of PETSc and petsc4py.
Thank you.
Nicola
--
_____________________________________________________________________
Nicola Creati
Istituto Nazionale di Oceanografia e di Geofisica Sperimentale - OGS
IRI (Ricerca Teconologica e Infrastrutture) Department
B.go Grotta Gigante - Brisciki 42/c
34010 Sgonico - Zgonik (TS) - Italy
Tel. +39-0402140213
Fax +39-040327307
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20160208/b54ffcb1/attachment.html>
More information about the petsc-users
mailing list