01 - Grids#

ReMKiT1D is spatially 1D with a Legendre harmonic representation of the effectively 2D velocity space, meaning that the angular dependence is captured in a harmonic dimension, while the magnitude of the velocity is discretised directly.

For a detailed explanation of the spatial and velocity grids, see section 3.1.4 of the code paper.

A brief recap is given here:

The spatial dimension is discretised into effectively two grids, the regular grid (representing cell centres) and the staggered/dual grid, representing cell edges.
The velocity space is discretised into an integer-values harmonic dimension h and a 1D velocity magnitude grid
Interpolation between spatial grids is provided, with the assumption that (in general) scalar fields of fluid quantities live in cell centres, and vector fields of fluid quantities live on cell edges. The electron distribution function harmonics are assumed to live in centres if they have an even l-number, and edges if odd. The user can circumvent this behaviour, however.

The Grid object#

Grid data is encapsulated in the Grid object provided by RMK_support module.

[30]:

from RMK_support import Grid,gridFromDict

import numpy as np

The signature of the Grid initialiser is:

[31]:

Grid?

Init signature:
Grid(
    xGrid: numpy.ndarray,
    vGrid: numpy.ndarray = array([1.]),
    lMax=0,
    mMax=0,
    interpretXGridAsWidths=False,
    interpretVGridAsWidths=False,
    isPeriodic=False,
    isLengthInMeters=False,
)
Docstring:      Class containing x and v-grid data
Init docstring:
Grid constructor

Args:
    xGrid (numpy.ndarray): x coordinates of each spatial grid cell centre or their widths. If using widths set interpretXGridAsWidths to True.
    vGrid (numpy.ndarray, optional): v coordinates of each velocity grid cell centres or their widths.  If using widths set interpretVGridAsWidths to True. Defaults to a single cell (effectively no v-grid)
    lMax (int, optional): Maximum l harmonic number. Defaults to 0.
    mMax (int, optional): Maximum m harmonic number. Defaults to 0.
    interpretXGridAsWidths (bool, optional): If True interprets xGrid as cell widths. Defaults to False.
    interpretVGridAsWidths (bool, optional): If True interprets vGrid as cell widths. Defaults to False.
    isPeriodic (bool, optional): If True the x grid is set to be periodic. This means that the right boundary of the rightmost cell is the left boundary of the leftmost cell. Defaults to False.
    isLengthInMeters (bool, optional): If True will instruct ReMKiT1D to use the built-in normalization to deduce the normalized coordinates of the spatial grid. CAUTION: This can lead to issues if the default normalization is not used in the rest of the simulation. Defaults to False.
File:           ~/.local/lib/python3.8/site-packages/RMK_support/grid.py
Type:           type
Subclasses:

While the m-numbers are included here, note that the current version of ReMKiT1D doesn’t support l,m-resolved kinetic simulations, and this is merely future-proofing.

We proceed to construct an example (small) grid.

[32]:

grid = Grid(xGrid = 0.1 * np.ones(16),
            interpretXGridAsWidths = True, # Together with the above,
                                           # this results in a uniform spatial grid
            vGrid = 0.1 * np.ones(8),
            interpretVGridAsWidths = True, # similarly for the velocity magnitude grid
            lMax = 3, # This will result in 4 harmonics - l=0,1,2,3
            )

We can then get the individual grid points from the object

[33]:

print("X: ")
print(grid.xGrid)

print("V: ")
print(grid.vGrid)

print("l: ")
print(grid.lGrid)

X:
[0.05 0.15 0.25 0.35 0.45 0.55 0.65 0.75 0.85 0.95 1.05 1.15 1.25 1.35
 1.45 1.55]
V:
[0.05 0.15 0.25 0.35 0.45 0.55 0.65 0.75]
l:
[0, 1, 2, 3]

For the spatial grid, we can also get the (right) cell edges/dual grid:

[34]:

print("X_dual: ")
print(grid.xGridDual)

X_dual:
[0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.  1.1 1.2 1.3 1.4 1.5 1.6]

We can set the cell face Jacobian (surface area) values of the grid

[35]:

grid.xJacobian = np.linspace(1.0,2.0,17) # 17 points because we are setting
                                         # the left cell face area of the fist cell as well

We can get the volumes of the spatial cells (see the paper for details):

[36]:

print("Volumes - regular grid:")
print(grid.xGridCellVolumes())

print("Volumes - dual grid:")
print(grid.xGridCellVolumesDual(extendedBoundaryCells=True)) # We extend the boundary cells - see the paper

