Hi Barry,
Thank you for your assistance . Since this is an ongoing development
project, I would prefer not to share the full code publicly at this stage
(public petsc-users). I’d be glad to share it with you privately via email,
if that's an option you’re open to.
Thanks again for your support
thanks
Art
El lun, 16 jun 2025 a las 20:51, Barry Smith (<bsm...@petsc.dev>) escribió:
>
> Can you please post the entire code so we can run it ourselves to
> reproduce the problem.
>
> I suspect the code maybe using a dense matrix to represent the
> Jacobian. This ill require a huge amount of memory. Sundials is likely
> using Jacobian free Newton Krylov to solve the linear system. You can do
> this also with PETSc but it might not be defaulting to it. With the code I
> can quickly check what is happening.
>
> Barry
>
>
> > On Jun 16, 2025, at 2:05 PM, Art <mac3...@gmail.com> wrote:
> >
> > Hi everyone,
> >
> > I’m porting a code from scikits.odes.sundials.cvode to PETSc, since it
> is compatible with FEniCSx and can run in parallel with MPI. First, I used
> Petsc with the "rk" solver, and it worked well, both serially and in
> parallel for a system with 14000 nodes (42000 dofs). However, when using
> an implicit solver like bdf, the solver takes up all the memory (16 gb),
> even on a small system. To do this, use this:
> >
> >
> > def ifunction(self, ts, t, y, ydot, f):
> >
> >
> y.ghostUpdate(PETSc.InsertMode.INSERT_VALUES,PETSc.ScatterMode.FORWARD)
> > y.copy(result=self.yv.x.petsc_vec)
> > self.yv.x.scatter_forward()
> >
> > dydt = self.rhs(self.yv)
> > dydt.x.scatter_forward()
> >
> > ydot.copy(result=f)
> > f.axpy(-1.0, dydt.x.petsc_vec)
> >
> > return 0
> >
> >
> > y = y0.petsc_vec.copy()
> > ts.setType(ts.Type.BDF)
> > ts.setIFunction(ifunction)
> > ts.setTime(0.0)
> > ts.setTimeStep(1e-14)
> > ts.setStepLimits(1e-17,1e-12)
> > ts.setMaxTime(1.0e-12)
> > ts.setExactFinalTime(PETSc.TS.ExactFinalTime.STEPOVER)
> > ts.setTolerances(rtol=1e-6, atol=1e-6)
> > snes = ts.getSNES()
> > ksp = snes.getKSP()
> > ksp.setType("gmres")
> > ksp.getPC().setType("none")
> > ksp.setFromOptions()
> >
> > For the scikits.odes.sundials.cvode library, in serial mode, I have
> used:
> >
> > solver = CVODE(rhs,
> > old_api=False,
> > linsolver='spgmr',
> > rtol=1e-6,
> > atol=1e-6,
> > max_steps=5000,
> > order=2)
> >
> > In this case, the solver worked perfectly and obtained similar results
> to the rk solver in PETSC. I suspect the issue might be related to the way
> the Jacobian is built in PETSC, but scikits.odes.sundials.cvode works
> perfectly without requiring the Jacobian. I would greatly appreciate any
> suggestions or examples on how to properly set up the BDF solver with PETSc.
> > Thanks
> > Art
>
>
