Add Two State-Variable Capability

[jrleeman]

Add the capability to run the rate and state model with two state variables, commonly used for more complex friction evolution.

[jrleeman]

Ok. I can think of a few ways to do this:

1) Have some state holder that sets a 1 state variable solution when b2 and dc2 are None, and sets a 2 state variable solution for when they are set. I don't like that we end up having basically two integration routines that are switched between in the solve method though. Seems like a lot of duplication, but it is easy to follow... sort of. Maybe something like:

        if self.n_state == 1:
            self.mu, self.theta, self.v = w
            self.v = self.vref * exp((self.mu - self.mu0 - self.b *
                                      log(self.vref * self.theta / self.dc)) / self.a)

            # Find the loadpoint_velocity corresponding to the most recent time
            # <= the current time.
            loadpoint_vel = self.loadpoint_velocity[self.model_time<=t][-1]
            dmu_dt = self.k * (loadpoint_vel - self.v)
            dtheta_dt = self.stateLaw(self)

            return [dmu_dt,dtheta_dt]

        elif self.n_state == 2:
            self.mu, self.theta, self.theta2, self.v = w
            self.v = self.vref * exp((self.mu - self.mu0 - self.b *
                                      log(self.vref * self.theta / self.dc) - self.b2 * log(self.vref * self.theta2 / self.dc2)) / self.a)

            # Find the loadpoint_velocity corresponding to the most recent time
            # <= the current time.
            loadpoint_vel = self.loadpoint_velocity[self.model_time<=t][-1]
            dmu_dt = self.k * (loadpoint_vel - self.v)
            dtheta1_dt = self.stateLaw(self)
            dtheta2_dt = self.stateLaw2(self)

2) Set b2 to zero and Dc2 to any arbitrary value. The zeroing of b2 will mean the second (now dummy) state variable will not contribute to the solution. We could even have the dtheta_2/dt return unity. The downside here is that we are integrating over an extra dummy variable and slow down the works in the process. Maybe not enough to worry about.

I'm leaning more towards the second approach currently, but an open to suggestion.

[dopplershift]

Would you be using two different state functions together (like Dietrich and Ruina) or would this be another standalone state function that uses two state variables?

[dopplershift]

I will say that as you generalize things, that's where the design questions start to be answered. While you can certainly go with the if block, as outlined above (and much scientific software does), it's not a pretty design. If you want to add a 3 state version, you'll have to come back and modify this if block. In general, what I'd shoot for is:

1. Be able to add more state laws/types later without modifying this block
2. Encapsulate fully the concept of a "state law"

The first thing that absolutely jumps out is when you start adding sets of things, where for every foo you now have a foo2. The first question that we need to decide is if the right abstraction is always a single state law (with 1 or more state variables) or is it better to model this as 1 or more state laws.

[jrleeman]

Each state variable would be solved using one of the common relations (Dietrich, Ruina, PRZ, etc), but with it's own respective Dc,b. Maybe some kind of encapsulation is what we need. I can imagine someone adding a third state variable. Not saying that I'd condone it. This means that we need to associate b and Dc with the state law formulation and then allow the integrator to take an arbitrary number of state laws.

[dopplershift]

Sounds like a StateLaw class initialized with b and Dc. In addition to some kind of dtheta_dt function, you might also want to capture the calculation of its contribution to the velocity function.

As an aside, why do you include v in your initial conditions (w0)? You don't return a derivative for it, nor do you ever use the value passed, nor the result of integrating it.

[jrleeman]

Ok. I'll make a pass at that. I just pushed a fix that removes that extra passing of v. Must have been left over junk from when I was first attempting this.

[jrleeman]

Okay, say we make our simple StateRelation class like this:

class StateRelation(object):
    """
    Encapsulates the state relation formulation.
    """

    def __init__(self,relation):
        self.b = None
        self.Dc = None
        self.state = None
        self.state_evolution = relation

We currently have the built in state relations as separate functions. Not sure that needs to/should change. The catch is that a user passed state relation could use any properties of the model (v, a, etc). Again, not sure why they would, but let's say they want to. Should we pass the RateState object to StateRelation or is there a more elegant OO way to do this?

[dopplershift]

You could make the case that the RateState object is just supposed to handle the integration and plotting, and that you should pull the velocity and dmu_dt calculations, as well as the parameters, into some other class (ExternalSystem ??).

I also think it makes sense to take the built-in state relations and make them classes:

class StateRelation(object):
    def __init__(self):
        self.b = None
        self.dc = None

    def velocity_component(self, system):
         return self.b * log(system.vref * system.theta / self.dc) 

class DieterichState(StateRelation):
    def evolve_state(self, system):
        return 1. - system.v * system.theta / self.dc

[jrleeman]

Something like this? I've played with it this afternoon, but I think I'm getting things too dependent on each other. From what we do to this, I can try to get the example running again.

    class StateRelation(object):
        def __init__(self, relation):
            self.b = None
            self.Dc = None
            self.state = None

        def velocity_componet(self, system):
            return self.b * log(system.vref * self.state / self.Dc)

    class DieterichState(StateRelation):
        def _set_steady_state(self, system):
            self.state = self.Dc/system.vref

        def evolve_state(self, system):
            if self.state == None:
                self.state = _steady_state(self, system)
            # return dtheta/dt
            return 1. - system.v * self.state / self.Dc

    class ExternalSystem(object):
        def __init__(self):
            # Rate and state model parameters
            self.mu0 = None
            self.a = None
            self.k = None
            self.v = None
            self.vref = None
            self.state_relations = []

        def velocity_evolution(self):
            for state in self.state_relations:
                v_contribution = state.velocity_componet(self)
                self.v = self.vref * exp((self.mu - self.mu0 - v_contribution) / self.a)
                print self.v

        def friction_evolution(self):
            return self.k * (loadpoint_vel - self.v)

[dopplershift]

Overall concept looks like what I thought--it doesn't look too coupled to me. I don't think velocity_evolution() is correct based on the original changes you made for the two state system above. How about:

        def velocity_evolution(self):
            v_contribution = 0.
            for state in self.state_relations:
                v_contribution += state.velocity_componet(self)
            self.v = self.vref * exp((self.mu - self.mu0 - v_contribution) / self.a)
            print self.v

And by all means, feel free to disregard anything I say. :)

[jrleeman]

I think these are all good suggestions! The coupling I was concerned about is when implementing this. There is a branch (StateLaws) that runs, but is ugly... also provides the wrong answer. I'll have to clean it up, but maybe it's in the right direction again? This is super educational for OO stuff!

[jrleeman]

Still a few things in plotting to fix, but I'm getting the correct answer now. Maybe this is getting close to a PR!

[dopplershift]

Don't let that hold up making a PR. Once you make a PR, and pushes you make to the branch will show up in the PR. The benefit of the PR is that I can comment on the code itself directly. You can also rebase your branch (to squash/remove commits), do a force push, and the PR will update.

That's my long-winded way of advising you to just make the PR.