Volumes - regular grid:
[0.103125 0.109375 0.115625 0.121875 0.128125 0.134375 0.140625 0.146875
 0.153125 0.159375 0.165625 0.171875 0.178125 0.184375 0.190625 0.196875]
Volumes - dual grid:
[0.1578125 0.1125    0.11875   0.125     0.13125   0.1375    0.14375
 0.15      0.15625   0.1625    0.16875   0.175     0.18125   0.1875
 0.2921875 1.       ]

We can see that the volumes increase due to the widening cross-section/increasing Jacobian.

Note that the final volume value on the dual grid is 1. This is simply a placeholder, since there is one fewer dual cell than regular (except when using periodic grids).

We can also ask for a multitude of other grid properties:

[37]:

print("Number of X cells:")
print(grid.numX)

print("Number of harmonics:")
print(grid.numH)

print("Number of velocity space cells:")
print(grid.numV)

print("Spatial grid widths:")
print(grid.xWidths)

print("Velocity grid widths:")
print(grid.vWidths)

Number of X cells:
16
Number of harmonics:
4
Number of velocity space cells:
8
Spatial grid widths:
[0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1]
Velocity grid widths:
[0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1]

Grids also provide a number of data manipulation methods, such as interpolation, velocity space moments, as well as spatial integrals. Ingoing data has to conform to expected data shapes.

[38]:

# A velocity space vector

vec = np.linspace(0.1,1.2,8)

# 4 * pi * integral of vec * v**3 dv - first moment of vec

print(grid.velocityMoment(vec,1)) # velocityMoment also supports arrays in x,h,v an x,v dimensions

1.3047034290358406

[39]:

# A spatial quantity

vec = np.linspace(0,15,16)

# Treated as a vector on cell centres and interpolated on cell edges

print(grid.gridToDual(vec))

# Treated as a vector on cell edges and interpolated to cell centres

print(grid.dualToGrid(vec))

[ 0.5  1.5  2.5  3.5  4.5  5.5  6.5  7.5  8.5  9.5 10.5 11.5 12.5 13.5
 14.5 15.5]
[-0.5  0.5  1.5  2.5  3.5  4.5  5.5  6.5  7.5  8.5  9.5 10.5 11.5 12.5
 13.5 14.5]

Let’s construct a full distribution function data on our grid. It should have a shape (16,4,8) - x,h,v

[40]:

f = np.ones((16,4,8))

for i in range(4):
    for j in range(8):
        f[:,i,j] = (i+1)*np.linspace(0,15,16)

# assuming this distribution has even l harmonics on cell centres and od ones on cell edges we can interpolate it

print("Original distribution at a velocity grid point:")
print(f[:,:,0].transpose())

Original distribution at a velocity grid point:
[[ 0.  1.  2.  3.  4.  5.  6.  7.  8.  9. 10. 11. 12. 13. 14. 15.]
 [ 0.  2.  4.  6.  8. 10. 12. 14. 16. 18. 20. 22. 24. 26. 28. 30.]
 [ 0.  3.  6.  9. 12. 15. 18. 21. 24. 27. 30. 33. 36. 39. 42. 45.]
 [ 0.  4.  8. 12. 16. 20. 24. 28. 32. 36. 40. 44. 48. 52. 56. 60.]]

[41]:

print("\n")
print("With all harmonics on cell edges - interpolating even l:")
print(grid.staggeredDistToDual(f)[:,:,0].transpose()) # This will interpolate only even harmonics



With all harmonics on cell edges - interpolating even l:
[[ 0.5  1.5  2.5  3.5  4.5  5.5  6.5  7.5  8.5  9.5 10.5 11.5 12.5 13.5
  14.5 15.5]
 [ 0.   2.   4.   6.   8.  10.  12.  14.  16.  18.  20.  22.  24.  26.
  28.  30. ]
 [ 1.5  4.5  7.5 10.5 13.5 16.5 19.5 22.5 25.5 28.5 31.5 34.5 37.5 40.5
  43.5 46.5]
 [ 0.   4.   8.  12.  16.  20.  24.  28.  32.  36.  40.  44.  48.  52.
  56.  60. ]]

[42]:

print("With all harmonics in cell centres - interpolating odd l:")
print(grid.staggeredDistToGrid(f)[:,:,0].transpose()) # This will interpolate only odd harmonics

