Pencil Mathematica Tutorials

Here you can find some tutorials on using Mathematica for post-processing.

Loading the package

We need to modify init.m so that the path to the package is automatically added to the $Path variable in Mathematica. First, type

FileNameJoin[{$UserBaseDirectory, "Kernel", "init.m"}]

in Mathematica to locate this init.m file. Then, add lines

AppendTo[$Path, "your/pencil/home/mathematica"]
AppendTo[$Path, "your/pencil/home/mathematica/special"]

in this file and save it. In general the path will be $PENCIL_HOME/mathematica/, but of course you may put it somewhere else. Mathematica will not search in subdirectories, so make sure the package in right in the folder.

After updating init.m, restart the Mathematica kernel (Evaluation -> Quit Kernel). To use the package, call Needs["pc`"] and then pcInitialize[] in a notebook or a script.

To use the package on subkernels, call pcParallelize[n]. This will launch n subkernels and load the package on each of them. Then you can do things like ParallelTable[readTS[...],...]. Only loading the package on the master kernel is not enough. See the discussions here, and the ‘Possible issues’ section here.

Each time you have updated the data, remember to do pcInitialize[] and pcParallelize[n] again. These two functions remove some persistent variables defined.

Pencil Code Commands in General

For a list of all Pencil Code commands, load the package and type pcFunction[]. To access the help of any command just type ‘?’ followed by the command, e.g. ?readTS. You can also check the full definition of the command by typing ‘??’ followed by the command.

Reading and Plotting Time Series

To read the time series, type

data = readTS[sim,var1,var2,...]

where var1, var2 etc. are entries in the time series. The return of the right side is a List object, with its elements corresponding to var1, var2 etc. You can then access, for example, the time series of var2 through data[[2]] (indexing in Mathematica starts from 1).

Alternatively, you may also put a List object on the left side so that var1, var2 etc. will be assigned to each of its elements. For example,

{t,urms} = readTS[sim,"t","urms"]

Make sure that Length of the left side is equal to the number of var; otherwise Mathematica will complain.

To plot the data, you can say

fig = ListPlot[Transpose[{t,urms}],Joined->True]

or, in a one-line command,

(** same as ListPlot[Transpose[readTS[sim,"t","urms"]]] **)
fig = readTS[sim,"t","urms"]//Transpose//ListPlot

A few options for some internal plotting functions have been reset by the package. For details check ??pcLabelStyle and ??pcPlotStyle.

To export the figure in .eps format,

Export["directory/to/export/figure.eps",fig]

Reading VAR files

VAR files can be read using

data = readVARN[sim,iVAR]

Here iVAR is the index of the VAR file and starts from 0.

By default, ghost zones will not be trimmed. You can do it using the option "ltrim"->True, or "ltrim"->More; the latter will trim 2*nghost cells on each boundary.

To compute “magic” variables, you can use

data = readVARN[sim,iVAR,{"oo","bb","jj"}]

Here “oo”, “bb”, “jj” refer to vorticity, magnetic, and current fields, respectively.

The return of readVARN is an Association object (i.e., Head[data]=Association). You can obtain all of its keys by Keys[data]. Here is an example of its return:

data = readVARN[sim,iVAR,{"oo","bb","jj"}];
Keys[data]
(* {"t", "dx", "dy", "dz", "deltay", "lx", "ly", "lz", "x", "y", "z",
   "uu1", "uu2", "uu3", "lnrho", "ooo1", "ooo2", "ooo3", "bbb1", "bbb2",
   "bbb3", "jjj1", "jjj2", "jjj3"} *)

Magic variables are named using triple characters, to avoid shadowing the auxilliary ones written by the code (which will be “oo1” etc.).

The x coordinates of the mesh points is then data["x"], which will have length (16+6)^3 if the resolutoin is 16^3 and nghost=3. One can form a three-dimensional map of uu1 using

uu1 = Transpose[ data/@{"x","y","z","uu1"} ];
(* {{x1,y1,z1,f1},{x2,y2,z2,f2},...} *)

Sometimes the following method is also useful:

Clear[uu1]
grid = Transpose[ data/@{"x","y","z"} ];
uu1 = Association[ Thread[ grid->data["uu1"] ] ];

Then uu1 becomes a “function” and its value at {x1,y1,z1} is simply uu1[{x1,y1,z1}].

Visualizing slices from VAR files

A quick way to make a density plot from data is

showSlice[data, "uu1", {"z", 0.2}]

Here {"z",0.2} instructs to plot the xy slice closest to z=0.2.

For vector fields one can also use

showSliceVector[data, "uu", {"z", 0.2}]

Notice the second argument is just "uu" with no index. The function then makes a density plot of the out-of-plane component of (here "uu3"), and a superposed vector plot of the in-plane components (here "uu1" and "uu2").

Reading video files

To read video or slice files, one uses

{slices,times,position}=readSlice[sim,"uu1","xy2"]

The returned slices variable is a List of all slices at different times, and can be visualized by, say, DensityPlot[ slices[[1]] ]. position tells you the spatial coordinate of the slices.

Here is an example to make a video:

Clear[makeFrame]
makeFrame[ slice_,time_ ] := DensityPlot[ slice, PlotLabel->"t="<>ToString@time]
frames = MapThread[ makeFrame, {slices,times} ];
(* to view the video in the notebook; can be slow if too many frames*)
ListAnimate[ frame, AnimationRunning->False ]
(* output to a movie file *)
Export[ "your/output/directory/video.mov", frames, FrameRate->24 ]

One can also visualize variables in a 3D box. For more information see the comments of makeBox and makeBoxes.

Running on supercomputers

First, make sure Mathematica is available on the machine. You can check this by saying which wolfram in the terminal. If it is not installed, contact your administrator to see if it can be loaded.

Once you have loaded the Mathematica module, try wolfram in the terminal. It should bring you to the text-based interface of Mathematica. You can then follow the steps in the previous sections to set up the package.

There is a sample script in the directory $PENCIL_HOME/mathematica/sample_script.wls. Modify its first line according to where your wolfram is. Remember to include the -script option.

To run a script, use wolframscript your_script.wls.

Revision history

  • introduced the folder mathematica/special, to place packages that are specific to projects by Hongzhe Zhou at 2022-06-16 15:55:55, 0ff0c1b

  • update mathematica tutorial by Hongzhe Zhou at 2022-05-03 06:32:56, 896ed01

  • tutorial on how to run scripts by Hongzhe Zhou at 2022-04-25 13:26:45, f7f9655

  • update the mathematica package for pcParallelize[] by Hongzhe Zhou at 2022-04-12 09:50:59, 7cec76d

  • update mathematica doc by Hongzhe Zhou at 2022-04-07 12:45:29, 5dc4c00

  • update mathematica tutorial, and minor change to the package by Hongzhe Zhou at 2021-12-13 10:36:35, 648d6f5

  • updates for the mathematica pacakge and tutorial by Hongzhe Zhou at 2021-10-14 12:13:22, 23d67e1

  • mathematica package changed default trimming to False in readVARN; improved some help messages by Hongzhe Zhou at 2021-08-17 08:49:06, 64ccc36

  • added the tutorial page for mathematica package by Hongzhe Zhou at 2021-08-16 14:57:31, 33d91b3