This tutorial introduces the PDE solver of Kadath
System_of_eqs is one of the main objects 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).
// Assume a space is known (see tutorial 1)
// Get the number of domains, if needed
// int ndom = space.get_nbr_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)
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).
// 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 referred to by the system.
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
// define a constant, being 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 referred to by the system.
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 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
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)") ;
}
At this point one 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 ;
// In this case only the inner boundary condition is not verified by our initial guess
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.
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 ;