[petsc-users] petsc4py - Convert c example code to Python
Nicola Creati
ncreati at inogs.it
Wed Apr 3 03:07:09 CDT 2019
Hello,
thanks to all, I used the wrong RHS and Jacobian function in my
conversion. I modified the code and I got the right solution. I'm
attaching the final code conversion. Someone could be interested in a
Python version of ex13 of the TS component tutorial.
"""
Python version of the ex13.c PETSc TS component at:
https://www.mcs.anl.gov/petsc/petsc-current/src/ts/examples/tutorials/ex13.c.html
"""
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, U, F, *args):
da = ts.getDM()
localU = da.getLocalVec()
da.globalToLocal(U, 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()
def Jacobian_func(ts, t, U, A, B, *args):
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()
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()
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.setRHSFunction(RHS_func, f)
da.setMatType('aij')
J = da.createMat()
ts.setRHSJacobian(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.setTimeStep(dt)
ts.setFromOptions()
ts.solve(u)
ftime = ts.getSolveTime()
steps = ts.getStepNumber()
Cheers.
Nicola
On 02/04/2019 16:54, Zhang, Hong wrote:
> 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 <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:
>> 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.
>>
--
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.
More information about the petsc-users
mailing list