KADATH Tutorials

1. Basic

Construct an Array

Arrays are constructed using templates. They can be of any dimension. For 1, 2 or 3 dimensions one can use direct constructors.

For arbitrary dimensions, one must pass the sizes of the array using the class Dim_array.

    // 1d array of 4 integers
    Array<int> onedint (4) ;
    // 2d array of (8x4) doubles :
    Array<double> twoddoubles (8, 4) ;

    // Array of dimension 5 (of booleans), constructed by a Dim_array :
    Dim_array dimensions(5) ;
    // The size of the various dimensions must be initialized
    for (int i=0 ; i<5 ; i++)
           dimensions.set(i) = i+1 ;
    Array<bool> multidbool (dimensions) ;

    // 1d array 
    for (int i=0 ; i<4 ; i++)
              onedint.set(i) = i ;
    cout << onedint << endl ; // Print the whole Array`enter code here`
    // 2d array, same value everywhere
    twoddoubles = 2.3 ;
    // Printing one particular element
    cout << twoddoubles(0, 2) << endl ;

0 1 2 3
2.3

Arrays are constructed uninitialized. They can be given a single value for the whole array or be accessed element by elements, using the set function.

If the dimension is arbitrary, they must be accessed via a class called Index.

    // Array of dimension 5  :
    Index pos(dimensions) ;
    // By default it is affected to the first element
    // Modify the index in dimension 2
    pos.set(2) = 1 ;
    cout << pos << endl ;

(0, 0, 1, 0, 0) in an array of (1, 2, 3, 4, 5) points.

    multidbool = false ; // Everything is false
    multidbool.set(pos) = true ; // but this one...
    // Sets the index to the first index
    pos.set_start() ;
    // Loop
    do {
          if (multidbool(pos)) {
                 cout << "True value found at " << endl ;
                 cout << pos << endl ;
        }
         }
    while (pos.inc()) ;

(0, 0, 1, 0, 0) in an array of (1, 2, 3, 4, 5) points.

2. Space

• the number of domains
• The type of basis (Chebyshev or Legendre)
• The number of points in each domain

Warning: points cannot be arbitrary and depends on the spectral basis. Rule of thumb : for phi variables use integers of the form $2^m 3^n 5^p$ and $2^m 3^n 5^p +1$ for other coordinates. For instance in the radial direction, 11, 13, 17, 21, 25, 33 are admissible numbers.

    // 3D :
    int dim = 3 ;
    // Number of points
    Dim_array res (dim) ;
    res.set(0) = 13 ; // Points in r 
    res.set(1) = 5 ; // Points in theta
    res.set(2) = 4 ; // Points in phi

    // Absolute coordinates of the center
    Point center (dim) ;
    for (int i=1 ; i<=dim ; i++)
         center.set(i) = 0 ;

    // Number of domains and boundaries :
    int ndom = 4 ;
    Array<double> bounds (ndom-1) ;
    bounds.set(0) = 1 ; bounds.set(1) = 2 ; bounds.set(2) = 4 ;
    // Chebyshev or Legendre :
    int type_coloc = CHEB_TYPE ;

    // Spherical space constructor
    Space_spheric space(type_coloc, center, res, bounds) ;
    cout << space << endl ;
-------------------------------------------------------------------
Space of 4 domains
Domain 1
Nucleus
Rmax    = 1
Center  = (0, 0, 0)
Nbr pts = (13, 5, 4)
********************
Domain 2
Shell
1 < r < 2
Center  = (0, 0, 0)
Nbr pts = (13, 5, 4)
********************
Domain 3
Shell
2 < r < 4
Center  = (0, 0, 0)
Nbr pts = (13, 5, 4)
********************
Domain 4
Compactified domain
Rmin    = 4
Center  = (0, 0, 0)
Nbr pts = (13, 5, 4)
********************
Chebyshev polynomials are used.

Individual domains can be accessed via the get_domain function that returns a pointer.

cout<<"-1- Prints the number of points in the domain 0"<<endl;
// Prints the number of points in the domain 0
cout << space.get_domain(0)->get_nbr_points() << endl ;

cout<<"-2- Prints the nulber of coefficients in the domain 1 "<<endl;
// Prints the nulber of coefficients in the domain 1 
cout << space.get_domain(1)->get_nbr_coefs() << endl ;

cout<<"-3- Returns the radius in the domain 0"<<endl;
// Returns the radius in the domain 0
cout << space.get_domain(0)->get_radius() << endl ;

cout<<"-4- Returns the X coordinate in the domain 1"<<endl;
// Returns the X coordinate in the domain 1
cout << space.get_domain(1)->get_cart(1) << endl ;
-------------------------------------------------------------------
-1- Prints the number of points in the domain 0
(13, 5, 4)
-2- Prints the nulber of coefficients in the domain 1
(13, 5, 6)
-3- Returns the radius in the domain 0
Configuration space :
Array of 3 dimension(s)
      0
0 : 0 0.130526 0.258819 0.382683 0.5 0.608761 0.707107 0.793353 0.866025 0.92388 0.965926 0.991445 1
1 : 0 0.130526 0.258819 0.382683 0.5 0.608761 0.707107 0.793353 0.866025 0.92388 0.965926 0.991445 1
2 : 0 0.130526 0.258819 0.382683 0.5 0.608761 0.707107 0.793353 0.866025 0.92388 0.965926 0.991445 1
3 : 0 0.130526 0.258819 0.382683 0.5 0.608761 0.707107 0.793353 0.866025 0.92388 0.965926 0.991445 1
4 : 0 0.130526 0.258819 0.382683 0.5 0.608761 0.707107 0.793353 0.866025 0.92388 0.965926 0.991445 1
      1
0 : 0 0.130526 0.258819 0.382683 0.5 0.608761 0.707107 0.793353 0.866025 0.92388 0.965926 0.991445 1
1 : 0 0.130526 0.258819 0.382683 0.5 0.608761 0.707107 0.793353 0.866025 0.92388 0.965926 0.991445 1
2 : 0 0.130526 0.258819 0.382683 0.5 0.608761 0.707107 0.793353 0.866025 0.92388 0.965926 0.991445 1
3 : 0 0.130526 0.258819 0.382683 0.5 0.608761 0.707107 0.793353 0.866025 0.92388 0.965926 0.991445 1
4 : 0 0.130526 0.258819 0.382683 0.5 0.608761 0.707107 0.793353 0.866025 0.92388 0.965926 0.991445 1
      2
0 : 0 0.130526 0.258819 0.382683 0.5 0.608761 0.707107 0.793353 0.866025 0.92388 0.965926 0.991445 1
1 : 0 0.130526 0.258819 0.382683 0.5 0.608761 0.707107 0.793353 0.866025 0.92388 0.965926 0.991445 1
2 : 0 0.130526 0.258819 0.382683 0.5 0.608761 0.707107 0.793353 0.866025 0.92388 0.965926 0.991445 1
3 : 0 0.130526 0.258819 0.382683 0.5 0.608761 0.707107 0.793353 0.866025 0.92388 0.965926 0.991445 1
4 : 0 0.130526 0.258819 0.382683 0.5 0.608761 0.707107 0.793353 0.866025 0.92388 0.965926 0.991445 1
      3
0 : 0 0.130526 0.258819 0.382683 0.5 0.608761 0.707107 0.793353 0.866025 0.92388 0.965926 0.991445 1
1 : 0 0.130526 0.258819 0.382683 0.5 0.608761 0.707107 0.793353 0.866025 0.92388 0.965926 0.991445 1
2 : 0 0.130526 0.258819 0.382683 0.5 0.608761 0.707107 0.793353 0.866025 0.92388 0.965926 0.991445 1

-4- Returns the X coordinate in the domain 1
Configuration space :
Array of 3 dimension(s)
      0
