# 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?

• 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