MWOperators¶
The MW operators discussed in this chapter is available to the application program by including:
#include "MRCPP/MWOperators"
ConvolutionOperator¶
Note
The convolution operators have separate precision parameters for their
construction and application. The build_prec
argument to the operator
constructors will affect e.g. the number of terms in the separated
representations of the Poisson/Helmholtz approximations, as well as the
operator bandwidth. The apply_prec
argument to the apply function relates
only to the adaptive construction of the output function, based on a wavelet
norm error estimate.
-
template<int
D
>
classmrcpp
::
IdentityConvolution
: public mrcpp::ConvolutionOperator<D>¶ Convolution with an identity kernel.
The identity kernel (Dirac’s delta function) is approximated by a narrow Gaussian function: \( I(r-r') = \delta(r-r') \approx \alpha e^{-\beta (r-r')^2} \)
Public Functions
-
IdentityConvolution
(const MultiResolutionAnalysis<D> &mra, double pr)¶ This will project a kernel of a single gaussian with exponent sqrt(10/build_prec).
- Return
New IdentityConvolution object
- Parameters
[in] mra
: Which MRA the operator is defined[in] pr
: Build precision, closeness to delta function
-
-
template<int
D
>
classmrcpp
::
DerivativeConvolution
: public mrcpp::ConvolutionOperator<D>¶ Convolution with a derivative kernel.
The derivative kernel (derivative of Dirac’s delta function) is approximated by the derivative of a narrow Gaussian function: \( D^x(r-r') = \frac{d}{dx}\delta(r-r') \approx \frac{d}{dx} \alpha e^{-\beta (r-r')^2} \)
Public Functions
-
DerivativeConvolution
(const MultiResolutionAnalysis<D> &mra, double pr)¶ This will project a kernel of a single differentiated gaussian with exponent sqrt(10/build_prec).
- Return
New DerivativeConvolution object
- Parameters
[in] mra
: Which MRA the operator is defined[in] pr
: Build precision, closeness to delta function
-
-
class
mrcpp
::
PoissonOperator
: public mrcpp::ConvolutionOperator<3>¶ Convolution with the Poisson Green’s function kernel.
The Poisson kernel is approximated as a sum of Gaussian functions in order to allow for separated application of the operator in the Cartesian directions: \( P(r-r') = \frac{1}{|r-r'|} \approx \sum_m^M \alpha_m e^{-\beta_m (r-r')^2} \)
Public Functions
-
PoissonOperator
(const MultiResolutionAnalysis<3> &mra, double pr = -1.0)¶ This will construct a gaussian expansion to approximate 1/r, and project each term into a one-dimensional MW operator. Subsequent application of this operator will apply each of the terms to the input function in all Cartesian directions.
- Return
New PoissonOperator object
- Parameters
[in] mra
: Which MRA the operator is defined[in] pr
: Build precision, closeness to 1/r
-
-
class
mrcpp
::
HelmholtzOperator
: public mrcpp::ConvolutionOperator<3>¶ Convolution with the Helmholtz Green’s function kernel.
The Helmholtz kernel is approximated as a sum of gaussian functions in order to allow for separated application of the operator in the Cartesian directions: \( H(r-r') = \frac{e^{-\mu|r-r'|}}{|r-r'|} \approx \sum_m^M \alpha_m e^{-\beta_m (r-r')^2} \)
Public Functions
-
HelmholtzOperator
(const MultiResolutionAnalysis<3> &mra, double m, double pr = -1.0)¶ This will construct a gaussian expansion to approximate exp(-mu*r)/r, and project each term into a one-dimensional MW operator. Subsequent application of this operator will apply each of the terms to the input function in all Cartesian directions.
- Return
New HelmholtzOperator object
- Parameters
[in] mra
: Which MRA the operator is defined[in] m
: Exponential parameter of the operator[in] pr
: Build precision, closeness to exp(-mu*r)/r
-
-
template<int
D
>
voidmrcpp
::
apply
(double prec, FunctionTree<D> &out, ConvolutionOperator<D> &oper, FunctionTree<D> &inp, int maxIter, bool absPrec)¶ Application of MW integral convolution operator.
The output function will be computed using the general algorithm:
Compute MW coefs on current grid
Refine grid where necessary based on
prec
Repeat until convergence or
maxIter
is reachedprec < 0
ormaxIter = 0
means NO refinementmaxIter < 0
means no bound
- Parameters
[in] prec
: Build precision of output function[out] out
: Output function to be built[in] oper
: Convolution operator to apply[in] inp
: Input function[in] maxIter
: Maximum number of refinement iterations in output tree, default -1[in] absPrec
: Build output tree based on absolute precision, default false
- Note
This algorithm will start at whatever grid is present in the
out
tree when the function is called (this grid should however be EMPTY, e.i. no coefs).
-
template<int
D
>
voidmrcpp
::
apply
(double prec, FunctionTree<D> &out, ConvolutionOperator<D> &oper, FunctionTree<D> &inp, FunctionTreeVector<D> &precTrees, int maxIter, bool absPrec)¶ Application of MW integral convolution operator.
The output function will be computed using the general algorithm:
Compute MW coefs on current grid
Refine grid where necessary based on scaled
prec
Repeat until convergence or
maxIter
is reachedprec < 0
ormaxIter = 0
means NO refinementmaxIter < 0
means no bound
- Parameters
[in] prec
: Build precision of output function[out] out
: Output function to be built[in] oper
: Convolution operator to apply[in] inp
: Input function[in] precTrees
: Precision trees[in] maxIter
: Maximum number of refinement iterations in output tree, default -1[in] absPrec
: Build output tree based on absolute precision, default false
The precision will be scaled locally by the maxNorms of the precTrees input vector.
- Note
This algorithm will start at whatever grid is present in the
out
tree when the function is called (this grid should however be EMPTY, e.i. no coefs).
DerivativeOperators¶
Note
The derivative operators have clearly defined requirements on the output grid structure, based on the grid of the input function. This means that there is no real grid adaptivity, and thus no precision parameter is needed for the application of such an operator.
-
template<int
D
>
classmrcpp
::
ABGVOperator
: public mrcpp::DerivativeOperator<D>¶ Derivative operator as defined by Alpert, Beylkin, Ginez and Vozovoi, J Comp Phys 182, 149-190 (2002).
Public Functions
-
ABGVOperator
(const MultiResolutionAnalysis<D> &mra, double a, double b)¶ Boundary parameters correspond to:
a=0.0
b=0.0
: Strictly local “center” differencea=0.5
b=0.5
: Semi-local central differencea=1.0
b=0.0
: Semi-local forward differencea=0.0
b=1.0
: Semi-local backward difference
- Return
New ABGVOperator object
- Parameters
[in] mra
: Which MRA the operator is defined[in] a
: Left boundary condition[in] b
: Right boundary condition
-
-
template<int
D
>
classmrcpp
::
PHOperator
: public mrcpp::DerivativeOperator<D>¶ Derivative operator based on the smoothing derivative of Pavel Holoborodko .
Public Functions
-
PHOperator
(const MultiResolutionAnalysis<D> &mra, int order)¶ - Return
New PHOperator object
- Parameters
[in] mra
: Which MRA the operator is defined[in] order
: Derivative order, defined for 1 and 2
-
-
template<int
D
>
classmrcpp
::
BSOperator
: public mrcpp::DerivativeOperator<D>¶ B-spline derivative operator as defined by Anderson etal, J Comp Phys X 4, 100033 (2019).
Public Functions
-
BSOperator
(const MultiResolutionAnalysis<D> &mra, int order)¶ - Return
New BSOperator object
- Parameters
[in] mra
: Which MRA the operator is defined[in] order
: Derivative order, defined for 1, 2 and 3
-
-
template<int
D
>
voidmrcpp
::
apply
(FunctionTree<D> &out, DerivativeOperator<D> &oper, FunctionTree<D> &inp, int dir)¶ Application of MW derivative operator.
The output function will be computed on a FIXED grid that is predetermined by the type of derivative operator. For a strictly local operator (ABGV_00), the grid is an exact copy of the input function. For operators that involve also neighboring nodes (ABGV_55, PH, BS) the base grid will be WIDENED by one node in the direction of application (on each side).
- Parameters
[out] out
: Output function to be built[in] oper
: Derivative operator to apply[in] inp
: Input function[in] dir
: Direction of derivative
- Note
The output function should contain only empty root nodes at entry.
-
template<int
D
>
voidmrcpp
::
divergence
(FunctionTree<D> &out, DerivativeOperator<D> &oper, FunctionTreeVector<D> &inp)¶ Calculation of divergence of a function vector.
The derivative operator is applied in each Cartesian direction to the corresponding components of the input vector and added up to the final output. The grid of the output is fixed as the union of the component grids (including any derivative widening, see derivative apply).
- Parameters
[out] out
: Output function[in] oper
: Derivative operator to apply[in] inp
: Input function vector
- Note
The length of the input vector must be the same as the template dimension D.
The output function should contain only empty root nodes at entry.
-
template<int
D
>
FunctionTreeVector<D>mrcpp
::
gradient
(DerivativeOperator<D> &oper, FunctionTree<D> &inp)¶ Calculation of gradient vector of a function.
The derivative operator is applied in each Cartesian direction to the input function and appended to the output vector.
- Return
FunctionTreeVector containing the gradient
- Parameters
[in] oper
: Derivative operator to apply[in] inp
: Input function
- Note
The length of the output vector will be the template dimension D.
Examples¶
PoissonOperator¶
The electrostatic potential \(g\) arising from a charge distribution \(f\) are related through the Poisson equation
This equation can be solved with respect to the potential by inverting the differential operator into the Green’s function integral convolution operator
This operator is available in the MW representation, and can be solved with
arbitrary (finite) precision in linear complexity with respect to system size.
Given an arbitrary charge dirtribution f_tree
in the MW representation, the
potential is computed in the following way:
double apply_prec; // Precision for operator application
double build_prec; // Precision for operator construction
mrcpp::PoissonOperator P(MRA, build_prec); // MW representation of Poisson operator
mrcpp::FunctionTree<3> f_tree(MRA); // Input function
mrcpp::FunctionTree<3> g_tree(MRA); // Output function
mrcpp::apply(apply_prec, g_tree, P, f_tree); // Apply operator adaptively
The Coulomb self-interaction energy can now be computed as the dot product:
double E = mrcpp::dot(g_tree, f_tree);
HelmholtzOperator¶
The Helmholtz operator is a generalization of the Poisson operator and is given as the integral convolution
The operator is the inverse of the shifted Laplacian
and appears e.g. when solving the SCF equations. The construction and application is similar to the Poisson operator, with an extra argument for the \(\mu\) parameter
double apply_prec; // Precision for operator application
double build_prec; // Precision for operator construction
double mu; // Must be a positive real number
mrcpp::HelmholtzOperator H(MRA, mu, build_prec);// MW representation of Helmholtz operator
mrcpp::FunctionTree<3> f_tree(MRA); // Input function
mrcpp::FunctionTree<3> g_tree(MRA); // Output function
mrcpp::apply(apply_prec, g_tree, H, f_tree); // Apply operator adaptively
ABGVOperator¶
The ABGV (Alpert, Beylkin, Gines, Vozovoi) derivative operator is initialized with two parameters \(a\) and \(b\) accounting for the boundary conditions between adjacent nodes, see Alpert et al.
double a = 0.0, b = 0.0; // Boundary conditions for operator
mrcpp::ABGVOperator<3> D(MRA, a, b); // MW derivative operator
mrcpp::FunctionTree<3> f(MRA); // Input function
mrcpp::FunctionTree<3> f_x(MRA); // Output function
mrcpp::FunctionTree<3> f_y(MRA); // Output function
mrcpp::FunctionTree<3> f_z(MRA); // Output function
mrcpp::apply(f_x, D, f, 0); // Operator application in x direction
mrcpp::apply(f_y, D, f, 1); // Operator application in y direction
mrcpp::apply(f_z, D, f, 2); // Operator application in z direction
The tree structure of the output function will depend on the choice of parameters \(a\) and \(b\): if both are zero, the output grid will be identical to the input grid; otherwise the grid will be widened by one node (on each side) in the direction of application.
PHOperator¶
The PH derivative operator is based on the noise reducing derivative of Pavel Holoborodko. This operator is also available as a direct second derivative.
mrcpp::PHOperator<3> D1(MRA, 1); // MW 1st derivative operator
mrcpp::PHOperator<3> D2(MRA, 2); // MW 2nd derivative operator
mrcpp::FunctionTree<3> f(MRA); // Input function
mrcpp::FunctionTree<3> f_x(MRA); // Output function
mrcpp::FunctionTree<3> f_xx(MRA); // Output function
mrcpp::apply(f_x, D1, f, 0); // Operator application in x direction
mrcpp::apply(f_xx, D2, f, 0); // Operator application in x direction
Special thanks to Prof. Robert J. Harrison (Stony Brook University) for sharing the operator coefficients.
BSOperator¶
The BS derivative operator is based on a pre-projection onto B-splines in order to remove the discontinuities in the MW basis, see Anderson et al. This operator is also available as a direct second and third derivative.
mrcpp::BSOperator<3> D1(MRA, 1); // MW 1st derivative operator
mrcpp::BSOperator<3> D2(MRA, 2); // MW 2nd derivative operator
mrcpp::BSOperator<3> D3(MRA, 3); // MW 3nd derivative operator
mrcpp::FunctionTree<3> f(MRA); // Input function
mrcpp::FunctionTree<3> f_x(MRA); // Output function
mrcpp::FunctionTree<3> f_yy(MRA); // Output function
mrcpp::FunctionTree<3> f_zzz(MRA); // Output function
mrcpp::apply(f_x, D1, f, 0); // Operator application in x direction
mrcpp::apply(f_yy, D2, f, 1); // Operator application in x direction
mrcpp::apply(f_zzz, D3, f, 2); // Operator application in x direction