[petsc-users] petsc4py - Convert c example code to Python
Matthew Knepley
knepley at gmail.com
Tue Apr 2 09:07:11 CDT 2019
On Tue, Apr 2, 2019 at 4:19 AM Nicola Creati via petsc-users <
petsc-users at mcs.anl.gov> 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:
>
In RHS_Function, you never use xdot, but this is an implicit TS. Thus, you
are solving the Laplace
equation, not the heat equation, and you get the exact answer after one
step.
Matt
> 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
> 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.
>
>
--
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
https://www.cse.buffalo.edu/~knepley/ <http://www.cse.buffalo.edu/~knepley/>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20190402/8b16e505/attachment-0001.html>
More information about the petsc-users
mailing list