0 : 0 0 0 0 0 0 0 0 0 0 0 0 0
1 : 0.382683 0.389203 0.408318 0.438726 0.478354 0.524502 0.574025 0.623548 0.669696 0.709324 0.739732 0.758847 0.765367
2 : 0.707107 0.719154 0.754474 0.81066 0.883883 0.969154 1.06066 1.15217 1.23744 1.31066 1.36685 1.40217 1.41421
3 : 0.92388 0.93962 0.985768 1.05918 1.15485 1.26626 1.38582 1.50538 1.61679 1.71246 1.78587 1.83202 1.84776
4 : 1 1.01704 1.06699 1.14645 1.25 1.37059 1.5 1.62941 1.75 1.85355 1.93301 1.98296 2
      1

3. Scalar Fields

Construction and assignment

Scalar fields are constructed from the Space

    // Assume a spherical space has been previously constructed
    // Recover the number of domains if need be
    //int ndom = space.get_nbr_domains() ;
    
    Scalar func (space) ;
    // Puts 1-r^2 in domain 0
    func.set_domain(0) = 1 - pow(space.get_domain(0)->get_radius(),2) ;
    // Puts 1/r in the other ones
    for (int d=1 ; d<ndom ; d++)
             func.set_domain(d)  = 1/space.get_domain(d)->get_radius() ;

    // Prints the values at each collocation point
    cout << "Scalar function func " << endl ;
    cout << func << endl ;
-------------------------------------------------------------------
Scalar function func
Scalar
Domain : 0
Configuration space :
Array of 3 dimension(s)
      0
0 : 1 0.982963 0.933013 0.853553 0.75 0.62941 0.5 0.37059 0.25 0.146447 0.0669873 0.0170371 0
1 : 1 0.982963 0.933013 0.853553 0.75 0.62941 0.5 0.37059 0.25 0.146447 0.0669873 0.0170371 0
2 : 1 0.982963 0.933013 0.853553 0.75 0.62941 0.5 0.37059 0.25 0.146447 0.0669873 0.0170371 0
3 : 1 0.982963 0.933013 0.853553 0.75 0.62941 0.5 0.37059 0.25 0.146447 0.0669873 0.0170371 0
4 : 1 0.982963 0.933013 0.853553 0.75 0.62941 0.5 0.37059 0.25 0.146447 0.0669873 0.0170371 0
      1

Dealing with infinite quantities

Kadath can perform standard operations using infinite quantities. For instance using 1/r up to infinity is fine, because the resulting value is zero.

However more complicated cases can lead to nan (not a number) quantities, when the library is not able to properly compute the result. In that case, one needs to pass the right value by hand, using functions like set_val_inf.

    // Put r exp(-r*r/10.) in the outer domains
    for (int d=1 ; d<ndom ; d++)
           func.set_domain(d)  = space.get_domain(d)->get_radius()* exp(- space.get_domain(d)->get_radius()*space.get_domain(d)->get_radius()/10.) ;

    // Nans appear at infinty
    cout << "Nans at infinity" << endl ;
    cout << func(ndom-1) << endl ;

    // Sets the right value at infinity by hand
    func.set_val_inf(0.) ;
    // No more Nans
    cout << "Nans are removed"  << endl ;
    cout << func(ndom-1) << endl ;
-------------------------------------------------------------------
Nans at infinity
Configuration space :
Array of 3 dimension(s)
      0
0 : 0.807586 0.776883 0.682245 0.521268 0.310222 0.111968 0.0132925 9.40958e-05 1.2195e-10 1.08726e-31 8.37833e-154 0 -nan
1 : 0.807586 0.776883 0.682245 0.521268 0.310222 0.111968 0.0132925 9.40958e-05 1.2195e-10 1.08726e-31 8.37833e-154 0 -nan
2 : 0.807586 0.776883 0.682245 0.521268 0.310222 0.111968 0.0132925 9.40958e-05 1.2195e-10 1.08726e-31 8.37833e-154 0 -nan
3 : 0.807586 0.776883 0.682245 0.521268 0.310222 0.111968 0.0132925 9.40958e-05 1.2195e-10 1.08726e-31 8.37833e-154 0 -nan
4 : 0.807586 0.776883 0.682245 0.521268 0.310222 0.111968 0.0132925 9.40958e-05 1.2195e-10 1.08726e-31 8.37833e-154 0 -nan
      1
      
Nans are removed
Configuration space :
Array of 3 dimension(s)
      0
0 : 0.807586 0.776883 0.682245 0.521268 0.310222 0.111968 0.0132925 9.40958e-05 1.2195e-10 1.08726e-31 8.37833e-154 0 0
1 : 0.807586 0.776883 0.682245 0.521268 0.310222 0.111968 0.0132925 9.40958e-05 1.2195e-10 1.08726e-31 8.37833e-154 0 0
2 : 0.807586 0.776883 0.682245 0.521268 0.310222 0.111968 0.0132925 9.40958e-05 1.2195e-10 1.08726e-31 8.37833e-154 0 0
3 : 0.807586 0.776883 0.682245 0.521268 0.310222 0.111968 0.0132925 9.40958e-05 1.2195e-10 1.08726e-31 8.37833e-154 0 0
4 : 0.807586 0.776883 0.682245 0.521268 0.310222 0.111968 0.0132925 9.40958e-05 1.2195e-10 1.08726e-31 8.37833e-154 0 0
      1

Setting the basis

the scalar field contains no information about the spectral basis. Then Another tutorial is devoted to the basis.

It is often useful to compute the value at a given point (i.e. not a collocation one). This is done by using the function val_point, which takes a Point as an argument. The class Point just contains the absolute (i.e. typically the X,Y,Z) coordinates of the point.

    // Sets the standard base for scalar fields
    func.std_base() ;    

    // Define a Point in a three dimensional space
    Point M (3) ;
    // Sets its coordinates
    M.set(1) = 2.93 ; // X coordinate
    M.set(2) = -3.19 ; // Y coordinate
    M.set(3) = 1.76 ; // Z coordinate
    // Prints the values of the scalar field at MM
    cout << "The point M is " << M << endl ;
    cout << "Field value at M " << func.val_point(M) << endl ;

The point M is (2.93, -3.19, 1.76)
Field value at M 0.525412

Operations on the scalar field

Using the spectral coefficients, Kadath can perform various operations on the scalar field, like taking its radial derivative or multiplying it by the radius. In this latter case, the operation is performed in the coefficient space which prevents the appearance of nan quantities and gives the right value. We can force the code to compute the values at the collocation points by invoking coef_i.

    // Compute the radial derivative
    Scalar derr (func.der_r()) ;
    cout << "Value of the radial derivative at M " << derr.val_point(M) << endl ;

    // Use the multiplication by r :
    Scalar multr (func.mult_r()) ;
    //Force the code to compute the values at the collocation points
    multr.coef_i() ;
    // No nans appear at infinity
    cout << "Value in the last domain" << endl ;
    cout << multr(ndom-1) << endl ;
-------------------------------------------------------------------
Value of the radial derivative at M -0.378456
    
Value in the last domain
Configuration space :
Array of 3 dimension(s)
      0
0 : 3.23034 3.16139 2.92491 2.44281 1.65452 0.711575 0.10634 0.00101563 1.9512e-09 1.08802e-14 3.15303e-14 7.4829e-14 0.0317112
1 : 3.23034 3.16139 2.92491 2.44281 1.65452 0.711575 0.10634 0.00101563 1.9512e-09 1.08802e-14 3.15303e-14 7.4829e-14 0.0317112
2 : 3.23034 3.16139 2.92491 2.44281 1.65452 0.711575 0.10634 0.00101563 1.9512e-09 1.08802e-14 3.15303e-14 7.4829e-14 0.0317112
3 : 3.23034 3.16139 2.92491 2.44281 1.65452 0.711575 0.10634 0.00101563 1.9512e-09 1.08802e-14 3.15303e-14 7.4829e-14 0.0317112
4 : 3.23034 3.16139 2.92491 2.44281 1.65452 0.711575 0.10634 0.00101563 1.9512e-09 1.08802e-14 3.15303e-14 7.4829e-14 0.0317112
      1   
    
