# Methods for solving discrete PDEs using algorithmic differentiation results

I'm looking for a method to solve a 20000 variable, 20000 residual non-linear PDE with a Galerkin method.

I have Fortran subroutines for:

• The residuals: $$\vec{r}(\vec{x})$$;

• Their Jacobian multiplied by a direct algorithmic derivative seed: $$\dot{r}(\dot{x})=\frac{\partial \vec{r}(\vec{x})}{\partial \vec{x}}\cdot\dot{x}$$; and

• Their transposed Jacobian multiplied by a reverse algorithmic derivative seed:

$$\bar{x}(\bar{r})=\left[\frac{\partial \vec{r}(\vec{x})}{\partial \vec{x}}\right]^t\cdot \bar{r}\, .$$

All these subroutines are run in no more than 0.05 s each. Due to the high number of variables, however, I cannot fully compute the Jacobian for a Newton-Raphson iteration because the run time would become impractical.

Is there any out-of-the-shelf method I can use in this case to solve the equation without fully calculating the Jacobian?

## 2 Answers

Is there any out-of-the-shelf method I can use in this case to solve the equation without fully calculating the Jacobian?

Using any integrator for stiff ODEs where the implicit equation is solved via Newton-Krylov methods will only require Jacobian-vector products. So something like DASKR, IDA, etc. can be made to do this.

So the size of the problem doesn't necessarily stop you from using the full Jacobian if you store it in a compressed/sparse storage format, you also should be able to know which residuals are dependent on which unknowns which should make it cheaper as well. That said, just use Jacobian free Newton-Krylov (JFNK) solvers. FGMRES is a really strong linear solver that uses only matrix vector products. Some simpler ones to implement could be BiCGStab (which isn't as strong though imo). Generalized Conjugate Residual (GCR) is also really strong and HANIM out of NASA and NIA has shown some great results with low storage GCR methods with strong preconditioners and line search technology. You won't be able to use HANIM itself, but you can find the papers on it and copy some parts pretty easily.

• The key piece is "you also should be able to know which residuals are dependent on which unknowns". The $i$th residual only depends on variable $j$ if shape functions $i$ and $j$ overlap. So the cost doesn't increase with the number of variables square, but just with the number of variables. – Wolfgang Bangerth Jul 3 '20 at 17:09