This tutorial shows how to use tensors in a System_of_eqs
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.
// Assume a spherical space has been defined (put 9 points in theta and 8 in phi to be safe)
// Assume one has two previously defined tensors vecU (vector) and tenT (rank-2, covariant one)
// The system of equations in domain 1 only
System_of_eqs syst (space, 1) ;
// Pass U as an unknown field
syst.add_var ("U", vecU) ;
// Pass T as a constant
syst.add_cst ("T", tenT) ;
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 index 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 ;
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).