Coefficient space :
Array of 3 dimension(s)
      0
0 : 1.0528 -1.65409 0.713435 0.00825485 -0.176361 0.0548712 0.042382 -0.01938 -0.00775406 0.0101449 0.00652708 0.000880517 -0
1 : -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0
2 : -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0
3 : -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0
4 : -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0
      1
0 : -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0
1 : -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0
2 : -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0
3 : -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0
4 : -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0
    
......5

4. Spectral basis

Standard basis

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 basis are stored as integers, the translation can be found in the file "spectral_basis.hpp"

    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 base, 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 ;
-------------------------------------------------------------------
Basis in domain 1
3-dimensional spectral base
Variable 0
Array of 2 dimension(s)
0 : 1 1 1 1 1
1 : 1 1 1 1 1
2 : 1 1 1 1 1
3 : 1 1 1 1 1
4 : 1 1 1 1 1
5 : 1 1 1 1 1
Variable 1
Array of 1 dimension(s)
6 6 10 10 6 6
Variable 2
Array of 1 dimension(s)
4
    
Basis in domain 0
3-dimensional spectral base
Variable 0
Array of 2 dimension(s)
0 : 2 2 2 2 2
1 : 2 2 2 2 2
2 : 3 3 3 3 3
3 : 3 3 3 3 3
4 : 2 2 2 2 2
5 : 2 2 2 2 2
Variable 1
Array of 1 dimension(s)
6 6 10 10 6 6
Variable 2
Array of 1 dimension(s)
4

Spectral basis and computations

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 basis

    // 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 basis ; 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 ;

Base for cos(f) is not defined
3-dimensional spectral base
Base not defined

Example of non standard basis

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 not in the nucleus
    cout << "Analytical value  = " << M(3) / rr / rr << endl ;
    cout << "Wrong numerical value = " << odd.val_point(M) << endl ;

    // Puts 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 ;

Analytical value = 0.179263
Wrong numerical value = 0.187244
True numerical value = 0.179263

5. The System Of Equations

This tutorial introduces the PDE solver of Kadath

System_of_eqs constructor

System_of_eqs is one of the main object of the Kadath library. It is used to describe a system of PDEs and to solve it. It is constructed from the space only and possibly the domains used. (no domain meaning the whole space, in one domain or in a range of domains).

    // Constructor 
    System_of_eqs syst (space, 1, ndom-1) ;
    // Here one will solve the equations in the domains 1 to ndom-1 (so everywhere but in the domain 0)

Defining the unknowns

One needs to say to the solver what are the unknowns of the system. They can be numbers, scalar fields or general tensors. First, one will consider a simple scalar field. This field has to be defined beforehand with the right spectral basis. Its value will be used as an initial guess by the solver. The function to be used is add_var (for add variable field, i.e. a field that the solver is allowed to modify)

Constants are fields or numbers that the solver will use but that it cannot change. As for the unknowns they must be defined beforehand. The function to use is add_cst

    // Construct a scalar field
    Scalar field (space) ;
    field = 1 ;
    field.std_base() ;

    // The field is an unknown
    syst.add_var ("F", field) ;
    // The string is the name by which the field will be refereed to by the system.


    // Define a constant, begin the radius of the first domain
    Index pos (space.get_domain(1)->get_nbr_points()) ; // Index on the domain 1
    double rad = space.get_domain(1)->get_radius()(pos) ;
    cout << "The inner radius of the first shell is " << rad << endl ;
    // This number is a constant for the system
    syst.add_cst ("a", rad) ;
    // The string is the name by which the constant will be refereed to by the system.

The inner radius of the first shell is 1

Bulk equations

One can now start implementing the equations that one wants to solve. In this particular tutorial, we want to solve Lap(F) = 0, in the whole computational domain. This is a second order equation and one must use the function add_eq_inside (inside because it needs to have two boundary conditions, an inner one and an outer one, being a second order equation).

    // The equation is the same in each domain
    for (int d=1 ; d<ndom ; d++) // Loop on the domains
            syst.add_eq_inside (d, "Lap(F)=0") ; 
    // Lap is a reserved word of Kadath that means the flat Laplacian and we have defined F before.

The boundary conditions

The bulk equation being second order, it needs an inner and an outer boundary conditions. One needs to use add_eq_bc. It differs from add_eq_inside because it needs not only the number of the domain but also a description of the boundary involved. Here one will use INNER_BC and OUTER_BC, two Kadath constants that describe, respectively the inner boundary and the outer one.

    // Inner boundary condition
    syst.add_eq_bc (1, INNER_BC, "dn(F) + 0.5*F/a = 0") ;
    // It is a Robin BC
    // dn stands for the normal derivative with respect to the boundary (d/dr in this case)

    // Outer boundary condition
    syst.add_eq_bc (ndom-1, OUTER_BC, "F=1") ;
    // The field is one at spatial infinity

Matching conditions

Last, one needs to state that the solution and its radial derivative are continuous across the numerical boundaries. This is done via the function add_eq_matching which is similar to add_eq_bc.

    // Loop on all the boundaries between the domains
    for (int d=1 ; d<ndom-1 ; d++) {
               // Matching of the field between domain d and d+1
            syst.add_eq_matching (d, OUTER_BC, "F") ;
            // Matching of the normal derivative of the field between domain d and d+1
            syst.add_eq_matching (d, OUTER_BC, "dn(F)") ;
    }

Checking the equations

At this point can check whether the equations are fulfilled or not. One can use check_equations. It will give an array containing the error on the various equations passed to the system

    Array<double> errors (syst.check_equations()) ;
    // Verify the number of equations, being the size of errors
    int neq = errors.get_size(0) ;
    cout << "The system has " << neq << " equations." << endl ;
    // Print the errors, when they are "big"
    for (int n=0 ; n<neq ; n++)
           if (fabs(errors(n)) > 1e-10)
                    cout << "eq " << n << " ; error " << errors(n) << endl ;

The system has 9 equations.
eq 3 ; error 0.5

Launching the solver

The solution is sought iteratively using the Newton-Raphson algorithm. The function to call is do_newton. It takes as an input the error threshold and gives as an output the current error. If the error is smaller than the threshold it just returns TRUE, if not, it does one iteration of the Newton-Raphson algorithm and returns FALSE.

Before doing the iteration the code will check that the number of unknowns (the coefficients of the unknown fields) is consistent with the number of equations. This is up to the user to give the right number of equations.

    // Define the threshold
    double prec = 1e-8 ;
    // The current error
    double error ;
    // End of the iteration ?
    bool endloop = false ;

    while (!endloop) {
            endloop = syst.do_newton(prec, error) ;
    }
    // At the first call the error is big and so the solver does one iteration and endloop is still false
    // At the second one, the error is smaller than the threshold and endloop is true thus ending the loop.
    // In that example the problem is linear and so converges in one iteration only.

Entering do_newton with error 0.5
Entering do_newton with error 6.04239e-13

Checking the solution

The scalar field has been modified by the code and should contain the right analytical solution. One will check that by taking a random point in the computational domain.

    // Random point
    Point MM(3) ;
    MM.set(1) = 1.3873 ;
    MM.set(2) = -0.827 ;
    MM.set(3) = 0.982 ;
    // The numerical solution
    double numerical = field.val_point(MM) ;
    // The analytic one
    // Get the radius
    double rr = sqrt(MM(1)*MM(1)+MM(2)*MM(2)+MM(3)*MM(3)) ;
    double analytic = 1 + rad  / rr ;
    cout << "Numerical = " << numerical << endl ;
    cout << "Analytic  = " << analytic << endl ;
    cout << "Relative difference = " << fabs(numerical-analytic) / numerical << endl ;

