|
The Open FUSION Toolkit 1.0.0-8905cc5
Modeling tools for plasma and fusion research and engineering
|
This page addresses general considerations with the implementation of finite element methods in the Open FUSION Toolkit (OFT). In particular basis set definition, field creation, and operator construction will be addressed. OFT currently support C0 Galerkin finite element methods with volume based interaction only, where the discretized system matrix has the form
\[ \int_{\Omega} g( u^T ) f( u ) dV + \int_{\partial \Omega} u^T h( u ) dS \]
where \( u \) and \( u^T \) are the basis and test functions respectively, which for Galerkin weighting are the same set of functions. Matrix elements exist for this system when the row and column correspond elements that share a common cell. This limitation is due to the linkage and seam construction methods, which currently do no support connectivity using internal faces.
This document explains the basic methods to define and implement a finite element representation using a scalar Lagrange representation as an example. For information about a specific finite element representation please view its documentation. These documents may be accessed through the finite element group.
A finite element representation is defined by the basis set on the unit tetrahedron. Each basis function is tied to a geometric entity (vertex, edge, etc.) according the logical coordinates it depends on. In the case of a quadratic Lagrange representation scalars are represented by a single basis function associated with each vertex and edge. These functions are
\( f^v_{i} (u_i) = 2 u_i ( u_i - 1/2 ) \)
\( f^e_{i} (u^1_i,u^2_i) = 4 u^1_i u^2_i \)
where \( u_i \) is the logical coordinate associated with the i-th vertex and \( u^1_i, u^2_i \) are the coordinates associated with the end points of the i-th edge.
In OFT a finite element representation is defined by the oft_fem_type class, which contains the structures required to define geometric connectivity, vector definition, and the underlying graph for discretized operators. Most of the setup and construction of this structure is automated and relies only on the mesh and association of the elements to geometric primitives.
Shown below is an example of the finite element setup procedure for the quadratic Lagrange set presented above. First the mesh is linked to the structure and the order of the representation is set, quadratic in this case. The number of basis functions per geometric primitive is then set in the gstruct field. Where gstruct is a 4 value integer array containing the number of basis functions associated with each vertex, edge, face and cell respectively. Finally, the remaining setup can be completed using the setup method, where the argument is the desired quadrature order.
Although the basic structure is setup automatically subroutines must be provided to evaluate the basis functions and their derivatives if desired. For Lagrange elements the basis functions do not require the Jacobian matrix for the spatial mapping, so instead the gradient evaluation routine will be presented to illustrate all features. The ordinary basis functions can be evaluated using the same structure, without the Jacobian terms.
Below is annotated source code for the oft_lag_geval subroutine, which evaluates the gradient in physical coordinates of a given basis function. The gradient is dependent both on the derivatives of the basis function in logical coordinates as well as the Jacobian of the logical to physical mapping for the current cell. Where physical derivatives are evaluated using the chain rule as
\( \frac{\partial f}{\partial x_i} = \frac{\partial f}{\partial u_j} \frac{\partial u_j}{\partial x_i} \)
Spatial derivatives of the logical coordinates are computed using the subroutine tetmesh_get_jacobian and passed to the routines through the gop parameter. In order to evaluate the function the position in the mesh must be specified, by the cell index cell and logical coordinates within the cell f, along with the desired basis function within the cell to evaluate dof. Basis functions are numbered by the fem_setup routine according to their geometric linkage. The resulting gradient is then returned using the val parameter. For the existing finite element representations element specific evaluation routines are called from this driver routine depending on the element linkage, specified in the cmap structure.
If the requested element is a vertex element no additional handling is needed and the evaluation routine may be called.
If the requested element is a edge element the edge must first be orient to be globally consistent. This is performed by first retrieving the orientation flag for the local edge lce and the end point vertices in the local orientation tet_ed. The end points can then be oriented to the global element orientation using orient_list2 and the evaluation routine called.
If the requested element is a face element the face must also be orient to be globally consistent. This is performed by first retrieving the orientation flag for the local face lcfo and the corner vertices in the local orientation tet_fc. The corner points can then be oriented to the global element orientation using orient_listn and the evaluation routine called.
If the requested element is a cell element no additional handling is needed and the evaluation routine may be called.
Each of the evaluation routines returns and array of values corresponding to \( \frac{\partial f}{\partial u_j} \). In order to reconstruct the gradient in physical coordinates these coefficients must be combined with the jacobian terms.
Once the finite element structure has been constructed it can be used to create a vector representation of the basis function weights. These vectors can then be used in the solution of linear systems formed using this finite element space. Vectors can be created from the finite element structure by using the vec_create method.
The oft_fem_type class also provides functionality for saving and reading FE weight vectors in binary restart files, using the vec_save and vec_load methods. Restart files utilize the HDF5 library and perform I/O in parallel using a common file for all processes.
The supporting structure presented so far is primarily used to support the construction of discretized operators defined by a finite element projection of the weak form of the desired PDE. The example below outlines the basics of this operation for the Laplacian operator, whose weak form is
\[ \nabla^2 u = - \int_{\Omega} \nabla \phi \cdot \nabla u dV + \int_{\Gamma} \phi \nabla u \cdot dS \]
For this example we will show the operator construction of the volume term using a Galerkin method where both the solution and test functions are expanded using a Lagrange basis. A Dirichlet boundary condition for all boundaries will also be imposed for this example, by eliminating the boundary rows and replacing the corresponding diagonal entry with 1.
\[ \phi = \sum \phi_i f_i \]
\[ u = \sum u_j f_j \]
\[ f \in F_{lagrange} \]
Before we can compute the entries of the matrix corresponding to this operator the matrix representation must be set up and space allocated for the all entries. This process is performed automatically for a standard FE representation by the mat_create method.
The mat_create method returns a bare matrix representation for an operator with full interaction between all FE weights that share a common cell. It contains all of the required mapping and structural information to support computing matrix entries and performing matrix-vector operations.
assemble called at least once before matrix vector products may be computed.Matrix entries for the desired operator correpsond to the result of the above volume integral for each basis/test function pair
\[ A_{ij} = \int \nabla f_i \cdot \nabla f_j dV \]
where \( f_i \) and \( f_j \) are independent Lagrange basis functions. As each function only has a non-zero value in cells which contain that particular element this integral simplfies to sum of contributions from each individual cell. Integration with in each cell is then carried out using a numerical integration with a given quadrature stencil. The desired order of the quadrature scheme for this case was set during finite element setup, so the appropriate quadrature object is already bound to the quadratic_lagrange structure.
The first section of the operator integration loop sets up local variables which are unique in each cell. These variables are declared private in the OpenMP clause and must be allocated after parallel execution has begun.
With the local environment setup the main integration begins, looping over all cells in the parent mesh. For each cell in the mesh a test is performed to see if any curved geometry exists for that cell. If there are no curved entities then the logical to physical Jacobian is constant over the entire cell and only needs to be computed once, otherwise it must be computed at each integration point. The list of elements in the current cell is the retrieved to allowing mapping of the local matrix contributions into the full matrix.
The integrand components are then computed at each quadrature point. The physical Jacobian is first computed if necessary at the current quadrature point. The weight value is then set by computing the product of the quadrature weight and the volume element Jacobian. Finally, the basis function gradients \( \left( \nabla f \right) \) are computed for each of the elements in the current cell.
The local matrix entries are now given as
\[ A_{ij} = \sum_k \nabla f_i (x_k) \cdot \nabla f_j (x_k) \omega_k \left| \mathbf{J}(x_k) \right| \]
where \( x_k \) and \( \omega_k \) are the k-th quadrature point and weight respectively and \( \mathbf{J} \) is the Jacobian matrix for the logical to physical mapping. These entries can be readily computed with the existing information by looping over the rows, columns, and quadrature points. Rows corresponding to boundary elements are skipped consistent with the desired boundary condition.
Finally, the local matrix contributions are added to the full matrix using the add_values method. This method takes the local matrix and list of local variables and uses them to add the local contributions into the full matrix. Depending on the linear algebra backed end being used the operations performed by this method are different but the calling structure is the same. This call must be surrounded by an !$omp critical clause to prevent data races when added contributions into the main matrix. This clause restricts execution of the given region to a single thread at a time.
Once the loop has completed local variables are deallocated and the matrix is assembled. This assembly call is needed for linear algebra backends such as PETSc. The diagonal entries corresponding to removed boundary rows are then set to 1 in order to produce a non-singular system and allowing imposing the boundary value by modifying the RHS vector. Diagonal entries are only set for rows that correspond to boundary elements and are only set on the domain with ownership of that element, as indicated by leo. Assemble is then called again to finalize the matrix and make it available for use in linear algebra operations.
1.9.8