We are going to implement a simple linear decay equation:
Where T is the temperature in Kelvin, a > 0 is the linear constant and F the forcing parameter.
To set up the model, some straight-forward steps are necessary. First, import the numericalmodel
module:
import numericalmodel
First initialize an object of NumericalModel
:
model = numericalmodel.numericalmodel.NumericalModel()
We may tell the model to start from time 0:
model.initial_time = 0
The StateVariable
, Parameter
, and ForcingValue
classes all derive from InterfaceValue
, which is a convenient class for time/value management. It also provides interpolation (InterfaceValue.__call__
). An InterfaceValue
has a (sensibly unique) id
for it to be referencable in a SetOfInterfaceValues
.
For convenience, let's import everything from the interfaces
submodule:
from numericalmodel.interfaces import *
Let's define our state variable. For the simple case of the linear decay equation, our only state variable is the temperature T:
temperature = StateVariable( id = "T", name = "temperature", unit = "K" )
Providing a name and a unit documents your model on the fly.
Tip
All classes in numericalmodel
are subclasses of ReprObject
. This makes them have a proper __repr__
method to provide an as-exact-as-possible representation. So at any time you might do a print(repr(temperature))
or just temperature<ENTER>
in an interactive python session to see a representation of this object:
- numericalmodel.interfaces.StateVariable(
time_function = numericalmodel.utils.utcnow, values = array([], dtype=float64), name = 'temperature', times = array([], dtype=float64), unit = 'K', id = 'T', interpolation = 'zero' )
The others - a and F - are created similarly:
parameter = Parameter( id = "a", name = "linear parameter", unit = "1/s" )
forcing = ForcingValue( id = "F", name = "forcing parameter", unit = "K/s" )
Now we add them to the model:
model.variables = SetOfStateVariables( [ temperature ] )
model.parameters = SetOfParameters( [ parameter ] )
model.forcing = SetOfForcingValues( [ forcing ] )
Tip
A lot of objects in numericalmodel
also have a sensible __str__
method, which enables them to print a summary of themselves. For example, if we do a print(model)
:
###
### "numerical model"
### - a numerical model -
### version 0.0.1
###
by: anonymous
This is a numerical model.
initial time: 1490648927.012074
"temperature"
--- T [K] ---currently: ? [K] interpolation: zero 0 total recorded values
"linear parameter"
--- a [1/s] ---currently: ? [1/s] interpolation: linear 0 total recorded values
"forcing parameter"
--- F [K/s] ---currently: ? [K/s] interpolation: linear 0 total recorded values
none
Note
When an InterfaceValue
's value
is set, a corresponding time is determined to record it. The default is to use the return value of the InterfaceValue.time_function
, which in turn defaults to the current utc timestamp. When the model was told to use temperature
, parameter
and forcing
, it automatically set the InterfaceValue.time_function
to its internal model_time
. That's why it makes sense to define initial values after adding the InterfaceValue
s to the model.
Now that we have defined our model and added the variables, parameters and forcing, we may set initial values:
temperature.value = 20 + 273.15
parameter.value = 0.1
forcing.value = 28
Tip
We could also have made use of SetOfInterfaceValues
' handy indexing features and have said:
model.variables["T"].value = 20 + 273.15 model.parameters["a"].value = 0.1 model.forcing["F"].value = 28
We proceed by defining our equation. In our case, we do this by subclassing PrognosticEquation
, since the linear decay equation is a prognostic equation:
class LinearDecayEquation(numericalmodel.equations.PrognosticEquation):
"""
Class for the linear decay equation
"""
def linear_factor(self, time = None ):
# take the "a" parameter from the input, interpolate it to the given
# "time" and return the negative value
return - self.input["a"](time)
def independent_addend(self, time = None ):
# take the "F" forcing parameter from the input, interpolate it to
# the given "time" and return it
return self.input["F"](time)
def nonlinear_addend(self, *args, **kwargs):
return 0 # nonlinear addend is always zero (LINEAR decay equation)
Now we initialize an object of this class:
decay_equation = LinearDecayEquation(
variable = temperature,
input = SetOfInterfaceValues( [parameter, forcing] ),
)
We can now calculate the derivative of the equation with the derivative
method:
>>> decay_equation.derivative()
-28.314999999999998
Alright, we have all input we need and an equation. Now everything that's missing is a numerical scheme to solve the equation. numericalmodel
ships with the most common numerical schemes. They reside in the submodule numericalmodel.numericalschemes
. For convenience, we import everything from there:
from numericalmodel.numericalschemes import *
For a linear decay equation whose parameters are independent of time, the EulerImplicit
scheme is a good choice:
implicit_scheme = numericalmodel.numericalschemes.EulerImplicit(
equation = decay_equation
)
We may now add the scheme to the model:
model.numericalschemes = SetOfNumericalSchemes( [ implicit_scheme ] )
That's it! The model is ready to run!
Running the model is as easy as telling it a final time:
model.integrate( final_time = model.model_time + 60 )
The model results are written directly into the StateVariable
's cache. You may either access the values directly via the values
property (a numpy.ndarray
) or interpolated via the InterfaceValue.__call__
method.
One may plot the results with matplotlib.pyplot
:
import matplotlib.pyplot as plt
plt.plot( temperature.times, temperature.values,
linewidth = 2,
label = temperature.name,
)
plt.xlabel( "time [seconds]" )
plt.ylabel( "{} [{}]".format( temperature.name, temperature.unit ) )
plt.legend()
plt.show()
The full code can be found in the examples
section.