Numerical = 1.52904
Analytic = 1.52904
Relative difference = 2.2685e-09

6. Using Definitions

Test problem

We will use the same physical problem as in tutorial 5, except we will work with the logarithm of the previous field. In order to do so one needs to do the following changes :

• The bulk equation becomes $\Delta P + \partial_i P \partial^i P = 0$
• The inner boundary condition is $dr(P) + \frac{1}{2a} = 0$
• The outer boundary condition is $P=0$

Doing so we have transformed a linear problem into a non-linear one.

Definitions in Kadath

Very often an expression appears several times in a system. Instead of rewritting it each time in terms of the variables and constants, one can pass it to the solver as a definition. This is done via the function add_def

    // Assume a system of equation syst and a scalar called field have been defined (as in tutorial 5).
    // Pass the term partial_i P partial^i P as a definition
    syst.add_var ("P", field) ;
    syst.add_def ("dP2 = scal(grad(P), grad(P))") ;
    //scal is a reserved word that denotes the flat scalar product
    //grad is a reserved word that denotes the flat gradient

After having pass a definition to the system, one can use its name in all Kadath equations

    // The bulk equations
    for (int d=1 ; d<ndom ; d++)
             syst.add_eq_inside (d, "Lap(P) + dP2 = 0") ;
    // Inner BC
    syst.add_eq_bc (1, INNER_BC, "dn(P) + 0.5/a = 0") ;
    // Outer BC 
    syst.add_eq_bc (ndom-1, OUTER_BC, "P=0") ;

You will notice that the system being non-linear the Newton-Raphson solver needs several steps to converge.

The system has 9 equations.
eq 3 ; error 0.5
eq 4 ; error 1
Entering do_newton with error 1
Entering do_newton with error 0.102199
Entering do_newton with error 0.00130704
Entering do_newton with error 1.78788e-06
Entering do_newton with error 4.82172e-12
Numerical = 0.424643
Analytic  = 0.424643
Relative difference = 8.53557e-09

Why use definitions ?

There are two main reasons :

• First it makes the writing of the system clearer because one does not need to write everything explicitly in terms of the variables and constants.
• Second, the value of the definition is computed only once per iteration and not each time it is encountered in the equations. This can lead to significant improvement in terms of computational time.

For instance, assume the radial derivative of a field appears many times in the equations. If you only write dr(F), Kadath will compute this derivative each time it is encountered. However, if you pass a definition that contains this derivative, it will be computed only once (per iteration), its value will be stored in the definition and just copied each time its value is needed.

Using definitions to do computations

Kadath allows the user to recover the value of a definition using the function give_val_def. It returns the current value of the definition, in the form of a Tensor (see Tutorial 7). Zeros will be put in the domains where the definition if not known.

    // Let us compute P^2 in the first shell, after the system resolution
    syst.add_def (1, "Psquare = P^2") ; 
    // The definition is only defined in the first shell.

    // Recover its value as a Tensor, that one will convert to a Scalar field
    Scalar P2 (syst.give_val_def("Psquare")) ;
    // Print P2 which is not zero only in the first shell
    cout << P2 << endl ;

    // Compare with the square computed directly
    Scalar P2direct (field*field) ;
    cout << "Direct computation" << endl ;
    cout << P2direct(1) << endl ;
-------------------------------------------------------------------
Scalar
Domain : 0
Null Val_domain
Domain : 1
Configuration space :
Array of 3 dimension(s)
      0
0 : 0.480453 0.468863 0.437256 0.393313 0.345493 0.300192 0.260943 0.229002 0.204291 0.186158 0.173848 0.166729 0.164402
1 : 0.480453 0.468863 0.437256 0.393313 0.345493 0.300192 0.260943 0.229002 0.204291 0.186158 0.173848 0.166729 0.164402
2 : 0.480453 0.468863 0.437256 0.393313 0.345493 0.300192 0.260943 0.229002 0.204291 0.186158 0.173848 0.166729 0.164402
3 : 0.480453 0.468863 0.437256 0.393313 0.345493 0.300192 0.260943 0.229002 0.204291 0.186158 0.173848 0.166729 0.164402
4 : 0.480453 0.468863 0.437256 0.393313 0.345493 0.300192 0.260943 0.229002 0.204291 0.186158 0.173848 0.166729 0.164402
      1
      2
      3
Domain : 2
Null Val_domain

Domain : 3
Null Val_domain

-------------------------------------------------------------------
Direct computation
Configuration space :
Array of 3 dimension(s)
      0
0 : 0.480453 0.468863 0.437256 0.393313 0.345493 0.300192 0.260943 0.229002 0.204291 0.186158 0.173848 0.166729 0.164402
1 : 0.480453 0.468863 0.437256 0.393313 0.345493 0.300192 0.260943 0.229002 0.204291 0.186158 0.173848 0.166729 0.164402
2 : 0.480453 0.468863 0.437256 0.393313 0.345493 0.300192 0.260943 0.229002 0.204291 0.186158 0.173848 0.166729 0.164402
3 : 0.480453 0.468863 0.437256 0.393313 0.345493 0.300192 0.260943 0.229002 0.204291 0.186158 0.173848 0.166729 0.164402
4 : 0.480453 0.468863 0.437256 0.393313 0.345493 0.300192 0.260943 0.229002 0.204291 0.186158 0.173848 0.166729 0.164402
      1

Of course in this example, it is easier to do a direct computation outside the System_of_eqs but this is not always the case.

7. Tensors

Tensorial basis

Before defining some tensors, one needs to describe the tensorial basis used. This object is different from the spectral basis and the two notions should not be confused. In each domains, the tensorial basis describes the basis on which the coordinates of a given tensor are known. In Kadath, there are two mains possibilities :

• Cartesian basis denoted by the reserved word CARTESIAN_BASIS
• Orthonormal spherical basis denoted by the reserved word SPHERICAL_BASIS

It is not because one is working in spherical coordinates than the Cartesian basis cannot be used. It just means that a vector V, will be described by its coordinates Vx, Vy and Vz, each on them being given as a function of r, theta and phi.

Different tensorial basis can also be used in different domains. Some of them are not usable for some spaces (for instance one cannot use a spherical basis in a bispherical space).

    // Assume a spherical space is known and ndom its number of domains
    Base_tensor basis_cart (space, CARTESIAN_BASIS) ;
    // Constructor with the same tensorial basis everywhere
    // Sets an orthonormal spherical basis in the first shell
    Base_tensor basis_mixed (space, CARTESIAN_BASIS) ;
    basis_mixed.set_basis(1) = SPHERICAL_BASIS ;
    cout << "The mixed tensorial basis is " << endl ;
    cout << basis_mixed << endl ;

The mixed tensorial basis is
Tensorial basis
Domain 0 : Cartesian
Domain 1 : Spherical
Domain 2 : Cartesian
Domain 3 : Cartesian

Constructing a vector field

After the scalar fields, vectors are the most simple tensorial fields. They are constructed from a Space. One must specify if it is in a covariant representation (reserved word COV) or a contravariant one (reserved word CON). So it can describe both true tensorial fields and 1-forms. The tensorial basis must also be specified.

    // Construction of a vector
    Vector vecU (space, CON, basis_cart) ;

Vector affectation and spectral basis.

The components of a vector are affected by the function set that returns a given component. Beware that the indices of the components range from 1 to the dimension (and so do not start at 0).

As for the scalar fields, the spectral basis must be provided. It depends on the space used, the tensorial basis and the component considered.

