This tutorial gives some hints about the spectral basis management
The function std_base, sets the standard basis to a scalar field. For most of the spaces, this basis assumes that the field can be expanded into polynomials of the form x^m y ^n z^2p (i.e. polynomials of the Cartesian coordinates).
Such form ensures regularity everywhere (for instance at the origin or on the axis of a spherical coordinate system). One can also notice that an even symmetry is assumed with respect to the plane z=0.
In each domain, the basis can be accessed. The result is given into the form of several arrays. For instance, in a spherical domain, the three coordinates are r, theta and phi and the coefficients are taken in the reverse order, first with respect to phi, then with respect to theta and then with respect to r. It follows that there is only one basis for phi, a one dimensional array of basis for theta (one for each index of phi) and a two dimensional array for r (one for each pair of angular indices). The various bases are stored as integers, the translation can be found in the file "spectral_basis.hpp"
// Assume a scalar field func has been previously computed (see tutorial 3)
// One sets the basis to the standard one
func.std_base() ;
// Print the basis in domain 1 (the first shell)
cout << "Basis in domain 1" << endl ;
cout << func(1).get_base() << endl ;
// The result is read as follows
// Variable 2 is phi, only one basis, labelled 4, which is a full sequence of sines and cosines
// Variable 1 is theta, an array of basis, with values 6 and 10.
// It corresponds to either even cosines or odd sines, depending on the index of phi.
// Variable 0 is r, a two dimensional array with values 1
// It corresponds to Chebyshev polynomials, not matter what the indices of the angles are.
// Print the basis in domain 0 (the nucleus)
cout << "Basis in domain 0" << endl ;
cout << func(0).get_base() << endl ;
// It differs from the previous case because of regularity conditions at the origin.
// It implies that the basis with respect to r are even or odd Chebyshev polynomials.
Kadath tries to keep track of the basis when fields are manipulated. For instance, when summing two fields, it will check that they have the same basis and assign it to the result. When multiplying fields the code is usually able to also guess the right basis. When calling operators like the derivative or the multiplication by r, the library is also able to compute the basis of the result.
However there are some cases, where the basis can not be obtained automatically. For instance if one takes the cosine of a quantity, it does not have a naturally define basis (the cosine of a Chebyshev polynomial is not a Chebyshev polynomial). In those cases, the basis has to be set by hand using functions like std_base.
Important rule : do not invoke std_base unless you really need to and try to understand why the code has lost track of the spectral basis.
// Combination of fields
Scalar f1 (2*func + func*func) ;
//Base is known and is the standard one ;
cout << "Base of 2*f + f^2 in the nucleus" << endl ;
cout << f1(0).get_base() << endl ;
// It is the standard base
// Taking the derivative
Scalar der (func.der_r()) ;
cout << "Base of of f' in the nucleus" << endl ;
cout << der(0).get_base() << endl ;
// It is not the standard base ; change in the r parity
// Operator cosine
Scalar test (cos(func)) ;
cout << "Base for cos(f) is not defined" << endl ;
cout << test(0).get_base() << endl ;
The standard basis, for scalars, assumes an even symmetry with respect to z=0. So if one defines an odd function, std_base will cause problems. There are a lot of other functions in Kadath that do set the spectral basis. It is a matter of choosing the right one.
For instance for a function of the form : x^m y ^n z^2p+1, one should use std_anti_base.
// Define a scalar field being z in the nucleus and z/r^2 everywhere else
Scalar odd (space) ;
odd.set_domain(0) = space.get_domain(0)->get_cart(3) ; // z coordinate
for (int d=1 ; d<ndom ; d++)
odd.set_domain(d) = space.get_domain(d)->get_cart(3)
/ space.get_domain(d)->get_radius() / space.get_domain(d)->get_radius() ;
// Sets the right value at infinity
odd.set_val_inf(0.) ;
// Try to use std_base()
odd.std_base() ;
// Compare the numerical value to the analytical one using val_point
Point M(3) ;
M.set(1) = 2.32 ;
M.set(2) = 1.32 ;
M.set(3) = 1.98 ;
double rr = sqrt(M(1)*M(1)+M(2)*M(2)+M(3)*M(3)) ; // Get the radius ; check it is not in the nucleus
cout << "Analytical value = " << M(3) / rr / rr << endl ;
cout << "Wrong numerical value = " << odd.val_point(M) << endl ;
// Put the right spectral basis :
odd.std_anti_base() ;
odd.set_in_conf() ; // One needs to kill the previously computed (wrong) coefficients.
cout << "True numerical value = " << odd.val_point(M) << endl ;