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 appCtx, DM dm, Vec V, 27 Vec Vloc, CeedData ceedData, Ceed ceed, 28 UserMult jacobianCtx) { 29 PetscErrorCode ierr; 30 31 PetscFunctionBeginUser; 32 33 // PETSc objects 34 jacobianCtx->comm = comm; 35 jacobianCtx->dm = dm; 36 37 // Work vectors 38 jacobianCtx->Xloc = Vloc; 39 ierr = VecDuplicate(Vloc, &jacobianCtx->Yloc); CHKERRQ(ierr); 40 jacobianCtx->Xceed = ceedData->xceed; 41 jacobianCtx->Yceed = ceedData->yceed; 42 43 // libCEED operator 44 jacobianCtx->op = ceedData->opJacob; 45 46 // Ceed 47 jacobianCtx->ceed = ceed; 48 49 PetscFunctionReturn(0); 50 }; 51 52 // Setup context data for prolongation and restriction operators 53 PetscErrorCode SetupProlongRestrictCtx(MPI_Comm comm, DM dmC, DM dmF, Vec VF, 54 Vec VlocC, Vec VlocF, CeedData ceedDataC, 55 CeedData ceedDataF, Ceed ceed, 56 UserMultProlongRestr prolongRestrCtx) { 57 PetscErrorCode ierr; 58 PetscScalar *multArray; 59 60 PetscFunctionBeginUser; 61 62 // PETSc objects 63 prolongRestrCtx->comm = comm; 64 prolongRestrCtx->dmC = dmC; 65 prolongRestrCtx->dmF = dmF; 66 67 // Work vectors 68 prolongRestrCtx->locVecC = VlocC; 69 prolongRestrCtx->locVecF = VlocF; 70 prolongRestrCtx->ceedVecC = ceedDataC->xceed; 71 prolongRestrCtx->ceedVecF = ceedDataF->xceed; 72 73 // libCEED operators 74 prolongRestrCtx->opProlong = ceedDataF->opProlong; 75 prolongRestrCtx->opRestrict = ceedDataF->opRestrict; 76 77 // Ceed 78 prolongRestrCtx->ceed = ceed; 79 80 // Multiplicity vector 81 // -- Set libCEED vector 82 ierr = VecZeroEntries(VlocF); 83 ierr = VecGetArray(VlocF, &multArray); CHKERRQ(ierr); 84 CeedVectorSetArray(ceedDataF->xceed, CEED_MEM_HOST, CEED_USE_POINTER, 85 multArray); 86 87 // -- Compute multiplicity 88 CeedElemRestrictionGetMultiplicity(ceedDataF->Erestrictu, ceedDataF->xceed); 89 90 // -- Restore PETSc vector 91 ierr = VecRestoreArray(VlocF, &multArray); CHKERRQ(ierr); 92 93 // -- Local-to-global 94 ierr = VecZeroEntries(VF); CHKERRQ(ierr); 95 ierr = DMLocalToGlobal(dmF, VlocF, ADD_VALUES, VF); CHKERRQ(ierr); 96 97 // -- Global-to-local 98 ierr = VecDuplicate(VlocF, &prolongRestrCtx->multVec); CHKERRQ(ierr); 99 ierr = DMGlobalToLocal(dmF, VF, INSERT_VALUES, prolongRestrCtx->multVec); 100 CHKERRQ(ierr); 101 102 // -- Reciprocal 103 ierr = VecReciprocal(prolongRestrCtx->multVec); CHKERRQ(ierr); 104 105 // -- Reset work arrays 106 ierr = VecZeroEntries(VF); CHKERRQ(ierr); 107 ierr = VecZeroEntries(VlocF); CHKERRQ(ierr); 108 109 PetscFunctionReturn(0); 110 }; 111 112 // ----------------------------------------------------------------------------- 113 // Jacobian setup 114 // ----------------------------------------------------------------------------- 115 PetscErrorCode FormJacobian(SNES snes, Vec U, Mat J, Mat Jpre, void *ctx) { 116 PetscErrorCode ierr; 117 118 PetscFunctionBeginUser; 119 120 // Context data 121 FormJacobCtx formJacobCtx = (FormJacobCtx)ctx; 122 PetscInt numLevels = formJacobCtx->numLevels; 123 Mat *jacobMat = formJacobCtx->jacobMat; 124 125 // Update Jacobian on each level 126 for (PetscInt level = 0; level < numLevels; level++) { 127 ierr = MatAssemblyBegin(jacobMat[level], MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 128 ierr = MatAssemblyEnd(jacobMat[level], MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 129 } 130 131 // Form coarse assembled matrix 132 ierr = VecZeroEntries(formJacobCtx->Ucoarse); CHKERRQ(ierr); 133 ierr = SNESComputeJacobianDefaultColor(formJacobCtx->snesCoarse, 134 formJacobCtx->Ucoarse, 135 formJacobCtx->jacobMat[0], 136 formJacobCtx->jacobMatCoarse, NULL); 137 CHKERRQ(ierr); 138 139 // Jpre might be AIJ (e.g., when using coloring), so we need to assemble it 140 ierr = MatAssemblyBegin(Jpre, MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 141 ierr = MatAssemblyEnd(Jpre, MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 142 if (J != Jpre) { 143 ierr = MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 144 ierr = MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 145 } 146 PetscFunctionReturn(0); 147 }; 148 149 // ----------------------------------------------------------------------------- 150 // Output solution for visualization 151 // ----------------------------------------------------------------------------- 152 PetscErrorCode ViewSolution(MPI_Comm comm, Vec U, PetscInt increment, 153 PetscScalar loadIncrement) { 154 PetscErrorCode ierr; 155 DM dm; 156 PetscViewer viewer; 157 char outputFilename[PETSC_MAX_PATH_LEN]; 158 159 PetscFunctionBeginUser; 160 161 // Build file name 162 ierr = PetscSNPrintf(outputFilename, sizeof outputFilename, 163 "solution-%03D.vtu", increment); CHKERRQ(ierr); 164 165 // Increment senquence 166 ierr = VecGetDM(U, &dm); CHKERRQ(ierr); 167 ierr = DMSetOutputSequenceNumber(dm, increment, loadIncrement); CHKERRQ(ierr); 168 169 // Output solution vector 170 ierr = PetscViewerVTKOpen(comm, outputFilename, FILE_MODE_WRITE, &viewer); 171 CHKERRQ(ierr); 172 ierr = VecView(U, viewer); CHKERRQ(ierr); 173 ierr = PetscViewerDestroy(&viewer); CHKERRQ(ierr); 174 175 PetscFunctionReturn(0); 176 }; 177