failure in calculating the solution

Good morning. i would like to report an issue i'm having with a non-newtonian flow and diffusion and convection system:
Specifically, when using either the variable Dof or the variable Dg it fails to provide a solution showing the following:

this happens also when i explicitly wrote the numerical value of the more complex constants. i have also tested all the variables to check for mistakes in writing them down, but they showed no such error. meaning that it should theoretically work.
the linear solver backslash in particular gives me the following errors:

Also, for reasons unknown to me, when i use other fixed diffusivities, regardless of the mesh used, it provides me negative values with the distribution shown in the following image:

and the solver backslash gives me the following errors:

As such, i would like to ask if i'm making any mistake in either defining the variables, a typo that escaped my check and tests or if it's any other issue i failed to notice and if those two issues may be related in some way.
thank you in advance
issue_non_calc_ready_to_present.fea

(Moving here for now as it is not clear it is an issue with the software, but seems more like a modeling and convergence issue.)

It is really hard to say anything with close to 60 parameters and variables. What you would need to do is significantly simplify and reduce the model (including geometry, mesh, parameters, number of physics modes etc.) to a "minimum reproducible example" where it is clear that a specific change or parameter causes an issue (for example when you change parameter "a" from "x" to "y").

good morning.
Thank you for the answer.
issue_non_calc_a_bit_more_presentable.fea

Here's the model with a minimized number of constants.
the error of non-solvability is given with the diffusivity already in the equation while the negative value in the plot shown in the previous post appear with the variable D_mol_liter.
i apologise in advance for not managing to further reduce the number of parameters.
thank you in advance.

Good morning.
I tried to simplify the system as much as possible.
issue_non_calc_as_presentable_as_possible.fea
of the 2 diffusivities present in the constants, the D_mol_liter is the one that gives me the concentration plot with the negative value shown in the first post and that the backslash solver is unable to solve. the other, Dof, is the one that does not give me a solution no matter what, resulting in the following error messages:

an analogous error message is given with any diffusivity (including arbitrary finite numbers) when multiplied by the factor "(1e-3/miu_nn)".
In case of the Dof parameter, i confirmed that the issue lays in the dG parameter as removing it from the equation makes everything work, as proved in the following model
proof_issue_is_dG.fea as removing dG from the parameters makes it technically work.

As such, i would like to ask which may be the issue in the "dG" parameter.
thank you in advance.

I have had a look at your model and can confirm that one problem seem to be in the "Dg" expression which is not valid for "c = 0" (which is the given initial condition).

For example try plotting "Dg" at "t = 0" (and clicking somewhere in the plot to directly evaluate it at specific points and print the expression in the command/log window). If there is no color in the plot this means the values are invalid. Exporting the "fea" struct and manually evaluating (with evalexpr) results in "nan" values (not a number) for "c = 0".

Such issues are typically due to division by or logarithm of zero. It might be helpful to try to use a symbolic tool to simplify the expressions further, as it is still very difficult to see what it is doing. For example in your expression:

Dg: 1-1.1514*(0.5091*1/3*(2*1+1*2^2))*sqrt(c*1000)/(1+(a*k*10^-8/sqrt(c*1000))*sqrt(c*1000))^2+2.303*10^(14*log10(a)-9.75)*c*1000

the factor "1.1514*(0.5091*1/3*(2*1+1*2^2))" can be simplified to just "1.1724", and so on.

After you find this issue, you might also want to have a look at the new section on solution strategy for highly-nonlinear problems https://featool.com/doc/solver#solver_nonlinear.

Thank you for the answer.
I have tried solving the issue by making use of the setnan function to compensate for the NaN value at null concentration and then following the suggestions for non-linear solving.
the results are the following.
corrected_as_resuggested.fea

as it can be seen, it still does not provide me an answer and, moreover, i've also noticed that it may not accept the P2P1 FEM. Since whenever i select it, when i go again to the "equation" settings, it does not show it and instead it shows as if i did not modify it.
May i ask if there are any mistakes in the procedure i followed?

Timmy098 wrote

I have tried solving the issue by making use of the setnan function to compensate for the NaN value at null concentration

