# Gridap.jl: Solving non-linear coupled PDEs

**URL:** <https://discourse.julialang.org/t/gridap-jl-solving-non-linear-coupled-pdes/106326>\
**Category:** Numerics\
**Tags:** gridap\
**Created:** [November 16, 2023, 4:35pm UTC](https://discourse.julialang.org/t/gridap-jl-solving-non-linear-coupled-pdes/106326 "2023-11-16T16:35:56Z")\
**Posts on this page:** 1\
**Page:** 1

<div class="post-metadata">

**Author:** ![charvey](https://avatars.discourse-cdn.com/v4/letter/c/7feea3/32.png) [@charvey](https://discourse.julialang.org/u/charvey)\
**Post date:** [November 16, 2023, 4:35pm UTC](https://discourse.julialang.org/t/gridap-jl-solving-non-linear-coupled-pdes/106326/1 "2023-11-16T16:35:56Z")

</div>

We are having some difficulties solving two coupled multi-variate non-linear PDEs using Gridap. The various tutorials haven’t quite filled the gaps in our knowledge.

We have two PDEs that we need to solve together. Their residuals are:

```julia
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Ω

```

Note that `uh_in` and `sh_in` are the fields at the end of the last increment.  
(We have Jacobians too, but with the automatic differentation, let’s not complicate matters yet!)

As per the tutorials, to solve one non-linear PDE, we can use

```julia
op_Disp = FEOperator(res_Disp, U_Disp, V0_Disp)
nls_Disp = NLSolver(show_trace=true, method=:newton, linesearch=BackTracking(), ftol = ftol, iterations = 4)
solver_Disp = FESolver(nls_Disp)
uh_out, = solve!(uh_in, solver_Disp, op_Disp)

```

where `U_Disp` and `V0_Disp` are the trial and test spaces, respectively.

Referring to [Tutorial 11](https://gridap.github.io/Tutorials/dev/pages/t011_fsi_tutorial/) on FSI, it suggests that we need to add our residuals, and solve like this:

```julia
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 = 5)
solver = FESolver(nls)
out = solve(solver, op)

```

This at least runs; however, it prompts two questions:

1. Is it correct to add residuals in this instance? I was expecting to solve the two residual functions in a 2D Newton’s method. If not, what is the correct way, please?

2. If it is the correct way, how can we specify the two initial fields (i.e. from the end of the last increment)? We can’t get any variation of the following to work:

```julia
out = solve([uh_in, sh_in], solver, op)

```

Many thanks!