For instance, for a vector in Cartesian coordinates, the function std_base assumes than the vector is the gradient of some scalar field. It follows that the x and y components are symmetric with respect to the plane z=0 and the z one is anti-symmetric, thus resulting in different basis.

    // Sets (x,y,z) in the nucleus :
    for (int cmp=1 ; cmp<=3 ; cmp++)
            vecU.set(cmp).set_domain(0) = space.get_domain(0)->get_cart(cmp) ;

    // Sets (x/r, y/r, z/r) for the other domains
    for (int d=1 ; d<ndom ; d++)
            for (int cmp=1 ; cmp<=3 ; cmp++)
                vecU.set(cmp).set_domain(d) = space.get_domain(d)->get_cart_surr(cmp) ;

    // Sets the appropriate spectral basis (here the standard one)
    vecU.std_base() ;

    // Prints the spectral basis in the nucleus
    for (int cmp=1 ; cmp<=3 ; cmp++) {
            cout << "Basis of cmp " << cmp << " in the nucleus: " << endl ;
            cout << vecU(cmp)(0).get_base() << endl ;
        }
-------------------------------------------------------------------
Basis of cmp 1 in the nucleus:
3-dimensional spectral base
Variable 0
Array of 2 dimension(s)
0 : 2 2 2 2 2
1 : 2 2 2 2 2
2 : 3 3 3 3 3
3 : 3 3 3 3 3
4 : 2 2 2 2 2
5 : 2 2 2 2 2

Variable 1
Array of 1 dimension(s)
6 6 10 10 6 6
Variable 2
Array of 1 dimension(s)
4

Basis of cmp 2 in the nucleus:
3-dimensional spectral base
Variable 0

Changing the tensorial basis

In some cases it is possible to change the tensorial basis by invoking :

• change_basis_spher_to_cart to go from the orthonormal spherical basis to the Cartesian one.
• change_basis_cart_to_spher to do the reverse.

The spectral basis are also changed in a consistent manner.

    // Change the tensorial basis
    vecU.change_basis_cart_to_spher() ;

    // Prints the new tensorial basis
    cout << vecU.get_basis() << endl ;


    // Check that the spectral basis has changed
    cout << "Spectral basis of component 3 in the nucleus: " << endl ;
    cout << vecU(3)(0).get_base() << endl ;
-------------------------------------------------------------------    
Spectral basis of component 3 in the nucleus:
3-dimensional spectral base
Variable 0
Array of 2 dimension(s)
0 : 3 3 3 3 3
1 : 3 3 3 3 3
2 : 2 2 2 2 2
3 : 2 2 2 2 2
4 : 3 3 3 3 3
5 : 3 3 3 3 3

Variable 1
Array of 1 dimension(s)
10 10 6 6 10 10
Variable 2
Array of 1 dimension(s)
4

General tensors

What is true for vectors, remains essentially true for higher order tensors. The only complication is that the valence must be passed to the constructor along with the type of all the indices.

    // Construct a rank-two tensor, both indices being contravariant
    Tensor Tone (space, 2, CON, basis_cart) ;
    cout << Tone << endl ; /// The values of the components are not defined at this point

    // A rank 3 tensor with COV, CON, CON indices
    int valence = 3 ;
    Array<int> type_indices (valence) ;
    type_indices.set(0) = COV ; 
    type_indices.set(1) = CON ;
    type_indices.set(2) = CON ;
    Tensor Ttwo (space, valence, type_indices, basis_cart) ;

Accessors

For low rank tensors (i.e. below a valence 4), accessors work by giving the indices one after the other.

For higher ranks, one needs to use the class Index (see tutorial 1). Beware that in this case, the indices should go from 0 to dim-1 (as one is using the Array description with Index).

    // Direct accessor
    Ttwo.set(1,2,2) = 1. ; // Sets the component (x,y,y) to one

    // Accessor using an Index
    Index indices (Ttwo) ;
    indices.set(0) = 2 ; indices.set(1) = 2 ; indices.set(2) = 0 ; // Corresponds to (z,z,x)
    Ttwo.set(indices) = 2. ;

    cout << Ttwo << endl ; // Only two components are defined
-------------------------------------------------------------------        
Class : N6Kadath6TensorE           Valence : 3
Tensorial basis
Domain 0 : Cartesian
Domain 1 : Cartesian
Domain 2 : Cartesian
Domain 3 : Cartesian

Type of the indices : index 0 :  covariant.
                      index 1 :  contravariant.
                      index 2 :  contravariant.

================ Component  1 1 1 ================

Scalar
Domain : 0

Domain : 1

Domain : 2

Domain : 3


================ Component  1 1 2 ================

Scalar
Domain : 0

Domain : 1

Domain : 2

Domain : 3


================ Component  1 1 3 ================

Scalar
Domain : 0

Domain : 1

Domain : 2

Domain : 3


================ Component  1 2 1 ================

Scalar
Domain : 0

Domain : 1

Domain : 2

Domain : 3


================ Component  1 2 2 ================

Scalar
Domain : 0
Configuration space :
Array of 3 dimension(s)
      0
0 : 1 1 1 1 1 1 1 1 1 1 1 1 1
1 : 1 1 1 1 1 1 1 1 1 1 1 1 1
2 : 1 1 1 1 1 1 1 1 1 1 1 1 1
3 : 1 1 1 1 1 1 1 1 1 1 1 1 1
4 : 1 1 1 1 1 1 1 1 1 1 1 1 1

      1
0 : 1 1 1 1 1 1 1 1 1 1 1 1 1
1 : 1 1 1 1 1 1 1 1 1 1 1 1 1
2 : 1 1 1 1 1 1 1 1 1 1 1 1 1
3 : 1 1 1 1 1 1 1 1 1 1 1 1 1
4 : 1 1 1 1 1 1 1 1 1 1 1 1 1

      2
0 : 1 1 1 1 1 1 1 1 1 1 1 1 1
1 : 1 1 1 1 1 1 1 1 1 1 1 1 1
2 : 1 1 1 1 1 1 1 1 1 1 1 1 1
3 : 1 1 1 1 1 1 1 1 1 1 1 1 1
4 : 1 1 1 1 1 1 1 1 1 1 1 1 1

      3
0 : 1 1 1 1 1 1 1 1 1 1 1 1 1
1 : 1 1 1 1 1 1 1 1 1 1 1 1 1
2 : 1 1 1 1 1 1 1 1 1 1 1 1 1
3 : 1 1 1 1 1 1 1 1 1 1 1 1 1
4 : 1 1 1 1 1 1 1 1 1 1 1 1 1



Domain : 1
Configuration space :
Array of 3 dimension(s)
      0
0 : 1 1 1 1 1 1 1 1 1 1 1 1 1
1 : 1 1 1 1 1 1 1 1 1 1 1 1 1
2 : 1 1 1 1 1 1 1 1 1 1 1 1 1
3 : 1 1 1 1 1 1 1 1 1 1 1 1 1
4 : 1 1 1 1 1 1 1 1 1 1 1 1 1

8. Tensors in systems

Tensors as variable or constants

When passing tensors as unknown or constant fields, there is no need to pass the indices. The library will deduce the proper valence and type of the indices from the tensor.

// The system of equation in domain 1 only
System_of_eqs syst (space, 1) ;
// Pass V as an unknown field
syst.add_var ("U", vecU) ;
// Pass T as a constant
syst.add_cst ("T", tenT) ;

Using tensors in definitions

When using tensors in definitions, one must provide a name for the indices. The covariant indices are denoted by "_" and the contravariant ones by "^". If a type is inconsistent with the definition of the tensorial field, the metric is used to manipulate the index, if it has been defined (see Tutorial 9)

Standard operations on tensors are available, like contraction, tensorial product and so on. When an indices name is repeated, Einstein summation is used (contraction).

// Tensorial product 
syst.add_def ("TV_ij^k = T_ij * U^k") ;

// Tensorial product with contraction
syst.add_def ("Contract_i = T_ij * U^j") ;

