1 static char help[] = "Demonstrates the use of the COO interface to PETSc matrices for finite element computations\n\n"; 2 3 /* 4 The COO interface for PETSc matrices provides a convenient way to provide finite element element stiffness matrices to PETSc matrix that should work 5 well on both CPUs and GPUs. It is an alternative to using MatSetValues() 6 7 This example is intended for people who are NOT using DMPLEX or libCEED or any other higher-level infrastructure for finite elements; 8 it is only to demonstrate the concepts in a simple way for those people who are interested and for those people who are using PETSc for 9 linear algebra solvers but are managing their own finite element process. 10 11 Please do NOT use this example as a starting point to writing your own finite element code from scratch! 12 13 Each element in this example has three vertices; hence the the usage below needs to be adjusted for elements of a different number of vertices. 14 */ 15 16 #include <petscmat.h> 17 #include "ex18.h" 18 19 static PetscErrorCode CreateFEStruct(FEStruct *fe) { 20 PetscFunctionBeginUser; 21 fe->Nv = 5; 22 fe->Ne = 3; 23 PetscCall(PetscMalloc1(3 * fe->Ne, &fe->vertices)); 24 /* the three vertices associated with each element in order of element */ 25 fe->vertices[0 + 0] = 0; 26 fe->vertices[0 + 1] = 1; 27 fe->vertices[0 + 2] = 2; 28 fe->vertices[3 + 0] = 2; 29 fe->vertices[3 + 1] = 1; 30 fe->vertices[3 + 2] = 3; 31 fe->vertices[6 + 0] = 2; 32 fe->vertices[6 + 1] = 4; 33 fe->vertices[6 + 2] = 3; 34 fe->n = 5; 35 PetscFunctionReturn(0); 36 } 37 38 static PetscErrorCode DestroyFEStruct(FEStruct *fe) { 39 PetscFunctionBeginUser; 40 PetscCall(PetscFree(fe->vertices)); 41 PetscCall(PetscFree(fe->coo)); 42 PetscFunctionReturn(0); 43 } 44 45 static PetscErrorCode CreateMatrix(FEStruct *fe, Mat *A) { 46 PetscInt *oor, *ooc, cnt = 0; 47 48 PetscFunctionBeginUser; 49 PetscCall(MatCreate(PETSC_COMM_WORLD, A)); 50 PetscCall(MatSetSizes(*A, fe->n, fe->n, PETSC_DECIDE, PETSC_DECIDE)); 51 PetscCall(MatSetFromOptions(*A)); 52 53 /* determine for each entry in each element stiffness matrix the global row and colum */ 54 /* since the element is triangular with piecewise linear basis functions there are three degrees of freedom per element, one for each vertex */ 55 PetscCall(PetscMalloc2(3 * 3 * fe->Ne, &oor, 3 * 3 * fe->Ne, &ooc)); 56 for (PetscInt e = 0; e < fe->Ne; e++) { 57 for (PetscInt vi = 0; vi < 3; vi++) { 58 for (PetscInt vj = 0; vj < 3; vj++) { 59 oor[cnt] = fe->vertices[3 * e + vi]; 60 ooc[cnt++] = fe->vertices[3 * e + vj]; 61 } 62 } 63 } 64 PetscCall(MatSetPreallocationCOO(*A, 3 * 3 * fe->Ne, oor, ooc)); 65 PetscCall(PetscFree2(oor, ooc)); 66 67 /* determine the offset into the COO value array the offset of each element stiffness; there are 9 = 3*3 entries for each element stiffness */ 68 /* for lists of elements with different numbers of degrees of freedom assocated with each element the offsets will not be uniform */ 69 PetscCall(PetscMalloc1(fe->Ne, &fe->coo)); 70 fe->coo[0] = 0; 71 for (PetscInt e = 1; e < fe->Ne; e++) { fe->coo[e] = fe->coo[e - 1] + 3 * 3; } 72 PetscFunctionReturn(0); 73 } 74 75 static PetscErrorCode FillMatrixCPU(FEStruct *fe, Mat A) { 76 PetscScalar s[9]; 77 78 PetscFunctionBeginUser; 79 /* simulation of traditional PETSc CPU based finite assembly process */ 80 for (PetscInt e = 0; e < fe->Ne; e++) { 81 for (PetscInt vi = 0; vi < 3; vi++) { 82 for (PetscInt vj = 0; vj < 3; vj++) { s[3 * vi + vj] = vi + 2 * vj; } 83 } 84 PetscCall(MatSetValues(A, 3, fe->vertices + 3 * e, 3, fe->vertices + 3 * e, s, ADD_VALUES)); 85 } 86 PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY)); 87 PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY)); 88 PetscFunctionReturn(0); 89 } 90 91 /* 92 Shows an example of tracking element offsets explicitly, which allows for 93 mixed-topology meshes and combining both volume and surface parts into the weak form. 94 */ 95 static PetscErrorCode FillMatrixCPUCOO(FEStruct *fe, Mat A) { 96 PetscScalar *v, *s; 97 98 PetscFunctionBeginUser; 99 /* simulation of CPU based finite assembly process with COO */ 100 PetscCall(PetscMalloc1(3 * 3 * fe->Ne, &v)); 101 for (PetscInt e = 0; e < fe->Ne; e++) { 102 s = v + fe->coo[e]; /* point to location in COO of current element stiffness */ 103 for (PetscInt vi = 0; vi < 3; vi++) { 104 for (PetscInt vj = 0; vj < 3; vj++) { s[3 * vi + vj] = vi + 2 * vj; } 105 } 106 } 107 PetscCall(MatSetValuesCOO(A, v, ADD_VALUES)); 108 PetscCall(PetscFree(v)); 109 PetscFunctionReturn(0); 110 } 111 112 /* 113 Uses a multi-dimensional indexing technique that works for homogeneous meshes 114 such as single-topology with volume integral only. 115 */ 116 static PetscErrorCode FillMatrixCPUCOO3d(FEStruct *fe, Mat A) { 117 PetscScalar(*s)[3][3]; 118 119 PetscFunctionBeginUser; 120 /* simulation of CPU based finite assembly process with COO */ 121 PetscCall(PetscMalloc1(fe->Ne, &s)); 122 for (PetscInt e = 0; e < fe->Ne; e++) { 123 for (PetscInt vi = 0; vi < 3; vi++) { 124 for (PetscInt vj = 0; vj < 3; vj++) { s[e][vi][vj] = vi + 2 * vj; } 125 } 126 } 127 PetscCall(MatSetValuesCOO(A, (PetscScalar *)s, INSERT_VALUES)); 128 PetscCall(PetscFree(s)); 129 PetscFunctionReturn(0); 130 } 131 132 int main(int argc, char **args) { 133 Mat A; 134 FEStruct fe; 135 PetscMPIInt size; 136 PetscBool is_kokkos, is_cuda; 137 138 PetscFunctionBeginUser; 139 PetscCall(PetscInitialize(&argc, &args, (char *)0, help)); 140 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size)); 141 PetscCheck(size <= 1, PETSC_COMM_WORLD, PETSC_ERR_WRONG_MPI_SIZE, "Demonstration is only for sequential runs"); 142 143 PetscCall(CreateFEStruct(&fe)); 144 PetscCall(CreateMatrix(&fe, &A)); 145 146 PetscCall(FillMatrixCPU(&fe, A)); 147 PetscCall(MatView(A, PETSC_VIEWER_STDOUT_WORLD)); 148 149 PetscCall(MatZeroEntries(A)); 150 PetscCall(FillMatrixCPUCOO(&fe, A)); 151 PetscCall(MatView(A, PETSC_VIEWER_STDOUT_WORLD)); 152 153 PetscCall(MatZeroEntries(A)); 154 PetscCall(FillMatrixCPUCOO3d(&fe, A)); 155 PetscCall(MatView(A, PETSC_VIEWER_STDOUT_WORLD)); 156 157 PetscCall(MatZeroEntries(A)); 158 PetscCall(PetscObjectBaseTypeCompare((PetscObject)A, MATSEQAIJKOKKOS, &is_kokkos)); 159 PetscCall(PetscObjectBaseTypeCompare((PetscObject)A, MATSEQAIJCUSPARSE, &is_cuda)); 160 #if defined(PETSC_HAVE_KOKKOS) 161 if (is_kokkos) PetscCall(FillMatrixKokkosCOO(&fe, A)); 162 #endif 163 #if defined(PETSC_HAVE_CUDA) 164 if (is_cuda) PetscCall(FillMatrixCUDACOO(&fe, A)); 165 #endif 166 PetscCall(MatView(A, PETSC_VIEWER_STDOUT_WORLD)); 167 168 PetscCall(MatDestroy(&A)); 169 PetscCall(DestroyFEStruct(&fe)); 170 PetscCall(PetscFinalize()); 171 return 0; 172 } 173 174 /*TEST 175 build: 176 requires: cuda kokkos_kernels 177 depends: ex18cu.cu ex18kok.kokkos.cxx 178 179 testset: 180 filter: grep -v "type" 181 output_file: output/ex18_1.out 182 183 test: 184 suffix: kok 185 requires: kokkos_kernels 186 args: -mat_type aijkokkos 187 188 test: 189 suffix: cuda 190 requires: cuda 191 args: -mat_type aijcusparse 192 193 TEST*/ 194