[petsc-users] petsc4py - Convert c example code to Python

Sat Mar 30 17:17:33 CDT 2019

   Nicola,

     What goes wrong with your code? Does it crash? Not converge in the same way as the C code? Produce wrong answers?

     If you think specifically the Jacobian is wrong you can run both codes with -mat_view to see any differences in the Jacobian
     For the computation of the right hand side you can call VecView() and compare it to the same output from the C code.

     I would debug by running each code next to each other for a very small problem looking at intermediate values computed and try to find the first location where the two are producing different results.

  Barry


> On Mar 28, 2019, at 4:17 AM, Nicola Creati via petsc-users <petsc-users at mcs.anl.gov> wrote:
> 
> Hello, I'm trying to write in Python one of the TS PETSc tutorial example:
> 
> https://www.mcs.anl.gov/petsc/petsc-current/src/ts/examples/tutorials/ex13.c.html 
> 
> I'm not able to fill the Jacobian matrix in the right way in Python. Maybe there are some other conversion errors.
> Might someone help, please?
> 
> This is my code:
> 
> # Example 13 petsc TS
> import sys, petsc4py
> petsc4py.init(sys.argv)
> 
> from petsc4py import PETSc
> from mpi4py import MPI
> import numpy as np
> import math
> 
> def RHS_func(ts, t, X, xdot, F):
>     da = ts.getDM()
> 
>     mx, my = da.getSizes()
> 
>     hx, hy = [1.0/(m-1) for m in [mx, my]]
>     sx = 1.0/(hx*hx)
>     sy = 1.0/(hy*hy)
> 
>     x = X.getArray(readonly=1).reshape(8, 8, order='C')
>     f = F.getArray(readonly=0).reshape(8, 8, order='C')
> 
>     (xs, xm), (ys, ym) = da.getRanges()
>     aa = np.zeros((8,8))
>     for j in range(ys, ym):
>         for i in range(xs, xm):
>             if i == 0 or j == 0 or i == (mx-1) or j == (my-1):
>                 f[i,j] = x[i,j]
>                 continue
>             u = x[i,j]
>             uxx = (-2.0 * u + x[i, j-1] + x[i, j+1]) * sx
>             uyy = (-2.0 * u + x[i-1, j] + x[i+1, j])* sy
>             f[i, j] = uxx + uyy
>     F.assemble()
> 
> def Jacobian_func(ts, t, x, xdot, a, A, B):
>     mx, my = da.getSizes()
>     hx = 1.0/(mx-1)
>     hy = 1.0/(my-1)
>     sx = 1.0/(hx*hx)
>     sy = 1.0/(hy*hy)
> 
>     B.zeroEntries()
>     (i0, i1), (j0, j1) = da.getRanges()
>     row = PETSc.Mat.Stencil()
>     col = PETSc.Mat.Stencil()
> 
>     for i in range(j0, j1):
>         for j in range(i0, i1):
>             row.index = (i,j)
>             row.field = 0
>             if i == 0 or j== 0 or i==(my-1) or j==(mx-1):
>                 B.setValueStencil(row, row, 1.0)
>             else:
>                 for index, value in [
>                     ((i-1, j),   sx),
>                     ((i+1, j),   sx),
>                     ((i,   j-1), sy),
>                     ((i-1, j+1), sy),
>                     ((i,   j),  -2*sx - 2*sy)]:
>                     col.index = index
>                     col.field = 0
>                     B.setValueStencil(row, col, value)
> 
>     B.assemble()
>     if A != B: B.assemble()
>     return PETSc.Mat.Structure.SAME_NONZERO_PATTERN
> 
> def make_initial_solution(da, U, c):
>     mx, my = da.getSizes()
>     hx = 1.0/(mx-1)
>     hy = 1.0/(my-1)
>     (xs, xm), (ys, ym) = da.getRanges()
> 
>     u = U.getArray(readonly=0).reshape(8, 8, order='C')
> 
>     for j in range(ys, ym):
>         y = j*hy
>         for i in range(xs, xm):
>             x = i*hx
>             r = math.sqrt((x-0.5)**2+(y-0.5)**2)
>             if r < (0.125):
>                 u[i, j] = math.exp(c*r*r*r)
>             else:
>                 u[i, j] = 0.0
>     U.assemble()
> 
> def monitor(ts, i, t, x):
>     xx = x[:].tolist()
>     history.append((i, t, xx))
> 
> history = []
> nx = 8
> ny = 8
> da = PETSc.DMDA().create([nx, ny], stencil_type= PETSc.DA.StencilType.STAR)
> da.setFromOptions()
> da.setUp()
> 
> u = da.createGlobalVec()
> f = u.duplicate()
> 
> c = -30.0
> 
> ts = PETSc.TS().create()
> ts.setDM(da)
> ts.setType(ts.Type.BEULER)
> 
> ts.setIFunction(RHS_func, f)
> ts.setIJacobian(Jacobian_func)
> 
> ftime = 1.0
> ts.setMaxTime(ftime)
> ts.setExactFinalTime(PETSc.TS.ExactFinalTime.STEPOVER)
> 
> make_initial_solution(da, u, c)
> dt = 0.01
> ts.setMonitor(monitor)
> ts.setTimeStep(dt)
> ts.setFromOptions()
> ts.solve(u)
> 
> ftime = ts.getSolveTime()
> steps = ts.getStepNumber()
> 
> Thanks.
> 
> Nicola
> 
> -- 
> 
> Nicola Creati
> Istituto Nazionale di Oceanografia e di Geofisica Sperimentale - OGS www.inogs.it
> Dipartimento di Geofisica della Litosfera Geophysics of Lithosphere Department
> CARS (Cartography and Remote Sensing) Research Group http://www.inogs.it/Cars/
> Borgo Grotta Gigante 42/c 34010 Sgonico - Trieste - ITALY ncreati at ogs.trieste.it
> off.   +39 040 2140 213
> fax.   +39 040 327307
> 
> _____________________________________________________________________
> This communication, that may contain confidential and/or legally privileged information, is intended solely for the use of the intended addressees. Opinions, conclusions and other information contained in this message, that do not relate to the official business of OGS, shall be considered as not given or endorsed by it. Every opinion or advice contained in this communication is subject to the terms and conditions provided by the agreement governing the engagement with such a client. Any use, disclosure, copying or distribution of the contents of this communication by a not-intended recipient or in violation of the purposes of this communication is strictly prohibited and may be unlawful. For Italy only: Ai sensi del D.Lgs.196/2003 - "T.U. sulla Privacy" si precisa che le informazioni contenute in questo messaggio sono riservate ed a uso esclusivo del destinatario.
>