// Transposition 
syst.add_def ("Transpose_ij = T_ji") ;

// Index manipulation produces an error as no metric has been defined
//syst.add_def ("Ucov_i = U_i") ;  

// The various definitions can be accessed by give_val_def
// For instance for printing
// cout << syst.give_val_def("Contract") << endl ;

Solving equations with tensors

Tensors are used in equations in the same way as in definitions.

Important point For now, it is not possible to use an orthonormal spherical basis in the nucleus (due to difficulties with the regularity conditions at the origin). Otherwise, all other cases should work (Cartesian basis everywhere and spherical one in domains avoiding the origin).

    // One will solve only in domain 1 using syst
    // Define a tensor for the inner BC (using twice the initial value of vecU)
    Vector vecInner (2*vecU) ;
    syst.add_cst ("I", vecInner) ;

    //Random inner BC
    syst.add_eq_bc (1, INNER_BC, "U^k= I^k") ; 
    // Bulk equation
    syst.add_eq_inside (1, "lap(U^k)= 0") ; 
    // Outer BC
    syst.add_eq_bc (1, OUTER_BC, "U^k= 0") ; 

    // Newton-Raphson solver
    double prec = 1e-8 ;
    double error ;
    bool endloop = false ;
    while (!endloop) {
           endloop = syst.do_newton(prec, error) ;
    }
    // The problem is linear and converges in one step.
    // One can print the solution
    for (int cmp=1 ; cmp<=3 ; cmp++) {
             cout << "Component " << cmp << endl ;
            cout << vecU(cmp)(1) << endl ;
    }
    // At this point vecU is no longer vecInner/2 (the first one has been modified by the solver).

9. Flat Metrics

Flat metric are constructed for the Space and the tensorial basis of decomposition only.

Those objects are designed to be used in System_of_eqs so there are not a lot of things that one can do at this point.

// Assume a spherical space has been defined (put 9 points in theta and 8 in phi to be safe)
// A Cartesian tensorial basis
Base_tensor basis_cart (space, CARTESIAN_BASIS) ;
Metric_flat fmet_cart (space, basis_cart) ; // The associated flat metric

// Same thing in orthonormal spherical basis
Base_tensor basis_spher (space, SPHERICAL_BASIS) ;
Metric_flat fmet_spher (space, basis_spher) ; // The associated flat metric

Associating the metric to the system

The link between a System_of_eqs and a Metric is done by invoking a function of the metric : set_system. Once this is done, it gives access to a lot of additional functions in the system like the manipulation of the indices of tensors and the covariant derivative.

    // The system only in the first shell
    System_of_eqs syst (space, 1) ;

    // Use the orthonormal spherical basis
    fmet_spher.set_system (syst, "f") ;
    // The string denotes the names by which the metric will be recognized by the system

    // Recover and print the inverse of the metric
    syst.add_def ("inv^ij = f^ij") ;
    //cout << syst.give_val_def("inv") << endl ;

    // Assume a contravariant vector vecU has been defined
    // Use it as a constant
    syst.add_cst ("U", vecU) ;
    // One can manipulate its indices :
    syst.add_def ("Ucov_i = U_i") ;
    // One can compute the scalar product (should be one)
    syst.add_def ("product = U_i * U^i") ;
    //cout << syst.give_val_def("product") << endl ;
    
    // The covariant derivative
    syst.add_def ("der_i^j = D_i U^j") ;
    cout << syst.give_val_def("der") << endl ;
-------------------------------------------------------------------    
Class : N6Kadath6TensorE           Valence : 2
Tensorial basis
Domain 0 : Undefined
Domain 1 : Spherical
Domain 2 : Undefined
Domain 3 : Undefined
Type of the indices : index 0 :  covariant.
                      index 1 :  contravariant.
================ Component  1 1 ================
Scalar
Domain : 0
Null Val_domain
Domain : 1
Coefficient space :
Array of 3 dimension(s)
      0
0 : 0 0 0 0 0 0 0 0 0 0 0 0 0
1 : 0 0 0 0 0 0 0 0 0 0 0 0 0
2 : 0 0 0 0 0 0 0 0 0 0 0 0 0
3 : 0 0 0 0 0 0 0 0 0 0 0 0 0
4 : 0 0 0 0 0 0 0 0 0 0 0 0 0
5 : 0 0 0 0 0 0 0 0 0 0 0 0 0
6 : 0 0 0 0 0 0 0 0 0 0 0 0 0
7 : 0 0 0 0 0 0 0 0 0 0 0 0 0
8 : 0 0 0 0 0 0 0 0 0 0 0 0 0
      1
      9
0 : 0 0 0 0 0 0 0 0 0 0 0 0 0
1 : 0 0 0 0 0 0 0 0 0 0 0 0 0
2 : 0 0 0 0 0 0 0 0 0 0 0 0 0
3 : 0 0 0 0 0 0 0 0 0 0 0 0 0
4 : 0 0 0 0 0 0 0 0 0 0 0 0 0
5 : 0 0 0 0 0 0 0 0 0 0 0 0 0
6 : 0 0 0 0 0 0 0 0 0 0 0 0 0
7 : 0 0 0 0 0 0 0 0 0 0 0 0 0
8 : 0 0 0 0 0 0 0 0 0 0 0 0 0

Domain : 2
Null Val_domain
Domain : 3
Null Val_domain
================ Component  1 2 ================
Scalar
Domain : 0
Null Val_domain

Domain : 1
Coefficient space :
Array of 3 dimension(s)
      0
0 : 0 0 0 0 0 0 0 0 0 0 0 0 0
1 : 0 0 0 0 0 0 0 0 0 0 0 0 0
2 : 0 0 0 0 0 0 0 0 0 0 0 0 0
3 : 0 0 0 0 0 0 0 0 0 0 0 0 0
4 : 0 0 0 0 0 0 0 0 0 0 0 0 0
5 : 0 0 0 0 0 0 0 0 0 0 0 0 0
6 : 0 0 0 0 0 0 0 0 0 0 0 0 0
7 : 0 0 0 0 0 0 0 0 0 0 0 0 0
8 : 0 0 0 0 0 0 0 0 0 0 0 0 0
      1
      9
0 : 0 0 0 0 0 0 0 0 0 0 0 0 0
1 : 0 0 0 0 0 0 0 0 0 0 0 0 0
2 : 0 0 0 0 0 0 0 0 0 0 0 0 0
3 : 0 0 0 0 0 0 0 0 0 0 0 0 0
4 : 0 0 0 0 0 0 0 0 0 0 0 0 0
5 : 0 0 0 0 0 0 0 0 0 0 0 0 0
6 : 0 0 0 0 0 0 0 0 0 0 0 0 0
7 : 0 0 0 0 0 0 0 0 0 0 0 0 0
8 : 0 0 0 0 0 0 0 0 0 0 0 0 0

Domain : 2
Null Val_domain
Domain : 3
Null Val_domain
================ Component  1 3 ================
Scalar
Domain : 0
Null Val_domain
Domain : 1
Coefficient space :
Array of 3 dimension(s)
      0
0 : 0 0 0 0 0 0 0 0 0 0 0 0 0
1 : 0 0 0 0 0 0 0 0 0 0 0 0 0
2 : 0 0 0 0 0 0 0 0 0 0 0 0 0
3 : 0 0 0 0 0 0 0 0 0 0 0 0 0
4 : 0 0 0 0 0 0 0 0 0 0 0 0 0
5 : 0 0 0 0 0 0 0 0 0 0 0 0 0
6 : 0 0 0 0 0 0 0 0 0 0 0 0 0
7 : 0 0 0 0 0 0 0 0 0 0 0 0 0
8 : 0 0 0 0 0 0 0 0 0 0 0 0 0

      1
      9
