Dear Jose,
Thankyou very much for your response.
On Sat, Oct 26, 2024 at 3:32 PM Jose E. Roman <jro...@dsic.upv.es> wrote:
> Your residue() function is computing something (f) that is not returned to
> the calling function. So evalFunction() returns a zero residual and the
> solver stops in one iteration.
>
Apologies, I am a little confused. The example bratu2d.py imports the
residue function from bratu2dnpy.py and it also does not return anything
(and it works). I was under the impression that
a pointer to `f` is passed to petsc, that overwrites it at each iteration.
Please correct me if I am wrong.
I also tried returning `f` in my residue function but it does not change
anything. Is this what you suggested?
> Mixing numpy arrays and petsc4py vectors is going to be tricky unless you
> know what you are doing. It is simpler to manipulate petsc4py vectors
> directly, see for instance bratu3d.py or other examples.
>
I have also tried this. I created a petsc vector for the input expected
solution and then converted it into a numpy array using the `getArray`
method. But noting changes.
I also printed some information about the convergence itself. I find
ConvergedReason is 2 (absolute tolerance)
getTolerances is [1e-8, 1e-50, 1e-8, 50] (default)
getErrorIfNotConverged is False
getIterationNumber is 0
In contrast, the convergence reason for the bratu2d example was 3 (relative
tolerance).
I also tried increasing the `max_it` but to no avail.
I appreciate your valuable suggestion and advice in this regard.
Thanks and Regards
Jose
>
>
> > El 26 oct 2024, a las 11:11, Vaishak Prasad <vaishakpra...@gmail.com>
> escribió:
> >
> > Hi,
> >
> > Thanks for pointing me to the updated example. The example itself works
> fine.
> >
> > However, on adapting the example to a simple problem ( recovering a
> given 1D vector), I find that there is only one function call happening and
> the solution has not converged. I am attaching the modified example here.
> Is there something silly I am doing here?
> >
> > Any help would be appreciated.
> >
> > Thanks and Regards
> >
> >
> > On Sat, Oct 26, 2024 at 12:42 PM Jose E. Roman <jro...@dsic.upv.es>
> wrote:
> > The examples in those old slides might be outdated. Try with the
> examples that come with the PETSc repo or tarball. In particular, that
> example can be found in src/binding/petsc4py/demo/legacy/bratu2d/bratu2d.py.
> >
> > Jose
> >
> > > El 25 oct 2024, a las 19:58, Vaishak Prasad <vaishakpra...@gmail.com>
> escribió:
> > >
> > > Dear All,
> > > Greetings from India.
> > >
> > > I am new to petsc4py. I was trying out the example on slides 18, and
> 19 at
> > >
> > > https://urldefense.us/v3/__https://www.bu.edu/pasi/files/2011/01/Lisandro-Dalcin-petsc4py.pdf__;!!G_uCfscf7eWS!bzL5mw0mqx-jt3BIhPqHNtxZUf469dbmQwxugUupECLT7CEMbk_vj-TZ_wEym_uN1O4CWKJKeODMn0a0tO72ky2PkUrV$
> > >
> > >
> > > However, I get a zero vector for the solution `x`.
> > >
> > > I am using petsc4py 3.21.3 on Linux.
> > >
> > > Any help would be appreciated.
> > >
> > > With regards
> > >
> > > --
> > > Vaishak Prasad
> > > Post Doctoral Fellow
> > > Astrophysical Relativity Group
> > > International Center for Theoretical Sciences
> > > Tata Institute of Fundamental Research
> > > Survey No. 151, Shivakote, Hesaraghatta Hobli,
> > > Bengaluru, Karnataka, India,
> > > Earth, Solar system,
> > > Milky Way Galaxy, The Local Group
> >
> >
> >
> > --
> > Vaishak Prasad
> > Post Doctoral Fellow
> > Astrophysical Relativity Group
> > International Center for Theoretical Sciences
> > Tata Institute of Fundamental Research
> > Survey No. 151, Shivakote, Hesaraghatta Hobli,
> > Bengaluru, Karnataka, India,
> > Earth, Solar system,
> > Milky Way Galaxy, The Local Group
> > <mwe.py>
>
>
--
Vaishak Prasad
Post Doctoral Fellow
Astrophysical Relativity Group
International Center for Theoretical Sciences
Tata Institute of Fundamental Research
Survey No. 151, Shivakote, Hesaraghatta Hobli,
Bengaluru, Karnataka, India,
Earth, Solar system,
Milky Way Galaxy, The Local Group
import sys, petsc4py
import numpy as np
petsc4py.init(sys.argv)
from petsc4py import PETSc
N=10
expected_soln_vec_npy = 100*np.random.randn(N)
def residue(expected_soln_vec_npy, x, f):
''' Residue calculator '''
print("Call")
f = (x - expected_soln_vec_npy)**2
print("Residue:", f)
return f
# this user class is an application
# context for the nonlinear problem
# at hand; it contains some parameters
# and knows how to compute residuals
class SimpleEqn:
def __init__(self, m, expected_soln_vec_petsc, impl='python'):
self.m = m
self.expected_soln_vec_petsc = expected_soln_vec_petsc
if impl == 'python':
order = 'c'
elif impl == 'fortran':
order = 'f'
else:
raise ValueError('invalid implementation')
self.compute = residue
self.order = order
def evalFunction(self, snes, X, F):
m = self.m
expected_soln_vec_petsc = self.expected_soln_vec_petsc
order = self.order
x = X.getArray(readonly=1).reshape(m, order=order)
f = F.getArray(readonly=0).reshape(m, order=order)
expected_soln_vec_npy = expected_soln_vec_petsc.getArray(readonly=1).reshape(m, order=order)
self.compute(expected_soln_vec_npy, x, f)
# convenience access to
# PETSc options database
OptDB = PETSc.Options()
m = OptDB.getInt('m', N)
impl = OptDB.getString('impl', 'python')
expected_soln_vec_petsc = PETSc.Vec().createSeq(m)
expected_soln_vec_petsc.setArray(expected_soln_vec_npy)
# create application context
# and PETSc nonlinear solver
appc = SimpleEqn(m, expected_soln_vec_petsc, impl)
snes = PETSc.SNES().create()
# register the function in charge of
# computing the nonlinear residual
f = PETSc.Vec().createSeq(m)
snes.setFunction(appc.evalFunction, f)
# configure the nonlinear solver
# to use a matrix-free Jacobian
snes.setUseMF(False)
snes.getKSP().setType('cg')
snes.setFromOptions()
snes.setTolerances(max_it=100)
# solve the nonlinear problem
b, x = None, f.duplicate()
x.set(0)
snes.solve(b, x)
cr = snes.getConvergedReason()
print('Converged Reason:', cr)
ncr = snes.getErrorIfNotConverged()
print("Not converged?", ncr)
tols = snes.getTolerances()
print("Tols:", tols)
snes.getIterationNumber()
