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 PC smoother on each level 126 for (int level = 0; level < numLevels; level++) { 127 // -- Update diagonal state counter 128 ierr = MatAssemblyBegin(jacobMat[level], MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 129 ierr = MatAssemblyEnd(jacobMat[level], MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 130 } 131 132 // Form coarse assembled matrix 133 ierr = VecZeroEntries(formJacobCtx->Ucoarse); CHKERRQ(ierr); 134 ierr = SNESComputeJacobianDefaultColor(formJacobCtx->snesCoarse, 135 formJacobCtx->Ucoarse, 136 formJacobCtx->jacobMat[0], 137 formJacobCtx->jacobMatCoarse, NULL); 138 CHKERRQ(ierr); 139 140 // Jpre might be AIJ (e.g., when using coloring), so we need to assemble it 141 ierr = MatAssemblyBegin(Jpre, MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 142 ierr = MatAssemblyEnd(Jpre, MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 143 144 PetscFunctionReturn(0); 145 }; 146 147 // ----------------------------------------------------------------------------- 148 // Output solution for visualization 149 // ----------------------------------------------------------------------------- 150 PetscErrorCode ViewSolution(MPI_Comm comm, Vec U, PetscInt increment, 151 PetscScalar loadIncrement) { 152 PetscErrorCode ierr; 153 DM dm; 154 PetscViewer viewer; 155 char outputFilename[PETSC_MAX_PATH_LEN]; 156 157 PetscFunctionBeginUser; 158 159 // Build file name 160 ierr = PetscSNPrintf(outputFilename, sizeof outputFilename, 161 "solution-%03D.vtu", increment); CHKERRQ(ierr); 162 163 // Increment senquence 164 ierr = VecGetDM(U, &dm); CHKERRQ(ierr); 165 ierr = DMSetOutputSequenceNumber(dm, increment, loadIncrement); CHKERRQ(ierr); 166 167 // Output solution vector 168 ierr = PetscViewerVTKOpen(comm, outputFilename, FILE_MODE_WRITE, &viewer); 169 CHKERRQ(ierr); 170 ierr = VecView(U, viewer); CHKERRQ(ierr); 171 ierr = PetscViewerDestroy(&viewer); CHKERRQ(ierr); 172 173 PetscFunctionReturn(0); 174 }; 175