I'm not sure this is quite the right approach. The "setnan" function was designed to visualize/plot fields with NaN values, not use in equations or expressions. What you are doing is effectively saying diffusion is zero everywhere, so you have a very difficult to solve pure convection problem, or an equation that just says "0 = 0". It would likely work better if you could (double check and) reformulate your equation and expression so that it is also valid for zero concentrations.

Timmy098 wrote

i've also noticed that it may not accept the P2P1 FEM. Since whenever i select it, when i go again to the "equation" settings, it does not show it and instead it shows as if i did not modify it.

The values in the popup menu always defaults to showing the first item first. It is the values in the edit field that are used in the equations, which correspond to the actual FEM shape functions https://featool.com/doc/dir_e6d29b8517d56e525a271ff5c6e3c994, "sflag2 sflag2 sflag1" for P2/P1 (sflag2 for velocities and sflag1 for pressure).

thank you for the answer.
the way i used the setnan property was to actually give a finite value to the "dG" function whenever it reaches the NaN value and then give a fixed viscosity for c=0 (which is when dG has the NaN value).
About the other issue, doing it this way, still seems to give the solver some issues, especially with the suggested "backslash" solver, while using the "mumps" solver gives me this loop.

 By using the "D_mol_liter" parameter for diffusivity, it results an erratic distribution of negative values of concentration at the initial times.
Moreover, this gets all exacerbated when any diffusivity parameter is multiplied by a "1e-3/miu_nn" factor.

Is there any correlation betwween these two issues?

Timmy098 wrote

Is there any correlation betwween these two issues?

Quite hard to say with such a complex model, without breaking it down and analyzing it in detail.

my assumption for the first one would be that, there being a dependence in dG to the square root of c, any negative value of it would result in a complex number which would break the solver (in fact, when receiving errors in the older simulations it often told me that variables with sqrt(c) in their formulas were complex values). However, i cannot think of anything that would then cause an issue in the specific case of D_mol_liter*1e-3/miu_nn that would prevent the solver from working properly.

Timmy098 wrote

dependence in dG to the square root of c, any negative value of it would result in a complex number

The "built-in" linear solver (UMFPACK) should be able to handle complex values, while "mumps" does not (or ignores the imaginary part). This cloud explain the discrepancy you observe between them. Again, you probably should try to reformulate the problem/expressions so complex values are avoided (in for example turbulence modeling sometimes one solves for the "log(dependent variable)" instead to circumvent issues with negative values).

thank you for the answer.
After checking for dimensional coherence, it does no longer give me negative values when i use D_mol_liter. However, in the diffusivity case dependent on Dg (which is based on the Onsager-Fuoss model with neglected electrophoretic term) it still ends up, in the case of the built-in solver, to give up the usual errors and, in the case of mumps, it gets stuck in an iteration loop.
there are also randomic cases with both solvers where they get stuck in the same iteration loops shown before.
I would also like to add that the main issue is that there are supposed to be no negative values in concentration at any given point in time, especially since, by the way they were shown in the images of my initial post, it appears to be due to computational imprecisions.Screenshot_2025-04-29_183343.png
Screenshot_2025-04-29_202428.png
I have also tried to elevate the order in the case of concentration to P4/Q4 but without success.
Nevertheless we can say that at least one issue has been solved as the D_mol_liter diffusivity no longer gives erratic negative values. Still, when multiplied by the variable "1e-3/miu_nn" still gives errors and/or iteration loops in both backslash and mumps for reasons beyond my current understanding.

a quick update: after rearranging the boundary condition for concentration in a time-dependent way, for reasons to me unknown it shows this peculiar problem where the absoulte error decreases while the relative one increases

this happens in all conditions with all diffusivities, with and without multiplying the diffusivities by the factor "1e-3/miu_nn", regardless of the solver and now also when i set it back to a constant value, despite it working before.
May i ask what might be the wrong input and more precisely the problem that causes this?
attempted_overhaul2.fea
the following file is the one with the boundary condition re-set to a constant value, however the issue appeared when i multiplied it by exp(-t).
thank you in advance