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 /// Setup libCEED for Navier-Stokes example using PETSc 19 20 #include "../navierstokes.h" 21 22 // Utility function - essential BC dofs are encoded in closure indices as -(i+1). 23 PetscInt Involute(PetscInt i) { 24 return i >= 0 ? i : -(i+1); 25 } 26 27 // Utility function to create local CEED restriction 28 PetscErrorCode CreateRestrictionFromPlex(Ceed ceed, DM dm, CeedInt height, 29 DMLabel domain_label, CeedInt value, CeedElemRestriction *elem_restr) { 30 PetscInt num_elem, elem_size, num_dof, num_comp, *elem_restr_offsets; 31 PetscErrorCode ierr; 32 33 PetscFunctionBeginUser; 34 ierr = DMPlexGetLocalOffsets(dm, domain_label, value, height, 0, &num_elem, 35 &elem_size, &num_comp, &num_dof, &elem_restr_offsets); 36 CHKERRQ(ierr); 37 38 CeedElemRestrictionCreate(ceed, num_elem, elem_size, num_comp, 39 1, num_dof, CEED_MEM_HOST, CEED_COPY_VALUES, 40 elem_restr_offsets, elem_restr); 41 ierr = PetscFree(elem_restr_offsets); CHKERRQ(ierr); 42 43 PetscFunctionReturn(0); 44 } 45 46 // Utility function to get Ceed Restriction for each domain 47 PetscErrorCode GetRestrictionForDomain(Ceed ceed, DM dm, CeedInt height, 48 DMLabel domain_label, PetscInt value, 49 CeedInt Q, CeedInt q_data_size, 50 CeedElemRestriction *elem_restr_q, 51 CeedElemRestriction *elem_restr_x, 52 CeedElemRestriction *elem_restr_qd_i) { 53 DM dm_coord; 54 CeedInt dim, loc_num_elem; 55 CeedInt Q_dim; 56 PetscErrorCode ierr; 57 PetscFunctionBeginUser; 58 59 ierr = DMGetDimension(dm, &dim); CHKERRQ(ierr); 60 dim -= height; 61 Q_dim = CeedIntPow(Q, dim); 62 ierr = DMGetCoordinateDM(dm, &dm_coord); CHKERRQ(ierr); 63 ierr = DMPlexSetClosurePermutationTensor(dm_coord, PETSC_DETERMINE, NULL); 64 CHKERRQ(ierr); 65 ierr = CreateRestrictionFromPlex(ceed, dm, height, domain_label, value, 66 elem_restr_q); 67 CHKERRQ(ierr); 68 ierr = CreateRestrictionFromPlex(ceed, dm_coord, height, domain_label, value, 69 elem_restr_x); 70 CHKERRQ(ierr); 71 CeedElemRestrictionGetNumElements(*elem_restr_q, &loc_num_elem); 72 CeedElemRestrictionCreateStrided(ceed, loc_num_elem, Q_dim, 73 q_data_size, q_data_size*loc_num_elem*Q_dim, 74 CEED_STRIDES_BACKEND, elem_restr_qd_i); 75 PetscFunctionReturn(0); 76 } 77 78 // Utility function to create CEED Composite Operator for the entire domain 79 PetscErrorCode CreateOperatorForDomain(Ceed ceed, DM dm, SimpleBC bc, 80 CeedData ceed_data, Physics phys, 81 CeedOperator op_apply_vol, CeedInt height, 82 CeedInt P_sur, CeedInt Q_sur, CeedInt q_data_size_sur, 83 CeedOperator *op_apply) { 84 CeedInt dim, num_face; 85 DMLabel domain_label; 86 PetscErrorCode ierr; 87 PetscFunctionBeginUser; 88 89 // Create Composite Operaters 90 CeedCompositeOperatorCreate(ceed, op_apply); 91 92 // --Apply Sub-Operator for the volume 93 CeedCompositeOperatorAddSub(*op_apply, op_apply_vol); 94 95 // -- Create Sub-Operator for in/outflow BCs 96 if (phys->has_neumann) { 97 // --- Setup 98 ierr = DMGetLabel(dm, "Face Sets", &domain_label); CHKERRQ(ierr); 99 ierr = DMGetDimension(dm, &dim); CHKERRQ(ierr); 100 if (dim == 2) num_face = 4; 101 if (dim == 3) num_face = 6; 102 103 // --- Get number of quadrature points for the boundaries 104 CeedInt num_qpts_sur; 105 CeedBasisGetNumQuadraturePoints(ceed_data->basis_q_sur, &num_qpts_sur); 106 107 // --- Create Sub-Operator for each face 108 for (CeedInt i=0; i<num_face; i++) { 109 CeedVector q_data_sur; 110 CeedOperator op_setup_sur, op_apply_sur; 111 CeedElemRestriction elem_restr_x_sur, elem_restr_q_sur, 112 elem_restr_qd_i_sur; 113 114 // ---- CEED Restriction 115 ierr = GetRestrictionForDomain(ceed, dm, height, domain_label, i+1, 116 Q_sur, q_data_size_sur, &elem_restr_q_sur, 117 &elem_restr_x_sur, &elem_restr_qd_i_sur); 118 CHKERRQ(ierr); 119 120 // ---- CEED Vector 121 PetscInt loc_num_elem_sur; 122 CeedElemRestrictionGetNumElements(elem_restr_q_sur, &loc_num_elem_sur); 123 CeedVectorCreate(ceed, q_data_size_sur*loc_num_elem_sur*num_qpts_sur, 124 &q_data_sur); 125 126 // ---- CEED Operator 127 // ----- CEED Operator for Setup (geometric factors) 128 CeedOperatorCreate(ceed, ceed_data->qf_setup_sur, NULL, NULL, &op_setup_sur); 129 CeedOperatorSetField(op_setup_sur, "dx", elem_restr_x_sur, 130 ceed_data->basis_x_sur, CEED_VECTOR_ACTIVE); 131 CeedOperatorSetField(op_setup_sur, "weight", CEED_ELEMRESTRICTION_NONE, 132 ceed_data->basis_x_sur, CEED_VECTOR_NONE); 133 CeedOperatorSetField(op_setup_sur, "surface qdata", elem_restr_qd_i_sur, 134 CEED_BASIS_COLLOCATED, CEED_VECTOR_ACTIVE); 135 136 // ----- CEED Operator for Physics 137 CeedOperatorCreate(ceed, ceed_data->qf_apply_sur, NULL, NULL, &op_apply_sur); 138 CeedOperatorSetField(op_apply_sur, "q", elem_restr_q_sur, 139 ceed_data->basis_q_sur, CEED_VECTOR_ACTIVE); 140 CeedOperatorSetField(op_apply_sur, "surface qdata", elem_restr_qd_i_sur, 141 CEED_BASIS_COLLOCATED, q_data_sur); 142 CeedOperatorSetField(op_apply_sur, "x", elem_restr_x_sur, 143 ceed_data->basis_x_sur, ceed_data->x_coord); 144 CeedOperatorSetField(op_apply_sur, "v", elem_restr_q_sur, 145 ceed_data->basis_q_sur, CEED_VECTOR_ACTIVE); 146 147 // ----- Apply CEED operator for Setup 148 CeedOperatorApply(op_setup_sur, ceed_data->x_coord, q_data_sur, 149 CEED_REQUEST_IMMEDIATE); 150 151 // ----- Apply Sub-Operator for the Boundary 152 CeedCompositeOperatorAddSub(*op_apply, op_apply_sur); 153 154 // ----- Cleanup 155 CeedVectorDestroy(&q_data_sur); 156 CeedElemRestrictionDestroy(&elem_restr_q_sur); 157 CeedElemRestrictionDestroy(&elem_restr_x_sur); 158 CeedElemRestrictionDestroy(&elem_restr_qd_i_sur); 159 CeedOperatorDestroy(&op_setup_sur); 160 CeedOperatorDestroy(&op_apply_sur); 161 } 162 } 163 164 PetscFunctionReturn(0); 165 } 166 167 PetscErrorCode SetupLibceed(Ceed ceed, CeedData ceed_data, DM dm, User user, 168 AppCtx app_ctx, ProblemData *problem, SimpleBC bc) { 169 PetscErrorCode ierr; 170 PetscFunctionBeginUser; 171 172 // ***************************************************************************** 173 // Set up CEED objects for the interior domain (volume) 174 // ***************************************************************************** 175 const PetscInt num_comp_q = 5; 176 const CeedInt dim = problem->dim, 177 num_comp_x = problem->dim, 178 q_data_size_vol = problem->q_data_size_vol, 179 P = app_ctx->degree + 1, 180 Q = P + app_ctx->q_extra; 181 182 // ----------------------------------------------------------------------------- 183 // CEED Bases 184 // ----------------------------------------------------------------------------- 185 CeedBasisCreateTensorH1Lagrange(ceed, dim, num_comp_q, P, Q, CEED_GAUSS, 186 &ceed_data->basis_q); 187 CeedBasisCreateTensorH1Lagrange(ceed, dim, num_comp_x, 2, Q, CEED_GAUSS, 188 &ceed_data->basis_x); 189 CeedBasisCreateTensorH1Lagrange(ceed, dim, num_comp_x, 2, P, 190 CEED_GAUSS_LOBATTO, &ceed_data->basis_xc); 191 192 // ----------------------------------------------------------------------------- 193 // CEED Restrictions 194 // ----------------------------------------------------------------------------- 195 // -- Create restriction 196 ierr = GetRestrictionForDomain(ceed, dm, 0, 0, 0, Q, q_data_size_vol, 197 &ceed_data->elem_restr_q, &ceed_data->elem_restr_x, 198 &ceed_data->elem_restr_qd_i); CHKERRQ(ierr); 199 // -- Create E vectors 200 CeedElemRestrictionCreateVector(ceed_data->elem_restr_q, &user->q_ceed, NULL); 201 CeedElemRestrictionCreateVector(ceed_data->elem_restr_q, &user->q_dot_ceed, 202 NULL); 203 CeedElemRestrictionCreateVector(ceed_data->elem_restr_q, &user->g_ceed, NULL); 204 205 // ----------------------------------------------------------------------------- 206 // CEED QFunctions 207 // ----------------------------------------------------------------------------- 208 // -- Create QFunction for quadrature data 209 CeedQFunctionCreateInterior(ceed, 1, problem->setup_vol, problem->setup_vol_loc, 210 &ceed_data->qf_setup_vol); 211 CeedQFunctionAddInput(ceed_data->qf_setup_vol, "dx", num_comp_x*dim, 212 CEED_EVAL_GRAD); 213 CeedQFunctionAddInput(ceed_data->qf_setup_vol, "weight", 1, CEED_EVAL_WEIGHT); 214 CeedQFunctionAddOutput(ceed_data->qf_setup_vol, "qdata", q_data_size_vol, 215 CEED_EVAL_NONE); 216 217 // -- Create QFunction for ICs 218 CeedQFunctionCreateInterior(ceed, 1, problem->ics, problem->ics_loc, 219 &ceed_data->qf_ics); 220 CeedQFunctionAddInput(ceed_data->qf_ics, "x", num_comp_x, CEED_EVAL_INTERP); 221 CeedQFunctionAddOutput(ceed_data->qf_ics, "q0", num_comp_q, CEED_EVAL_NONE); 222 223 // -- Create QFunction for RHS 224 if (problem->apply_vol_rhs) { 225 CeedQFunctionCreateInterior(ceed, 1, problem->apply_vol_rhs, 226 problem->apply_vol_rhs_loc, &ceed_data->qf_rhs_vol); 227 CeedQFunctionAddInput(ceed_data->qf_rhs_vol, "q", num_comp_q, CEED_EVAL_INTERP); 228 CeedQFunctionAddInput(ceed_data->qf_rhs_vol, "dq", num_comp_q*dim, 229 CEED_EVAL_GRAD); 230 CeedQFunctionAddInput(ceed_data->qf_rhs_vol, "qdata", q_data_size_vol, 231 CEED_EVAL_NONE); 232 CeedQFunctionAddInput(ceed_data->qf_rhs_vol, "x", num_comp_x, CEED_EVAL_INTERP); 233 CeedQFunctionAddOutput(ceed_data->qf_rhs_vol, "v", num_comp_q, 234 CEED_EVAL_INTERP); 235 CeedQFunctionAddOutput(ceed_data->qf_rhs_vol, "dv", num_comp_q*dim, 236 CEED_EVAL_GRAD); 237 } 238 239 // -- Create QFunction for IFunction 240 if (problem->apply_vol_ifunction) { 241 CeedQFunctionCreateInterior(ceed, 1, problem->apply_vol_ifunction, 242 problem->apply_vol_ifunction_loc, &ceed_data->qf_ifunction_vol); 243 CeedQFunctionAddInput(ceed_data->qf_ifunction_vol, "q", num_comp_q, 244 CEED_EVAL_INTERP); 245 CeedQFunctionAddInput(ceed_data->qf_ifunction_vol, "dq", num_comp_q*dim, 246 CEED_EVAL_GRAD); 247 CeedQFunctionAddInput(ceed_data->qf_ifunction_vol, "q dot", num_comp_q, 248 CEED_EVAL_INTERP); 249 CeedQFunctionAddInput(ceed_data->qf_ifunction_vol, "qdata", q_data_size_vol, 250 CEED_EVAL_NONE); 251 CeedQFunctionAddInput(ceed_data->qf_ifunction_vol, "x", num_comp_x, 252 CEED_EVAL_INTERP); 253 CeedQFunctionAddOutput(ceed_data->qf_ifunction_vol, "v", num_comp_q, 254 CEED_EVAL_INTERP); 255 CeedQFunctionAddOutput(ceed_data->qf_ifunction_vol, "dv", num_comp_q*dim, 256 CEED_EVAL_GRAD); 257 } 258 259 // --------------------------------------------------------------------------- 260 // Element coordinates 261 // --------------------------------------------------------------------------- 262 // -- Create CEED vector 263 CeedElemRestrictionCreateVector(ceed_data->elem_restr_x, &ceed_data->x_coord, 264 NULL); 265 266 // -- Copy PETSc vector in CEED vector 267 Vec X_loc; 268 const PetscScalar *X_loc_array; 269 ierr = DMGetCoordinatesLocal(dm, &X_loc); CHKERRQ(ierr); 270 ierr = VecScale(X_loc, problem->dm_scale); CHKERRQ(ierr); 271 ierr = VecGetArrayRead(X_loc, &X_loc_array); CHKERRQ(ierr); 272 CeedVectorSetArray(ceed_data->x_coord, CEED_MEM_HOST, CEED_COPY_VALUES, 273 (PetscScalar *)X_loc_array); 274 ierr = VecRestoreArrayRead(X_loc, &X_loc_array); CHKERRQ(ierr); 275 276 // ----------------------------------------------------------------------------- 277 // CEED vectors 278 // ----------------------------------------------------------------------------- 279 // -- Create CEED vector for geometric data 280 CeedInt num_qpts_vol; 281 PetscInt loc_num_elem_vol; 282 CeedBasisGetNumQuadraturePoints(ceed_data->basis_q, &num_qpts_vol); 283 CeedElemRestrictionGetNumElements(ceed_data->elem_restr_q, &loc_num_elem_vol); 284 CeedVectorCreate(ceed, q_data_size_vol*loc_num_elem_vol*num_qpts_vol, 285 &ceed_data->q_data); 286 287 // ----------------------------------------------------------------------------- 288 // CEED Operators 289 // ----------------------------------------------------------------------------- 290 // -- Create CEED operator for quadrature data 291 CeedOperatorCreate(ceed, ceed_data->qf_setup_vol, NULL, NULL, 292 &ceed_data->op_setup_vol); 293 CeedOperatorSetField(ceed_data->op_setup_vol, "dx", ceed_data->elem_restr_x, 294 ceed_data->basis_x, CEED_VECTOR_ACTIVE); 295 CeedOperatorSetField(ceed_data->op_setup_vol, "weight", 296 CEED_ELEMRESTRICTION_NONE, ceed_data->basis_x, CEED_VECTOR_NONE); 297 CeedOperatorSetField(ceed_data->op_setup_vol, "qdata", 298 ceed_data->elem_restr_qd_i, CEED_BASIS_COLLOCATED, CEED_VECTOR_ACTIVE); 299 300 // -- Create CEED operator for ICs 301 CeedOperatorCreate(ceed, ceed_data->qf_ics, NULL, NULL, &ceed_data->op_ics); 302 CeedOperatorSetField(ceed_data->op_ics, "x", ceed_data->elem_restr_x, 303 ceed_data->basis_xc, CEED_VECTOR_ACTIVE); 304 CeedOperatorSetField(ceed_data->op_ics, "q0", ceed_data->elem_restr_q, 305 CEED_BASIS_COLLOCATED, CEED_VECTOR_ACTIVE); 306 307 // Create CEED operator for RHS 308 if (ceed_data->qf_rhs_vol) { 309 CeedOperator op; 310 CeedOperatorCreate(ceed, ceed_data->qf_rhs_vol, NULL, NULL, &op); 311 CeedOperatorSetField(op, "q", ceed_data->elem_restr_q, ceed_data->basis_q, 312 CEED_VECTOR_ACTIVE); 313 CeedOperatorSetField(op, "dq", ceed_data->elem_restr_q, ceed_data->basis_q, 314 CEED_VECTOR_ACTIVE); 315 CeedOperatorSetField(op, "qdata", ceed_data->elem_restr_qd_i, 316 CEED_BASIS_COLLOCATED, ceed_data->q_data); 317 CeedOperatorSetField(op, "x", ceed_data->elem_restr_x, ceed_data->basis_x, 318 ceed_data->x_coord); 319 CeedOperatorSetField(op, "v", ceed_data->elem_restr_q, ceed_data->basis_q, 320 CEED_VECTOR_ACTIVE); 321 CeedOperatorSetField(op, "dv", ceed_data->elem_restr_q, ceed_data->basis_q, 322 CEED_VECTOR_ACTIVE); 323 user->op_rhs_vol = op; 324 } 325 326 // -- CEED operator for IFunction 327 if (ceed_data->qf_ifunction_vol) { 328 CeedOperator op; 329 CeedOperatorCreate(ceed, ceed_data->qf_ifunction_vol, NULL, NULL, &op); 330 CeedOperatorSetField(op, "q", ceed_data->elem_restr_q, ceed_data->basis_q, 331 CEED_VECTOR_ACTIVE); 332 CeedOperatorSetField(op, "dq", ceed_data->elem_restr_q, ceed_data->basis_q, 333 CEED_VECTOR_ACTIVE); 334 CeedOperatorSetField(op, "q dot", ceed_data->elem_restr_q, ceed_data->basis_q, 335 user->q_dot_ceed); 336 CeedOperatorSetField(op, "qdata", ceed_data->elem_restr_qd_i, 337 CEED_BASIS_COLLOCATED, ceed_data->q_data); 338 CeedOperatorSetField(op, "x", ceed_data->elem_restr_x, ceed_data->basis_x, 339 ceed_data->x_coord); 340 CeedOperatorSetField(op, "v", ceed_data->elem_restr_q, ceed_data->basis_q, 341 CEED_VECTOR_ACTIVE); 342 CeedOperatorSetField(op, "dv", ceed_data->elem_restr_q, ceed_data->basis_q, 343 CEED_VECTOR_ACTIVE); 344 user->op_ifunction_vol = op; 345 } 346 347 // ***************************************************************************** 348 // Set up CEED objects for the exterior domain (surface) 349 // ***************************************************************************** 350 CeedInt height = 1, 351 dim_sur = dim - height, 352 P_sur = app_ctx->degree + 1, 353 Q_sur = P_sur + app_ctx->q_extra; 354 const CeedInt q_data_size_sur = problem->q_data_size_sur; 355 356 // ----------------------------------------------------------------------------- 357 // CEED Bases 358 // ----------------------------------------------------------------------------- 359 CeedBasisCreateTensorH1Lagrange(ceed, dim_sur, num_comp_q, P_sur, Q_sur, 360 CEED_GAUSS, &ceed_data->basis_q_sur); 361 CeedBasisCreateTensorH1Lagrange(ceed, dim_sur, num_comp_x, 2, Q_sur, CEED_GAUSS, 362 &ceed_data->basis_x_sur); 363 364 // ----------------------------------------------------------------------------- 365 // CEED QFunctions 366 // ----------------------------------------------------------------------------- 367 // -- Create QFunction for quadrature data 368 CeedQFunctionCreateInterior(ceed, 1, problem->setup_sur, problem->setup_sur_loc, 369 &ceed_data->qf_setup_sur); 370 CeedQFunctionAddInput(ceed_data->qf_setup_sur, "dx", num_comp_x*dim_sur, 371 CEED_EVAL_GRAD); 372 CeedQFunctionAddInput(ceed_data->qf_setup_sur, "weight", 1, CEED_EVAL_WEIGHT); 373 CeedQFunctionAddOutput(ceed_data->qf_setup_sur, "surface qdata", 374 q_data_size_sur, 375 CEED_EVAL_NONE); 376 377 // -- Creat QFunction for the physics on the boundaries 378 if (problem->apply_sur) { 379 CeedQFunctionCreateInterior(ceed, 1, problem->apply_sur, problem->apply_sur_loc, 380 &ceed_data->qf_apply_sur); 381 CeedQFunctionAddInput(ceed_data->qf_apply_sur, "q", num_comp_q, 382 CEED_EVAL_INTERP); 383 CeedQFunctionAddInput(ceed_data->qf_apply_sur, "surface qdata", q_data_size_sur, 384 CEED_EVAL_NONE); 385 CeedQFunctionAddInput(ceed_data->qf_apply_sur, "x", num_comp_x, 386 CEED_EVAL_INTERP); 387 CeedQFunctionAddOutput(ceed_data->qf_apply_sur, "v", num_comp_q, 388 CEED_EVAL_INTERP); 389 } 390 391 // ***************************************************************************** 392 // CEED Operator Apply 393 // ***************************************************************************** 394 // -- Apply CEED Operator for the geometric data 395 CeedOperatorApply(ceed_data->op_setup_vol, ceed_data->x_coord, 396 ceed_data->q_data, CEED_REQUEST_IMMEDIATE); 397 398 // -- Create and apply CEED Composite Operator for the entire domain 399 if (!user->phys->implicit) { // RHS 400 ierr = CreateOperatorForDomain(ceed, dm, bc, ceed_data, user->phys, 401 user->op_rhs_vol, height, P_sur, Q_sur, 402 q_data_size_sur, &user->op_rhs); CHKERRQ(ierr); 403 } else { // IFunction 404 ierr = CreateOperatorForDomain(ceed, dm, bc, ceed_data, user->phys, 405 user->op_ifunction_vol, height, P_sur, Q_sur, 406 q_data_size_sur, &user->op_ifunction); CHKERRQ(ierr); 407 } 408 409 PetscFunctionReturn(0); 410 } 411