Solving non-linear coupled PDEs

We are having some difficulties solving two coupled multi-variate non-linear PDEs using Gridap.

We have two PDEs that we need to solve together in increments for the vector field `uh` and the scalar field `sh`. Their residuals are:
````
res_PF(s,ϕ) = ∫( Gc(tags)*ls*∇(ϕ)⋅ ∇(s) + 2*ψPlusPrev_in*s*ϕ + (Gc(tags)/ls)*s*ϕ )*dΩ - ∫( (Gc(tags)/ls)*ϕ )*dΩ
res_Disp(u,v) = ∫( (ε(v) ⊙ (σ_mod∘(ε(u),ε(uh_in),sh_in)) ) )*dΩ
````

The following function solves one increment starting from the initial conditions, `uh_in` and `sh_in`.

````
function step(uh_in,sh_in,vApp,ψPlusPrev_in,ftol)
uApp1(x) = VectorValue(0.0,0.0)
uApp2(x) = VectorValue(0.0,vApp)
uApp3(x) = VectorValue(0.0,-vApp)
U_Disp = TrialFESpace(V0_Disp,[uApp1, uApp1, uApp2, uApp3])

residual((u,s),(v,ϕ)) = ∫( (ε(v) ⊙ (σ_mod∘(ε(u),ε(uh_in),sh_in)) ) )*dΩ + ∫( Gc(tags)*ls*∇(ϕ)⋅ ∇(s) + 2*ψPlusPrev_in*s*ϕ + (Gc(tags)/ls)*s*ϕ )*dΩ - ∫( (Gc(tags)/ls)*ϕ )*dΩ
op = FEOperator(residual, MultiFieldFESpace([U_Disp; U_PF]), MultiFieldFESpace([V0_Disp; V0_PF]))
nls = NLSolver(show_trace=true, method=:newton, linesearch=BackTracking(), ftol = ftol, iterations = 10)
solver = FESolver(nls)
initial = Gridap.MultiField.MultiFieldFEFunction([get_free_dof_values(uh_in); get_free_dof_values(sh_in)], MultiFieldFESpace([V0_Disp; V0_PF]), [uh_in; sh_in])
out, = solve!(initial, solver, op)

return out, norm(Gridap.Algebra.residual(op, out), Inf)
end
````

As you can see, our approach is to add the residuals from each PDE, combine the trial spaces and test spaces into `MultiFieldFESpace`s, and to combine `uh_in` and `sh_in` into `MultiFieldFEFunction`s.

Unfortunately, although the model runs, the results are not correct against benchmarks analytical solutions.

I should point out that our previous approach used the "splitting method" to solve the system of PDEs, in which one equation is solved at a time, and the solution fed into the next equation, and cycled until the residual of both is below the tolerance. We found that this method was not very efficient and hard to converge, and hoped that this new method based on the fully implicit (aka monolithic residual) approach would work better. The only difference between the codes is in the `step` function above, so we are confident the error is somewhere here.

Can anyone see what is going wrong?

Start with the step function and its residual, MultiFieldFESpace construction, and NLSolver setup. Reproduce the coupled solve against the analytical benchmarks, then compare its results and residuals with the previous splitting-method implementation. Done means identifying and correcting the cause of the incorrect monolithic results while preserving convergence below the requested tolerance.

Written by the indexing model from the issue text.