I'm working on a project where I have two adv-diff coupled domains through their respective source terms (one domain adds mass, the other subtracts mass). For brevity, I'm modeling them in steady state. The equations are your standard advection-diffusion transport equation with a source term look like this:
$$ \frac{\partial c_1}{\partial t} = 0 = \mathcal{F}_1 + \mathcal{Q}_1(c_1,c_2) \\ \frac{\partial c_2}{\partial t} = 0 = \mathcal{F}_2 + \mathcal{Q}_2(c_1,c_2) $$
Where $\mathcal{F}_i$ is diffusive and advective flux for species $i$, and $\mathcal{Q}_i$ is the source term for species $i$.
I have been able to write a solver for my problem using the Newton-Raphson method, and have completely coupled the two domains using a block mass matrix, ie:
$$ F_{coupled} = \left[\begin{array}{c c} A_1 & 0 \\ 0 & A_2 \\ \end{array}\right]\underbrace{ \left[\begin{array}{c} c_{1,i} \\ c_{2,i} \\ \end{array}\right] }_{x_i} - \left[\begin{array}{c} b_1(c_{1,i}, c_{2,i}) \\ b_2(c_{1,i}, c_{2,i}) \\ \end{array}\right] $$
The term $F_{coupled}$ is used to determine the Jacobian matrix and update both $c_1$ and $c_2$:
$$ \mathcal{J}(x_i) \left[x_{i+1} - x_i \right] = -\mathcal{F}_{coupled} $$
or
$$ x_{i+1} = x_i - \left(\mathcal{J}(x_i)\right)^{-1} \mathcal{F}_{coupled} $$
To speed things up, I don't calculate the Jacobian every iteration - right now I'm playing with every five iterations, which has seemed to work well enough and keep the solution steady.
The problem is: I'm going to be moving to a larger system where both domains are in 2D/2.5D, and calculating the Jacobian matrix is going to quickly deplete my available computer resources. I'm building this model to be used in an optimization setting later on, so I also can't be behind the wheel at every iteration tuning the damping factor, etc..
Am I right to be looking elsewhere for a more robust and algorithm for my problem, or is this is as good as it gets? I've looked a bit into Quasi-linearization, but am not sure how applicable it is w.r.t. to my system.
Are there any other slick algorithms that I may have missed that can solve a system of nonlinear equations without resorting to re-calculating the Jacobian as offen?