Hi everyone,
I am currently transitioning our lab’s electrochemical modeling workflow from COMSOL to Gridap.jl. My goal is to move away from proprietary licenses and leverage the flexibility of the Julia ecosystem.
The system I am modeling is a 1D stationary system with 6 coupled nonlinear PDEs (poisson equations). The unknowns are potentials and species concentrations. In the future, I would need to dive into transient terms as well. For now, stationnary problem would do.
The coupling is highly nonlinear due to Butler-Volmer reaction rates acting as source terms in the porous electrodes.
The Problem: When calling the solver, the REPL hangs indefinitely (10+ minutes) with no output of any sorts. Even with autodiff=false and a very low maxiter, I am seeing a “black screen” during the solving phase.
Current Implementation: I am using a MultiFieldFESpace. To troubleshoot, I have attempted to bypass Automatic Differentiation to rule out compilation hangs:
op = FEOperator(res, U, V)
ls = LUSolver()
nls = NLSolver(ls; show_trace=true, method=:newton, linesearch=BackTracking(), atol=1e-8, rtol=1e-10, maxiter=1, autodiff=false)
solver = FESolver(nls)
println("\nsolving phase")
u_solution = solve(solver, op)
where res is the residual function that I defined, U and V are of following types.
typeof(U) MultiFieldFESpace{Gridap.MultiField.ConsecutiveMultiFieldStyle, Gridap.FESpaces.UnConstrained, Vector{Float64}}
typeof(V) MultiFieldFESpace{Gridap.MultiField.ConsecutiveMultiFieldStyle, Gridap.FESpaces.UnConstrained, Vector{Float64}}
I verfied that the res function is properly implemented as well, using the following code, which didn’t cause any errors and returns DomainContribution()
u_test = interpolate_everywhere([1.0, 1.0, 1.0, 1.0, 1.0, 1.0], U)
v_test = get_fe_basis(V)
println("running residual for testing: ", res(u_test, v_test))
And I am using this particular solve(), where I am giving only the nonlinear solver (FESolver) and the FEoperator. I would need a method to be able to define the initial values manually that are currently defaulted to zero in this method. I have tried solve!(u_test, solver, op). But the aforementioned problem persists.
Any pointers to what I am missing out, to any documentation are welcome.
I am currently using Julia 1.10.10 and Gridap 0.19.17. I had to use julia 1.10 because I ran into some compatibility issues installing Gridap and ForwardDiff in the latest julia. I had to downgrade to 1.10. Kindly, let me know if you require any further information.