Refactor functions in modular_temp_timestepping.py?

[andreww]

I've now looked at the code in a bit more detail and think the main place to look to make things more flexible / reusable are the large functions in modular_temp_timestepping.py. There are several functions in there with quite similar (but not identical) code but also which do things like define internal functions for K(T) and dK(T)/dT. I'm pretty sure that some moving around of code will lead to more code reuse and give users the ability feed in their own functions (functions can be assigned to variables in python). It's also probably worth pulling out the core time-stepping code as that should be optimisable if needed.

Anyway, I think this needs some discussion and planning. I don't think this is a particularly big job, but getting it right should mean that it's possible to do the next project with (much) of the same code. That's probably important!

[murphyqm]

Quick sketch of code structure

[murphyqm]

Quick draft of potential new core code:

Will test later today

class Core:
    """ Core class represents and manipulates core temp and latent heat """

    def __init__(self, temp, melt, r, rho, cp, k, maxlh, lat=0):
        """ Create a new core with temperature and latent heat """
        self.temperature = temp
        self.latent = lat
        self.melting = melt
        self.radius = r
        self.density = rho
        self.heatcap = cp
        self.cmbk = k
        self.maxlatent = maxlh

    def __str__(self):
        return"Core at {0} K. Latent heat extracted: {1}".format(self.temperature, self.latent)

    def cool(self):
        """Test cooling method"""
        self.temperature -= 100.0
        self.latent += 10.0

    def cooling(self, mantletemps, timestep, dr):
        if (self.temperature > self.melting) or (
        self.latent >= self.maxlatent):
            self.temperature = self.temperature - (
            3.0
            * self.cmbk
            * ((mantletemps[0, i] - mantletemps[1, i]) / dr)
            * timestep) / (self.density * self.heatcap * self.radius)
        else:
            self.latent = (
                self.latent
                + (4.0 * np.pi * self.radius ** 2)
                * self.cmbk
                * ((mantletemps[0, i] - mantletemps[1, i]) / dr)
                * timestep)

[andreww]

The overall structure looks pretty good to me.

You may want to add additional objects for the physical properties of the mantle. These would be passed into "temp timestepping" and would have methods to return the properties given a temperature (and maybe a pressure and time for future extensibility). What I'm not sure about is if you should have a single "mantle physical properties" object, or different objects for the conductivity, density, and heat capacity. I would be tempted to bundle them together so you would need four (I think) methods called as something like properties.k(temp), properties.dkdt(temp), properties.density(temp) and properties.heat_capacity(temp). I think it would be useful if these take scalars or arrays as arguments (I can show you how to do this with numba) and you would need two implementations - one for constant values and one (or more) for the variable models. Probably best to do this via inheritance.

Thinking about modules, the three columns in the sketch give (I think) a starting point. You can probably have a module that covers much of the material on the left hand side (a set up function, a timestepping function, and an analysis function along with something that manages the overall flow), a module to deal with analysis of output (so some of the middle column) and a module to deal with plotting (the right column). Physical properties probably belong in another module, as does the core object.

[andreww]

Some thoughts on the core object below. In general I think it's a good idea to think about how some other code would "use" the object. Mostly this is going to be inside the mantle time stepper. This needs to do two things (in general)

1. Find the temperature at the top of the core. This should probably look like an attribute so you would do something like: temp[i,-1] = core.boundary_temperature.
2. Extract some heat from the core. This would be a method (coolingbelow) but I think it's best to pass in a amount of heat removed (probably in W along with the timestep duration) rather than expect the core to know how to calculate this. This would move half of the calculation in cooling out of the core object but it would keep core stuff in the core and mantle stuff in the mantle (which seems like the right place for them). For what it's worth I would call the method core.extract_heat(power, time). You'll need to document the direction of heat flow very carefully (what is negative).

Anyway, the way to decide if this is sensible is to think about another possible core object that does something different. You should be able to think about this with exactly the same attributes and methods (used by the mantle) so that all the changes are internal to the new core object. From memory Hack had a slightly different core (and Claire has something more complex in the preprint too).

Beyond the core time stepping we should decide if the core object needs to interact with the analysis stuff or if the mantle is responsible for keeping track of things to plot. I would put that stuff in the core (and make the core keep its own history) but that is something to discuss.

A few more detailed comments below.

class Core:

Have you decided if this is supposed to run under Python 2.7 or are you Python 3.x only? I ask because of the difference in the default way classes work has changed (a little bit). In 2.7 class Core is an "old style" class and to make a "new style" class you need to do class Core(object). In python 3 all objects behave like new style classes and (by default) inherent from object. class Core(object) is thus portable (and, I think, idiomatic). It's what I tend to do but the cool kids may have moved on and accepted the default is the default and python 2.7 is dead. It's probably worth thinking about future more complex cores. I would call this IsothermalEutecticCore or IsothermalSimpleCore I think.

""" Core class represents and manipulates core temp and latent heat """

def __init__(self, temp, melt, r, rho, cp, k, maxlh, lat=0):
    """ Create a new core with temperature and latent heat """
    self.temperature = temp
    self.latent = lat
    self.melting = melt
    self.radius = r

Is it worth having an inner and outer radius (for now, even if the inner radius is always 0). self.density = rho self.heatcap = cp self.cmbk = k This belongs "in" the mantle self.maxlatent = maxlh This can be calculated when the object is set up (given the volume etc.

I would store a 1D array of times (record the time whenever heat is extracted) and temperatures. For an isothermal core this makes it possible to dump out all temperature radius information at the end. def str(self): return"Core at {0} K. Latent heat extracted: {1}".format(self.temperature, self.latent)

def cool(self):
    """Test cooling method"""
    self.temperature -= 100.0
    self.latent += 10.0

Tests are good, but I would make this kind of thing a unit test (set up a model, extract heat, check the temperature). def cooling(self, mantletemps, timestep, dr): See above. I would avoid passing mantletemps and dr in here. if (self.temperature > self.melting) or ( self.latent >= self.maxlatent): self.temperature = self.temperature - ( 3.0

• self.cmbk
• ((mantletemps[0, i] - mantletemps[1, i]) / dr)
• timestep) / (self.density self.heatcap self.radius) else: self.latent = ( self.latent
  • (4.0 np.pi self.radius ** 2)
  • self.cmbk
  • ((mantletemps[0, i] - mantletemps[1, i]) / dr)
  • timestep)

Anyway, a very solid start.

[murphyqm]

Have uploaded the bits I was working on from earlier today (written before you added these notes, so these aren't incorporated yet) - it's reproducing the actual unit tests I had applied to the current core function and works with a slightly modified mantle timestepping function. These are really useful pointers, thanks - will incorporate them on Monday.

[murphyqm]

Functions have been reformatted.