With all harmonics in cell centres - interpolating odd l:
[[ 0.  1.  2.  3.  4.  5.  6.  7.  8.  9. 10. 11. 12. 13. 14. 15.]
 [-1.  1.  3.  5.  7.  9. 11. 13. 15. 17. 19. 21. 23. 25. 27. 29.]
 [ 0.  3.  6.  9. 12. 15. 18. 21. 24. 27. 30. 33. 36. 39. 42. 45.]
 [-2.  2.  6. 10. 14. 18. 22. 26. 30. 34. 38. 42. 46. 50. 54. 58.]]

[43]:

print("With even harmonics on cell edges and odd in cell centres:")
print(grid.distFullInterp(f)[:,:,0].transpose()) # This will interpolate all harmonics to their opposite grids

With even harmonics on cell edges and odd in cell centres:
[[ 0.5  1.5  2.5  3.5  4.5  5.5  6.5  7.5  8.5  9.5 10.5 11.5 12.5 13.5
  14.5 15.5]
 [-1.   1.   3.   5.   7.   9.  11.  13.  15.  17.  19.  21.  23.  25.
  27.  29. ]
 [ 1.5  4.5  7.5 10.5 13.5 16.5 19.5 22.5 25.5 28.5 31.5 34.5 37.5 40.5
  43.5 46.5]
 [-2.   2.   6.  10.  14.  18.  22.  26.  30.  34.  38.  42.  46.  50.
  54.  58. ]]

Finally, it is possible to serialise and de-serialise grid data

[44]:

gridSerial = grid.dict()

print("Old grid:")
print(gridSerial)

newGrid = gridFromDict(gridSerial)

print("New grid:")
print(newGrid.dict())

Old grid:
{'xGrid': {'isPeriodic': False, 'isLengthInMeters': False, 'cellCentreCoords': [0.05, 0.15000000000000002, 0.25, 0.35, 0.44999999999999996, 0.5499999999999999, 0.6499999999999999, 0.7499999999999999, 0.8499999999999999, 0.9499999999999998, 1.0499999999999998, 1.15, 1.25, 1.35, 1.4500000000000002, 1.5500000000000003], 'faceJacobians': [1.0, 1.0625, 1.125, 1.1875, 1.25, 1.3125, 1.375, 1.4375, 1.5, 1.5625, 1.625, 1.6875, 1.75, 1.8125, 1.875, 1.9375, 2.0]}, 'vGrid': {'cellCentreCoords': [0.05, 0.15000000000000002, 0.25, 0.35, 0.44999999999999996, 0.5499999999999999, 0.6499999999999999, 0.7499999999999999], 'maxL': 3, 'maxM': 0}}
New grid:
{'xGrid': {'isPeriodic': False, 'isLengthInMeters': False, 'cellCentreCoords': [0.05, 0.15000000000000002, 0.25, 0.35, 0.44999999999999996, 0.5499999999999999, 0.6499999999999999, 0.7499999999999999, 0.8499999999999999, 0.9499999999999998, 1.0499999999999998, 1.15, 1.25, 1.35, 1.4500000000000002, 1.5500000000000003], 'faceJacobians': [1.0, 1.0625, 1.125, 1.1875, 1.25, 1.3125, 1.375, 1.4375, 1.5, 1.5625, 1.625, 1.6875, 1.75, 1.8125, 1.875, 1.9375, 2.0]}, 'vGrid': {'cellCentreCoords': [0.05, 0.15000000000000002, 0.25, 0.35, 0.44999999999999996, 0.5499999999999999, 0.6499999999999999, 0.7499999999999999], 'maxL': 3, 'maxM': 0}}

Profile objects#

Profile objects are wrappers for 1D data corresponding to one of the three ReMKiT1D grid dimensions X,H, or V.

These can be constructed by directly invoking the Profile constructor, but this is not recommended. Instead, a grid can be used to wrap a numpy array as a Profile while performing bounds checking.

[45]:

xProfile = grid.profile(np.ones(16),dim="X")
vProfile = grid.profile(np.ones(8),dim="V")
hProfile = grid.profile(np.ones(4),dim="H")

print("X profile:")
print(xProfile.data)
print(xProfile.dim)

print("\n")
print("V profile:")
print(vProfile.data)
print(vProfile.dim)

print("\n")
print("H profile:")
print(hProfile.data)
print(hProfile.dim)

X profile:
[1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1.]
X


V profile:
[1. 1. 1. 1. 1. 1. 1. 1.]
V


H profile:
[1. 1. 1. 1.]
H

Profile wrappers are the preferred way of specifying fixed coordinate dependence, and are used in the construction of stencils and terms (see in later tutorials).

Table of Contents

01 - Grids#

The Grid object#

Profile objects#

This Page