[petsc-users] petsc4py snes behavior
Matthew Knepley
knepley at gmail.com
Mon Feb 8 05:02:44 CST 2016
On Mon, Feb 8, 2016 at 2:41 AM, Nicola Creati <ncreati at inogs.it> wrote:
> 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:
>
You did not provide the Jacobian function, resulting in a finite difference
approximation that uses more function evaluations.
Matt
> ## # 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:
>
> 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):
>
> 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
>
>
--
What most experimenters take for granted before they begin their
experiments is infinitely more interesting than any results to which their
experiments lead.
-- Norbert Wiener
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20160208/ce6c4f81/attachment-0001.html>
More information about the petsc-users
mailing list