Parameters#
The computational and physical parameters are handled by xcompact3d_toolbox.Parameters. It is built on top of Traitlets, which aims to make the parameters compatible with what XCompact3d expects, and also brings some advantages:
Attributes are type-checked;
Default values, restrictions and connections between related parameters are applied when necessary;
‘On change’ callbacks for validation and observation;
Two-way linking with ipywidgets.
import numpy as np
import xcompact3d_toolbox as x3d
The first step is to establish numerical precision. Use np.float64 if Xcompact3d was compiled with the flag -DDOUBLE_PREC (check the Makefile), use np.float32 otherwise:
x3d.param["mytype"] = np.float32
Initialization#
There are a few ways to initialize the class. First, calling it with no arguments initializes all variables with default value:
prm = x3d.Parameters()
You can access a list with all the available variables at the Api reference.
Let’s see how it looks like:
print(prm)
! -*- mode: f90 -*-
!===================
&BasicParam
!===================
C_filter = 0.49 !
beta = 1.0 ! Refinement parameter
dt = 0.001 ! Time step
gravx = 0.0 ! Gravity unitary vector in x-direction
gravy = 0.0 ! Gravity unitary vector in y-direction
gravz = 0.0 ! Gravity unitary vector in z-direction
ifilter = 0 !
ifirst = 0 ! The number for the first iteration
iibm = 0 ! Flag for immersed boundary method (0: No, 1: Yes)
iin = 0 ! Defines perturbation at initial condition
ilast = 0 ! The number for the last iteration
ilesmod = 0 ! Enables Large-Eddy methodologies (0: No, 1: Yes)
ilmn = .false. !
inflow_noise = 0.0 ! Turbulence intensity (1=100%) !! Inflow condition
init_noise = 0.0 ! Turbulence intensity (1=100%) !! Initial condition
ipost = 0 ! Enables online postprocessing at a frequency iprocessing (0: No, 1: Yes)
iscalar = 0 !
istret = 0 ! y mesh refinement (0:no, 1:center, 2:both sides, 3:bottom)
iturbine = 0 !
itype = 12 ! Flow configuration (1:Lock-exchange, 2:TGV, 3:Channel, and others)
ivisu = 1 ! Enable store snapshots at a frequency ioutput (0: No, 1: Yes)
nclx1 = 2 ! Velocity boundary condition where x=0
nclxn = 2 ! Velocity boundary condition where x=xlx
ncly1 = 2 ! Velocity boundary condition where y=0
nclyn = 2 ! Velocity boundary condition where y=yly
nclz1 = 2 ! Velocity boundary condition where z=0
nclzn = 2 ! Velocity boundary condition where z=zlz
numscalar = 0 ! Number of scalar fractions
nx = 17 ! X-direction nodes
ny = 17 ! Y-direction nodes
nz = 17 ! Z-direction nodes
p_col = 0 ! Column partition for domain decomposition and parallel computation
p_row = 0 ! Row partition for domain decomposition and parallel computation
re = 1000.0 ! Reynolds number
u1 = 2.0 !
u2 = 1.0 !
xlx = 1.0 ! Size of the box in x-direction
yly = 1.0 ! Size of the box in y-direction
zlz = 1.0 ! Size of the box in z-direction
/End
!===================
&NumOptions
!===================
cnu = 0.44 ! Ratio between hyperviscosity at km=2/3π and kc=π (dissipation factor range)
ifirstder = 4 !
iimplicit = 0 !
isecondder = 4 ! Scheme for second order derivative
itimescheme = 3 ! Time integration scheme (1: Euler, 2: AB2, 3: AB3, 5: RK3)
nu0nu = 4.0 ! Ratio between hyperviscosity/viscosity at nu (dissipation factor intensity)
/End
!===================
&InOutParam
!===================
icheckpoint = 1000 ! Frequency for writing backup file
ioutflow = 0 !
ioutput = 1000 ! Frequency for visualization file
iprocessing = 1000 ! Frequency for online postprocessing
irestart = 0 ! Read initial flow field (0: No, 1: Yes)
ninflows = 1 !
nprobes = 0 !
ntimesteps = 1 !
nvisu = 1 ! Size for visualization collection
output2D = 0 !
/End
!===================
&Statistics
!===================
/End
!===================
&CASE
!===================
/End
It is possible to access and/or set values afterwards:
# Reynolds Number
print(prm.re)
# attribute new value
prm.re = 1e6
print(prm.re)
1000.0
1000000.0
Second, we can specify some values, and let the missing ones be initialized with default value:
prm = x3d.Parameters(
filename="example.i3d",
itype=10,
nx=129,
ny=65,
nz=32,
xlx=15.0,
yly=10.0,
zlz=3.0,
nclx1=2,
nclxn=2,
ncly1=1,
nclyn=1,
nclz1=0,
nclzn=0,
iin=1,
istret=2,
re=300.0,
init_noise=0.0125,
inflow_noise=0.0125,
dt=0.0025,
ifirst=1,
ilast=45000,
irestart=0,
icheckpoint=45000,
ioutput=200,
iprocessing=50,
)
It is easy to write example.i3d to disc, just type:
prm.write()
And finally, it is possible to read the parameters from the disc:
prm = x3d.Parameters(filename="example.i3d")
prm.load()
The same result is obtained in a more concise way:
prm = x3d.Parameters(loadfile="example.i3d")
The class can also read the previous parameters format (se more information here):
prm = x3d.Parameters(loadfile="incompact3d.prm")
There are extra objects to read and write the raw binary files from XCompact3d on-demand.
Read a binary field from the disc:
ux = prm.dataset.load_array("ux-0000.bin")
Read the entire time series for a given variable:
ux = prm.dataset.load_time_series("ux")
Read all variables for a given snapshot:
snapshot = prm.dataset.load_snapshot(10)
Write
xdmffiles, so the binary files can be open in any external visualization tool:prm.dataset.write_xdmf()
Compute the coordinates, including support for mesh refinement in y:
prm.get_mesh()
{'x': array([ 0. , 0.1171875, 0.234375 , 0.3515625, 0.46875 ,
0.5859375, 0.703125 , 0.8203125, 0.9375 , 1.0546875,
1.171875 , 1.2890625, 1.40625 , 1.5234375, 1.640625 ,
1.7578125, 1.875 , 1.9921875, 2.109375 , 2.2265625,
2.34375 , 2.4609375, 2.578125 , 2.6953125, 2.8125 ,
2.9296875, 3.046875 , 3.1640625, 3.28125 , 3.3984375,
3.515625 , 3.6328125, 3.75 , 3.8671875, 3.984375 ,
4.1015625, 4.21875 , 4.3359375, 4.453125 , 4.5703125,
4.6875 , 4.8046875, 4.921875 , 5.0390625, 5.15625 ,
5.2734375, 5.390625 , 5.5078125, 5.625 , 5.7421875,
5.859375 , 5.9765625, 6.09375 , 6.2109375, 6.328125 ,
6.4453125, 6.5625 , 6.6796875, 6.796875 , 6.9140625,
7.03125 , 7.1484375, 7.265625 , 7.3828125, 7.5 ,
7.6171875, 7.734375 , 7.8515625, 7.96875 , 8.0859375,
8.203125 , 8.3203125, 8.4375 , 8.5546875, 8.671875 ,
8.7890625, 8.90625 , 9.0234375, 9.140625 , 9.2578125,
9.375 , 9.4921875, 9.609375 , 9.7265625, 9.84375 ,
9.9609375, 10.078125 , 10.1953125, 10.3125 , 10.4296875,
10.546875 , 10.6640625, 10.78125 , 10.8984375, 11.015625 ,
11.1328125, 11.25 , 11.3671875, 11.484375 , 11.6015625,
11.71875 , 11.8359375, 11.953125 , 12.0703125, 12.1875 ,
12.3046875, 12.421875 , 12.5390625, 12.65625 , 12.7734375,
12.890625 , 13.0078125, 13.125 , 13.2421875, 13.359375 ,
13.4765625, 13.59375 , 13.7109375, 13.828125 , 13.9453125,
14.0625 , 14.1796875, 14.296875 , 14.4140625, 14.53125 ,
14.6484375, 14.765625 , 14.8828125, 15. ], dtype=float32),
'y': array([ 0. , 0.04504655, 0.09029294, 0.13594207, 0.18220428,
0.22930034, 0.27746603, 0.32695502, 0.3780484 , 0.4310567 ,
0.48632994, 0.5442658 , 0.60532016, 0.6700244 , 0.7389975 ,
0.8129757 , 0.89283013, 0.9796148 , 1.0746001 , 1.1793398 ,
1.2957419 , 1.426162 , 1.5735211 , 1.7414374 , 1.9343703 ,
2.157731 , 2.417885 , 2.7218597 , 3.076467 , 3.4864438 ,
3.9514093 , 4.462363 , 5. , 5.537637 , 6.0485907 ,
6.5135565 , 6.923533 , 7.27814 , 7.5821147 , 7.842269 ,
8.06563 , 8.258562 , 8.426478 , 8.573838 , 8.704258 ,
8.820661 , 8.9254 , 9.020385 , 9.10717 , 9.187024 ,
9.261003 , 9.329975 , 9.39468 , 9.455734 , 9.51367 ,
9.568943 , 9.621952 , 9.673045 , 9.722534 , 9.7706995 ,
9.817796 , 9.864058 , 9.909707 , 9.954953 , 10. ],
dtype=float32),
'z': array([0. , 0.09375, 0.1875 , 0.28125, 0.375 , 0.46875, 0.5625 ,
0.65625, 0.75 , 0.84375, 0.9375 , 1.03125, 1.125 , 1.21875,
1.3125 , 1.40625, 1.5 , 1.59375, 1.6875 , 1.78125, 1.875 ,
1.96875, 2.0625 , 2.15625, 2.25 , 2.34375, 2.4375 , 2.53125,
2.625 , 2.71875, 2.8125 , 2.90625], dtype=float32)}
More details about I/O and array manipulations with xarray are available at:
Traitlets#
Type-checking#
All parameters are type-checked, to make sure that they are what XCompact3d expects. Use the cellcode below to see how a TraitError pops out when we try:
prm.itype = 10.5
prm.itype = -5
prm.itype = 20
prm.itype = 'sandbox'
Validation#
Some parameters, like mesh points (nx, ny and nz), trigger a validation operation when a new value is attributed to them.
Due to restrictions at the FFT library, they must be equal to:
where \(a\), \(b\) and \(c\) are non negative integers. In addition, the derivatives stencil imposes that:
Again, give it a try at the cellcode below:
prm.nx = 129
prm.nx = 4
prm.nx = 60
prm.nx = 61
Observation#
Other parameters, like mesh resolution (dx, dy and dz), are automatically updated when any new attribution occurs to mesh points and/or domain size. Let’s create a quick print functions to play with:
def show_param():
for var in "nclx1 nclxn nx xlx dx".split():
print(f"{var:>5} = {getattr(prm, var)}")
We are starting with:
show_param()
nclx1 = 2
nclxn = 2
nx = 129
xlx = 15.0
dx = 0.1171875
Let’s change just the domain’s length:
prm.xlx = 50.0
show_param()
nclx1 = 2
nclxn = 2
nx = 129
xlx = 50.0
dx = 0.390625
The resolution was updated as well. Now the number of mesh points:
prm.nx = 121
show_param()
nclx1 = 2
nclxn = 2
nx = 121
xlx = 50.0
dx = 0.4166666666666667
Again, the resolution was updated. Now we set a new mesh resolution, this time, xlx will be updated in order to satisfy the new resolution:
prm.dx = 1e-2
show_param()
nclx1 = 2
nclxn = 2
nx = 121
xlx = 1.2
dx = 0.01
Boundary conditions are observed as well. Xcompact3d allows three different BC for velocity:
Periodic
0;Free-slip
1;Dirichlet
2.
They can be assigned individually for each of the six boundaries:
nclx1andnclxn, where \(x=0\) and \(x=xlx\);ncly1andnclyn, where \(y=0\) and \(y=yly\);nclz1andnclzn, where \(z=0\) and \(z=zlz\).
It leads to 5 possibilities (00, 11, 12, 21 and 22), because both boundary must be periodic, or not, so 0 cannot be combined.
Let’s check it out, we are starting with:
show_param()
nclx1 = 2
nclxn = 2
nx = 121
xlx = 1.2
dx = 0.01
We will change just one side to periodic (nclx1 = 0), for consistence, the other side should be periodic too. Let’s see:
prm.nclx1 = 0
show_param()
nclx1 = 0
nclxn = 0
nx = 120
xlx = 1.2
dx = 0.01
Now free-slip in one side (nclx1 = 1), and the other should be non-periodic:
prm.nclx1 = 1
show_param()
nclx1 = 1
nclxn = 1
nx = 121
xlx = 1.2
dx = 0.01
Setting the other boundary to periodic:
prm.nclxn = 0
show_param()
nclx1 = 0
nclxn = 0
nx = 120
xlx = 1.2
dx = 0.01
and now back to Dirichlet:
prm.nclxn = 2
show_param()
nclx1 = 2
nclxn = 2
nx = 121
xlx = 1.2
dx = 0.01
This time, free-slip:
prm.nclxn = 1
show_param()
nclx1 = 2
nclxn = 1
nx = 121
xlx = 1.2
dx = 0.01
There was no need to update nclx1, because 1 and 2 can be freely combined. Notice that nx was modified properly from 121 to 120 and then back, according to the possible values, dx and xlx stayed untouched.
Metadata#
Traitlets types constructors have a tag method to store metadata in a dictionary. In the case of Xcompact3d-toolbox, two are especially useful:
groupdefines to what namespace a given parameter belongs when the class is written to.i3dfile (.write()method) or read from.i3dor.prmfiles (.load()method), parameters without a group are ignored for both methods;desccontains a brief description of each parameter that is shown on screen as we saw above, and also printed with the.write()method.
Declaring new parameters#
You probably would like to add more parameters for your own flow configuration, or because some of them were not implemented yet (it is a work in progress).
To do so, any Auxiliary variable can be included after initialization, like:
prm.my_variable = 0 # or any other datatype
It was called Auxiliary variable because, in this way, it will be available only for the Python application.
In order to include it at the .i3d file and make it available for XCompact3d, we can create a subclass that inherits all the functionality from xcompact3d_toolbox.Parameters:
import traitlets
# Create a class named MyParameters, which inherits the properties all properties and methods
class MyParameters(x3d.Parameters):
# .tag with group and description guarantees that the new variable will
# be compatible with all functionalities (like .write() and .load())
my_variable = traitlets.Int(default_value=0, min=0).tag(
group="BasicParam", desc="An example at the Tutorial <------"
)
# And a custom method, for instance
def my_method(self):
return self.my_variable * 2
prm = MyParameters(nx=257, ny=129, nz=31, my_variable=10) # and here we go
# Testing the method
print(prm.my_method())
# Show all parameters on screen
print(prm)
20
! -*- mode: f90 -*-
!===================
&BasicParam
!===================
C_filter = 0.49 !
beta = 1.0 ! Refinement parameter
dt = 0.001 ! Time step
gravx = 0.0 ! Gravity unitary vector in x-direction
gravy = 0.0 ! Gravity unitary vector in y-direction
gravz = 0.0 ! Gravity unitary vector in z-direction
ifilter = 0 !
ifirst = 0 ! The number for the first iteration
iibm = 0 ! Flag for immersed boundary method (0: No, 1: Yes)
iin = 0 ! Defines perturbation at initial condition
ilast = 0 ! The number for the last iteration
ilesmod = 0 ! Enables Large-Eddy methodologies (0: No, 1: Yes)
ilmn = .false. !
inflow_noise = 0.0 ! Turbulence intensity (1=100%) !! Inflow condition
init_noise = 0.0 ! Turbulence intensity (1=100%) !! Initial condition
ipost = 0 ! Enables online postprocessing at a frequency iprocessing (0: No, 1: Yes)
iscalar = 0 !
istret = 0 ! y mesh refinement (0:no, 1:center, 2:both sides, 3:bottom)
iturbine = 0 !
itype = 12 ! Flow configuration (1:Lock-exchange, 2:TGV, 3:Channel, and others)
ivisu = 1 ! Enable store snapshots at a frequency ioutput (0: No, 1: Yes)
my_variable = 10 ! An example at the Tutorial <------
nclx1 = 2 ! Velocity boundary condition where x=0
nclxn = 2 ! Velocity boundary condition where x=xlx
ncly1 = 2 ! Velocity boundary condition where y=0
nclyn = 2 ! Velocity boundary condition where y=yly
nclz1 = 2 ! Velocity boundary condition where z=0
nclzn = 2 ! Velocity boundary condition where z=zlz
numscalar = 0 ! Number of scalar fractions
nx = 257 ! X-direction nodes
ny = 129 ! Y-direction nodes
nz = 31 ! Z-direction nodes
p_col = 0 ! Column partition for domain decomposition and parallel computation
p_row = 0 ! Row partition for domain decomposition and parallel computation
re = 1000.0 ! Reynolds number
u1 = 2.0 !
u2 = 1.0 !
xlx = 1.0 ! Size of the box in x-direction
yly = 1.0 ! Size of the box in y-direction
zlz = 1.0 ! Size of the box in z-direction
/End
!===================
&NumOptions
!===================
cnu = 0.44 ! Ratio between hyperviscosity at km=2/3π and kc=π (dissipation factor range)
ifirstder = 4 !
iimplicit = 0 !
isecondder = 4 ! Scheme for second order derivative
itimescheme = 3 ! Time integration scheme (1: Euler, 2: AB2, 3: AB3, 5: RK3)
nu0nu = 4.0 ! Ratio between hyperviscosity/viscosity at nu (dissipation factor intensity)
/End
!===================
&InOutParam
!===================
icheckpoint = 1000 ! Frequency for writing backup file
ioutflow = 0 !
ioutput = 1000 ! Frequency for visualization file
iprocessing = 1000 ! Frequency for online postprocessing
irestart = 0 ! Read initial flow field (0: No, 1: Yes)
ninflows = 1 !
nprobes = 0 !
ntimesteps = 1 !
nvisu = 1 ! Size for visualization collection
output2D = 0 !
/End
!===================
&Statistics
!===================
/End
!===================
&CASE
!===================
/End
Take a look at the source code of parameters.py if you need more examples for different data types.
Graphical User Interface#
For an interactive experience launch this tutorial on Binder, the widgets are not so responsive when disconnected from a Python application.
To conclude this part of the tutorial, let’s see another option to handle the parameters. The class ParametersGui is a subclass of Parameters, and includes all the features described above. In addition, ParametersGui offers an user interface with IPywidgets.
It is still under development, more parameters and features are going to be included soon, as well as more widgets.
Just like before, we start with:
prm = x3d.ParametersGui()
Widgets are returned on demand when any instance of class ParametersGui is called, let’s see:
prm("nx", "xlx", "dx", "nclx1", "nclxn")
You can play around with the widgets above and see the effect of the observations made previously.
Notice that the Traitlets parameters are related to the value at their widgets in a two-way link, in this way, a print will show the actual value on the widgets:
show_param()
nclx1 = 2
nclxn = 2
nx = 17
xlx = 1.0
dx = 0.0625
Give it a try, modify the values at the widgets and print them again.
It also works on the other way, set a new value to a parameters will change its widget, try it:
# prm.nclx1 = 0
And of course, different widgets for the same parameter are always synchronized, change the widget below and see what happens with the widget above:
prm("nx")
A last example is about the domain decomposition for parallel computation, Xcompact3d uses 2DECOMP&FFT.
The available options for p_row and p_col are presented as functions of the number of computational cores ncores, notice that p_row * p_col = ncores should be respected and p_row * p_col = 0 activates the auto-tunning mode. The widgets are prepared to respect these restrictions:
prm("ncores", "p_row", "p_col")
To conclude this part of the tutorial, let’s see what happens when class ParametersGui is presented on screen, hover the mouse over some variable to see its description:
prm