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!crzDSnMPvGJ19hJ6wP30I5DQnV7XMkJmxebYeLV9MFjpkLL2d5VtXHhHE2si93mzUmxzKxeVo0dkoyx-Cd04hRPUqrAd$
> >  
> >
> > 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
import sys, petsc4py
petsc4py.init(sys.argv)
from petsc4py import PETSc


# Number of components of the exact solution
N=10

# Create the exact solution vector
expected_soln_vec = 100*np.random.randn(N)



def residue(expected_soln_vec,  x, f):
    ''' Residue equation implementation '''
    print("Call")
    f = (x - expected_soln_vec)**2
    print("Residue:", 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, impl='python'):
        self.m = m
        self.expected_soln_vec = expected_soln_vec
        
        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 = self.expected_soln_vec
        order = self.order
        x = X.getArray(readonly=1).reshape(m, order=order)
        f = F.getArray(readonly=0).reshape(m, order=order)

        self.compute(expected_soln_vec, x, f)


    
# convenience access to
# PETSc options database
OptDB = PETSc.Options()
m = OptDB.getInt('m', N)
impl  = OptDB.getString('impl', 'python')

# create application context
# and PETSc nonlinear solver
appc = SimpleEqn(m, expected_soln_vec, 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(True)
snes.getKSP().setType('cg')
snes.setFromOptions()


# solve the nonlinear problem
b, x = None, f.duplicate()
x.set(0) # zero initial guess
snes.solve(b, x)

Reply via email to