05 - Model construction#
In ReMKiT1D, terms live in Models, which enclose connected contributions to various equations being solved, as well as any model-bound data - data accessible by default only to the terms enclosed in the relevant model.
In this tutorial we examine how models are constructed, and in particular:
Adding terms to models,
ddtcomponents andTermCollectionsModelpropertiesModelboundDataand using modelbound variables inMatrixTermsModelCollectionsand applying filters to models
[1]:
from RMK_support import Grid, Variable, Model, VarlikeModelboundData, DiagonalStencil, ModelCollection, node, varFromNode
import numpy as np
[2]:
grid = Grid(xGrid = 0.1 * np.ones(16),
interpretXGridAsWidths = True,
vGrid = 0.1 * np.ones(8),
interpretVGridAsWidths = True,
lMax = 3,
)
a,b,c = (Variable(name,grid,data=(i+1)*np.ones(16)) for i,name in enumerate(["a","b","c"]))
The main component of a Model are the ddt entries, representing additive collections of terms that the model contributes to the evolution of individual variables.
For user convenience, ddt is indexed using the evolved variable, and this enables the quick construction of models by defining terms inline and adding them.
Quick note about the diagonal stencil#
Here we will use the built-in DiagonalStencil, which we’ve already used in the previous tutorial. However, we will look at slightly more involved terms, so it’s useful to go over what the stencil does exactly.
If the column variable (the implicit variable of the term) lives on the same grid as the evolved (row) variable of the term (the one on which we use to index
ddt) it simply uses the values of the column variables as they are.If the row and column variables live on different grids, the column variables will be interpolated onto the grid that the evolved (row) variable lives on.
[3]:
newModel = Model("model") # Name of the model
diag = DiagonalStencil()
# inline addition of terms to the term collection ddt[a] and renaming them
newModel.ddt[a] += 2*diag(b).rename("term1") - 0.1*diag(b*c).rename("term2")
newModel.ddt[b] += 2*diag(c).rename("term3") - diag(a*b).rename("term4")
Following the above points about the DiagonalStencil, in term2 above b*c will be projected onto the grid a lives on, and are assumed to both live on the same grid.
NOTE: Terms in tutorials are only meant to be examples of the interface, and do not produce actual relevant code. For Fortran-runnable configurations see the examples.
The above has automatically added 4 terms to the model, and set their evolved variables.
Each element of ddt is a TermCollection, and we can access the terms inside it.
[4]:
print(newModel.evolvedVars) # Evolved variables in the model
print(newModel.ddt[a].termNames)
print(newModel.ddt[b].termNames)
print(newModel.ddt[a]["term1"].evolvedVar.name) # The correct evolved variable was set automatically
print(newModel.ddt[b]["term3"].evolvedVar.name) # The correct evolved variable was set automatically
['a', 'b']
['term1', 'term2']
['term3', 'term4']
a
b
Modelbound data#
ReMKiT1D allows for binding data to Models, giving direct data access only to terms within that model.
MatrixTerms can then include those variables in their column and row variables multiplicatively.
We will showcase modelbound data using the VarlikeModelboundData class, which allows for adding regular derived Variables as modelbound data
[5]:
mbData = VarlikeModelboundData()
d,f,g = (varFromNode(name,grid,(i+1)*node(a)) for i,name in enumerate(["d","f","g"])) # Some dummy variables, normally they would have more complicated derivation rules associated with them
mbData.addVar(d,f,g)
print(mbData.varNames)
print(mbData["d"].name) # Accessors like for VariableContainer
newModel.setModelboundData(mbData) # Now the data is associated with the model
['d', 'f', 'g']
d
To add modelbound variables to column variables they should simply be the second argument to the stencil. This then assumes they live on the same grid as the column/implicit variable of the term.
For row modelbound variables we left multiply the term with modelbound Variables or MultiplicativeArguments made from them using the @ operator.
[6]:
newModel.ddt[a] += diag(c,d*f).rename("mb_term1") # d*f will now be modelbound column variables, so the term will effectively be a*d*f
print(newModel.ddt[a]["mb_term1"].__modelboundC__.argMultiplicity) # modelbound column components
newModel.ddt[b] += (d*g) @ diag(c).rename("mb_term2") # d*f will now be modelbound row variables, so the term will effectively be d*g*c
print(newModel.ddt[b]["mb_term2"].__modelboundR__.argMultiplicity) # modelbound row components
{'f': 1.0, 'd': 1.0}
{'g': 1.0, 'd': 1.0}
NOTE: For the example above, since we’re using the DiagonalStencil, if all variables live on the same grid, it doesn’t matter whether the variables are put in as the column variables or row variables. However, if variables live on different grids, the above examples translate into the following:
mb_term1assumes thatd*flive on the same grid asc, which could be the opposite gridalives on (one of them could be the dual grid and the other the regular grid), and will interpolatec*d*fonto the gridalives on.mb_term2assumes thatd*glive on the same grid asb- the row variable. This means that onlycmight get interpolated if it doesn’t live on the same grid asb, whiled*gwill be used as is.
Filtering terms/models and ModelCollections#
Terms can be grouped into different implicit (can only contain MatrixTerms) and general groups to allow for filtering and selective integration. This can be done inline using the regroup function when adding terms to the model
[7]:
newModel.ddt[a] += c*diag(a).rename("group2_term").regroup(implicitGroups=[2]) # The new term has been added to implicit group 2 only
print(newModel.activeImplicitGroups) # Check that both group 1 and 2 are now active
[1, 2]
We can filter a model to only have terms in some of the groups.
[8]:
filteredModel = newModel.filterByGroup(groups=[2]).rename("filteredModel")
print(filteredModel.evolvedVars) # Only a
print(filteredModel.ddt[a].termNames) # Only the term we regrouped into group 2
['a']
['group2_term']
Similarly, we can also filter models based on evolved variables.
[9]:
filteredModel2 = newModel.onlyEvolving(b).rename("filteredModel2")
print(filteredModel2.evolvedVars) # Only b
print(filteredModel2.ddt[b].termNames) # Only the terms evolving b
['b']
['term3', 'term4', 'mb_term2']
Finally models can be grouped further into ModelCollections, providing accessors for models by name, as well as the same filtering (applied to all contained models).
ModelCollections also let us know the total number of implicit and general groups across all models in the collection, as well as providing utilities for getting all models and terms that evolved any given variable.
[10]:
newCollection = ModelCollection(newModel) # We can initialise the collection with any number of models
# or add them
newCollection.add(filteredModel,filteredModel2)
print(newCollection.modelNames)
print(newCollection["model"].ddt[a].termNames)
print(newCollection.numGroups()) # implicit,general groups
print(newCollection.getTermsThatEvolveVar(a)) # model,term tuples evolving variable a
print(newCollection.getTermsThatEvolveVar(b))
print("\n")
filteredCollection = newCollection.onlyEvolving(a) # Will lose filteredModel2 because that doesn't evolve a
print(filteredCollection.modelNames)
print(filteredCollection["model"].ddt[a].termNames)
print(filteredCollection.numGroups()) # implicit,general groups
print(filteredCollection.getTermsThatEvolveVar(a))
print(filteredCollection.getTermsThatEvolveVar(b))
['model', 'filteredModel', 'filteredModel2']
['term1', 'term2', 'mb_term1', 'group2_term']
(2, 1)
[('model', 'term1'), ('model', 'term2'), ('model', 'mb_term1'), ('model', 'group2_term'), ('filteredModel', 'group2_term')]
[('model', 'term3'), ('model', 'term4'), ('model', 'mb_term2'), ('filteredModel2', 'term3'), ('filteredModel2', 'term4'), ('filteredModel2', 'mb_term2')]
['model', 'filteredModel']
['term1', 'term2', 'mb_term1', 'group2_term']
(2, 1)
[('model', 'term1'), ('model', 'term2'), ('model', 'mb_term1'), ('model', 'group2_term'), ('filteredModel', 'group2_term')]
[]