C++ and C for Python¶
Expectation from Python¶
Linear Wave¶
The governing equation is
We will see a propagating wave from left to right with phase velocity of unity.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 | import numpy as np
from matplotlib import pyplot as plt
import libst
grid = libst.Grid(0, 4*2*np.pi, 4*64)
cfl = 1
dx = (grid.xmax - grid.xmin) / grid.ncelm
dt = dx * cfl
svr = libst.LinearScalarSolver(grid=grid, time_increment=dt)
# Initialize
for e in svr.selms(odd_plane=False):
if e.xctr < 2*np.pi or e.xctr > 2*2*np.pi:
v = 0
dv = 0
else:
v = np.sin(e.xctr)
dv = np.cos(e.xctr)
e.set_so0(0, v)
e.set_so1(0, dv)
plt.figure(figsize=(15,10))
plt.xlim((0, 8))
plt.xlabel('$x$ $(\pi)$')
plt.ylabel('$u$')
plt.grid()
plt.plot(svr.xctr() / np.pi, svr.get_so0(0).ndarray, '-', label='begin')
svr.setup_march()
svr.march_alpha2(50)
plt.plot(svr.xctr() / np.pi, svr.get_so0(0).ndarray, '-', label='end')
plt.legend()
|
Inviscid Burgers Equation¶
The second example is a non-linear equation:
The wave propagates in a way that is not predictable.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 | import numpy as np
from matplotlib import pyplot as plt
import libst
res = 32
xcrd = np.arange(res+1) / res * 2 * np.pi
time_stop = 2*np.pi
grid = libst.Grid(xcrd)
cfl_max = 1.0
dx = (grid.xmax - grid.xmin) / grid.ncelm
dt_max = dx * cfl_max
nstep = int(np.ceil(time_stop / dt_max))
dt = time_stop / nstep
cfl = dt / dx
svr = libst.InviscidBurgersSolver(grid=grid, time_increment=dt)
# Initialize
svr.set_so0(0, np.sin(svr.xctr()))
svr.set_so1(0, np.cos(svr.xctr()))
plt.figure(figsize=(15,10))
plt.xlim((0, 2))
plt.xlabel('$x$ $(\pi)$')
plt.ylabel('$u$')
plt.grid()
plt.plot(svr.xctr() / np.pi, svr.get_so0(0).ndarray, '-', label='begin')
svr.setup_march()
svr.march_alpha2(50)
plt.plot(svr.xctr() / np.pi, svr.get_so0(0).ndarray, '-', label='end')
plt.legend()
|
Pybind11 Build System¶
Pybind11 is a header-only C++ template library, that allows calling CPython API and provides C++ friendly semantics to allow Python to call C++ constructs and vise versa.
A header-only library doesn’t have anything to be built. When we say “building” pybind11, we mean to build the project that uses pybind11.
To build pybind11, we need CPython. It optionally depends on numpy and eigen. There are several suggested ways to build. Here list those I think important:
Setuptools¶
Setuptools is an enhancement to Python built-in distutils. Because pybind11 is released to PyPI as a Python package (https://pypi.org/project/pybind11/), both setuptools and distutils can get the header files and use them to build C++ file that use pybind11.
There is an example for using setuptools to build pybind11 (https://github.com/pybind/python_example/blob/master/setup.py):
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 | from setuptools import setup, Extension
from setuptools.command.build_ext import build_ext
import sys
import setuptools
__version__ = '0.0.1'
class get_pybind_include(object):
"""Helper class to determine the pybind11 include path
The purpose of this class is to postpone importing pybind11
until it is actually installed, so that the ``get_include()``
method can be invoked. """
def __init__(self, user=False):
self.user = user
def __str__(self):
import pybind11
return pybind11.get_include(self.user)
ext_modules = [
Extension(
'python_example',
['src/main.cpp'],
include_dirs=[
# Path to pybind11 headers
get_pybind_include(),
get_pybind_include(user=True)
],
language='c++'
),
]
# As of Python 3.6, CCompiler has a `has_flag` method.
# cf http://bugs.python.org/issue26689
def has_flag(compiler, flagname):
"""Return a boolean indicating whether a flag name is supported on
the specified compiler.
"""
import tempfile
with tempfile.NamedTemporaryFile('w', suffix='.cpp') as f:
f.write('int main (int argc, char **argv) { return 0; }')
try:
compiler.compile([f.name], extra_postargs=[flagname])
except setuptools.distutils.errors.CompileError:
return False
return True
def cpp_flag(compiler):
"""Return the -std=c++[11/14/17] compiler flag.
The newer version is prefered over c++11 (when it is available).
"""
flags = ['-std=c++17', '-std=c++14', '-std=c++11']
for flag in flags:
if has_flag(compiler, flag): return flag
raise RuntimeError('Unsupported compiler -- at least C++11 support '
'is needed!')
class BuildExt(build_ext):
"""A custom build extension for adding compiler-specific options."""
c_opts = {
'msvc': ['/EHsc'],
'unix': [],
}
l_opts = {
'msvc': [],
'unix': [],
}
if sys.platform == 'darwin':
darwin_opts = ['-stdlib=libc++', '-mmacosx-version-min=10.7']
c_opts['unix'] += darwin_opts
l_opts['unix'] += darwin_opts
def build_extensions(self):
ct = self.compiler.compiler_type
opts = self.c_opts.get(ct, [])
link_opts = self.l_opts.get(ct, [])
if ct == 'unix':
opts.append('-DVERSION_INFO="%s"' % self.distribution.get_version())
opts.append(cpp_flag(self.compiler))
if has_flag(self.compiler, '-fvisibility=hidden'):
opts.append('-fvisibility=hidden')
elif ct == 'msvc':
opts.append('/DVERSION_INFO=\\"%s\\"' % self.distribution.get_version())
for ext in self.extensions:
ext.extra_compile_args = opts
ext.extra_link_args = link_opts
build_ext.build_extensions(self)
setup(
name='python_example',
version=__version__,
author='Sylvain Corlay',
author_email='sylvain.corlay@gmail.com',
url='https://github.com/pybind/python_example',
description='A test project using pybind11',
long_description='',
ext_modules=ext_modules,
install_requires=['pybind11>=2.4'],
setup_requires=['pybind11>=2.4'],
cmdclass={'build_ext': BuildExt},
zip_safe=False,
)
|
Cmake with a Sub-Directory¶
When the source tree is put in a sub-directory in your project, as mentioned in
the document,
you can use cmake add_subdirectory
to include the pybind11:
cmake_minimum_required(VERSION 2.8.12)
project(example)
add_subdirectory(pybind11)
pybind11_add_module(example example.cpp)
Pybind11 provides the cmake command pybind11_add_module
. It set various
flags to build your C++ code as an extension module.
Cmake with installed pybind11¶
If pybind11 is installed using cmake itself, the *.cmake
files that
pybind11 supplies are installed to the specified location. It’s not needed to
write add_subdirectory
in the CMakeLists.txt
in your project.
Additional Wrapping Layer for Customization¶
Wrapper needs to take care of the differences between the dynamic behaviors in Python and the staticity in C++. You can directly call pybind11 API. But a better way is to create another wrapping layer between the pybind11 and your library code. It allows us to insert additional code in a systematic way. Since it is difficult to see the point in a small example, I pull the code for a bigger project ‘turgon’ for demonstration.
Here is one way to implement the additional wrapping layer:
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 | #pragma once
/*
* Copyright (c) 2017, Yung-Yu Chen <yyc@solvcon.net>
* BSD 3-Clause License, see COPYING
*/
#include <pybind11/pybind11.h>
#include <pybind11/numpy.h>
#define NPY_NO_DEPRECATED_API NPY_1_7_API_VERSION
#include <numpy/arrayobject.h>
#include <memory>
#include <type_traits>
PYBIND11_DECLARE_HOLDER_TYPE(T, std::unique_ptr<T>);
PYBIND11_DECLARE_HOLDER_TYPE(T, std::shared_ptr<T>);
#ifdef __GNUG__
# define SPACETIME_PYTHON_WRAPPER_VISIBILITY __attribute__((visibility("hidden")))
#else
# define SPACETIME_PYTHON_WRAPPER_VISIBILITY
#endif
namespace spacetime
{
namespace python
{
/**
* Helper template for pybind11 class wrappers.
*/
template< class Wrapper, class Wrapped, class Holder = std::unique_ptr<Wrapped>, class WrappedBase = Wrapped >
class
SPACETIME_PYTHON_WRAPPER_VISIBILITY
WrapBase {
public:
using wrapper_type = Wrapper;
using wrapped_type = Wrapped;
using wrapped_base_type = WrappedBase;
using holder_type = Holder;
using base_type = WrapBase< wrapper_type, wrapped_type, holder_type, wrapped_base_type >;
using class_ = typename std::conditional<
std::is_same< Wrapped, WrappedBase >::value
, pybind11::class_< wrapped_type, holder_type >
, pybind11::class_< wrapped_type, wrapped_base_type, holder_type >
>::type;
static wrapper_type & commit(pybind11::module * mod, const char * pyname, const char * clsdoc) {
static wrapper_type derived(mod, pyname, clsdoc);
return derived;
}
WrapBase() = delete;
WrapBase(WrapBase const & ) = default;
WrapBase(WrapBase &&) = delete;
WrapBase & operator=(WrapBase const & ) = default;
WrapBase & operator=(WrapBase &&) = delete;
~WrapBase() = default;
#define DECL_ST_PYBIND_CLASS_METHOD(METHOD) \
template< class... Args > \
/* NOLINTNEXTLINE(bugprone-macro-parentheses) */ \
wrapper_type & METHOD(Args&&... args) { \
m_cls.METHOD(std::forward<Args>(args)...); \
return *static_cast<wrapper_type*>(this); \
}
DECL_ST_PYBIND_CLASS_METHOD(def)
DECL_ST_PYBIND_CLASS_METHOD(def_readwrite)
DECL_ST_PYBIND_CLASS_METHOD(def_property)
DECL_ST_PYBIND_CLASS_METHOD(def_property_readonly)
DECL_ST_PYBIND_CLASS_METHOD(def_property_readonly_static)
#undef DECL_ST_PYBIND_CLASS_METHOD
protected:
WrapBase(pybind11::module * mod, const char * pyname, const char * clsdoc)
: m_cls(*mod, pyname, clsdoc)
{}
private:
class_ m_cls;
}; /* end class WrapBase */
} /* end namespace python */
} /* end namespace spacetime */
|
The following classes WrapGrid
, WrapField
, WrapSolver
,
WrapCelm
, and WrapSelm
show how WrapBase
is used:
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 | #pragma once
/*
* Copyright (c) 2018, Yung-Yu Chen <yyc@solvcon.net>
* BSD 3-Clause License, see COPYING
*/
#include "spacetime/python/common.hpp"
namespace spacetime
{
namespace python
{
class
SPACETIME_PYTHON_WRAPPER_VISIBILITY
WrapGrid
: public WrapBase< WrapGrid, Grid, std::shared_ptr<Grid> >
{
friend base_type;
WrapGrid(pybind11::module * mod, const char * pyname, const char * clsdoc)
: base_type(mod, pyname, clsdoc)
{
namespace py = pybind11;
(*this)
.def(
py::init([](real_type xmin, real_type xmax, size_t nelm) {
return Grid::construct(xmin, xmax, nelm);
}),
py::arg("xmin"), py::arg("xmax"), py::arg("nelm")
)
.def(
py::init([](xt::pyarray<wrapped_type::value_type> & xloc) {
return Grid::construct(xloc);
}),
py::arg("xloc")
)
.def("__str__", &detail::to_str<wrapped_type>)
.def_property_readonly("xmin", &wrapped_type::xmin)
.def_property_readonly("xmax", &wrapped_type::xmax)
.def_property_readonly("ncelm", &wrapped_type::ncelm)
.def_property_readonly("nselm", &wrapped_type::nselm)
.def_property_readonly(
"xcoord",
static_cast<wrapped_type::array_type & (wrapped_type::*)()>(&wrapped_type::xcoord)
)
.def_property_readonly_static("BOUND_COUNT", [](py::object const &){ return Grid::BOUND_COUNT; })
;
}
}; /* end class WrapGrid */
class
SPACETIME_PYTHON_WRAPPER_VISIBILITY
WrapField
: public WrapBase< WrapField, Field, std::shared_ptr<Field> >
{
friend base_type;
WrapField(pybind11::module * mod, const char * pyname, const char * clsdoc)
: base_type(mod, pyname, clsdoc)
{
namespace py = pybind11;
(*this)
.def("__str__", &detail::to_str<wrapped_type>)
.def_property_readonly("grid", [](wrapped_type & self){ return self.grid().shared_from_this(); })
.def_property_readonly("nvar", &wrapped_type::nvar)
.def_property(
"time_increment",
&wrapped_type::time_increment,
&wrapped_type::set_time_increment
)
.def_property_readonly("dt", &wrapped_type::dt)
.def_property_readonly("hdt", &wrapped_type::hdt)
.def_property_readonly("qdt", &wrapped_type::qdt)
.def(
"celm",
static_cast<Celm (wrapped_type::*)(sindex_type, bool)>(&wrapped_type::celm_at<Celm>),
py::arg("ielm"), py::arg("odd_plane")=false
)
.def(
"selm",
static_cast<Selm (wrapped_type::*)(sindex_type, bool)>(&wrapped_type::selm_at<Selm>),
py::arg("ielm"), py::arg("odd_plane")=false
)
;
}
}; /* end class WrapField */
class
SPACETIME_PYTHON_WRAPPER_VISIBILITY
WrapSolver
: public WrapSolverBase< WrapSolver, Solver >
{
using base_type = WrapSolverBase< WrapSolver, Solver >;
using wrapper_type = typename base_type::wrapper_type;
using wrapped_type = typename base_type::wrapped_type;
friend base_type;
friend base_type::base_type;
WrapSolver(pybind11::module * mod, const char * pyname, const char * clsdoc)
: base_type(mod, pyname, clsdoc)
{
namespace py = pybind11;
(*this)
.def
(
py::init(static_cast<std::shared_ptr<wrapped_type> (*) (
std::shared_ptr<Grid> const &, typename wrapped_type::value_type, size_t
)>(&wrapped_type::construct))
, py::arg("grid"), py::arg("time_increment"), py::arg("nvar")
)
;
}
}; /* end class WrapSolver */
class
SPACETIME_PYTHON_WRAPPER_VISIBILITY
WrapCelm
: public WrapCelmBase< WrapCelm, Celm >
{
using base_type = WrapCelmBase< WrapCelm, Celm >;
friend base_type::base_type::base_type;
WrapCelm(pybind11::module * mod, const char * pyname, const char * clsdoc)
: base_type(mod, pyname, clsdoc)
{}
}; /* end class WrapCelm */
class
SPACETIME_PYTHON_WRAPPER_VISIBILITY
WrapSelm
: public WrapSelmBase< WrapSelm, Selm >
{
using base_type = WrapSelmBase< WrapSelm, Selm >;
friend base_type::base_type::base_type;
WrapSelm(pybind11::module * mod, const char * pyname, const char * clsdoc)
: base_type(mod, pyname, clsdoc)
{}
}; /* end class WrapSelm */
} /* end namespace python */
} /* end namespace spacetime */
// vim: set et sw=4 ts=4:
|
Note that the use of WrapperBase
saves little duplicated code, if any. But
it facilitates to keep the .cpp
file for building the Python extension to
be very short:
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 | /*
* Copyright (c) 2018, Yung-Yu Chen <yyc@solvcon.net>
* BSD 3-Clause License, see COPYING
*/
#include "spacetime/python.hpp" // must be first
#include "spacetime.hpp"
#include <algorithm>
#include <cstring>
#include <memory>
#include <utility>
#include <vector>
namespace
{
PyObject * initialize_spacetime(pybind11::module * mod)
{
namespace spy = spacetime::python;
xt::import_numpy(); // otherwise numpy c api segfault.
mod->doc() = "_libst: One-dimensional space-time CESE method code";
spy::WrapGrid::commit(mod, "Grid", "Spatial grid data");
spy::WrapField::commit(mod, "Field", "Solution data");
return mod->ptr();
}
} /* end namespace */
PYBIND11_MODULE(_libst, mod) // NOLINT
{
namespace spy = spacetime::python;
spy::ModuleInitializer::get_instance()
.add(initialize_spacetime)
.add_solver<spy::WrapSolver, spy::WrapCelm, spy::WrapSelm>
(&mod, "", "no equation")
.add_solver<spy::WrapLinearScalarSolver, spy::WrapLinearScalarCelm, spy::WrapLinearScalarSelm>
(&mod, "LinearScalar", "a linear scalar equation")
.add_solver<spy::WrapInviscidBurgersSolver, spy::WrapInviscidBurgersCelm, spy::WrapInviscidBurgersSelm>
(&mod, "InviscidBurgers", "the inviscid Burgers equation")
.initialize(&mod)
;
}
|
Wrapping API¶
Let’s take a look at the helpers that pybind11 provides.
Function and Property¶
Let’s use the Grid class as an example to demonstrate how to expose functions and properties. We have a constructor:
.def
(
py::init
(
[](real_type xmin, real_type xmax, size_t nelm)
{
return Grid::construct(xmin, xmax, nelm);
}
)
, py::arg("xmin"), py::arg("xmax"), py::arg("nelm")
)
and a special-purpose function:
.def("__str__", &detail::to_str<wrapped_type>)
>>> grid = libst.Grid(0, 8, 4*64)
>>> print('directly call Grid.__str__():', grid.__str__())
directly call Grid.__str__(): Grid(xmin=0, xmax=8, ncelm=256)
>>> print('call str(Grid):', str(grid))
call str(Grid): Grid(xmin=0, xmax=8, ncelm=256)
Sometimes properties are more suitable for certain getters:
.def_property_readonly("xmin", &wrapped_type::xmin)
.def_property_readonly("xmax", &wrapped_type::xmax)
.def_property_readonly_static("BOUND_COUNT", [](py::object const &){ return Grid::BOUND_COUNT; })
As shown above, pybind11 supports static properties that are associated on the class instead of the instance.
>>> print('grid.BOUND_COUNT =', grid.BOUND_COUNT)
grid.BOUND_COUNT = 2
>>> print('grid.xmin =', grid.xmin)
grid.xmin = 0.0
>>> print('grid.xmax =', grid.xmax)
grid.xmax = 8.0
>>> print('Grid.BOUND_COUNT (number of points beyond spatial boundary) =', libst.Grid.BOUND_COUNT)
Grid.BOUND_COUNT (number of points beyond spatial boundary) = 2
>>> print('Grid.xmin =', libst.Grid.xmin)
Grid.xmin = <property object at 0x110e9ffb0>
>>> print('Grid.xmax =', libst.Grid.xmax)
Grid.xmax = <property object at 0x110ea60b0>
Compare to pure Python object:
class PythonGrid:
BOUND_COUNT = 2
@property
def xmin(self):
return 0
@property
def xmax(self):
return 8
>>> print(PythonGrid.BOUND_COUNT)
2
>>> print(PythonGrid.xmin)
<property object at 0x1112daad0>
>>> print(PythonGrid.xmax)
<property object at 0x1112dab30>
In addition to def_property_readonly
and def_property_readonly_static
,
pybind11 also provides:
def_property
anddef_property_static
for read/write properties with C++ accessors.def_readonly
anddef_readonly_static
for read-only access to C++ data members.def_readwrite
anddef_readwrite_static
for read/write access to C++ data members.
Named and Keyword Arguments¶
Pybind11 allows named arguments. I take the advantage for wrapping the
constructor of Grid
:
.def
(
py::init
(
[](real_type xmin, real_type xmax, size_t nelm)
{
return Grid::construct(xmin, xmax, nelm);
}
)
, py::arg("xmin"), py::arg("xmax"), py::arg("nelm")
)
and a special-purpose function:
.def("__str__", &detail::to_str<wrapped_type>)
>>> grid = libst.Grid(xmin=0, xmax=8, nelm=4*64)
>>> print('directly call Grid.__str__():', grid.__str__())
directly call Grid.__str__(): Grid(xmin=0, xmax=8, ncelm=256)
>>> print('call str(Grid):', str(grid))
call str(Grid): Grid(xmin=0, xmax=8, ncelm=256)
pybind11::arg
also allows default value to the arguments (keyword
arguments). For example, see what I have for all the solver classes:
.def
(
"selms"
, [](wrapped_type & self, bool odd_plane)
{ return elm_iter_type(self.shared_from_this(), odd_plane, 0, true); }
, py::arg("odd_plane")=false
)
grid = libst.Grid(0, 4*2*np.pi, 4*64)
cfl = 1
dx = (grid.xmax - grid.xmin) / grid.ncelm
dt = dx * cfl
svr = libst.LinearScalarSolver(grid=grid, time_increment=dt)
>>> print(svr.selms())
SolverElementIterator(selm, on_even_plane, current=0, nelem=257)
>>> print(svr.selms(False))
SolverElementIterator(selm, on_even_plane, current=0, nelem=257)
>>> print(svr.selms(True))
SolverElementIterator(selm, on_odd_plane, current=0, nelem=256)
>>> print(svr.selms(odd_plane=False))
SolverElementIterator(selm, on_even_plane, current=0, nelem=257)
>>> print(svr.selms(odd_plane=True))
SolverElementIterator(selm, on_odd_plane, current=0, nelem=256)
What Happens in Python Stays in Python (or Pybind11)¶
When wrapping from C++ to Python, there are constructs only available in the scripting language but not the low-level implementation. When it happens, write the adapting code in the pybind11 layer and do not pollute the low-level implementation.
One example is Python iterator protocol. The adapting code is:
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 | template< typename ST >
class SolverElementIterator
{
public:
using solver_type = ST;
SolverElementIterator() = delete;
SolverElementIterator(std::shared_ptr<ST> sol, bool odd_plane, size_t starting, bool selm)
: m_solver(std::move(sol)), m_odd_plane(odd_plane), m_current(starting), m_selm(selm)
{}
typename ST::celm_type next_celm()
{
size_t ncelm = m_solver->grid().ncelm();
if (m_odd_plane) { --ncelm; }
if (m_current >= ncelm) { throw pybind11::stop_iteration(); }
typename ST::celm_type ret = m_solver->celm(m_current, m_odd_plane);
++m_current;
return ret;
}
typename ST::selm_type next_selm()
{
size_t nselm = m_solver->grid().nselm();
if (m_odd_plane) { --nselm; }
if (m_current >= nselm) { throw pybind11::stop_iteration(); }
typename ST::selm_type ret = m_solver->selm(m_current, m_odd_plane);
++m_current;
return ret;
}
bool is_selm() const { return m_selm; }
bool on_odd_plane() const { return m_odd_plane; }
size_t current() const { return m_current; }
size_t nelem() const
{
size_t ret = is_selm() ? m_solver->grid().nselm() : m_solver->grid().ncelm();
if (m_odd_plane) { --ret; }
return ret;
}
private:
std::shared_ptr<solver_type> m_solver;
bool m_odd_plane;
size_t m_current = 0;
bool m_selm = false;
}; /* end class SolverElementIterator */
|
The wrapping code is:
using elm_iter_type = SolverElementIterator<wrapped_type>;
std::string elm_pyname = std::string(pyname) + "ElementIterator";
pybind11::class_< elm_iter_type >(*mod, elm_pyname.c_str())
.def("__str__", &detail::to_str<elm_iter_type>)
.def("__iter__", [](elm_iter_type & self){ return self; })
.def(
"__next__"
, [](elm_iter_type & self)
{
py::object ret;
if (self.is_selm()) { ret = py::cast(self.next_selm()); }
else { ret = py::cast(self.next_celm()); }
return ret;
}
)
;
Let’s use a concrete solver of linear wave (governing equation is \(u_t + u_x = 0\)) to demonstrate how it works in Python:
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 | import numpy as np
from matplotlib import pyplot as plt
import libst
grid = libst.Grid(0, 4*2*np.pi, 4*64)
cfl = 1
dx = (grid.xmax - grid.xmin) / grid.ncelm
dt = dx * cfl
svr = libst.LinearScalarSolver(grid=grid, time_increment=dt)
# Initialize
# NOTE: 'selms' returns a template instance of SolverElementIterator
for e in svr.selms(odd_plane=False):
if e.xctr < 2*np.pi or e.xctr > 2*2*np.pi:
v = 0
dv = 0
else:
v = np.sin(e.xctr)
dv = np.cos(e.xctr)
e.set_so0(0, v)
e.set_so1(0, dv)
# Plot it
plt.figure(figsize=(15,10))
plt.xlim((0, 8))
plt.xlabel('$x$ $(\pi)$')
plt.ylabel('$u$')
plt.grid()
plt.plot(svr.xctr() / np.pi, svr.get_so0(0).ndarray, '-')
|
Manipulate Python Objects in C++¶
Pybind11 provides C++ API for manipulating Python object (the C struct
PyObject
), so that we don’t need to dig into the Python C API
and worry about the reference counting by hand.
The first example is to create a None
object from C++.
#include "pybind11/pybind11.h"
PYBIND11_MODULE(code_none, m)
{
namespace py = pybind11;
m
.def
(
"create_none", []() { return py::none(); }
)
;
}
>>> print(type(create_none()))
<class 'NoneType'>
>>> assert None is create_none()
>>> print(create_none())
None
pybind11::object
is the C++ counterpart of PyObject
in C,
but automatically does reference counting for us.
#include "pybind11/pybind11.h"
PYBIND11_MODULE(code_object, m)
{
namespace py = pybind11;
m
.def
(
"return_none"
, []()
{
py::object ret = py::none();
return ret;
}
)
;
}
>>> print(return_none, return_none())
<built-in method return_none of PyCapsule object at 0x1111b4300> None
Pybind11 allows to use pybind11::object::attr
to assign attribute to a
Python object.
#include "pybind11/pybind11.h"
PYBIND11_MODULE(code_attr, m)
{
namespace py = pybind11;
m.attr("string_name") = "string_content";
}
>>> print(type(string_name), string_name)
<class 'str'> string_content
Pybind11 provides helpers to import Python module and access attibutes of every Python object, including a Python module.
#include "pybind11/pybind11.h"
PYBIND11_MODULE(code1, m)
{
namespace py = pybind11;
m
.def
(
"return_numpy_version"
, []()
{
py::object numpy = py::module::import("numpy");
return numpy.attr("__version__");
}
)
;
m.attr("alias_to_return_numpy_version") = m.attr("return_numpy_version");
}
>>> print(return_numpy_version())
1.17.0.dev0+3c3ba10
>>> import numpy as np
>>> print(np.__version__)
1.17.0.dev0+3c3ba10
>>> assert np.__version__ is return_numpy_version()
>>> print(return_numpy_version)
<built-in method return_numpy_version of PyCapsule object at 0x1111b4060>
>>> print(alias_to_return_numpy_version)
<built-in method return_numpy_version of PyCapsule object at 0x1111b4060>
Python Containers¶
Pybind11 provides C++ API for creating and manipulating the most important
Python containers: tuple
, list
, and
dict
. See
https://pybind11.readthedocs.io/en/stable/advanced/pycpp/object.html and the
unit tests for more information.
tuple
¶
#include "pybind11/pybind11.h"
PYBIND11_MODULE(code_tuple, m)
{
namespace py = pybind11;
py::tuple my_tuple = py::make_tuple("string_data_in_tuple", 10, 3.1415926);
m.attr("my_tuple") = my_tuple;
}
>>> print(type(my_tuple), my_tuple)
<class 'tuple'> ('string_data_in_tuple', 10, 3.1415926)
list
¶
#include "pybind11/pybind11.h"
PYBIND11_MODULE(code2, m)
{
namespace py = pybind11;
py::list my_list = py::list();
my_list.append("string_data_in_list");
my_list.append(11);
my_list.append(2.71828);
py::list my_list2 = py::make_tuple("string_data_in_list2", 12);
m.attr("my_list") = my_list;
m.attr("my_list2") = my_list2;
}
>>> print(type(my_list), my_list)
<class 'list'> ['string_data_in_list', 11, 2.71828]
>>> print(type(my_list2), my_list2)
<class 'list'> ['string_data_in_list2', 12]
dict
¶
#include "pybind11/pybind11.h"
PYBIND11_MODULE(code2, m)
{
namespace py = pybind11;
py::dict my_dict;
my_dict["key_string"] = "string_data_in_dict";
my_dict["key_int"] = 13;
my_dict["key_real"] = 1.414;
m.attr("my_dict") = my_dict;
}
>>> print(type(my_dict), my_dict)
<class 'dict'> {'key_string': 'string_data_in_dict', 'key_int': 13, 'key_real': 1.414}
Use CPython API with Pybind11¶
We can use any Python C API with pybind11.
When we import pybind11/pybind11.h
, we don’t need to import Python.h
,
becuase the former does it for us. But please note that
pybind11/pybind11.h
or Python.h
should be included before every other
inclusion.
#include "pybind11/pybind11.h"
#include "Python.h" // Unnecessary
using namespace pybind11;
PYBIND11_MODULE(ex_long, m)
{
PyObject * v = PyLong_FromLong(2000000);
m.attr("integer_value") = v;
Py_DECREF(v);
}
>>> print(type(integer_value), integer_value)
<class 'int'> 2000000
PyObject
reference counting¶
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 | #include "pybind11/pybind11.h"
using namespace pybind11;
static PyObject * s;
PYBIND11_MODULE(ex_str, m)
{
s = PyUnicode_FromString("string_from_c");
m.attr("string_value") = s;
Py_DECREF(s);
m
.def
(
"show_ref_count_with_handle"
, [](handle const & h)
{
return Py_REFCNT(h.ptr());
}
)
.def
(
"show_ref_count_with_object"
, [](object const & o)
{
return Py_REFCNT(o.ptr());
}
)
.def
(
"show_string_value_ref_count"
, [&]()
{
return Py_REFCNT(s);
}
)
;
}
|
1 2 3 4 5 6 7 8 9 10 11 12 13 | def check_string_value():
print(type(string_value), string_value)
print('before aliasing')
print(show_ref_count_with_object(string_value), 'refcnt by object')
print(show_ref_count_with_handle(string_value), 'refcnt by handle')
print(sys.getrefcount(string_value), 'refcnt by sys')
print(show_string_value_ref_count(), 'refcnt from c++')
string_value_aliasing = string_value
print('after aliasing')
print(show_ref_count_with_object(string_value), 'refcnt by object')
print(show_ref_count_with_handle(string_value), 'refcnt by handle')
print(sys.getrefcount(string_value), 'refcnt by sys')
print(show_string_value_ref_count(), 'refcnt from c++')
|
>>> check_string_value()
<class 'str'> string_from_c
before aliasing
7 refcnt by object
6 refcnt by handle
5 refcnt by sys
4 refcnt from c++
after aliasing
8 refcnt by object
7 refcnt by handle
6 refcnt by sys
5 refcnt from c++
Pybind11 offers two low-level shorthands for reference counting:
handle::inc_ref()
and handle::dec_ref()
. If we don’t want to go so
low-level, it provides reinterpret_borrow
and reinterpret_steal
function templates.
Built-in Types¶
It is possible to use Python C API along with pybind11. I am going to should how to do it. Please keep in mind that the examples in the notes omit a lot of error checking code that is necessary for a system to run correctly. When you need to use the C API, consult the manual: Python/C API Reference Manual.
Cached Value¶
Python caches small (-5 to 256) integers (see the code). Don’t get surprised when you see a large reference count for some of them integers:
>>> print('ref counts of 0:', sys.getrefcount(0))
ref counts of 0: 10198
>>> print('ref counts of 257:', sys.getrefcount(257))
ref counts of 257: 3
Real number doesn’t have that cache:
>>> print(sys.getrefcount(0.0))
3
Python interns strings consisting of alphanumerical and underscore characters.
>>> print('' is '')
True
>>> print(sys.getrefcount(''))
5552
def check_string_intern():
s1 = 'numerical'
print(sys.getrefcount('numerical'))
print(s1 is 'numerical')
s2 = 'num' + 'erical'
print(s1 is s2)
print(sys.getrefcount('numerical'))
>>> check_string_intern()
4
True
True
5
Attribute Access¶
The Python object protocol defines a set of API for accessing object attributes. Here is a simple example that sets and gets an attribute of an object using the API:
int PyObject_SetAttr(PyObject *o, PyObject *attr_name, PyObject *v);
PyObject* PyObject_GetAttr(PyObject *o, PyObject *attr_name);
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 | #include "pybind11/pybind11.h"
using namespace pybind11;
void attach_attr(PyObject * o, PyObject * name, PyObject * attr)
{
/*int ret =*/
PyObject_SetAttr(o, name, attr);
}
PyObject * retrieve_attr(PyObject * o, PyObject * name)
{
PyObject * ret = PyObject_GetAttr(o, name);
return ret;
}
PYBIND11_MODULE(ex_attr, m)
{
m
.def
(
"attach_attr"
, [](object & o, object & name, object & attr)
{
attach_attr(o.ptr(), name.ptr(), attr.ptr());
}
)
.def
(
"retrieve_attr"
, [](object & o, object & name)
{
handle(retrieve_attr(o.ptr(), name.ptr()));
}
)
;
}
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 | >>> class Cls():
>>> pass
>>> obj = Cls()
>>> val = 'attached value'
>>> print(sys.getrefcount(val))
3
>>>
>>> attach_attr(obj, 'name', val)
>>> print(sys.getrefcount(val))
4
>>>
>>> print(obj.name is val)
True
>>> print(sys.getrefcount(val))
4
>>>
>>> val2 = retrieve_attr(obj, 'name')
>>> print(sys.getrefcount(val))
5
|
There are shorthand versions of the API that takes C string for the attribute name:
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 | #include "pybind11/pybind11.h"
#include <string>
using namespace pybind11;
void attach_attr(PyObject * o, char const * name, PyObject * attr)
{
/*int ret =*/
PyObject_SetAttrString(o, name, attr);
}
PyObject * retrieve_attr(PyObject * o, char const * name)
{
PyObject * ret = PyObject_GetAttrString(o, name);
return ret;
}
PYBIND11_MODULE(ex_attr_by_string, m)
{
m
.def
(
"attach_attr_by_string"
, [](object & o, object & name, object & attr)
{
std::string name_str = cast<std::string>(name);
attach_attr(o.ptr(), name_str.c_str(), attr.ptr());
}
)
.def
(
"retrieve_attr_by_string"
, [](object & o, object & name)
{
std::string name_str = cast<std::string>(name);
handle(retrieve_attr(o.ptr(), name_str.c_str()));
}
)
;
}
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 | >>> class Cls():
>>> pass
>>> obj = Cls()
>>> val = 'attached value'
>>> print(sys.getrefcount(val))
3
>>>
>>> attach_attr_by_string(obj, 'name', val)
>>> print(sys.getrefcount(val))
4
>>>
>>> print(obj.name is val)
True
>>> print(sys.getrefcount(val))
4
>>>
>>> val2 = retrieve_attr_by_string(obj, 'name')
>>> print(sys.getrefcount(val))
5
|
See also the documentation of Object Protocol.
Function Call¶
This section shows how to make Python function call from C.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 | #include "pybind11/pybind11.h"
#include <string>
using namespace pybind11;
PyObject * function_call(PyObject * callable, PyObject * args, PyObject * kw)
{
PyObject * ret = PyObject_Call(callable, args, kw);
return ret;
}
PYBIND11_MODULE(ex_call, m)
{
m
.def
(
"function_call"
, [](object & o, tuple & t, dict & kw)
{
return handle(function_call(o.ptr(), t.ptr(), kw.ptr()));
}
)
;
}
|
1 2 3 4 5 6 7 8 9 10 11 12 | >>> def my_func(arg1, kw1='default'):
>>> return 'results: {}, {}'.format(arg1, kw1)
>>>
>>> print('(direct call) ', my_func('first argument'))
(direct call) results: first argument, default
>>> print('(function_call)', function_call(my_func, ('first argument',), {}))
(function_call) results: first argument, default
>>>
>>> print('(direct call) ', my_func('first argument', kw1='non default'))
(direct call) results: first argument, non default
>>> print('(function_call)', function_call(my_func, ('first argument',), {'kw1': 'non default'}))
(function_call) results: first argument, non default
|
See the documentation of Object Protocol for other variants of the API.
Python C API for tuple
¶
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 | #include "pybind11/pybind11.h"
#include <string>
using namespace pybind11;
PyObject * reverse_tuple(PyObject * tup)
{
PyObject * ret = PyTuple_New(PyTuple_Size(tup));
for (Py_ssize_t i = 0 ; i < PyTuple_Size(tup) ; ++i)
{
PyObject * item = PyTuple_GetItem(tup, i);
Py_INCREF(item);
PyTuple_SetItem(ret, i, item); // This only works when 1 == Py_REFCNT(ret)
}
return ret;
}
PYBIND11_MODULE(ex_tuple, m)
{
m
.def
(
"reverse_tuple"
, [](tuple & t)
{
return handle(reverse_tuple(t.ptr()));
}
)
;
}
|
>>> tv0 = "value0"
>>> tv1 = object()
>>> tup = (tv0, tv1)
>>> print(sys.getrefcount(tv1))
3
>>> rtup = reverse_tuple(tup)
>>> print(sys.getrefcount(tv1))
4
See the C API documentation for the tuple protocol
and the code implementing
PyTuple_SetItem()
.
Python C API for dict
¶
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 | #include "pybind11/pybind11.h"
#include <string>
using namespace pybind11;
PyObject * make_dict()
{
PyObject * ret = PyDict_New();
return ret;
}
void add_dict_item(PyObject * d, PyObject * k, PyObject * v)
{
/*int ret =*/
PyDict_SetItem(d, k, v);
}
PYBIND11_MODULE(ex_dict, m)
{
m
.def
(
"make_dict"
, []()
{
return handle(make_dict());
}
)
.def
(
"add_dict_item"
, [](dict & d, object & k, object & v)
{
add_dict_item(d.ptr(), k.ptr(), v.ptr());
}
)
;
}
|
>>> d0 = {}
>>> d1 = make_dict()
>>> print(d0 is d1)
False
>>> print(d0 == d1)
True
>>> d0['k1'] = 'value1'
>>> print(d0)
{'k1': 'value1'}
>>> add_dict_item(d1, 'k1', 'value1')
>>> print(d1)
{'k1': 'value1'}
>>> print(d0 == d1)
True
Python C API for list
¶
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 | #include "pybind11/pybind11.h"
#include <string>
using namespace pybind11;
PyObject * make_list_from_iterator(PyObject * o)
{
PyObject * iter = PyObject_GetIter(o);
PyObject * ret = PyList_New(0);
PyObject * item = nullptr;
while (nullptr != (item = PyIter_Next(iter)))
{
PyList_Append(ret, item);
Py_DECREF(item);
}
Py_DECREF(iter);
return ret;
}
PYBIND11_MODULE(ex_list, m)
{
m
.def
(
"make_list_from_iterator"
, [](object & o)
{
PyObject * ret = make_list_from_iterator(o.ptr());
return handle(ret);
}
)
;
}
|
1 2 3 4 5 6 7 8 9 10 | >>> v0 = 'first value'
>>> v1 = 'second value'
>>> tup = (v0, v1)
>>> print(sys.getrefcount(v1))
4
>>> lst = make_list_from_iterator(tup)
>>> print(sys.getrefcount(v1))
5
>>> print(lst)
['first value', 'second value']
|
See the C API documentation for the list protocol and the C API documentation for the iterator protocol.
Useful Operations¶
Import¶
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 | #include "pybind11/pybind11.h"
#include <string>
using namespace pybind11;
PyObject * get_modules()
{
PyObject * sysmod = PyImport_ImportModule("sys");
PyObject * modules = PyObject_GetAttrString(sysmod, "modules");
Py_DECREF(sysmod);
return modules;
}
PYBIND11_MODULE(ex_import, m)
{
m
.def
(
"get_modules"
, []()
{
PyObject * ret = get_modules();
return handle(ret);
}
)
;
}
|
>>> modules = get_modules();
>>> print(type(modules), len(modules))
<class 'dict'> 1146
Exception¶
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 | #include "pybind11/pybind11.h"
#include <string>
using namespace pybind11;
PyObject * function_with_exception(PyObject * o)
{
// Normal processing code.
PyObject * ret = nullptr;
PyObject * item = nullptr;
PyObject * iter = PyObject_GetIter(o);
if (nullptr == iter) { goto error; }
ret = PyList_New(0);
if (nullptr == ret) { goto error; }
while (nullptr != (item = PyIter_Next(iter)))
{
int status = PyList_Append(ret, item);
Py_DECREF(item);
if (0 != status) { goto error; }
}
Py_DECREF(iter);
// Exception.
PyErr_SetString(PyExc_RuntimeError, "intentional exception");
error: // A good use of goto: clean up.
Py_XDECREF(iter);
Py_XDECREF(ret);
return nullptr;
}
PYBIND11_MODULE(ex_except, m)
{
m
.def
(
"function_with_exception"
, [](object & o)
{
PyObject * ret = function_with_exception(o.ptr());
if (nullptr == ret) { throw error_already_set(); }
return handle(ret);
}
)
;
}
|
>>> try:
>>> function_with_exception(1)
>>> except TypeError as e:
>>> print(e)
>>> else:
>>> print('error not raised')
'int' object is not iterable
>>> tup = ('first value', 'second value')
>>> try:
>>> function_with_exception(('first value', 'second value'))
>>> except RuntimeError as e:
>>> print(e)
>>> else:
>>> print('error not raised')
intentional exception
See also Exceptions and Exception Handling.
Python Memory Management¶
Python has its own memory manager. When writing Python extension, they should
be used for PyObject
. The memory managing system has three levels:
- Raw memory interface: wrapper to the C standard memory managers. It allows
distinct addressed returned when requesting 0 byte. GIL is not involved.
- Normal memory interface: ‘pymalloc’ with small memory optimization. GIL is
required when calling.
- Object memory interface: allocate for
PyObject
. GIL is
required when calling.
The public API are:
void * PyMem_RawMalloc(size_t size);
void * PyMem_RawCalloc(size_t nelem, size_t elsize);
void * PyMem_RawRealloc(void *ptr, size_t new_size);
void PyMem_RawFree(void *ptr);
void * PyMem_Malloc(size_t size);
void * PyMem_Calloc(size_t nelem, size_t elsize);
void * PyMem_Realloc(void *ptr, size_t new_size);
void PyMem_Free(void *ptr);
void * PyObject_Malloc(size_t size);
void * PyObject_Calloc(size_t nelem, size_t elsize);
void * PyObject_Realloc(void *ptr, size_t new_size);
void PyObject_Free(void *ptr);
In [Include/cpython/pymem.h](https://github.com/python/cpython/blob/v3.8.0/Include/cpython/pymem.h#L53), Python provides a struct and a set of API to switch to custom memory managers:
typedef struct {
/* user context passed as the first argument to the 4 functions */
void *ctx;
/* allocate a memory block */
void* (*malloc) (void *ctx, size_t size);
/* allocate a memory block initialized by zeros */
void* (*calloc) (void *ctx, size_t nelem, size_t elsize);
/* allocate or resize a memory block */
void* (*realloc) (void *ctx, void *ptr, size_t new_size);
/* release a memory block */
void (*free) (void *ctx, void *ptr);
} PyMemAllocatorEx;
/* Get the memory block allocator of the specified domain. */
PyAPI_FUNC(void) PyMem_GetAllocator(PyMemAllocatorDomain domain,
PyMemAllocatorEx *allocator);
/* Set the memory block allocator of the specified domain.
The new allocator must return a distinct non-NULL pointer when requesting
zero bytes.
For the PYMEM_DOMAIN_RAW domain, the allocator must be thread-safe: the GIL
is not held when the allocator is called.
If the new allocator is not a hook (don't call the previous allocator), the
PyMem_SetupDebugHooks() function must be called to reinstall the debug hooks
on top on the new allocator. */
PyAPI_FUNC(void) PyMem_SetAllocator(PyMemAllocatorDomain domain,
PyMemAllocatorEx *allocator);
See the official documentation customize-memory-allocators. The public API is wrappers to the functions populated in the struct, e.g.:
void *
PyMem_RawMalloc(size_t size)
{
/*
* Limit ourselves to PY_SSIZE_T_MAX bytes to prevent security holes.
* Most python internals blindly use a signed Py_ssize_t to track
* things without checking for overflows or negatives.
* As size_t is unsigned, checking for size < 0 is not required.
*/
if (size > (size_t)PY_SSIZE_T_MAX)
return NULL;
return _PyMem_Raw.malloc(_PyMem_Raw.ctx, size);
}
Also see the code.
Small Memory Optimization¶
Take a look at the documentation in `the code<https://github.com/python/cpython/blob/v3.8.0/Objects/obmalloc.c#L766>`__. This is the ‘pymalloc’, and it uses 256 KB for allocation not greater than 512 bytes.
Tracemalloc¶
Tracemalloc is a standard library that uses the custom memory manager to profile and debug Python memory use: tracemalloc — Trace memory allocations. We can follow the implementation to create more specific analysis.
Exercises¶
- Use pybind11 to expose a memory buffer to Python as a numpy ndarray.
References¶
[1] | Project turgon (work in progress): https://github.com/yungyuc/turgon. |
[2] | S.C. Chang, “The Method of Space-Time Conservation Element and Solution Element – A New Approach for Solving the Navier-Stokes and Euler Equations,” J. Comput. Phys., 119, pp. 295-324, (1995). DOI: 10.1006/jcph.1995.1137 |