Hi @Francyrad,
I will anwer your questions below, but let me start with a general comment. I understand implementing a numerical model can be stressful, in particular if it is tied to a proposal or thesis and there is a lot of time pressure. However, please refrain from using ALL CAPS in your messages and demand people’s attention to solve a problem. No one in this forum is paid to answer your questions, the people reading your comments are volunteers and take a part of their time that they could have used for their own research problems to reply to your messages. Demanding an answer will only lead to others not replying to your messages.
That being said, here is how I understand your problem at the moment:
I think your pull request 6168 correctly fixes a bug in the dynamic core plugin for the case of a fully molten or fully frozen core. These cases were apparently not tested by the original author. Unfortunately the changes I made in 6177 and 6187 to improve that plugin led to some conflicts with your PR. The version of the PR you updated is still correct, but also still contains some conflicts with the main branch. I have opened a new PR 6224 in which I took your original commit and modified it so it works with the current ASPECT and also added two new tests for fully molten and fully frozen cores. The reason I waited so long to do this was that you didnt reply to my question here for how to reproduce the problem you are having. Without being able to reproduce the problem it can easily take days to even understand the problem of another scientist and spending my time this way is simply not efficient.
While creating the test cases for the new PR I also noticed that there are two more bugs for the fully frozen case (a division by zero). I fixed those too. Note, that I am not working on core models, and have not compared the implementation of these cases with the original Nimmo paper. Therefore, as is the case for most open-source software: There is no guaranty that this implementation correctly reproduces the paper. It is the responsibility of every user to make sure that the software works for the application case you are considering. I think you are aware of this, since you already tried to reproduce the figure by the paper above.
Now to your question about input parameters: This is indeed confusing and I am not sure what the correct solution is. The input parameter is documented as being in the unit of 1/TPa, and 1/TPa^2 respectively, but if I check the code (e.g. in the get_solidus function) the value is actually multiplied with pressure in GPa, so maybe it should be converted into 1/GPa? Can you give this a try, and if it is indeed incorrectly documented we should fix the documentation (and possibly wrong values).
Best,
Rene