4.2.2. Internals: What is happening in pyoomph
When the problem is initialized - which either happens by an explicit call of initialise() or implicitly by the first call of e.g. run(), solve(), output(), etc., the define_problem() method is invoked. Depending on the present problem, meshes and equations are added to the problem.
The equations are assembled in an equation tree. This tree has the bulk domains as base nodes and each bulk domain can have further sub-nodes by equations defined on the boundary interfaces of this domain. If further equations are added to the boundaries of interfaces, the tree can go deeper. A text file containing the equation tree can be found relative to the output directory in _ccode/_equation_tree.txt . For each entry of the tree, all equations defined on this domain/interfaces are merged and a corresponding C code is generated. This C code performs the assembly of the discretized residual vector and its Jacobian on the basis of a single element of the domain mesh. It furthermore instructs the core of pyoomph to allocate all fields defined in the define_fields() methods of all equations defined on this particular domain. Furthermore, initial conditions and the values of the Dirichlet conditions are incorporated in the C code as well.
For the C code generation, all fields occurring in the weak expressions of the residuals added to this domain are expanded by sums over the shape functions of the selected basis functions. Likewise, the corresponding test functions are expanded to sums of the same basis functions. Spatial derivatives are carried out on the basis functions within these expansions (considering the selected coordinate system) and potential temporal derivative schemes are applied on the nodal values. During this expansion, also physical dimensions are plugged in (if used), i.e. the residual form is nondimensionalized based on the set scales.
Whenever the user is solving the system, for each element in each domain, the C code is called to assemble the residual (and the Jacobian) based on the very same steps as discussed in the previous section, but of course more general, i.e. allowing also for non-linear systems. Instead of the stiffness matrix \(\mathbf{K}\), the Jacobian will enter the system (in case of the Poisson problem, both matrices are the same). The Jacobian \(\mathbf{J}\) can be constructed from the residual vector (entries \(R_i\)) by the entries \(J_{ij}=\partial_{u_j} R_i\). For nonlinear problems, the assembly of the residual and the Jacobian must be repeated iteratively (Newton’s method), since \(J_{ij}\) may change with the vector of unknowns \(\vec{u}\). Therefore, \(\vec{u}\) is updated according to
until the residual vector vanishes (i.e. its maximum entry by modulus falls below a specific threshold, controlled by the newton_solver_tolerance property of the Problem class).
To assemble the Jacobian and the residual vector, the contribution of each element is added to the global matrix/vector. To that end, each element has a map of the local degrees of freedom to the global equation numbers. The generated C code just gets passed the information of the values of each field stored at the nodes of the current element and the equation mapping from the pyoomph core and calculates the contribution to the Jacobian and residuals for a single element. The core iterates over all elements and accumulates the individual contributions of all elements, which usually overlap due to the shared nodes between different elements. This gives finally the global Jacobian and residual vector. The solution of the matrix equation (4.10) is done by a linear solver back-end, which can be selected by the set_linear_solver() method of the Problem.