1 // Copyright (c) 2017, Lawrence Livermore National Security, LLC. Produced at 2 // the Lawrence Livermore National Laboratory. LLNL-CODE-734707. All Rights 3 // reserved. See files LICENSE and NOTICE for details. 4 // 5 // This file is part of CEED, a collection of benchmarks, miniapps, software 6 // libraries and APIs for efficient high-order finite element and spectral 7 // element discretizations for exascale applications. For more information and 8 // source code availability see http://github.com/ceed. 9 // 10 // The CEED research is supported by the Exascale Computing Project 17-SC-20-SC, 11 // a collaborative effort of two U.S. Department of Energy organizations (Office 12 // of Science and the National Nuclear Security Administration) responsible for 13 // the planning and preparation of a capable exascale ecosystem, including 14 // software, applications, hardware, advanced system engineering and early 15 // testbed platforms, in support of the nation's exascale computing imperative. 16 17 /// @file 18 /// Helper functions for solid mechanics example using PETSc 19 20 #include "../elasticity.h" 21 22 // ----------------------------------------------------------------------------- 23 // Create libCEED operator context 24 // ----------------------------------------------------------------------------- 25 // Setup context data for Jacobian evaluation 26 PetscErrorCode SetupJacobianCtx(MPI_Comm comm, AppCtx app_ctx, DM dm, Vec V, 27 Vec V_loc, CeedData ceed_data, Ceed ceed, 28 CeedQFunctionContext ctx_phys, 29 CeedQFunctionContext ctx_phys_smoother, 30 UserMult jacobian_ctx) { 31 PetscErrorCode ierr; 32 33 PetscFunctionBeginUser; 34 35 // PETSc objects 36 jacobian_ctx->comm = comm; 37 jacobian_ctx->dm = dm; 38 39 // Work vectors 40 jacobian_ctx->X_loc = V_loc; 41 ierr = VecDuplicate(V_loc, &jacobian_ctx->Y_loc); CHKERRQ(ierr); 42 jacobian_ctx->x_ceed = ceed_data->x_ceed; 43 jacobian_ctx->y_ceed = ceed_data->y_ceed; 44 45 // libCEED operator 46 jacobian_ctx->op = ceed_data->op_jacob; 47 jacobian_ctx->qf = ceed_data->qf_jacob; 48 49 // Ceed 50 jacobian_ctx->ceed = ceed; 51 52 // Physics 53 jacobian_ctx->ctx_phys = ctx_phys; 54 jacobian_ctx->ctx_phys_smoother = ctx_phys_smoother; 55 PetscFunctionReturn(0); 56 }; 57 58 // Setup context data for prolongation and restriction operators 59 PetscErrorCode SetupProlongRestrictCtx(MPI_Comm comm, AppCtx app_ctx, DM dm_c, 60 DM dm_f, Vec V_f, Vec V_loc_c, Vec V_loc_f, 61 CeedData ceed_data_c, CeedData ceed_data_f, 62 Ceed ceed, UserMultProlongRestr prolong_restr_ctx) { 63 PetscFunctionBeginUser; 64 65 // PETSc objects 66 prolong_restr_ctx->comm = comm; 67 prolong_restr_ctx->dm_c = dm_c; 68 prolong_restr_ctx->dm_f = dm_f; 69 70 // Work vectors 71 prolong_restr_ctx->loc_vec_c = V_loc_c; 72 prolong_restr_ctx->loc_vec_f = V_loc_f; 73 prolong_restr_ctx->ceed_vec_c = ceed_data_c->x_ceed; 74 prolong_restr_ctx->ceed_vec_f = ceed_data_f->x_ceed; 75 76 // libCEED operators 77 prolong_restr_ctx->op_prolong = ceed_data_f->op_prolong; 78 prolong_restr_ctx->op_restrict = ceed_data_f->op_restrict; 79 80 // Ceed 81 prolong_restr_ctx->ceed = ceed; 82 PetscFunctionReturn(0); 83 }; 84 85 // ----------------------------------------------------------------------------- 86 // Jacobian setup 87 // ----------------------------------------------------------------------------- 88 PetscErrorCode FormJacobian(SNES snes, Vec U, Mat J, Mat J_pre, void *ctx) { 89 PetscErrorCode ierr; 90 91 PetscFunctionBeginUser; 92 93 // Context data 94 FormJacobCtx form_jacob_ctx = (FormJacobCtx)ctx; 95 PetscInt num_levels = form_jacob_ctx->num_levels; 96 Mat *jacob_mat = form_jacob_ctx->jacob_mat; 97 98 // Update Jacobian on each level 99 for (PetscInt level = 0; level < num_levels; level++) { 100 ierr = MatAssemblyBegin(jacob_mat[level], MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 101 ierr = MatAssemblyEnd(jacob_mat[level], MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 102 } 103 104 // Form coarse assembled matrix 105 ierr = VecZeroEntries(form_jacob_ctx->u_coarse); CHKERRQ(ierr); 106 ierr = SNESComputeJacobianDefaultColor(form_jacob_ctx->snes_coarse, 107 form_jacob_ctx->u_coarse, 108 form_jacob_ctx->jacob_mat[0], 109 form_jacob_ctx->jacob_mat_coarse, NULL); 110 CHKERRQ(ierr); 111 112 // J_pre might be AIJ (e.g., when using coloring), so we need to assemble it 113 ierr = MatAssemblyBegin(J_pre, MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 114 ierr = MatAssemblyEnd(J_pre, MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 115 if (J != J_pre) { 116 ierr = MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 117 ierr = MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 118 } 119 PetscFunctionReturn(0); 120 }; 121 122 // ----------------------------------------------------------------------------- 123 // Output solution for visualization 124 // ----------------------------------------------------------------------------- 125 PetscErrorCode ViewSolution(MPI_Comm comm, AppCtx app_ctx, Vec U, 126 PetscInt increment, PetscScalar load_increment) { 127 PetscErrorCode ierr; 128 DM dm; 129 PetscViewer viewer; 130 char output_filename[PETSC_MAX_PATH_LEN]; 131 PetscMPIInt rank; 132 133 PetscFunctionBeginUser; 134 135 // Create output directory 136 MPI_Comm_rank(comm, &rank); 137 if (!rank) {ierr = PetscMkdir(app_ctx->output_dir); CHKERRQ(ierr);} 138 139 // Build file name 140 ierr = PetscSNPrintf(output_filename, sizeof output_filename, 141 "%s/solution-%03D.vtu", app_ctx->output_dir, 142 increment); CHKERRQ(ierr); 143 144 // Increment sequence 145 ierr = VecGetDM(U, &dm); CHKERRQ(ierr); 146 ierr = DMSetOutputSequenceNumber(dm, increment, load_increment); CHKERRQ(ierr); 147 148 // Output solution vector 149 ierr = PetscViewerVTKOpen(comm, output_filename, FILE_MODE_WRITE, &viewer); 150 CHKERRQ(ierr); 151 ierr = VecView(U, viewer); CHKERRQ(ierr); 152 ierr = PetscViewerDestroy(&viewer); CHKERRQ(ierr); 153 154 PetscFunctionReturn(0); 155 }; 156 157 // ----------------------------------------------------------------------------- 158 // Output diagnostic quantities for visualization 159 // ----------------------------------------------------------------------------- 160 PetscErrorCode ViewDiagnosticQuantities(MPI_Comm comm, DM dmU, 161 UserMult user, AppCtx app_ctx, Vec U, 162 CeedElemRestriction elem_restr_diagnostic) { 163 PetscErrorCode ierr; 164 Vec Diagnostic, Y_loc, mult_vec; 165 CeedVector y_ceed; 166 CeedScalar *x, *y; 167 PetscMemType x_mem_type, y_mem_type; 168 PetscInt loc_size; 169 PetscViewer viewer; 170 char output_filename[PETSC_MAX_PATH_LEN]; 171 172 PetscFunctionBeginUser; 173 174 // --------------------------------------------------------------------------- 175 // PETSc and libCEED vectors 176 // --------------------------------------------------------------------------- 177 ierr = DMCreateGlobalVector(user->dm, &Diagnostic); CHKERRQ(ierr); 178 ierr = PetscObjectSetName((PetscObject)Diagnostic, ""); CHKERRQ(ierr); 179 ierr = DMCreateLocalVector(user->dm, &Y_loc); CHKERRQ(ierr); 180 ierr = VecGetSize(Y_loc, &loc_size); CHKERRQ(ierr); 181 CeedVectorCreate(user->ceed, loc_size, &y_ceed); 182 183 // --------------------------------------------------------------------------- 184 // Compute quantities 185 // --------------------------------------------------------------------------- 186 // -- Global-to-local 187 ierr = VecZeroEntries(user->X_loc); CHKERRQ(ierr); 188 ierr = DMPlexInsertBoundaryValues(dmU, PETSC_TRUE, user->X_loc, 189 user->load_increment, NULL, NULL, NULL); 190 CHKERRQ(ierr); 191 ierr = DMGlobalToLocal(dmU, U, INSERT_VALUES, user->X_loc); CHKERRQ(ierr); 192 ierr = VecZeroEntries(Y_loc); CHKERRQ(ierr); 193 194 // -- Setup CEED vectors 195 ierr = VecGetArrayReadAndMemType(user->X_loc, (const PetscScalar **)&x, 196 &x_mem_type); CHKERRQ(ierr); 197 ierr = VecGetArrayAndMemType(Y_loc, &y, &y_mem_type); CHKERRQ(ierr); 198 CeedVectorSetArray(user->x_ceed, MemTypeP2C(x_mem_type), CEED_USE_POINTER, x); 199 CeedVectorSetArray(y_ceed, MemTypeP2C(y_mem_type), CEED_USE_POINTER, y); 200 201 // -- Apply CEED operator 202 CeedOperatorApply(user->op, user->x_ceed, y_ceed, CEED_REQUEST_IMMEDIATE); 203 204 // -- Restore PETSc vector; keep y_ceed viewing memory of Y_loc for use below 205 CeedVectorTakeArray(user->x_ceed, MemTypeP2C(x_mem_type), NULL); 206 ierr = VecRestoreArrayReadAndMemType(user->X_loc, (const PetscScalar **)&x); 207 CHKERRQ(ierr); 208 209 // -- Local-to-global 210 ierr = VecZeroEntries(Diagnostic); CHKERRQ(ierr); 211 ierr = DMLocalToGlobal(user->dm, Y_loc, ADD_VALUES, Diagnostic); 212 CHKERRQ(ierr); 213 214 // --------------------------------------------------------------------------- 215 // Scale for multiplicity 216 // --------------------------------------------------------------------------- 217 // -- Setup vectors 218 ierr = VecDuplicate(Diagnostic, &mult_vec); CHKERRQ(ierr); 219 ierr = VecZeroEntries(Y_loc); CHKERRQ(ierr); 220 221 // -- Compute multiplicity 222 CeedElemRestrictionGetMultiplicity(elem_restr_diagnostic, y_ceed); 223 224 // -- Restore vectors 225 CeedVectorTakeArray(y_ceed, MemTypeP2C(y_mem_type), NULL); 226 ierr = VecRestoreArrayAndMemType(Y_loc, &y); CHKERRQ(ierr); 227 228 // -- Local-to-global 229 ierr = VecZeroEntries(mult_vec); CHKERRQ(ierr); 230 ierr = DMLocalToGlobal(user->dm, Y_loc, ADD_VALUES, mult_vec); 231 CHKERRQ(ierr); 232 233 // -- Scale 234 ierr = VecReciprocal(mult_vec); CHKERRQ(ierr); 235 ierr = VecPointwiseMult(Diagnostic, Diagnostic, mult_vec); 236 237 // --------------------------------------------------------------------------- 238 // Output solution vector 239 // --------------------------------------------------------------------------- 240 ierr = PetscSNPrintf(output_filename, sizeof output_filename, 241 "%s/diagnostic_quantities.vtu", 242 app_ctx->output_dir); CHKERRQ(ierr); 243 ierr = PetscViewerVTKOpen(comm, output_filename, FILE_MODE_WRITE, &viewer); 244 CHKERRQ(ierr); 245 ierr = VecView(Diagnostic, viewer); CHKERRQ(ierr); 246 ierr = PetscViewerDestroy(&viewer); CHKERRQ(ierr); 247 248 // --------------------------------------------------------------------------- 249 // Cleanup 250 // --------------------------------------------------------------------------- 251 ierr = VecDestroy(&Diagnostic); CHKERRQ(ierr); 252 ierr = VecDestroy(&mult_vec); CHKERRQ(ierr); 253 ierr = VecDestroy(&Y_loc); CHKERRQ(ierr); 254 CeedVectorDestroy(&y_ceed); 255 256 PetscFunctionReturn(0); 257 }; 258