Converting data tables from entropy to pressure

[djbower]

In scripts/ there is a convert_datatables.py prototype that eventually should convert the original C Spider lookup tables by entropy to lookup by temperature. The script currently maps entropy to temperature and plots the data (note that the plot is misleading because it suggests a regular grid, but in fact it interpolates outside the actual data bounds which are revealed by the contours). Depending how the lookup is implemented, it may be necessary to interpolate the mapped pressure-temperature lookup data to a regular grid to allow fast and consistent indexing. Prototype output in data/1TPa-dk09-elec-free/temperature shows that temperature is not regularly spaced.

It would probably be most sensible to somehow mask the data so that regions outside data have NaNs, whilst still allowing a regular grid format. Hence this script requires a few modifications to output the T-P data. Also note that the prototype temperature output is not currently generated by the script, but rather older code that is below the new class structure.

[timlichtenberg]

Also relevant: https://github.com/ExPlanetology/aragog/issues/7

[planetmariana]

@lsoucasse @timlichtenberg @djbower

The translation from entropy to temperature into a regular grid is finished. The points where the interpolation is out of the bounds were filled with NaNs. You can check the comparison plots below. The only thing is I couldn't resolve the small features near the surface when projecting the interpolation into the regular grid. The functions that receives any pressure array and returns the solidus temperature according different parametrizations and the fit the original P-T anhydrous solidus from Aragog are given in the file solidus_functions.py in the temperature_regular_grid folder in the OSF repository.

[planetmariana]

@lsoucasse @djbower @timlichtenberg

Using the branch "add-output" with the last changes I was able to plot the different quantities with constant properties, but when introducing the lookup tables with regular grids and NaNs as showed above, at least one by one for testing purposes. It complains for arrays with NaNs or infs, I was debugging to find where was punctually the problem. Turns out, when computing dTdt, is calling the quantity capacitance that is computed through the -staggered density- (density_melt.dat was the first lookup table I tried) and somehow it's computing these properties at the nodes with out-of bounds values where the NaNs are. That's my idea about the problem. Do you have any insight on how to avoid this?

[djbower]

Depending on how the boundary conditions are implemented, it may be that the code is trying to calculate values at the innermost and outermost boundaries of the mesh. Due to numerical noise or zero values, the values at these points might not be evaluated correctly. Notably, P is zero or maybe even slightly negative at the uppermost boundary, so perhaps this is propagating an error into the calculation (NaN). It's often convenient to evaluate all properties at all grid points and then apply boundary conditions afterwards, but maybe it will be necessary to only evaluate properties at the innermost nodes of the basic mesh (excluding the top and bottom boundaries). More debugging might be required to identify exactly which boundary is throwing the NaN.

[planetmariana]

@djbower @timlichtenberg @lsoucasse.

I managed to put all together the lookup tables only for P<100 GPa. I changed the NaNs discussed above for values 1e-15. The problem is, when setting a final end_time, let's say 1000 the simulation stops at 300 (image below), like it crashes but still no error is showed in the terminal, so it's a bit difficult to know where is crashing punctually. Is there some way to see a .log while running? I included these tables into PROTEUS with the latest update of a dummy model with Aragog but takes so long to complete the first iteration. PROTEUS .log shows that is running and that's all. I changed the tolerances but also didn't worked. This simulation took around 40mins, and it's just the beginning of the solidification. There's only one warning that shows:

_RuntimeWarning: invalid value encountered in sqrt self._inviscid_velocity[self._isconvective] = np.sqrt(

Which could be related with the whole problem. But, I still cannot find where those variables are directly related to the lookup tables.

[djbower]

@planetmariana Test the latest main and let me know if everything seems in order. Then we can look to remove the old entropy-pressure data tables which are no longer required.

[planetmariana]

@djbower I tested the latest main but when running mixed phase it complains on output.py when declaring:

    def pressure_GPa_basic(self) -> npt.NDArray:
        return self.evaluator.mesh.basic.eos.pressure * self.parameters.scalings.pressure * 1.0e-9

AttributeError: 'FixedMesh' object has no attribute 'eos'

[djbower]

Fixed in 8650910 (I believe)

[planetmariana]

@djbower I tested for solid and liquid, it's working and they look different but still negative values in the thermal expansivity for the mixed phase producing that the code crashes in the first time step:

Solid phase:

Liquid phase:

Mixed phase: