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 164 DM dm_coord; 165 CeedInt dim, loc_num_elem; 166 CeedInt Q_dim; 167 PetscErrorCode ierr; 168 PetscFunctionBeginUser; 169 170 ierr = DMGetDimension(dm, &dim); CHKERRQ(ierr); 171 dim -= height; 172 Q_dim = CeedIntPow(Q, dim); 173 ierr = DMGetCoordinateDM(dm, &dm_coord); CHKERRQ(ierr); 174 ierr = DMPlexSetClosurePermutationTensor(dm_coord, PETSC_DETERMINE, NULL); 175 CHKERRQ(ierr); 176 ierr = CreateRestrictionFromPlex(ceed, dm, P, height, domain_label, value, 177 elem_restr_q); 178 CHKERRQ(ierr); 179 ierr = CreateRestrictionFromPlex(ceed, dm_coord, 2, height, domain_label, value, 180 elem_restr_x); 181 CHKERRQ(ierr); 182 CeedElemRestrictionGetNumElements(*elem_restr_q, &loc_num_elem); 183 CeedElemRestrictionCreateStrided(ceed, loc_num_elem, Q_dim, 184 q_data_size, q_data_size*loc_num_elem*Q_dim, 185 CEED_STRIDES_BACKEND, elem_restr_qd_i); 186 PetscFunctionReturn(0); 187 } 188 189 // Utility function to create CEED Composite Operator for the entire domain 190 PetscErrorCode CreateOperatorForDomain(Ceed ceed, DM dm, SimpleBC bc, 191 CeedData ceed_data, Physics phys, 192 CeedOperator op_apply_vol, CeedInt height, 193 CeedInt P_sur, CeedInt Q_sur, CeedInt q_data_size_sur, 194 CeedOperator *op_apply) { 195 CeedInt dim, num_face; 196 DMLabel domain_label; 197 PetscErrorCode ierr; 198 PetscFunctionBeginUser; 199 200 // Create Composite Operaters 201 CeedCompositeOperatorCreate(ceed, op_apply); 202 203 // --Apply Sub-Operator for the volume 204 CeedCompositeOperatorAddSub(*op_apply, op_apply_vol); 205 206 // -- Create Sub-Operator for in/outflow BCs 207 if (phys->has_neumann == PETSC_TRUE) { 208 // --- Setup 209 ierr = DMGetLabel(dm, "Face Sets", &domain_label); CHKERRQ(ierr); 210 ierr = DMGetDimension(dm, &dim); CHKERRQ(ierr); 211 if (dim == 2) num_face = 4; 212 if (dim == 3) num_face = 6; 213 214 // --- Get number of quadrature points for the boundaries 215 CeedInt num_qpts_sur; 216 CeedBasisGetNumQuadraturePoints(ceed_data->basis_q_sur, &num_qpts_sur); 217 218 // --- Create Sub-Operator for each face 219 for (CeedInt i=0; i<num_face; i++) { 220 CeedVector q_data_sur; 221 CeedOperator op_setup_sur, op_apply_sur; 222 CeedElemRestriction elem_restr_x_sur, elem_restr_q_sur, 223 elem_restr_qd_i_sur; 224 225 // ---- CEED Restriction 226 ierr = GetRestrictionForDomain(ceed, dm, height, domain_label, i+1, P_sur, 227 Q_sur, q_data_size_sur, &elem_restr_q_sur, 228 &elem_restr_x_sur, &elem_restr_qd_i_sur); 229 CHKERRQ(ierr); 230 231 // ---- CEED Vector 232 PetscInt loc_num_elem_sur; 233 CeedElemRestrictionGetNumElements(elem_restr_q_sur, &loc_num_elem_sur); 234 CeedVectorCreate(ceed, q_data_size_sur*loc_num_elem_sur*num_qpts_sur, 235 &q_data_sur); 236 237 // ---- CEED Operator 238 // ----- CEED Operator for Setup (geometric factors) 239 CeedOperatorCreate(ceed, ceed_data->qf_setup_sur, NULL, NULL, &op_setup_sur); 240 CeedOperatorSetField(op_setup_sur, "dx", elem_restr_x_sur, 241 ceed_data->basis_x_sur, 242 CEED_VECTOR_ACTIVE); 243 CeedOperatorSetField(op_setup_sur, "weight", CEED_ELEMRESTRICTION_NONE, 244 ceed_data->basis_x_sur, CEED_VECTOR_NONE); 245 CeedOperatorSetField(op_setup_sur, "q_data_sur", elem_restr_qd_i_sur, 246 CEED_BASIS_COLLOCATED, CEED_VECTOR_ACTIVE); 247 248 // ----- CEED Operator for Physics 249 CeedOperatorCreate(ceed, ceed_data->qf_apply_sur, NULL, NULL, &op_apply_sur); 250 CeedOperatorSetField(op_apply_sur, "q", elem_restr_q_sur, 251 ceed_data->basis_q_sur, CEED_VECTOR_ACTIVE); 252 CeedOperatorSetField(op_apply_sur, "q_data_sur", elem_restr_qd_i_sur, 253 CEED_BASIS_COLLOCATED, q_data_sur); 254 CeedOperatorSetField(op_apply_sur, "x", elem_restr_x_sur, 255 ceed_data->basis_x_sur, ceed_data->x_coord); 256 CeedOperatorSetField(op_apply_sur, "v", elem_restr_q_sur, 257 ceed_data->basis_q_sur, CEED_VECTOR_ACTIVE); 258 259 // ----- Apply CEED operator for Setup 260 CeedOperatorApply(op_setup_sur, ceed_data->x_coord, q_data_sur, 261 CEED_REQUEST_IMMEDIATE); 262 263 // ----- Apply Sub-Operator for the Boundary 264 CeedCompositeOperatorAddSub(*op_apply, op_apply_sur); 265 266 // ----- Cleanup 267 CeedVectorDestroy(&q_data_sur); 268 CeedElemRestrictionDestroy(&elem_restr_q_sur); 269 CeedElemRestrictionDestroy(&elem_restr_x_sur); 270 CeedElemRestrictionDestroy(&elem_restr_qd_i_sur); 271 CeedOperatorDestroy(&op_setup_sur); 272 CeedOperatorDestroy(&op_apply_sur); 273 } 274 } 275 276 PetscFunctionReturn(0); 277 } 278 279 PetscErrorCode SetupLibceed(Ceed ceed, CeedData ceed_data, DM dm, User user, 280 AppCtx app_ctx, ProblemData *problem, SimpleBC bc) { 281 PetscErrorCode ierr; 282 PetscFunctionBeginUser; 283 284 // ***************************************************************************** 285 // Set up CEED objects for the interior domain (volume) 286 // ***************************************************************************** 287 const PetscInt num_comp_q = 5; 288 const CeedInt dim = problem->dim, 289 num_comp_x = problem->dim, 290 q_data_size_vol = problem->q_data_size_vol, 291 P = app_ctx->degree + 1, 292 Q = P + app_ctx->q_extra; 293 294 // ----------------------------------------------------------------------------- 295 // CEED Bases 296 // ----------------------------------------------------------------------------- 297 CeedBasisCreateTensorH1Lagrange(ceed, dim, num_comp_q, P, Q, CEED_GAUSS, 298 &ceed_data->basis_q); 299 300 CeedBasisCreateTensorH1Lagrange(ceed, dim, num_comp_x, 2, Q, CEED_GAUSS, 301 &ceed_data->basis_x); 302 303 CeedBasisCreateTensorH1Lagrange(ceed, dim, num_comp_x, 2, P, 304 CEED_GAUSS_LOBATTO, &ceed_data->basis_xc); 305 306 // ----------------------------------------------------------------------------- 307 // CEED Restrictions 308 // ----------------------------------------------------------------------------- 309 // -- Create restriction 310 ierr = GetRestrictionForDomain(ceed, dm, 0, 0, 0, P, Q, 311 q_data_size_vol, &ceed_data->elem_restr_q, &ceed_data->elem_restr_x, 312 &ceed_data->elem_restr_qd_i); CHKERRQ(ierr); 313 // -- Create E vectors 314 CeedElemRestrictionCreateVector(ceed_data->elem_restr_q, &user->q_ceed, NULL); 315 CeedElemRestrictionCreateVector(ceed_data->elem_restr_q, &user->q_dot_ceed, 316 NULL); 317 CeedElemRestrictionCreateVector(ceed_data->elem_restr_q, &user->g_ceed, NULL); 318 319 // ----------------------------------------------------------------------------- 320 // CEED QFunctions 321 // ----------------------------------------------------------------------------- 322 // -- Create QFunction for quadrature data 323 CeedQFunctionCreateInterior(ceed, 1, problem->setup_vol, problem->setup_vol_loc, 324 &ceed_data->qf_setup_vol); 325 CeedQFunctionAddInput(ceed_data->qf_setup_vol, "dx", num_comp_x*dim, 326 CEED_EVAL_GRAD); 327 CeedQFunctionAddInput(ceed_data->qf_setup_vol, "weight", 1, CEED_EVAL_WEIGHT); 328 CeedQFunctionAddOutput(ceed_data->qf_setup_vol, "q_data", q_data_size_vol, 329 CEED_EVAL_NONE); 330 331 // -- Create QFunction for ICs 332 CeedQFunctionCreateInterior(ceed, 1, problem->ics, problem->ics_loc, 333 &ceed_data->qf_ics); 334 CeedQFunctionAddInput(ceed_data->qf_ics, "x", num_comp_x, CEED_EVAL_INTERP); 335 CeedQFunctionAddOutput(ceed_data->qf_ics, "q0", num_comp_q, CEED_EVAL_NONE); 336 337 // -- Create QFunction for RHS 338 if (problem->apply_vol_rhs) { 339 CeedQFunctionCreateInterior(ceed, 1, problem->apply_vol_rhs, 340 problem->apply_vol_rhs_loc, &ceed_data->qf_rhs_vol); 341 CeedQFunctionAddInput(ceed_data->qf_rhs_vol, "q", num_comp_q, CEED_EVAL_INTERP); 342 CeedQFunctionAddInput(ceed_data->qf_rhs_vol, "dq", num_comp_q*dim, 343 CEED_EVAL_GRAD); 344 CeedQFunctionAddInput(ceed_data->qf_rhs_vol, "q_data", q_data_size_vol, 345 CEED_EVAL_NONE); 346 CeedQFunctionAddInput(ceed_data->qf_rhs_vol, "x", num_comp_x, CEED_EVAL_INTERP); 347 CeedQFunctionAddOutput(ceed_data->qf_rhs_vol, "v", num_comp_q, 348 CEED_EVAL_INTERP); 349 CeedQFunctionAddOutput(ceed_data->qf_rhs_vol, "dv", num_comp_q*dim, 350 CEED_EVAL_GRAD); 351 } 352 353 // -- Create QFunction for IFunction 354 if (problem->apply_vol_ifunction) { 355 CeedQFunctionCreateInterior(ceed, 1, problem->apply_vol_ifunction, 356 problem->apply_vol_ifunction_loc, &ceed_data->qf_ifunction_vol); 357 CeedQFunctionAddInput(ceed_data->qf_ifunction_vol, "q", num_comp_q, 358 CEED_EVAL_INTERP); 359 CeedQFunctionAddInput(ceed_data->qf_ifunction_vol, "dq", num_comp_q*dim, 360 CEED_EVAL_GRAD); 361 CeedQFunctionAddInput(ceed_data->qf_ifunction_vol, "q_dot", num_comp_q, 362 CEED_EVAL_INTERP); 363 CeedQFunctionAddInput(ceed_data->qf_ifunction_vol, "q_data", q_data_size_vol, 364 CEED_EVAL_NONE); 365 CeedQFunctionAddInput(ceed_data->qf_ifunction_vol, "x", num_comp_x, 366 CEED_EVAL_INTERP); 367 CeedQFunctionAddOutput(ceed_data->qf_ifunction_vol, "v", num_comp_q, 368 CEED_EVAL_INTERP); 369 CeedQFunctionAddOutput(ceed_data->qf_ifunction_vol, "dv", num_comp_q*dim, 370 CEED_EVAL_GRAD); 371 } 372 373 // --------------------------------------------------------------------------- 374 // Element coordinates 375 // --------------------------------------------------------------------------- 376 // -- Create CEED vector 377 CeedElemRestrictionCreateVector(ceed_data->elem_restr_x, &ceed_data->x_coord, 378 NULL); 379 380 // -- Copy PETSc vector in CEED vector 381 Vec X_loc; 382 const PetscScalar *X_loc_array; 383 ierr = DMGetCoordinatesLocal(dm, &X_loc); CHKERRQ(ierr); 384 ierr = VecGetArrayRead(X_loc, &X_loc_array); CHKERRQ(ierr); 385 CeedVectorSetArray(ceed_data->x_coord, CEED_MEM_HOST, CEED_COPY_VALUES, 386 (PetscScalar *)X_loc_array); 387 ierr = VecRestoreArrayRead(X_loc, &X_loc_array); CHKERRQ(ierr); 388 389 // ----------------------------------------------------------------------------- 390 // CEED vectors 391 // ----------------------------------------------------------------------------- 392 // -- Create CEED vector for geometric data 393 CeedInt num_qpts_vol; 394 PetscInt loc_num_elem_vol; 395 CeedBasisGetNumQuadraturePoints(ceed_data->basis_q, &num_qpts_vol); 396 CeedElemRestrictionGetNumElements(ceed_data->elem_restr_q, &loc_num_elem_vol); 397 CeedVectorCreate(ceed, q_data_size_vol*loc_num_elem_vol*num_qpts_vol, 398 &ceed_data->q_data); 399 400 // ----------------------------------------------------------------------------- 401 // CEED Operators 402 // ----------------------------------------------------------------------------- 403 // -- Create CEED operator for quadrature data 404 CeedOperatorCreate(ceed, ceed_data->qf_setup_vol, NULL, NULL, 405 &ceed_data->op_setup_vol); 406 CeedOperatorSetField(ceed_data->op_setup_vol, "dx", ceed_data->elem_restr_x, 407 ceed_data->basis_x, CEED_VECTOR_ACTIVE); 408 CeedOperatorSetField(ceed_data->op_setup_vol, "weight", 409 CEED_ELEMRESTRICTION_NONE, 410 ceed_data->basis_x, CEED_VECTOR_NONE); 411 CeedOperatorSetField(ceed_data->op_setup_vol, "q_data", 412 ceed_data->elem_restr_qd_i, 413 CEED_BASIS_COLLOCATED, CEED_VECTOR_ACTIVE); 414 415 // -- Create CEED operator for ICs 416 CeedOperatorCreate(ceed, ceed_data->qf_ics, NULL, NULL, &ceed_data->op_ics); 417 CeedOperatorSetField(ceed_data->op_ics, "x", ceed_data->elem_restr_x, 418 ceed_data->basis_xc, CEED_VECTOR_ACTIVE); 419 CeedOperatorSetField(ceed_data->op_ics, "q0", ceed_data->elem_restr_q, 420 CEED_BASIS_COLLOCATED, CEED_VECTOR_ACTIVE); 421 422 // Create CEED operator for RHS 423 if (ceed_data->qf_rhs_vol) { 424 CeedOperator op; 425 CeedOperatorCreate(ceed, ceed_data->qf_rhs_vol, NULL, NULL, &op); 426 CeedOperatorSetField(op, "q", ceed_data->elem_restr_q, ceed_data->basis_q, 427 CEED_VECTOR_ACTIVE); 428 CeedOperatorSetField(op, "dq", ceed_data->elem_restr_q, ceed_data->basis_q, 429 CEED_VECTOR_ACTIVE); 430 CeedOperatorSetField(op, "q_data", ceed_data->elem_restr_qd_i, 431 CEED_BASIS_COLLOCATED, 432 ceed_data->q_data); 433 CeedOperatorSetField(op, "x", ceed_data->elem_restr_x, ceed_data->basis_x, 434 ceed_data->x_coord); 435 CeedOperatorSetField(op, "v", ceed_data->elem_restr_q, ceed_data->basis_q, 436 CEED_VECTOR_ACTIVE); 437 CeedOperatorSetField(op, "dv", ceed_data->elem_restr_q, ceed_data->basis_q, 438 CEED_VECTOR_ACTIVE); 439 user->op_rhs_vol = op; 440 } 441 442 // -- CEED operator for IFunction 443 if (ceed_data->qf_ifunction_vol) { 444 CeedOperator op; 445 CeedOperatorCreate(ceed, ceed_data->qf_ifunction_vol, NULL, NULL, &op); 446 CeedOperatorSetField(op, "q", ceed_data->elem_restr_q, ceed_data->basis_q, 447 CEED_VECTOR_ACTIVE); 448 CeedOperatorSetField(op, "dq", ceed_data->elem_restr_q, ceed_data->basis_q, 449 CEED_VECTOR_ACTIVE); 450 CeedOperatorSetField(op, "q_dot", ceed_data->elem_restr_q, ceed_data->basis_q, 451 user->q_dot_ceed); 452 CeedOperatorSetField(op, "q_data", ceed_data->elem_restr_qd_i, 453 CEED_BASIS_COLLOCATED, 454 ceed_data->q_data); 455 CeedOperatorSetField(op, "x", ceed_data->elem_restr_x, ceed_data->basis_x, 456 ceed_data->x_coord); 457 CeedOperatorSetField(op, "v", ceed_data->elem_restr_q, ceed_data->basis_q, 458 CEED_VECTOR_ACTIVE); 459 CeedOperatorSetField(op, "dv", ceed_data->elem_restr_q, ceed_data->basis_q, 460 CEED_VECTOR_ACTIVE); 461 user->op_ifunction_vol = op; 462 } 463 464 // ***************************************************************************** 465 // Set up CEED objects for the exterior domain (surface) 466 // ***************************************************************************** 467 CeedInt height = 1, 468 dim_sur = dim - height, 469 P_sur = app_ctx->degree + 1, 470 Q_sur = P_sur + app_ctx->q_extra; 471 const CeedInt q_data_size_sur = problem->q_data_size_sur; 472 473 // ----------------------------------------------------------------------------- 474 // CEED Bases 475 // ----------------------------------------------------------------------------- 476 CeedBasisCreateTensorH1Lagrange(ceed, dim_sur, num_comp_q, P_sur, Q_sur, 477 CEED_GAUSS, &ceed_data->basis_q_sur); 478 479 CeedBasisCreateTensorH1Lagrange(ceed, dim_sur, num_comp_x, 2, Q_sur, CEED_GAUSS, 480 &ceed_data->basis_x_sur); 481 482 CeedBasisCreateTensorH1Lagrange(ceed, dim_sur, num_comp_x, 2, P_sur, 483 CEED_GAUSS_LOBATTO, &ceed_data->basis_xc_sur); 484 485 // ----------------------------------------------------------------------------- 486 // CEED QFunctions 487 // ----------------------------------------------------------------------------- 488 // -- Create QFunction for quadrature data 489 CeedQFunctionCreateInterior(ceed, 1, problem->setup_sur, problem->setup_sur_loc, 490 &ceed_data->qf_setup_sur); 491 CeedQFunctionAddInput(ceed_data->qf_setup_sur, "dx", num_comp_x*dim_sur, 492 CEED_EVAL_GRAD); 493 CeedQFunctionAddInput(ceed_data->qf_setup_sur, "weight", 1, CEED_EVAL_WEIGHT); 494 CeedQFunctionAddOutput(ceed_data->qf_setup_sur, "q_data_sur", q_data_size_sur, 495 CEED_EVAL_NONE); 496 497 // -- Creat QFunction for the physics on the boundaries 498 if (problem->apply_sur) { 499 CeedQFunctionCreateInterior(ceed, 1, problem->apply_sur, problem->apply_sur_loc, 500 &ceed_data->qf_apply_sur); 501 CeedQFunctionAddInput(ceed_data->qf_apply_sur, "q", num_comp_q, 502 CEED_EVAL_INTERP); 503 CeedQFunctionAddInput(ceed_data->qf_apply_sur, "q_data_sur", q_data_size_sur, 504 CEED_EVAL_NONE); 505 CeedQFunctionAddInput(ceed_data->qf_apply_sur, "x", num_comp_x, 506 CEED_EVAL_INTERP); 507 CeedQFunctionAddOutput(ceed_data->qf_apply_sur, "v", num_comp_q, 508 CEED_EVAL_INTERP); 509 } 510 511 // ***************************************************************************** 512 // CEED Operator Apply 513 // ***************************************************************************** 514 // -- Apply CEED Operator for the geometric data 515 CeedOperatorApply(ceed_data->op_setup_vol, ceed_data->x_coord, 516 ceed_data->q_data, CEED_REQUEST_IMMEDIATE); 517 518 // -- Create and apply CEED Composite Operator for the entire domain 519 if (!user->phys->implicit) { // RHS 520 ierr = CreateOperatorForDomain(ceed, dm, bc, ceed_data, user->phys, 521 user->op_rhs_vol, height, P_sur, Q_sur, 522 q_data_size_sur, &user->op_rhs); CHKERRQ(ierr); 523 } else { // IFunction 524 ierr = CreateOperatorForDomain(ceed, dm, bc, ceed_data, user->phys, 525 user->op_ifunction_vol, height, P_sur, Q_sur, 526 q_data_size_sur, &user->op_ifunction); CHKERRQ(ierr); 527 } 528 529 PetscFunctionReturn(0); 530 } 531 532 // Set up contex for QFunctions 533 PetscErrorCode SetupContextForProblems(Ceed ceed, CeedData ceed_data, 534 AppCtx app_ctx, SetupContext setup_ctx, Physics phys) { 535 PetscFunctionBeginUser; 536 537 // ICs 538 CeedQFunctionContextCreate(ceed, &ceed_data->setup_context); 539 CeedQFunctionContextSetData(ceed_data->setup_context, CEED_MEM_HOST, 540 CEED_USE_POINTER, 541 sizeof(*setup_ctx), setup_ctx); 542 if (ceed_data->qf_ics && strcmp(app_ctx->problem_name, "euler_vortex") != 0) 543 CeedQFunctionSetContext(ceed_data->qf_ics, ceed_data->setup_context); 544 545 // DENSITY_CURRENT 546 if (strcmp(app_ctx->problem_name, "density_current") == 0) { 547 CeedQFunctionContextCreate(ceed, &ceed_data->dc_context); 548 CeedQFunctionContextSetData(ceed_data->dc_context, CEED_MEM_HOST, 549 CEED_USE_POINTER, 550 sizeof(*phys->dc_ctx), phys->dc_ctx); 551 if (ceed_data->qf_rhs_vol) 552 CeedQFunctionSetContext(ceed_data->qf_rhs_vol, ceed_data->dc_context); 553 if (ceed_data->qf_ifunction_vol) 554 CeedQFunctionSetContext(ceed_data->qf_ifunction_vol, ceed_data->dc_context); 555 556 // EULER_VORTEX 557 } else if (strcmp(app_ctx->problem_name, "euler_vortex") == 0) { 558 CeedQFunctionContextCreate(ceed, &ceed_data->euler_context); 559 CeedQFunctionContextSetData(ceed_data->euler_context, CEED_MEM_HOST, 560 CEED_USE_POINTER, 561 sizeof(*phys->euler_ctx), phys->euler_ctx); 562 if (ceed_data->qf_ics) 563 CeedQFunctionSetContext(ceed_data->qf_ics, ceed_data->euler_context); 564 if (ceed_data->qf_apply_sur) 565 CeedQFunctionSetContext(ceed_data->qf_apply_sur, ceed_data->euler_context); 566 567 // ADVECTION and ADVECTION2D 568 } else { 569 CeedQFunctionContextCreate(ceed, &ceed_data->advection_context); 570 CeedQFunctionContextSetData(ceed_data->advection_context, CEED_MEM_HOST, 571 CEED_USE_POINTER, 572 sizeof(*phys->advection_ctx), phys->advection_ctx); 573 if (ceed_data->qf_rhs_vol) 574 CeedQFunctionSetContext(ceed_data->qf_rhs_vol, ceed_data->advection_context); 575 if (ceed_data->qf_ifunction_vol) 576 CeedQFunctionSetContext(ceed_data->qf_ifunction_vol, 577 ceed_data->advection_context); 578 if (ceed_data->qf_apply_sur) 579 CeedQFunctionSetContext(ceed_data->qf_apply_sur, ceed_data->advection_context); 580 } 581 582 PetscFunctionReturn(0); 583 } 584