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 P, 29 CeedInt height, DMLabel domain_label, 30 CeedInt value, CeedElemRestriction *elem_restr) { 31 PetscSection section; 32 PetscInt p, num_elem, num_dofs, *elem_restrict, elem_offset, num_fields, 33 dim, depth; 34 DMLabel depth_label; 35 IS depth_IS, iter_IS; 36 Vec U_loc; 37 const PetscInt *iter_indices; 38 PetscErrorCode ierr; 39 PetscFunctionBeginUser; 40 41 ierr = DMGetDimension(dm, &dim); CHKERRQ(ierr); 42 dim -= height; 43 ierr = DMGetLocalSection(dm, §ion); CHKERRQ(ierr); 44 ierr = PetscSectionGetNumFields(section, &num_fields); CHKERRQ(ierr); 45 PetscInt num_comp[num_fields], field_off[num_fields+1]; 46 field_off[0] = 0; 47 for (PetscInt f=0; f<num_fields; f++) { 48 ierr = PetscSectionGetFieldComponents(section, f, &num_comp[f]); CHKERRQ(ierr); 49 field_off[f+1] = field_off[f] + num_comp[f]; 50 } 51 52 ierr = DMPlexGetDepth(dm, &depth); CHKERRQ(ierr); 53 ierr = DMPlexGetDepthLabel(dm, &depth_label); CHKERRQ(ierr); 54 ierr = DMLabelGetStratumIS(depth_label, depth - height, &depth_IS); 55 CHKERRQ(ierr); 56 if (domain_label) { 57 IS domain_IS; 58 ierr = DMLabelGetStratumIS(domain_label, value, &domain_IS); CHKERRQ(ierr); 59 if (domain_IS) { // domain_IS is non-empty 60 ierr = ISIntersect(depth_IS, domain_IS, &iter_IS); CHKERRQ(ierr); 61 ierr = ISDestroy(&domain_IS); CHKERRQ(ierr); 62 } else { // domain_IS is NULL (empty) 63 iter_IS = NULL; 64 } 65 ierr = ISDestroy(&depth_IS); CHKERRQ(ierr); 66 } else { 67 iter_IS = depth_IS; 68 } 69 if (iter_IS) { 70 ierr = ISGetLocalSize(iter_IS, &num_elem); CHKERRQ(ierr); 71 ierr = ISGetIndices(iter_IS, &iter_indices); CHKERRQ(ierr); 72 } else { 73 num_elem = 0; 74 iter_indices = NULL; 75 } 76 ierr = PetscMalloc1(num_elem*PetscPowInt(P, dim), &elem_restrict); 77 CHKERRQ(ierr); 78 for (p=0, elem_offset=0; p<num_elem; p++) { 79 PetscInt c = iter_indices[p]; 80 PetscInt num_indices, *indices, num_nodes; 81 ierr = DMPlexGetClosureIndices(dm, section, section, c, PETSC_TRUE, 82 &num_indices, &indices, NULL, NULL); 83 CHKERRQ(ierr); 84 bool flip = false; 85 if (height > 0) { 86 PetscInt num_cells, num_faces, start = -1; 87 const PetscInt *orients, *faces, *cells; 88 ierr = DMPlexGetSupport(dm, c, &cells); CHKERRQ(ierr); 89 ierr = DMPlexGetSupportSize(dm, c, &num_cells); CHKERRQ(ierr); 90 if (num_cells != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, 91 "Expected one cell in support of exterior face, but got %D cells", 92 num_cells); 93 ierr = DMPlexGetCone(dm, cells[0], &faces); CHKERRQ(ierr); 94 ierr = DMPlexGetConeSize(dm, cells[0], &num_faces); CHKERRQ(ierr); 95 for (PetscInt i=0; i<num_faces; i++) {if (faces[i] == c) start = i;} 96 if (start < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_CORRUPT, 97 "Could not find face %D in cone of its support", 98 c); 99 ierr = DMPlexGetConeOrientation(dm, cells[0], &orients); CHKERRQ(ierr); 100 if (orients[start] < 0) flip = true; 101 } 102 if (num_indices % field_off[num_fields]) SETERRQ1(PETSC_COMM_SELF, 103 PETSC_ERR_ARG_INCOMP, "Number of closure indices not compatible with Cell %D", 104 c); 105 num_nodes = num_indices / field_off[num_fields]; 106 for (PetscInt i=0; i<num_nodes; i++) { 107 PetscInt ii = i; 108 if (flip) { 109 if (P == num_nodes) ii = num_nodes - 1 - i; 110 else if (P*P == num_nodes) { 111 PetscInt row = i / P, col = i % P; 112 ii = row + col * P; 113 } else SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_SUP, 114 "No support for flipping point with %D nodes != P (%D) or P^2", 115 num_nodes, P); 116 } 117 // Check that indices are blocked by node and thus can be coalesced as a single field with 118 // field_off[num_fields] = sum(num_comp) components. 119 for (PetscInt f=0; f<num_fields; f++) { 120 for (PetscInt j=0; j<num_comp[f]; j++) { 121 if (Involute(indices[field_off[f]*num_nodes + ii*num_comp[f] + j]) 122 != Involute(indices[ii*num_comp[0]]) + field_off[f] + j) 123 SETERRQ4(PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, 124 "Cell %D closure indices not interlaced for node %D field %D component %D", 125 c, ii, f, j); 126 } 127 } 128 // Essential boundary conditions are encoded as -(loc+1), but we don't care so we decode. 129 PetscInt loc = Involute(indices[ii*num_comp[0]]); 130 elem_restrict[elem_offset++] = loc; 131 } 132 ierr = DMPlexRestoreClosureIndices(dm, section, section, c, PETSC_TRUE, 133 &num_indices, &indices, NULL, NULL); 134 CHKERRQ(ierr); 135 } 136 if (elem_offset != num_elem*PetscPowInt(P, dim)) 137 SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_LIB, 138 "ElemRestriction of size (%D,%D) initialized %D nodes", num_elem, 139 PetscPowInt(P, dim),elem_offset); 140 if (iter_IS) { 141 ierr = ISRestoreIndices(iter_IS, &iter_indices); CHKERRQ(ierr); 142 } 143 ierr = ISDestroy(&iter_IS); CHKERRQ(ierr); 144 145 ierr = DMGetLocalVector(dm, &U_loc); CHKERRQ(ierr); 146 ierr = VecGetLocalSize(U_loc, &num_dofs); CHKERRQ(ierr); 147 ierr = DMRestoreLocalVector(dm, &U_loc); CHKERRQ(ierr); 148 CeedElemRestrictionCreate(ceed, num_elem, PetscPowInt(P, dim), 149 field_off[num_fields], 150 1, num_dofs, CEED_MEM_HOST, CEED_COPY_VALUES, elem_restrict, 151 elem_restr); 152 ierr = PetscFree(elem_restrict); CHKERRQ(ierr); 153 PetscFunctionReturn(0); 154 } 155 156 // Utility function to get Ceed Restriction for each domain 157 PetscErrorCode GetRestrictionForDomain(Ceed ceed, DM dm, CeedInt height, 158 DMLabel domain_label, PetscInt value, 159 CeedInt P, CeedInt Q, CeedInt q_data_size, 160 CeedElemRestriction *elem_restr_q, 161 CeedElemRestriction *elem_restr_x, 162 CeedElemRestriction *elem_restr_qd_i) { 163 DM dm_coord; 164 CeedInt dim, loc_num_elem; 165 CeedInt Q_dim; 166 PetscErrorCode ierr; 167 PetscFunctionBeginUser; 168 169 ierr = DMGetDimension(dm, &dim); CHKERRQ(ierr); 170 dim -= height; 171 Q_dim = CeedIntPow(Q, dim); 172 ierr = DMGetCoordinateDM(dm, &dm_coord); CHKERRQ(ierr); 173 ierr = DMPlexSetClosurePermutationTensor(dm_coord, PETSC_DETERMINE, NULL); 174 CHKERRQ(ierr); 175 ierr = CreateRestrictionFromPlex(ceed, dm, P, height, domain_label, value, 176 elem_restr_q); CHKERRQ(ierr); 177 ierr = CreateRestrictionFromPlex(ceed, dm_coord, 2, height, domain_label, value, 178 elem_restr_x); CHKERRQ(ierr); 179 CeedElemRestrictionGetNumElements(*elem_restr_q, &loc_num_elem); 180 CeedElemRestrictionCreateStrided(ceed, loc_num_elem, Q_dim, 181 q_data_size, q_data_size*loc_num_elem*Q_dim, 182 CEED_STRIDES_BACKEND, elem_restr_qd_i); 183 PetscFunctionReturn(0); 184 } 185 186 // Utility function to create CEED Composite Operator for the entire domain 187 PetscErrorCode CreateOperatorForDomain(Ceed ceed, DM dm, SimpleBC bc, 188 CeedData ceed_data, Physics phys, 189 CeedOperator op_apply_vol, CeedInt height, 190 CeedInt P_sur, CeedInt Q_sur, CeedInt q_data_size_sur, 191 CeedOperator *op_apply) { 192 CeedInt dim, num_face; 193 DMLabel domain_label; 194 PetscErrorCode ierr; 195 PetscFunctionBeginUser; 196 197 // Create Composite Operaters 198 CeedCompositeOperatorCreate(ceed, op_apply); 199 200 // --Apply Sub-Operator for the volume 201 CeedCompositeOperatorAddSub(*op_apply, op_apply_vol); 202 203 // -- Create Sub-Operator for in/outflow BCs 204 if (phys->has_neumann) { 205 // --- Setup 206 ierr = DMGetLabel(dm, "Face Sets", &domain_label); CHKERRQ(ierr); 207 ierr = DMGetDimension(dm, &dim); CHKERRQ(ierr); 208 if (dim == 2) num_face = 4; 209 if (dim == 3) num_face = 6; 210 211 // --- Get number of quadrature points for the boundaries 212 CeedInt num_qpts_sur; 213 CeedBasisGetNumQuadraturePoints(ceed_data->basis_q_sur, &num_qpts_sur); 214 215 // --- Create Sub-Operator for each face 216 for (CeedInt i=0; i<num_face; i++) { 217 CeedVector q_data_sur; 218 CeedOperator op_setup_sur, op_apply_sur; 219 CeedElemRestriction elem_restr_x_sur, elem_restr_q_sur, 220 elem_restr_qd_i_sur; 221 222 // ---- CEED Restriction 223 ierr = GetRestrictionForDomain(ceed, dm, height, domain_label, i+1, P_sur, 224 Q_sur, q_data_size_sur, &elem_restr_q_sur, 225 &elem_restr_x_sur, &elem_restr_qd_i_sur); 226 CHKERRQ(ierr); 227 228 // ---- CEED Vector 229 PetscInt loc_num_elem_sur; 230 CeedElemRestrictionGetNumElements(elem_restr_q_sur, &loc_num_elem_sur); 231 CeedVectorCreate(ceed, q_data_size_sur*loc_num_elem_sur*num_qpts_sur, 232 &q_data_sur); 233 234 // ---- CEED Operator 235 // ----- CEED Operator for Setup (geometric factors) 236 CeedOperatorCreate(ceed, ceed_data->qf_setup_sur, NULL, NULL, &op_setup_sur); 237 CeedOperatorSetField(op_setup_sur, "dx", elem_restr_x_sur, 238 ceed_data->basis_x_sur, CEED_VECTOR_ACTIVE); 239 CeedOperatorSetField(op_setup_sur, "weight", CEED_ELEMRESTRICTION_NONE, 240 ceed_data->basis_x_sur, CEED_VECTOR_NONE); 241 CeedOperatorSetField(op_setup_sur, "q_data_sur", elem_restr_qd_i_sur, 242 CEED_BASIS_COLLOCATED, CEED_VECTOR_ACTIVE); 243 244 // ----- CEED Operator for Physics 245 CeedOperatorCreate(ceed, ceed_data->qf_apply_sur, NULL, NULL, &op_apply_sur); 246 CeedOperatorSetField(op_apply_sur, "q", elem_restr_q_sur, 247 ceed_data->basis_q_sur, CEED_VECTOR_ACTIVE); 248 CeedOperatorSetField(op_apply_sur, "q_data_sur", elem_restr_qd_i_sur, 249 CEED_BASIS_COLLOCATED, q_data_sur); 250 CeedOperatorSetField(op_apply_sur, "x", elem_restr_x_sur, 251 ceed_data->basis_x_sur, ceed_data->x_coord); 252 CeedOperatorSetField(op_apply_sur, "v", elem_restr_q_sur, 253 ceed_data->basis_q_sur, CEED_VECTOR_ACTIVE); 254 255 // ----- Apply CEED operator for Setup 256 CeedOperatorApply(op_setup_sur, ceed_data->x_coord, q_data_sur, 257 CEED_REQUEST_IMMEDIATE); 258 259 // ----- Apply Sub-Operator for the Boundary 260 CeedCompositeOperatorAddSub(*op_apply, op_apply_sur); 261 262 // ----- Cleanup 263 CeedVectorDestroy(&q_data_sur); 264 CeedElemRestrictionDestroy(&elem_restr_q_sur); 265 CeedElemRestrictionDestroy(&elem_restr_x_sur); 266 CeedElemRestrictionDestroy(&elem_restr_qd_i_sur); 267 CeedOperatorDestroy(&op_setup_sur); 268 CeedOperatorDestroy(&op_apply_sur); 269 } 270 } 271 272 PetscFunctionReturn(0); 273 } 274 275 PetscErrorCode SetupLibceed(Ceed ceed, CeedData ceed_data, DM dm, User user, 276 AppCtx app_ctx, ProblemData *problem, SimpleBC bc) { 277 PetscErrorCode ierr; 278 PetscFunctionBeginUser; 279 280 // ***************************************************************************** 281 // Set up CEED objects for the interior domain (volume) 282 // ***************************************************************************** 283 const PetscInt num_comp_q = 5; 284 const CeedInt dim = problem->dim, 285 num_comp_x = problem->dim, 286 q_data_size_vol = problem->q_data_size_vol, 287 P = app_ctx->degree + 1, 288 Q = P + app_ctx->q_extra; 289 290 // ----------------------------------------------------------------------------- 291 // CEED Bases 292 // ----------------------------------------------------------------------------- 293 CeedBasisCreateTensorH1Lagrange(ceed, dim, num_comp_q, P, Q, CEED_GAUSS, 294 &ceed_data->basis_q); 295 CeedBasisCreateTensorH1Lagrange(ceed, dim, num_comp_x, 2, Q, CEED_GAUSS, 296 &ceed_data->basis_x); 297 CeedBasisCreateTensorH1Lagrange(ceed, dim, num_comp_x, 2, P, 298 CEED_GAUSS_LOBATTO, &ceed_data->basis_xc); 299 300 // ----------------------------------------------------------------------------- 301 // CEED Restrictions 302 // ----------------------------------------------------------------------------- 303 // -- Create restriction 304 ierr = GetRestrictionForDomain(ceed, dm, 0, 0, 0, P, Q, 305 q_data_size_vol, &ceed_data->elem_restr_q, &ceed_data->elem_restr_x, 306 &ceed_data->elem_restr_qd_i); CHKERRQ(ierr); 307 // -- Create E vectors 308 CeedElemRestrictionCreateVector(ceed_data->elem_restr_q, &user->q_ceed, NULL); 309 CeedElemRestrictionCreateVector(ceed_data->elem_restr_q, &user->q_dot_ceed, 310 NULL); 311 CeedElemRestrictionCreateVector(ceed_data->elem_restr_q, &user->g_ceed, NULL); 312 313 // ----------------------------------------------------------------------------- 314 // CEED QFunctions 315 // ----------------------------------------------------------------------------- 316 // -- Create QFunction for quadrature data 317 CeedQFunctionCreateInterior(ceed, 1, problem->setup_vol, problem->setup_vol_loc, 318 &ceed_data->qf_setup_vol); 319 CeedQFunctionAddInput(ceed_data->qf_setup_vol, "dx", num_comp_x*dim, 320 CEED_EVAL_GRAD); 321 CeedQFunctionAddInput(ceed_data->qf_setup_vol, "weight", 1, CEED_EVAL_WEIGHT); 322 CeedQFunctionAddOutput(ceed_data->qf_setup_vol, "q_data", q_data_size_vol, 323 CEED_EVAL_NONE); 324 325 // -- Create QFunction for ICs 326 CeedQFunctionCreateInterior(ceed, 1, problem->ics, problem->ics_loc, 327 &ceed_data->qf_ics); 328 CeedQFunctionAddInput(ceed_data->qf_ics, "x", num_comp_x, CEED_EVAL_INTERP); 329 CeedQFunctionAddOutput(ceed_data->qf_ics, "q0", num_comp_q, CEED_EVAL_NONE); 330 331 // -- Create QFunction for RHS 332 if (problem->apply_vol_rhs) { 333 CeedQFunctionCreateInterior(ceed, 1, problem->apply_vol_rhs, 334 problem->apply_vol_rhs_loc, &ceed_data->qf_rhs_vol); 335 CeedQFunctionAddInput(ceed_data->qf_rhs_vol, "q", num_comp_q, CEED_EVAL_INTERP); 336 CeedQFunctionAddInput(ceed_data->qf_rhs_vol, "dq", num_comp_q*dim, 337 CEED_EVAL_GRAD); 338 CeedQFunctionAddInput(ceed_data->qf_rhs_vol, "q_data", q_data_size_vol, 339 CEED_EVAL_NONE); 340 CeedQFunctionAddInput(ceed_data->qf_rhs_vol, "x", num_comp_x, CEED_EVAL_INTERP); 341 CeedQFunctionAddOutput(ceed_data->qf_rhs_vol, "v", num_comp_q, 342 CEED_EVAL_INTERP); 343 CeedQFunctionAddOutput(ceed_data->qf_rhs_vol, "dv", num_comp_q*dim, 344 CEED_EVAL_GRAD); 345 } 346 347 // -- Create QFunction for IFunction 348 if (problem->apply_vol_ifunction) { 349 CeedQFunctionCreateInterior(ceed, 1, problem->apply_vol_ifunction, 350 problem->apply_vol_ifunction_loc, &ceed_data->qf_ifunction_vol); 351 CeedQFunctionAddInput(ceed_data->qf_ifunction_vol, "q", num_comp_q, 352 CEED_EVAL_INTERP); 353 CeedQFunctionAddInput(ceed_data->qf_ifunction_vol, "dq", num_comp_q*dim, 354 CEED_EVAL_GRAD); 355 CeedQFunctionAddInput(ceed_data->qf_ifunction_vol, "q_dot", num_comp_q, 356 CEED_EVAL_INTERP); 357 CeedQFunctionAddInput(ceed_data->qf_ifunction_vol, "q_data", q_data_size_vol, 358 CEED_EVAL_NONE); 359 CeedQFunctionAddInput(ceed_data->qf_ifunction_vol, "x", num_comp_x, 360 CEED_EVAL_INTERP); 361 CeedQFunctionAddOutput(ceed_data->qf_ifunction_vol, "v", num_comp_q, 362 CEED_EVAL_INTERP); 363 CeedQFunctionAddOutput(ceed_data->qf_ifunction_vol, "dv", num_comp_q*dim, 364 CEED_EVAL_GRAD); 365 } 366 367 // --------------------------------------------------------------------------- 368 // Element coordinates 369 // --------------------------------------------------------------------------- 370 // -- Create CEED vector 371 CeedElemRestrictionCreateVector(ceed_data->elem_restr_x, &ceed_data->x_coord, 372 NULL); 373 374 // -- Copy PETSc vector in CEED vector 375 Vec X_loc; 376 const PetscScalar *X_loc_array; 377 ierr = DMGetCoordinatesLocal(dm, &X_loc); CHKERRQ(ierr); 378 ierr = VecGetArrayRead(X_loc, &X_loc_array); CHKERRQ(ierr); 379 CeedVectorSetArray(ceed_data->x_coord, CEED_MEM_HOST, CEED_COPY_VALUES, 380 (PetscScalar *)X_loc_array); 381 ierr = VecRestoreArrayRead(X_loc, &X_loc_array); CHKERRQ(ierr); 382 383 // ----------------------------------------------------------------------------- 384 // CEED vectors 385 // ----------------------------------------------------------------------------- 386 // -- Create CEED vector for geometric data 387 CeedInt num_qpts_vol; 388 PetscInt loc_num_elem_vol; 389 CeedBasisGetNumQuadraturePoints(ceed_data->basis_q, &num_qpts_vol); 390 CeedElemRestrictionGetNumElements(ceed_data->elem_restr_q, &loc_num_elem_vol); 391 CeedVectorCreate(ceed, q_data_size_vol*loc_num_elem_vol*num_qpts_vol, 392 &ceed_data->q_data); 393 394 // ----------------------------------------------------------------------------- 395 // CEED Operators 396 // ----------------------------------------------------------------------------- 397 // -- Create CEED operator for quadrature data 398 CeedOperatorCreate(ceed, ceed_data->qf_setup_vol, NULL, NULL, 399 &ceed_data->op_setup_vol); 400 CeedOperatorSetField(ceed_data->op_setup_vol, "dx", ceed_data->elem_restr_x, 401 ceed_data->basis_x, CEED_VECTOR_ACTIVE); 402 CeedOperatorSetField(ceed_data->op_setup_vol, "weight", 403 CEED_ELEMRESTRICTION_NONE, ceed_data->basis_x, CEED_VECTOR_NONE); 404 CeedOperatorSetField(ceed_data->op_setup_vol, "q_data", 405 ceed_data->elem_restr_qd_i, CEED_BASIS_COLLOCATED, CEED_VECTOR_ACTIVE); 406 407 // -- Create CEED operator for ICs 408 CeedOperatorCreate(ceed, ceed_data->qf_ics, NULL, NULL, &ceed_data->op_ics); 409 CeedOperatorSetField(ceed_data->op_ics, "x", ceed_data->elem_restr_x, 410 ceed_data->basis_xc, CEED_VECTOR_ACTIVE); 411 CeedOperatorSetField(ceed_data->op_ics, "q0", ceed_data->elem_restr_q, 412 CEED_BASIS_COLLOCATED, CEED_VECTOR_ACTIVE); 413 414 // Create CEED operator for RHS 415 if (ceed_data->qf_rhs_vol) { 416 CeedOperator op; 417 CeedOperatorCreate(ceed, ceed_data->qf_rhs_vol, NULL, NULL, &op); 418 CeedOperatorSetField(op, "q", ceed_data->elem_restr_q, ceed_data->basis_q, 419 CEED_VECTOR_ACTIVE); 420 CeedOperatorSetField(op, "dq", ceed_data->elem_restr_q, ceed_data->basis_q, 421 CEED_VECTOR_ACTIVE); 422 CeedOperatorSetField(op, "q_data", ceed_data->elem_restr_qd_i, 423 CEED_BASIS_COLLOCATED, ceed_data->q_data); 424 CeedOperatorSetField(op, "x", ceed_data->elem_restr_x, ceed_data->basis_x, 425 ceed_data->x_coord); 426 CeedOperatorSetField(op, "v", ceed_data->elem_restr_q, ceed_data->basis_q, 427 CEED_VECTOR_ACTIVE); 428 CeedOperatorSetField(op, "dv", ceed_data->elem_restr_q, ceed_data->basis_q, 429 CEED_VECTOR_ACTIVE); 430 user->op_rhs_vol = op; 431 } 432 433 // -- CEED operator for IFunction 434 if (ceed_data->qf_ifunction_vol) { 435 CeedOperator op; 436 CeedOperatorCreate(ceed, ceed_data->qf_ifunction_vol, NULL, NULL, &op); 437 CeedOperatorSetField(op, "q", ceed_data->elem_restr_q, ceed_data->basis_q, 438 CEED_VECTOR_ACTIVE); 439 CeedOperatorSetField(op, "dq", ceed_data->elem_restr_q, ceed_data->basis_q, 440 CEED_VECTOR_ACTIVE); 441 CeedOperatorSetField(op, "q_dot", ceed_data->elem_restr_q, ceed_data->basis_q, 442 user->q_dot_ceed); 443 CeedOperatorSetField(op, "q_data", ceed_data->elem_restr_qd_i, 444 CEED_BASIS_COLLOCATED, ceed_data->q_data); 445 CeedOperatorSetField(op, "x", ceed_data->elem_restr_x, ceed_data->basis_x, 446 ceed_data->x_coord); 447 CeedOperatorSetField(op, "v", ceed_data->elem_restr_q, ceed_data->basis_q, 448 CEED_VECTOR_ACTIVE); 449 CeedOperatorSetField(op, "dv", ceed_data->elem_restr_q, ceed_data->basis_q, 450 CEED_VECTOR_ACTIVE); 451 user->op_ifunction_vol = op; 452 } 453 454 // ***************************************************************************** 455 // Set up CEED objects for the exterior domain (surface) 456 // ***************************************************************************** 457 CeedInt height = 1, 458 dim_sur = dim - height, 459 P_sur = app_ctx->degree + 1, 460 Q_sur = P_sur + app_ctx->q_extra; 461 const CeedInt q_data_size_sur = problem->q_data_size_sur; 462 463 // ----------------------------------------------------------------------------- 464 // CEED Bases 465 // ----------------------------------------------------------------------------- 466 CeedBasisCreateTensorH1Lagrange(ceed, dim_sur, num_comp_q, P_sur, Q_sur, 467 CEED_GAUSS, &ceed_data->basis_q_sur); 468 CeedBasisCreateTensorH1Lagrange(ceed, dim_sur, num_comp_x, 2, Q_sur, CEED_GAUSS, 469 &ceed_data->basis_x_sur); 470 471 // ----------------------------------------------------------------------------- 472 // CEED QFunctions 473 // ----------------------------------------------------------------------------- 474 // -- Create QFunction for quadrature data 475 CeedQFunctionCreateInterior(ceed, 1, problem->setup_sur, problem->setup_sur_loc, 476 &ceed_data->qf_setup_sur); 477 CeedQFunctionAddInput(ceed_data->qf_setup_sur, "dx", num_comp_x*dim_sur, 478 CEED_EVAL_GRAD); 479 CeedQFunctionAddInput(ceed_data->qf_setup_sur, "weight", 1, CEED_EVAL_WEIGHT); 480 CeedQFunctionAddOutput(ceed_data->qf_setup_sur, "q_data_sur", q_data_size_sur, 481 CEED_EVAL_NONE); 482 483 // -- Creat QFunction for the physics on the boundaries 484 if (problem->apply_sur) { 485 CeedQFunctionCreateInterior(ceed, 1, problem->apply_sur, problem->apply_sur_loc, 486 &ceed_data->qf_apply_sur); 487 CeedQFunctionAddInput(ceed_data->qf_apply_sur, "q", num_comp_q, 488 CEED_EVAL_INTERP); 489 CeedQFunctionAddInput(ceed_data->qf_apply_sur, "q_data_sur", q_data_size_sur, 490 CEED_EVAL_NONE); 491 CeedQFunctionAddInput(ceed_data->qf_apply_sur, "x", num_comp_x, 492 CEED_EVAL_INTERP); 493 CeedQFunctionAddOutput(ceed_data->qf_apply_sur, "v", num_comp_q, 494 CEED_EVAL_INTERP); 495 } 496 497 // ***************************************************************************** 498 // CEED Operator Apply 499 // ***************************************************************************** 500 // -- Apply CEED Operator for the geometric data 501 CeedOperatorApply(ceed_data->op_setup_vol, ceed_data->x_coord, 502 ceed_data->q_data, CEED_REQUEST_IMMEDIATE); 503 504 // -- Create and apply CEED Composite Operator for the entire domain 505 if (!user->phys->implicit) { // RHS 506 ierr = CreateOperatorForDomain(ceed, dm, bc, ceed_data, user->phys, 507 user->op_rhs_vol, height, P_sur, Q_sur, 508 q_data_size_sur, &user->op_rhs); CHKERRQ(ierr); 509 } else { // IFunction 510 ierr = CreateOperatorForDomain(ceed, dm, bc, ceed_data, user->phys, 511 user->op_ifunction_vol, height, P_sur, Q_sur, 512 q_data_size_sur, &user->op_ifunction); CHKERRQ(ierr); 513 } 514 515 PetscFunctionReturn(0); 516 } 517