Solving non-linear coupled PDEs

[christopherfear]

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 MultiFieldFESpaces, and to combine uh_in and sh_in into MultiFieldFEFunctions.

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?

[JordiManyer]

@christopherfear When you are creating your initial condition initial, you are using

MultiFieldFESpace([V0_Disp; V0_PF])

which is your test FESpace, i.e it has zero dirichlet boundary conditions. I would probably try to change that to your trial FESpace, i.e

initial = Gridap.MultiField.MultiFieldFEFunction([get_free_dof_values(uh_in); get_free_dof_values(sh_in)], MultiFieldFESpace([U_Disp; U_PF]), [uh_in; sh_in])

I also think it would be safer (at least to debug) to interpolate instead of creating your own object:

initial = interpolate([uh_in,sh_in],MultiFieldFESpace([U_Disp; U_PF]))

[christopherfear]

@JordiManyer

Thanks for the suggestion, I had tried this before, and checked again after you commented but still no luck