0 : 0 0 0 0 0 0 0 0 0 0 0 0 0
1 : 0 0 0 0 0 0 0 0 0 0 0 0 0
2 : 0 0 0 0 0 0 0 0 0 0 0 0 0
3 : 0 0 0 0 0 0 0 0 0 0 0 0 0
4 : 0 0 0 0 0 0 0 0 0 0 0 0 0
5 : 0 0 0 0 0 0 0 0 0 0 0 0 0
6 : 0 0 0 0 0 0 0 0 0 0 0 0 0
7 : 0 0 0 0 0 0 0 0 0 0 0 0 0
8 : 0 0 0 0 0 0 0 0 0 0 0 0 0

Domain : 2
Null Val_domain
Domain : 3
Null Val_domain
================ Component  2 1 ================
Scalar
Domain : 0
Null Val_domain
Domain : 1
Configuration space :
Array of 3 dimension(s)
      0
0 : -3.09252e-18 -3.04071e-18 -2.89836e-18 -2.69748e-18 -2.47401e-18 -2.25634e-18 -2.06168e-18 -1.89794e-18 -1.76715e-18 -1.66843e-18 -1.59984e-18 -1.55954e-18 -1.54626e-18
1 : -9.0852e-16 -8.93301e-16 -8.51482e-16 -7.92466e-16 -7.26816e-16 -6.62868e-16 -6.0568e-16 -5.57576e-16 -5.19154e-16 -4.90151e-16 -4.70002e-16 -4.58163e-16 -4.5426e-16
2 : 1.62754e-15 1.60028e-15 1.52536e-15 1.41964e-15 1.30203e-15 1.18747e-15 1.08503e-15 9.98853e-16 9.30023e-16 8.78065e-16 8.41971e-16 8.20762e-16 8.1377e-16
3 : -2.23128e-15 -2.1939e-15 -2.0912e-15 -1.94626e-15 -1.78502e-15 -1.62797e-15 -1.48752e-15 -1.36938e-15 -1.27502e-15 -1.20378e-15 -1.1543e-15 -1.12522e-15 -1.11564e-15
4 : 1.84389e-15 1.813e-15 1.72812e-15 1.60835e-15 1.47511e-15 1.34532e-15 1.22926e-15 1.13163e-15 1.05365e-15 9.94784e-16 9.53892e-16 9.29864e-16 9.21943e-16
5 : -1.14004e-15 -1.12094e-15 -1.06847e-15 -9.94414e-16 -9.12034e-16 -8.31789e-16 -7.60028e-16 -6.99666e-16 -6.51453e-16 -6.15058e-16 -5.89775e-16 -5.74919e-16 -5.70021e-16
6 : 1.05748e-15 1.03977e-15 9.91091e-16 9.22399e-16 8.45985e-16 7.71552e-16 7.04988e-16 6.48997e-16 6.04275e-16 5.70516e-16 5.47064e-16 5.33283e-16 5.28741e-16
7 : -1.05038e-15 -1.03278e-15 -9.84433e-16 -9.16203e-16 -8.40302e-16 -7.66369e-16 -7.00252e-16 -6.44637e-16 -6.00216e-16 -5.66683e-16 -5.43389e-16 -5.29701e-16 -5.25189e-16
8 : 1.84889e-32 1.81792e-32 1.73282e-32 1.61272e-32 1.47911e-32 1.34898e-32 1.2326e-32 1.1347e-32 1.05651e-32 9.97486e-33 9.56482e-33 9.32389e-33 9.24446e-33
      1

10. General Metrics

The Metric_tensor class

In order to deal with metrics, one must use a special type of tensor Metric_tensor. It is a rank-2 symmetric tensor. Both indices are of the same type but can be covariant or contravariant, so that this class can describe the metric or its inverse.

// Assume a spherical space has been defined (put 9 points in theta and 8 in phi to be safe)
// The tensorial basis
Base_tensor basis (space, CARTESIAN_BASIS) ; 
Metric_tensor gmet (space, COV, basis) ;

// Affectation to flat metric everywhere
for (int i=1 ; i<=3 ; i++)  {
     gmet.set(i,i) = 1. ;
    for (int j=i+1 ; j<=3 ; j++)
         gmet.set(i,j) = 0 ;
}
// Non trivial stuff in domain 1
for (int i=1 ; i<=3 ; i++)
     for (int j=i ; j<=3 ; j++)
           gmet.set(i,j).set_domain(1) = 
              gmet(i,j)(1) + space.get_domain(1)->get_cart(i)*space.get_domain(1)->get_cart(j) ;

// Standard base in that case
gmet.std_base() ;
//cout << gmet << endl ;

Using constant metrics

Constant metrics are described by the class Metric_const. By constant, one means that it is given once and for all and that the system is not allowed to change it. In that respect it works as the constant fields defined using add_cst. (see tutorial 5). Obviously it does not mean that the metric is constant in space.

It then works in the same manner as the flat metrics of tutorial 9.

// Construction of a constant metric
// The numerical value is contained in a Metric_tensor that can be covariant or contravariant.
Metric_const met (gmet) ;

// Passing it to a system 
System_of_eqs syst (space, 1) ;
met.set_system (syst, "g") ;

// One can get the connection coefficients (reserved word Gam)
syst.add_def ("Christo_ij^k = Gam_ij^k") ;
//cout << syst.give_val_def("Christo") << endl ;

// One can also get the Ricci tensor (reserved word R)
syst.add_def ("Ricci_ij = R_ij") ;
//cout << syst.give_val_def("Ricci") << endl ;

Metrics as unknown

Often the metric is given but is part of the physical problem, meaning it is an unknown field that the solver is trying to find. When this is the case, one must use the class Metric_general. It works as the Metric_const except that the metric is now treated as a variable field (see tutorial 5).

Solving the system of equations, will change the Metric_tensor used to define the metric.

// Construction of the fields for inner and outer BC
Metric_tensor ginner (gmet) ; // Contains the initial value of gmet
// No need to call std_base 
Metric_tensor gouter (space, COV, basis) ;
for (int i=1 ; i<=3 ; i++) {
    gouter.set(i,i) = 1. ;
    for (int j=i+1 ; j<=3 ; j++)
        gouter.set(i,j) = 0 ;
}
gouter.std_base() ;

// Construction of a unknown metric
// The initial value numerical value is contained in a Metric_tensor.
Metric_general metgen (gmet) ;

// Passing it to a new system 
System_of_eqs systtwo  (space, 1) ;
metgen.set_system (systtwo, "g") ;

// Pass the fields ginner and gouter as constants
systtwo.add_cst ("gin", ginner) ;
systtwo.add_cst ("gout", gouter) ;

// Random equations
systtwo.add_eq_bc (1, INNER_BC, "g_ij = gin_ij") ;
systtwo.add_eq_inside (1, "lap(g_ij) = 0") ;
systtwo.add_eq_bc (1, OUTER_BC, "g_ij = gout_ij") ;

// Newton-Raphson solver
double prec = 1e-8 ;
double error ;
bool endloop = false ;
while (!endloop) {
    endloop = systtwo.do_newton(prec, error) ;
}

// The solution is contained in gmet
// One can check that it is different from ginner that contains the initial value of gmet
for (int i=1 ; i<=3 ; i++)
    for (int j=i ; j<=3 ; j++)
        cout << "Diff in comp " << i << ", " << j << " : " << diffmax(gmet(i,j)(1), ginner(i,j)(1)) << endl ;
-------------------------------------------------------------------    
Diff in comp 1, 1 : 4
Diff in comp 1, 2 : 2
Diff in comp 1, 3 : 2
Diff in comp 2, 2 : 4
Diff in comp 2, 3 : 2
Diff in comp 3, 3 : 4

11. Components in systems

Position of the problem

Sometime ones needs to use different equation for different components of the fields. This can be done by specifying the components to be used, when passing the concerned equation to the System_of_eqs.

