Here you can find some tutorials on using Mathematica for post-processing.
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.
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.
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]
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}]
.
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"
).
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
.
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
.