<div dir="ltr">Hi Barry, <br><br>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.<br><br>Thanks again for your support<div><br></div><div>thanks </div><div><br></div><div>Art</div></div><br><div class="gmail_quote gmail_quote_container"><div dir="ltr" class="gmail_attr">El lun, 16 jun 2025 a las 20:51, Barry Smith (<<a href="mailto:bsmith@petsc.dev">bsmith@petsc.dev</a>>) escribió:<br></div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><br>
   Can you please post the entire code so we can run it ourselves to reproduce the problem.<br>
<br>
    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.<br>
<br>
  Barry<br>
<br>
<br>
> On Jun 16, 2025, at 2:05 PM, Art <<a href="mailto:mac3bar@gmail.com" target="_blank">mac3bar@gmail.com</a>> wrote:<br>
> <br>
> Hi everyone,<br>
> <br>
> 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:<br>
> <br>
> <br>
> def ifunction(self, ts, t, y, ydot, f):<br>
> <br>
>       y.ghostUpdate(PETSc.InsertMode.INSERT_VALUES,PETSc.ScatterMode.FORWARD)<br>
>       y.copy(result=self.yv.x.petsc_vec)<br>
>       self.yv.x.scatter_forward()<br>
> <br>
>       dydt = self.rhs(self.yv)    <br>
>       dydt.x.scatter_forward()<br>
> <br>
>       ydot.copy(result=f) <br>
>       f.axpy(-1.0, dydt.x.petsc_vec) <br>
> <br>
>       return 0<br>
> <br>
> <br>
> y = y0.petsc_vec.copy()<br>
> ts.setType(ts.Type.BDF)<br>
> ts.setIFunction(ifunction)    <br>
> ts.setTime(0.0)<br>
> ts.setTimeStep(1e-14)<br>
> ts.setStepLimits(1e-17,1e-12)<br>
> ts.setMaxTime(1.0e-12)<br>
> ts.setExactFinalTime(PETSc.TS.ExactFinalTime.STEPOVER)<br>
> ts.setTolerances(rtol=1e-6, atol=1e-6) <br>
> snes = ts.getSNES()<br>
> ksp = snes.getKSP()<br>
> ksp.setType("gmres")       <br>
> ksp.getPC().setType("none")    <br>
> ksp.setFromOptions()<br>
> <br>
> For the scikits.odes.sundials.cvode library, in serial mode, I  have used:<br>
> <br>
> solver = CVODE(rhs,<br>
>                old_api=False,<br>
>                linsolver='spgmr',<br>
>                rtol=1e-6,<br>
>                atol=1e-6,<br>
>                max_steps=5000,<br>
>                order=2)<br>
> <br>
> 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.<br>
> Thanks<br>
> Art<br>
<br>
</blockquote></div>