It also happens that the library lost track of some symmetries. For instance, assume the unknown is a metric. It is a symmetric tensor, resulting into a set of 6 unknown components. However, when the equations are complicated, Kadath can lose track of the symmetry and find itself considering a general tensor (hence 9 components). This would result in an inconsistency. In order to cope with that, one can force the code to consider only the "right" components.

// Assume a spherical space has been defined (put 9 points in theta and 8 in phi to be safe)
// Assume a Metric_tensor gmet has been defined
// Assume two Metric_tensor ginner and gouter used for the BC (as in tutorial 10).

// General metric
Metric_general met (gmet) ;

// The system
System_of_eqs syst (space, 1) ;
met.set_system (syst, "g") ;

// One can print the number full number of unknowns
// This corresponds roughly to the number of coefficients of 6 components.
cout << "Number of unknowns : " << syst.get_nbr_unknowns() << endl ;

// Pass the constants
syst.add_cst ("gin", ginner) ;
syst.add_cst ("gout", gouter) ;

// Now the equations
// Let us use a "complicated" but symmetric bulk equation
syst.add_eq_bc (1, INNER_BC, "g_ij = gin_ij") ;
syst.add_eq_inside (1, "lap(g_ij) = g_ik * R_j^k + g_jk * R_i^k") ;
syst.add_eq_bc (1, OUTER_BC, "g_ij = gout_ij") ;

// Compute the initial errors
// This ensure that the code properly determines the number of conditions
Array<double> errors (syst.sec_member()) ;

// One can print the number of conditions araising from the equations
cout << "Wrong number of conditions : " << syst.get_nbr_conditions() << endl ;

Number of unknowns : 4706
Wrong number of conditions : 6609
Right number of conditions : 4706

What happens is that the code has lost track that the bulk equation was symmetric and so it is considered as a 9-components tensor, thus leading to too many conditions.

Using List_comp

In order to solve the issues raised above, one needs to specify, when defining a tensorial equation, which components to consider. This is done be explicitly describing the list of components by means of the class List_comp.

// Constructor with two arguments.
// First : the number of components to be considered.
// Second : the rank (i.e. the number of indices)
List_comp used (6, 2) ; // Here 6 components of a rank-2 tensor.

// Passing of the components by hand
used.set(0)->set(0) = 1 ; used.set(0)->set(1) = 1 ; // Component xx
used.set(1)->set(0) = 1 ; used.set(1)->set(1) = 2 ; // Component xy
used.set(2)->set(0) = 1 ; used.set(2)->set(1) = 3 ; // Component xz
used.set(3)->set(0) = 2 ; used.set(3)->set(1) = 2 ; // Component yy
used.set(4)->set(0) = 2 ; used.set(4)->set(1) = 3 ; // Component yz
used.set(5)->set(0) = 3 ; used.set(5)->set(1) = 3 ; // Component zz

Now when defining the system, one can pas this List_comp to the bulk equation that will now take into account only the right components. This will lead to a number of conditions consistent with the number of unknowns.

// Define a new system, systtwo, like the first one and just use
systtwo.add_eq_inside (1, "lap(g_ij) = g_ik * R_j^k + g_jk * R_i^k", used) ;

12. Using Files

Writing files

The file must be open via the fopen function.

For integers and doubles one needs to use the function fwrite_be, that uses the big endian convention.

Most of Kadath classes have a function save that must be invoked. Please not that the various fields need to be constructed not only from the file but also from the Space so that this one needs to be saved first.

// Opening the file in write mode.
FILE* fout = fopen ("file.dat", "w") ;

// Two values to be saved
int valint = 2 ;
double valdouble = 1.31732 ;

// Write using fwrite_be
fwrite_be (&valint, sizeof(int), 1, fout) ;
fwrite_be (&valdouble, sizeof(double), 1, fout) ;

// Assume a space, a scalar and a vector have been previously defined
// They are saved by invoking save
space.save(fout) ;
field.save(fout) ;
vecU.save(fout) ;

// Dont forget to close the file in the end
fclose(fout) ;

Reading files

The objects in the files must obviously be read in the same order as they have been saved.

For integers and doubles, this in done using fread_be.

The various classes have constructors from a file. As said above, the constructors also needs a Space so that this one must be saved first.

// Opening the file in read mode.
FILE* fin = fopen ("file.dat", "r") ;

// Two values to be read
int valint  ;
double valdouble  ;

// Read using fread_be
fread_be (&valint, sizeof(int), 1, fin) ;
fread_be (&valdouble, sizeof(double), 1, fin) ;

cout << "Integer read " << valint << endl ;
cout << "Double read " << valdouble << endl ;

// Reading Kadath classes
// Space first, its exact type must be known
Space_spheric space (fin) ;
cout << space << endl ;

// Then the fields
Scalar field (space, fin) ;
Vector vecU (space, fin) ;

// Dont forget to close the file in the end
fclose(fin) ;

13. Schwarz solution test

#include "kadath_spheric.hpp"

using namespace Kadath ;
int main() {

    // 3D :
    int dim = 3 ;

    // Number of points
    int nbr  = 13 ;
    Dim_array res (dim) ;
    res.set(0) = nbr ; res.set(1) = 5 ; res.set(2) = 4 ;

    // Center of the coordinates
    Point center (dim) ;
    for (int i=1 ; i<=dim ; i++)
         center.set(i) = 0 ;

    // Number of domains and boundaries :
    int ndom = 4 ;
    Array<double> bounds (ndom-1) ;
    // Radius of the BH !
    double aa = 1.323 ;
    bounds.set(0) = aa ; bounds.set(1) = 1.7557*aa ; bounds.set(2) = 2.9861*aa ;

    // Chebyshev or Legendre :
    int type_coloc = CHEB_TYPE ;

    // Sherical space :
    Space_spheric space(type_coloc, center, res, bounds) ;

    // Initial guess for the conformal factor :
    Scalar conf (space) ;    
    conf = 1. ;
    conf.std_base() ;

    // Solve the equation in space outside the nucleus
    System_of_eqs syst (space, 1, ndom-1) ;
    // Only one unknown
    syst.add_var ("P", conf) ;
    // One user defined constant
    syst.add_cst ("a", aa) ;    

    // Inner BC
    space.add_inner_bc (syst, "dn(P) + 0.5 / a * P = 0") ;
    // Equation
    space.add_eq (syst, "Lap(P) = 0", "P", "dn(P)") ;
    // Outer BC
    space.add_outer_bc (syst, "P=1") ;

    // Newton-Raphson
    double conv ;
    bool endloop = false ;
    int ite = 1 ;
    while (!endloop) {
        endloop = syst.do_newton(1e-8, conv) ;
        ite++ ;
    }

    // Check of the solution
    int resol = 100 ;
    double xxmin = bounds(0)*1.01 ;
    double xxmax = bounds(2)*5 ;
    double step = (xxmax-xxmin)/resol ;
    double xx = xxmin+step ;
    double error_max = 0 ;

    double tet = M_PI/2. ;
    double phi = -2.349 ;

    double xunit = sin(tet)*cos(phi) ;
    double yunit = sin(tet)*sin(phi) ;
    double zunit = cos(tet) ;

    Point M (3) ;
    for (int i=0 ; i<resol-1 ; i++) {

        M.set(1)=xx*xunit ;
        M.set(2)=xx*yunit ;
        M.set(3)=xx*zunit ;

        double ana = 1. + aa/xx ;
        double error = fabs (ana - conf.val_point(M)) ;
        if (error > error_max)
            error_max = error ;
        xx+=step ;
        }

    cout << "Error max " << error_max << endl ;
    return EXIT_SUCCESS ;
}

Entering do_newton with error 0.377929
Entering do_newton with error 5.17635e-13
Error max 2.95179e-10