1 // Copyright (c) 2017-2022, Lawrence Livermore National Security, LLC and other CEED contributors. 2 // All Rights Reserved. See the top-level LICENSE and NOTICE files for details. 3 // 4 // SPDX-License-Identifier: BSD-2-Clause 5 // 6 // This file is part of CEED: http://github.com/ceed 7 8 // libCEED + PETSc Example: Surface Area 9 // 10 // This example demonstrates a simple usage of libCEED with PETSc to calculate 11 // the surface area of a simple closed surface, such as the one of a cube or a 12 // tensor-product discrete sphere via the mass operator. 13 // 14 // The code uses higher level communication protocols in DMPlex. 15 // 16 // Build with: 17 // 18 // make area [PETSC_DIR=</path/to/petsc>] [CEED_DIR=</path/to/libceed>] 19 // 20 // Sample runs: 21 // Sequential: 22 // 23 // ./area -problem cube -degree 3 -dm_refine 2 24 // ./area -problem sphere -degree 3 -dm_refine 2 25 // 26 // In parallel: 27 // 28 // mpiexec -n 4 ./area -problem cube -degree 3 -dm_refine 2 29 // mpiexec -n 4 ./area -problem sphere -degree 3 -dm_refine 2 30 // 31 // The above example runs use 2 levels of refinement for the mesh. 32 // Use -dm_refine k, for k levels of uniform refinement. 33 // 34 //TESTARGS -ceed {ceed_resource} -test -degree 3 -dm_refine 1 35 36 /// @file 37 /// libCEED example using the mass operator to compute a cube or a cubed-sphere surface area using PETSc with DMPlex 38 static const char help[] = "Compute surface area of a cube or a cubed-sphere using DMPlex in PETSc\n"; 39 40 #include "area.h" 41 42 #include <ceed.h> 43 #include <petsc.h> 44 #include <petscdmplex.h> 45 #include <stdbool.h> 46 #include <string.h> 47 48 #include "include/areaproblemdata.h" 49 #include "include/libceedsetup.h" 50 #include "include/matops.h" 51 #include "include/petscutils.h" 52 #include "include/petscversion.h" 53 #include "include/structs.h" 54 55 #if PETSC_VERSION_LT(3, 12, 0) 56 #ifdef PETSC_HAVE_CUDA 57 #include <petsccuda.h> 58 // Note: With PETSc prior to version 3.12.0, providing the source path to 59 // include 'cublas_v2.h' will be needed to use 'petsccuda.h'. 60 #endif 61 #endif 62 63 #ifndef M_PI 64 #define M_PI 3.14159265358979323846 65 #endif 66 67 int main(int argc, char **argv) { 68 MPI_Comm comm; 69 char filename[PETSC_MAX_PATH_LEN], ceed_resource[PETSC_MAX_PATH_LEN] = "/cpu/self"; 70 PetscInt l_size, g_size, xl_size, 71 q_extra = 1, // default number of extra quadrature points 72 num_comp_x = 3, // number of components of 3D physical coordinates 73 num_comp_u = 1, // dimension of field to which apply mass operator 74 topo_dim = 2, // topological dimension of manifold 75 degree = 3; // default degree for finite element bases 76 PetscBool read_mesh = PETSC_FALSE, test_mode = PETSC_FALSE, simplex = PETSC_FALSE; 77 Vec U, U_loc, V, V_loc; 78 DM dm; 79 OperatorApplyContext op_apply_ctx; 80 Ceed ceed; 81 CeedData ceed_data; 82 ProblemType problem_choice; 83 VecType vec_type; 84 PetscMemType mem_type; 85 86 PetscCall(PetscInitialize(&argc, &argv, NULL, help)); 87 comm = PETSC_COMM_WORLD; 88 89 // Read command line options 90 PetscOptionsBegin(comm, NULL, "CEED surface area problem with PETSc", NULL); 91 problem_choice = SPHERE; 92 PetscCall(PetscOptionsEnum("-problem", "Problem to solve", NULL, problem_types, (PetscEnum)problem_choice, (PetscEnum *)&problem_choice, NULL)); 93 PetscCall(PetscOptionsInt("-q_extra", "Number of extra quadrature points", NULL, q_extra, &q_extra, NULL)); 94 PetscCall(PetscOptionsString("-ceed", "CEED resource specifier", NULL, ceed_resource, ceed_resource, sizeof(ceed_resource), NULL)); 95 PetscCall(PetscOptionsBool("-test", "Testing mode (do not print unless error is large)", NULL, test_mode, &test_mode, NULL)); 96 PetscCall(PetscOptionsString("-mesh", "Read mesh from file", NULL, filename, filename, sizeof(filename), &read_mesh)); 97 PetscCall(PetscOptionsBool("-simplex", "Use simplices, or tensor product cells", NULL, simplex, &simplex, NULL)); 98 PetscCall(PetscOptionsInt("-degree", "Polynomial degree of tensor product basis", NULL, degree, °ree, NULL)); 99 PetscOptionsEnd(); 100 101 // Setup DM 102 if (read_mesh) { 103 PetscCall(DMPlexCreateFromFile(PETSC_COMM_WORLD, filename, NULL, PETSC_TRUE, &dm)); 104 } else { 105 // Create the mesh as a 0-refined sphere. This will create a cubic surface, not a box 106 PetscCall(DMPlexCreateSphereMesh(PETSC_COMM_WORLD, topo_dim, simplex, 1., &dm)); 107 if (problem_choice == CUBE) { 108 PetscCall(DMPlexCreateCoordinateSpace(dm, 1, NULL)); 109 } 110 // Set the object name 111 PetscCall(PetscObjectSetName((PetscObject)dm, problem_types[problem_choice])); 112 // Refine DMPlex with uniform refinement using runtime option -dm_refine 113 PetscCall(DMPlexSetRefinementUniform(dm, PETSC_TRUE)); 114 PetscCall(DMSetFromOptions(dm)); 115 // View DMPlex via runtime option 116 PetscCall(DMViewFromOptions(dm, NULL, "-dm_view")); 117 } 118 119 // Create DM 120 PetscCall(SetupDMByDegree(dm, degree, q_extra, num_comp_u, topo_dim, false)); 121 122 // Create vectors 123 PetscCall(DMCreateGlobalVector(dm, &U)); 124 PetscCall(VecGetLocalSize(U, &l_size)); 125 PetscCall(VecGetSize(U, &g_size)); 126 PetscCall(DMCreateLocalVector(dm, &U_loc)); 127 PetscCall(VecGetSize(U_loc, &xl_size)); 128 PetscCall(VecDuplicate(U, &V)); 129 PetscCall(VecDuplicate(U_loc, &V_loc)); 130 131 // Setup op_apply_ctx structure 132 PetscCall(PetscMalloc1(1, &op_apply_ctx)); 133 134 // Set up libCEED 135 CeedInit(ceed_resource, &ceed); 136 CeedMemType mem_type_backend; 137 CeedGetPreferredMemType(ceed, &mem_type_backend); 138 139 PetscCall(DMGetVecType(dm, &vec_type)); 140 if (!vec_type) { // Not yet set by op_apply_ctx -dm_vec_type 141 switch (mem_type_backend) { 142 case CEED_MEM_HOST: 143 vec_type = VECSTANDARD; 144 break; 145 case CEED_MEM_DEVICE: { 146 const char *resolved; 147 CeedGetResource(ceed, &resolved); 148 if (strstr(resolved, "/gpu/cuda")) vec_type = VECCUDA; 149 else if (strstr(resolved, "/gpu/hip/occa")) vec_type = VECSTANDARD; // https://github.com/CEED/libCEED/issues/678 150 else if (strstr(resolved, "/gpu/hip")) vec_type = VECHIP; 151 else vec_type = VECSTANDARD; 152 } 153 } 154 PetscCall(DMSetVecType(dm, vec_type)); 155 } 156 157 // Print summary 158 if (!test_mode) { 159 PetscInt P = degree + 1, Q = P + q_extra; 160 const char *used_resource; 161 CeedGetResource(ceed, &used_resource); 162 PetscCall(PetscPrintf(comm, 163 "\n-- libCEED + PETSc Surface Area of a Manifold --\n" 164 " libCEED:\n" 165 " libCEED Backend : %s\n" 166 " libCEED Backend MemType : %s\n" 167 " Mesh:\n" 168 " Solution Order (P) : %" CeedInt_FMT "\n" 169 " Quadrature Order (Q) : %" CeedInt_FMT "\n" 170 " Additional quadrature points (q_extra) : %" CeedInt_FMT "\n" 171 " Global nodes : %" PetscInt_FMT "\n" 172 " DoF per node : %" PetscInt_FMT "\n" 173 " Global DoFs : %" PetscInt_FMT "\n", 174 used_resource, CeedMemTypes[mem_type_backend], P, Q, q_extra, g_size / num_comp_u, num_comp_u, g_size)); 175 } 176 177 // Setup libCEED's objects and apply setup operator 178 PetscCall(PetscMalloc1(1, &ceed_data)); 179 PetscCall(SetupLibceedByDegree(dm, ceed, degree, topo_dim, q_extra, num_comp_x, num_comp_u, g_size, xl_size, problem_options[problem_choice], 180 ceed_data, false, (CeedVector)NULL, (CeedVector *)NULL)); 181 182 // Setup output vector 183 PetscScalar *v; 184 PetscCall(VecZeroEntries(V_loc)); 185 PetscCall(VecGetArrayAndMemType(V_loc, &v, &mem_type)); 186 CeedVectorSetArray(ceed_data->y_ceed, MemTypeP2C(mem_type), CEED_USE_POINTER, v); 187 188 // Compute the mesh volume using the mass operator: area = 1^T \cdot M \cdot 1 189 if (!test_mode) { 190 PetscCall(PetscPrintf(comm, "Computing the mesh area using the formula: area = 1^T M 1\n")); 191 } 192 193 // Initialize u with ones 194 CeedVectorSetValue(ceed_data->x_ceed, 1.0); 195 196 // Apply the mass operator: 'u' -> 'v' 197 CeedOperatorApply(ceed_data->op_apply, ceed_data->x_ceed, ceed_data->y_ceed, CEED_REQUEST_IMMEDIATE); 198 199 // Gather output vector 200 CeedVectorTakeArray(ceed_data->y_ceed, CEED_MEM_HOST, NULL); 201 PetscCall(VecRestoreArrayAndMemType(V_loc, &v)); 202 PetscCall(VecZeroEntries(V)); 203 PetscCall(DMLocalToGlobalBegin(dm, V_loc, ADD_VALUES, V)); 204 PetscCall(DMLocalToGlobalEnd(dm, V_loc, ADD_VALUES, V)); 205 206 // Compute and print the sum of the entries of 'v' giving the mesh surface area 207 PetscScalar area; 208 PetscCall(VecSum(V, &area)); 209 210 // Compute the exact surface area and print the result 211 CeedScalar exact_surface_area = 4 * M_PI; 212 if (problem_choice == CUBE) { 213 exact_surface_area = 6 * 2 * 2; // surface of [-1, 1]^3 214 } 215 216 PetscReal error = fabs(area - exact_surface_area); 217 PetscReal tol = 5e-6; 218 if (!test_mode || error > tol) { 219 PetscCall(PetscPrintf(comm, "Exact mesh surface area : % .14g\n", exact_surface_area)); 220 PetscCall(PetscPrintf(comm, "Computed mesh surface area : % .14g\n", area)); 221 PetscCall(PetscPrintf(comm, "Area error : % .14g\n", error)); 222 } 223 224 // Cleanup 225 PetscCall(DMDestroy(&dm)); 226 PetscCall(VecDestroy(&U)); 227 PetscCall(VecDestroy(&U_loc)); 228 PetscCall(VecDestroy(&V)); 229 PetscCall(VecDestroy(&V_loc)); 230 PetscCall(PetscFree(op_apply_ctx)); 231 PetscCall(CeedDataDestroy(0, ceed_data)); 232 CeedDestroy(&ceed); 233 return PetscFinalize(); 234 } 235