New to Julia gridap

**Personal context**

I am a phd student who works in numerical characterisation of electrochemical performance of fuel cells/electrolysis cells. This involves resolution of a system of coupled, non linear, time dependant system of PDE using finite element method. The lab I work in has a tradition of using COMSOL and matlab with Livelink (the latter to write some scripts to automate certain things such as resolving a system of pde for one given parameter set). Lately this workflow proves to be tedious and less effective.

And therefore, I am trying to transition to julia to make the workflow less cumbersome.

**Reason for this post**

I have trouble finding right gridap framework/tutorial that suits my needs. I have an impression that the available tutorials are not suitable to my problem. I will describe the problem in a simplified version, without going into the details in the next section and in the ensuing paragraph, I try to explain why I get the above written impressions.

**Problem simplified**

The physical problem in itself is complex and is founded on a hypothetical reaction mechanism occurring at the electrodes. I focus on the reactions occuring at the oxygen electrode (doi: 10.1149/1945-7111/abf40a). This hypothetical reaction mechanism implies the presence of several species. The transport and mass conservation equations of these species result in a system of poisson equations that has to be solved simultaneously. The physical quantities we would like to calculate as a function of electrode geometry and time are the following: normalised surface concentrations \theta_i of assumed surface species i, ionic potential \varphi_{io}, electronic potential \varphi_{el}, oxygen concentration in the bulk of the material C_O, oxygen gas concentration y_{O_2}. The poisson equations of all these physical quantities will be written as follows (combine table I and table II in the above mentioned DOI link):

\Delta \varphi_{io} = f(C_O, \theta_1, t)

\Delta C_O = g(C_O, \theta_1) + \frac{\partial C_O}{\partial t}

\Delta \theta_2 = h(\theta_1, \theta_2, \theta_3) + \frac{\partial \theta_2}{\partial t}

\Delta \theta_1 = k(\theta_1, C_O) + \frac{\partial \theta_1}{\partial t}

\Delta \theta_3 = l(\theta_2, \theta_3)+ \frac{\partial \theta_3}{\partial t}

\Delta y_{O_2} = m(\theta_3) + \frac{\partial y_{O_2}}{\partial t}

\Delta \varphi_{el} = n((\theta_1, C_O) + \frac{\partial \varphi_{el}}{\partial t}

I have the boundary conditions of each of these variables. For simplicity I will leave this out.

I have looked into the tutorials that are available in gridap framework. All problems that involves poisson equations have an independant function (of space or both space and time) on the right hand side and itâ€™s mostly user defined. The tutorial that fits the closest to what I want is in tutorial 18. Transient non linear equation and it follows the same pattern of having an user defined, known functions on the other side of the equality. Whereas in my case, all the terms are coupled to one another.

I would like to know whether I am missing something in my understanding. If anyone could provide a direction to pursue, it would be very helpful. Any other suggestions for libraries to explore are also welcome!