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
> <[email protected]> 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
> [email protected]
> 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.
>
