PETSc input format for linear solvers

Question

I’m going through some considerable effort to translate one of my codes from MATLAB. It’s a type of finite element code and I haven’t implemented the solver yet but comparing CPU times for simply assembling the global system between MATLAB and C, the latter is more than two orders of magnitude faster...

Anyways, I have no idea how to use PETSc but will attempt to learn it in the next few weeks. Just reading over some code and documentation, it seems like it might be pretty hard to translate my code for PETSc.

I noticed that there’s an option to read a matrix from a “file” (see here). I’m just wondering if I go through with translating my code as I’m doing now, can I just save matrix/rhs to a .dat file or something then read from there? I don’t have all that much experience with C so I don’t know if this would be difficult to do or not, especially for larger systems like I’ll have.

In my code, I’ll have to solve a linear system at every Newton iteration (it’s a nonlinear system), and combine this with the fact that the problem is time-dependent. The only time in my entire code that I’ll need PETSc is for the linear solve (admittedly, it’s a lot of linear solves...).

If I’m going to have to learn the syntax for PETSc now, I’d just like to know so I’m not frustrated later.

On another note, I’m really surprised how much faster the code is in C so far. Maybe my MATLAB coding isn’t that efficient, but I’ve been keeping the for loops at a minimum where possible.

Answer 1

I would recommend looking at a time-dependent example because PETSc can provide a lot more diagnostics if you formulate at that level. For example, you could use a Rosenbrock method, an additive Runge-Kutta IMEX method, or others. Some involve Newton iteration, but that is not required and is not always the best approach. To use those methods, you can start by just providing your residual function, $f(t,u)$, as in $\dot u = f(t,u)$, or alternatively, $g$ in $g(t,u,\dot u) = 0$. PETSc can automatically compute a Jacobian via coloring, and if you later decide to assemble the Jacobian for efficiency, PETSc can check whether what you assembled was correct and show you the errors. If you roll your own time integration and your own Newton, these algorithmic variants and debugging features will not be available unless you write them yourself.

The Users Manual explains the interfaces. Also see the manual pages for specific API documentation. The man pages are cross-referenced to examples. A good example to start with would be the Brusselator src/ts/examples/tutorials/ex25.c. You may also find slides from tutorials useful. If you get stuck, feel free to ask here or on the petsc-users mailing list.

Your proposal to use the file system will be a huge bottleneck even in serial. It complicates the workflow greatly and will limit the methods that you can use. But if you just want to test a matrix once, you can use PetscBinaryWrite() from MATLAB and load the system and solve using src/ksp/ksp/examples/tutorials/ex10.c.