04 - Matrix Terms#

ReMKiT1D equations are constructed from terms, and the main kind of term used is the MatrixTerm.

A matrix term represents contributions to equations of the form

\[\frac{\partial n}{\partial t} = Mu\]

where \(n\) is the evolved variable and both \(n\) and \(u\) are implicit variables (NOTE: if the evolved variable is stationary, the LHS of the above equation is set to 0). The matrix elements of \(M\) are decomposed as

\[M_{ij} = cTF_iR_iS_{ij}C_j\]

where \(c\) is a constant scalar component, \(T\) is a TimeSignalData component encoding periodic time dependence, \(F_i\) is the Profile components, encoding explicit \(X\),\(V\), and \(H\) dependence, \(S_{ij}\) is the Stencil, which contains the structure of the operator, and \(R_i\) and \(C_i\) are MultiplicativeArgument objects, i.e. products of variables raised to powers (similar to SimpleDerivation from the previous tutorial).

Except for the Stencil all other MatrixTerm components are optional.

While MatrixTerm objects can be constructed directly using a constructor, from v2.0.0, ReMKiT1D offers a more intuitive way to construct these terms, which will be demonstrated in this tutorial.

[1]:

from RMK_support import Grid, Variable, DiagonalStencil
import RMK_support.stencils as st

import numpy as np

We start by preparing some variables

[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"]))

To start creating a MatrixTerm we need a stencil. Many stencils are offered, and the user is encouraged to browse the examples and documentation. Here we focus on the DiagonalStencil and StaggeredDivStencil as examples.

[3]:

diag = DiagonalStencil()
div = st.StaggeredDivStencil()

We can create a term by acting on a Variable or a MultiplicativeArgument. If acting only on a Variable, that variable is taken as \(u\) in the above equation (i.e. the implicit variable). If acting on a MultiplicativeArgument, the rightmost argument with power 1 or greater is the implicit variable, and the remainder become the column variables \(C_j\).

[4]:

# divergence of a (assuming the evolved variable - not yet set - and a live on opposite grids)
term1 = div(a).rename("div_term") #inline renaming is the recommended approach

# Diagonal stencil (performs interpolation if the variables it acts on are not on the same grid as the evolved variable)
term2 = diag(b*a**2/c) # All variables the stencil acts on must live on the same grid

[5]:

print(term1.name)
print(term1.implicitVar.name)
print(type(term1.stencil))

print(term2.name) # Default name
print(term2.implicitVar.name)
print(term2.__C__.argMultiplicity) # The column variables and their powers, note that a's power was reduced by 1 since it's used as the implicit variable

div_term
a
<class 'RMK_support.stencils.StaggeredDivStencil'>
unnamed_term
a
{'a': 1.0, 'b': 1.0, 'c': -1.0}

Evolved variables are not yet set! We can confirm this

[6]:

print(term1.evolvedVar)

None

Evolved variables can be set manually, but we will cover the preferred method in the next tutorial.

Terms can be left multiplied by any of the following:

  • Variables or MultiplicativeArguments - this will multiply row variables \(R_i\) (see the next tutorial for treating model-bound variables). If the MultiplicativeArgument has a scalar component it will multiply the existing scalar component.

  • Scalars (float or int) - this will multiply the scalar component of the term

  • Profiles - this will multiply the corresponding profile dependence by a given profile in one of the coordinates - unity by default (note that the evolved variable must live on the corresponding grids, e.g. v-grid if Profile is a V-profile, etc.)

  • TimeSignalData - this will add the corresponding periodic time-dependence (multiplicatively)

We demonstrate some of these here (for TimeSignalData see SOL-KiT-like examples)

[7]:

term1 = 2*b**2*term1 # multiplication my MultiplicativeArgument (the scalar is folded into the argument)

xProfile = grid.profile(2*np.ones(16))
term1 = xProfile * term1

# Note that all multiplication must be explicitly left-multiplication, so sometimes awkward nested brackets are unavoidable

term2 = 3*b*(xProfile *(xProfile * term2)).rename("term2") # we can rename the term at any point in this process

We can interrogate some of the private components to see that the terms are correctly constructed

[9]:

print(term1.__profiles__["X"].data) # Spatial profile dependence
print(term1.__R__.argMultiplicity) # row variables and their powers
print(term1.multConst) # scalar multiplicative constant

print("\n")
print(term2.__profiles__["X"].data) # Spatial profile dependence (note that the profile was squared)
print(term2.__R__.argMultiplicity) # row variables and their powers
print(term2.multConst) # scalar multiplicative constant

[2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2.]
{'b': 2.0}
2


[4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4.]
{'b': 1.0}
3

Terms have their own LaTeX representation (more on this in a later tutorial)

[11]:

print(term1.latex())

# We can change the constant component's latex

term1.constLatex = "c"

print(term1.latex())

 2 X  \text{b}^{2}  \nabla\cdot\left( \text{a}\right)
 c X  \text{b}^{2}  \nabla\cdot\left( \text{a}\right)

Finally, all of the special MatrixTerm properties that can be set at construction can also be set either inline or through setters. Examples include letting ReMKiT1D know that the exact sparsity pattern of a term already exists (skipPattern - saving on startup time) or that a matrix will never be updated (fixedMatrix - saving on solver time). See MatrixTerm documentation for more details.

This Page