Using SymPy, NRPy+ is able to provide a lot of different functions out of the box, covering many different use cases. However, there are some more obscure functions that we can't access; we may also want to use a different method than the SymPy function uses or use data from multiple different points.
The code below first tells SymPy that nrpyMyFunc
should be treated as a SymPy function. This name is what should be used in the python code. Then, a dictionary is imported from outputC
so we can add an entry to it. The key for the entry should be identical to the function; the entry (myFunc_Ccode
in the example below) should be identical to the name of the function in the C code.
import sympy as sp
nrpyMyFunc = sp.Function('nrpyMyFunc')
from outputC import custom_functions_for_SymPy_ccode
custom_functions_for_SymPy_ccode["nrpyMyFunc"] = "myFunc_Ccode"
The above method is not restricted to functions; macros can be used as well. Additionally, SymPy interprets the argument of the function quite generously; the argument of the SymPy function can be more SymPy code, a string, or even nothing (possibly more!). Consider the following examples; the Python code will be on the left and the resulting C code on the right:
x = nrpyMyFunc(y)
$\rightarrow$ x = myFunc_Ccode(y)
x = nrpyMyFunc()
$\rightarrow$ x = myFunc_Ccode()
x = nrpyMyFunc
$\rightarrow$ x = myFunc_Ccode
x = nrpyMyFunc("SOME_GF")
$\rightarrow$ x = myFunc_Ccode(SOME_GF)
As you can see, we have considerable flexibility in how the C code turns out.
Most functions that we may need to use (such as trigonometric or logarithmic functions) are already included in SymPy. However, there are myriad other types of functions that can show up in various problems, such as the dilogarithm function, which is needed for the Split Monopole initial data in GiRaFFE_NRPy
. In this case, we were able to use a preexisting library, but this is not necessary in general; it is just as possible to add a function that you have written yourself.
Sometimes, SymPy may have the function you need, but implemented in a way that is not ideal. This was the case when we wrote Tutorial-Min_Max_and_Piecewise_Expressions. The SymPy implementation sp.Abs
assumed a complex input (even when we specified that the input was real), resulting in a function that was much more computationally expensive than it needed to be and introducing errors into the produced C code. We solved this by adding nrpyAbs
to the dictionary, pointing to the basic C implementation of the function.
Operations like interpolation are critical to some applications (GiRaFFE_NRPy
in particular depends on interpolation for its staggered grids). However, the only context in which data is read from multiple gridpoints in NRPy+ is finite-difference derivatives. For the interpolation of metric gridfunctions required in GiRaFFE_NRPy
, we handwrote an entire function to replicate a macro in the original GiRaFFE
, requiring several extra gridfunctions of storage compared to the original code. An alternative would be to port the macro to the GiRaFFE_NRPy
C code and use the above method to allow SymPy-generated kernels to use it.
A basic interpolator could be written as follows:
#define A0 1.0
#define A1 1.0
#define METRIC auxevol_gfs[IDX4S(GF_TO_INTERP, i0,i1,i2)]
#define METRICp1 auxevol_gfs[IDX4S(GF_TO_INTERP, i0+(flux_dirn==0),i1+(flux_dirn==1),i2+(flux_dirn==2))]
#define METRIC_INTERPED(GF_TO_INTERP) (A0*(METRIC) + A1*(METRICp1))
This computes the average of a values of the some gridfunction GF_TO_INTERP
at the current point and the next one in a given flux_dirn
, which is a basic interpolation of the value at the cell face.
Naturally, one must take care to use the correct gridfunction macros as defined in the C code. A function could be written to automate this because SymPy expressions can be easily cast to strings and NRPy+'s gridfunction #define
s are straightforwrard to generate. For instance, the following could be used:
import os,sys
nrpy_dir_path = os.path.join("..")
if nrpy_dir_path not in sys.path:
sys.path.append(nrpy_dir_path)
from outputC import outputC,custom_functions_for_SymPy_ccode # NRPy+: Core C code output module
import sympy as sp # SymPy: The Python computer algebra package upon which NRPy+ depends
import grid as gri # NRPy+: Functions having to do with numerical grids
import indexedexp as ixp # NRPy+: Symbolic indexed expression (e.g., tensors, vectors, etc.) support
gammaDD = ixp.register_gridfunctions_for_single_rank2("AUXEVOL","gammaDD","sym01")
alpha = gri.register_gridfunctions("AUXEVOL","alpha")
nrpyMetricInterp = sp.Function('nrpyMetricInterp')
def gf_id_of(gridfunction):
return (str(gridfunction).upper() + "GF")
custom_functions_for_SymPy_ccode["nrpyMetricInterp"] = "METRIC_INTERPED"
x = nrpyMetricInterp(gf_id_of(alpha)) * nrpyMetricInterp(gf_id_of(gammaDD[0][0]))
outputC(x,"x")
/* * Original SymPy expression: * "x = nrpyMetricInterp(ALPHAGF)*nrpyMetricInterp(GAMMADD00GF)" */ { x = METRIC_INTERPED(ALPHAGF)*METRIC_INTERPED(GAMMADD00GF); }