It would be nice, wouldn't it. I spent a bit of time seeing if I could reproduce the 'reference' model of Haack et al. (1990). It looks like this is non-trivial but could be done (if we allow some API breaking changes). (Not so) brief notes below.
The 'reference' model is quite complex. One issue is that it has source terms in the mantle. In principle these are easy to add (it's just another term to add to the temperature change here:
|
temperatures[j, i] = temperatures[j, i - 1] + A_1 + B_1 + C_1 |
Assuming this defaults to zero and we can add to the API without loosing backwards compatibility this would be okay.
The values to feed into the source term are a bit more complex. There are two parts. A radiogenic component (which is different in the 'crust' and the 'mantle' - the only difference between the two regions) and a latent heat component (the model starts off above the solids for the mantle). Ultimately this is just a function of time, depth and temperature (which would return a heat generation rate, that we would divide by the heat capacity and multiply by the time step before adding this to the change in temperature).
A second issue is that we'll need a more complex core model. It is not totally clear to me from looking at the paper if the core is isothermal or if it is part of a single conduction problem. I think it is the latter. This wouldn't be a problem (it's just a different set of material properties and a different source term for the latent heat). This time there is no radiogenic heating, just a different latent heat term. I'm not sure how to handle this once it cools to the eutectic temperature (the model is cooling through a two phase region). Maybe this would need a separate core (with its own conductive solution and a zero flux boundary condition in the middle). Grrr.
Then the third issue (that actually made me stop) is that as far as I can see the density is no specified anywhere in the paper.
It would be nice, wouldn't it. I spent a bit of time seeing if I could reproduce the 'reference' model of Haack et al. (1990). It looks like this is non-trivial but could be done (if we allow some API breaking changes). (Not so) brief notes below.
The 'reference' model is quite complex. One issue is that it has source terms in the mantle. In principle these are easy to add (it's just another term to add to the temperature change here:
pytesimal/pytesimal/numerical_methods.py
Line 380 in d166aba
Assuming this defaults to zero and we can add to the API without loosing backwards compatibility this would be okay.
The values to feed into the source term are a bit more complex. There are two parts. A radiogenic component (which is different in the 'crust' and the 'mantle' - the only difference between the two regions) and a latent heat component (the model starts off above the solids for the mantle). Ultimately this is just a function of time, depth and temperature (which would return a heat generation rate, that we would divide by the heat capacity and multiply by the time step before adding this to the change in temperature).
A second issue is that we'll need a more complex core model. It is not totally clear to me from looking at the paper if the core is isothermal or if it is part of a single conduction problem. I think it is the latter. This wouldn't be a problem (it's just a different set of material properties and a different source term for the latent heat). This time there is no radiogenic heating, just a different latent heat term. I'm not sure how to handle this once it cools to the eutectic temperature (the model is cooling through a two phase region). Maybe this would need a separate core (with its own conductive solution and a zero flux boundary condition in the middle). Grrr.
Then the third issue (that actually made me stop) is that as far as I can see the density is no specified anywhere in the paper.