Your python implementation of the residual function and the Jacobian function is intended for solving ODEs in the form x' = f(t,x). So you should use ts.setRHSFunction() and ts.setRHSJacobian() instead of the implicit version.
Hong > On Apr 2, 2019, at 3:18 AM, Nicola Creati via petsc-users > <[email protected]> wrote: > > Hello, I have run the two codes using: -pc_type lu -ts_monitor -snes_monitor > -ksp_monitor -ts_max_steps 5. The codes start to diverge after the first time > step: > Python run : > > 0 TS dt 0.01 time 0. > 0 SNES Function norm 2.327405179696e+02 > 0 KSP Residual norm 1.939100379240e+00 > 1 KSP Residual norm 1.013760235597e-15 > 1 SNES Function norm 9.538686413612e-14 > 1 TS dt 0.01 time 0.01 > 0 SNES Function norm 8.565210784140e-14 > 0 KSP Residual norm 1.266678408555e-15 > 1 KSP Residual norm 2.663404219322e-31 > 1 SNES Function norm 5.528402779439e-29 > 2 TS dt 0.01 time 0.02 > 0 SNES Function norm 5.829656072531e-29 > 0 KSP Residual norm 3.805606380832e-31 > 1 KSP Residual norm 3.283104975114e-46 > 1 SNES Function norm 1.867771204856e-44 > 3 TS dt 0.01 time 0.03 > 0 SNES Function norm 1.644541571708e-44 > 0 KSP Residual norm 1.969573456719e-46 > 1 KSP Residual norm 1.297976290541e-61 > 1 SNES Function norm 1.780215821230e-59 > 4 TS dt 0.01 time 0.04 > 0 SNES Function norm 1.700894451833e-59 > 5 TS dt 0.01 time 0.05 > > ex13.c run: > > 0 TS dt 0.01 time 0. > 0 SNES Function norm 2.327405179696e+02 > 0 KSP Residual norm 9.158498731719e-01 > 1 KSP Residual norm 5.129005976821e-16 > 1 SNES Function norm 1.141297389434e-13 > 1 TS dt 0.01 time 0.01 > 0 SNES Function norm 9.158498731719e+01 > 0 KSP Residual norm 4.156347391374e-01 > 1 KSP Residual norm 1.996081836894e-16 > 1 SNES Function norm 3.609761734594e-14 > 2 TS dt 0.01 time 0.02 > 0 SNES Function norm 4.156347391374e+01 > 0 KSP Residual norm 2.203298795311e-01 > 1 KSP Residual norm 9.047356740098e-17 > 1 SNES Function norm 3.292648977280e-14 > 3 TS dt 0.01 time 0.03 > 0 SNES Function norm 2.203298795311e+01 > 0 KSP Residual norm 1.367458644503e-01 > 1 KSP Residual norm 5.055030661919e-17 > 1 SNES Function norm 1.272304808486e-14 > 4 TS dt 0.01 time 0.04 > 0 SNES Function norm 1.367458644503e+01 > 0 KSP Residual norm 9.695393464180e-02 > 1 KSP Residual norm 1.918401734795e-17 > 1 SNES Function norm 1.130247628824e-14 > 5 TS dt 0.01 time 0.05 > > I rechecked the Python code against the C one but I did not find the source > of the issue. > > # 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() > localU = da.getLocalVec() > > da.globalToLocal(X, localU) > > 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) > > uarray = localU.getArray(readonly=1).reshape(8, 8, order='C') > f = F.getArray(readonly=0).reshape(8, 8, order='C') > > (xs, xm), (ys, ym) = da.getRanges() > 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] = uarray[i,j] > continue > u = uarray[i,j] > uxx = (-2.0 * u + uarray[i, j-1] + uarray[i, j+1]) * sx > uyy = (-2.0 * u + uarray[i-1, j] + uarray[i+1, j])* sy > f[i, j] = uxx + uyy > F.assemble() > #PETSc.Log.logFlops(11.0*xm*ym) > > > 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.setOption(PETSc.Mat.Option.NEW_NONZERO_ALLOCATION_ERR, False) > 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, j+1), sy), > ((i, j), -2*sx - 2*sy)]: > col.index = index > col.field = 0 > B.setValueStencil(row, col, value) > > B.assemble() > #B.view() > 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() > #U.view() > > > 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, > stencil_width=1, dof=1) > > 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) > > da.setMatType('aij') > J = da.createMat() > ts.setIJacobian(Jacobian_func, J, J) > > 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. >
