1 // Copyright (c) 2017-2022, Lawrence Livermore National Security, LLC and other CEED contributors. 2 // All Rights Reserved. See the top-level LICENSE and NOTICE files for details. 3 // 4 // SPDX-License-Identifier: BSD-2-Clause 5 // 6 // This file is part of CEED: http://github.com/ceed 7 8 /// @file 9 /// Setup libCEED for Navier-Stokes example using PETSc 10 11 #include "../navierstokes.h" 12 13 // Utility function - essential BC dofs are encoded in closure indices as -(i+1). 14 PetscInt Involute(PetscInt i) { 15 return i >= 0 ? i : -(i+1); 16 } 17 18 // Utility function to create local CEED restriction 19 PetscErrorCode CreateRestrictionFromPlex(Ceed ceed, DM dm, CeedInt height, 20 DMLabel domain_label, CeedInt value, CeedElemRestriction *elem_restr) { 21 PetscInt num_elem, elem_size, num_dof, num_comp, *elem_restr_offsets; 22 PetscErrorCode ierr; 23 24 PetscFunctionBeginUser; 25 ierr = DMPlexGetLocalOffsets(dm, domain_label, value, height, 0, &num_elem, 26 &elem_size, &num_comp, &num_dof, &elem_restr_offsets); 27 CHKERRQ(ierr); 28 29 CeedElemRestrictionCreate(ceed, num_elem, elem_size, num_comp, 30 1, num_dof, CEED_MEM_HOST, CEED_COPY_VALUES, 31 elem_restr_offsets, elem_restr); 32 ierr = PetscFree(elem_restr_offsets); CHKERRQ(ierr); 33 34 PetscFunctionReturn(0); 35 } 36 37 // Utility function to get Ceed Restriction for each domain 38 PetscErrorCode GetRestrictionForDomain(Ceed ceed, DM dm, CeedInt height, 39 DMLabel domain_label, PetscInt value, 40 CeedInt Q, CeedInt q_data_size, 41 CeedElemRestriction *elem_restr_q, 42 CeedElemRestriction *elem_restr_x, 43 CeedElemRestriction *elem_restr_qd_i) { 44 DM dm_coord; 45 CeedInt dim, loc_num_elem; 46 CeedInt Q_dim; 47 CeedElemRestriction elem_restr_tmp; 48 PetscErrorCode ierr; 49 PetscFunctionBeginUser; 50 51 ierr = DMGetDimension(dm, &dim); CHKERRQ(ierr); 52 dim -= height; 53 Q_dim = CeedIntPow(Q, dim); 54 ierr = CreateRestrictionFromPlex(ceed, dm, height, domain_label, value, 55 &elem_restr_tmp); 56 CHKERRQ(ierr); 57 if (elem_restr_q) *elem_restr_q = elem_restr_tmp; 58 if (elem_restr_x) { 59 ierr = DMGetCoordinateDM(dm, &dm_coord); CHKERRQ(ierr); 60 ierr = DMPlexSetClosurePermutationTensor(dm_coord, PETSC_DETERMINE, NULL); 61 CHKERRQ(ierr); 62 ierr = CreateRestrictionFromPlex(ceed, dm_coord, height, domain_label, value, 63 elem_restr_x); 64 CHKERRQ(ierr); 65 } 66 if (elem_restr_qd_i) { 67 CeedElemRestrictionGetNumElements(elem_restr_tmp, &loc_num_elem); 68 CeedElemRestrictionCreateStrided(ceed, loc_num_elem, Q_dim, 69 q_data_size, q_data_size*loc_num_elem*Q_dim, 70 CEED_STRIDES_BACKEND, elem_restr_qd_i); 71 } 72 if (!elem_restr_q) CeedElemRestrictionDestroy(&elem_restr_tmp); 73 PetscFunctionReturn(0); 74 } 75 76 // Utility function to create CEED Composite Operator for the entire domain 77 PetscErrorCode CreateOperatorForDomain(Ceed ceed, DM dm, SimpleBC bc, 78 CeedData ceed_data, Physics phys, 79 CeedOperator op_apply_vol, CeedInt height, 80 CeedInt P_sur, CeedInt Q_sur, CeedInt q_data_size_sur, 81 CeedOperator *op_apply) { 82 //CeedInt dim; 83 DMLabel domain_label; 84 PetscErrorCode ierr; 85 PetscFunctionBeginUser; 86 87 // Create Composite Operaters 88 CeedCompositeOperatorCreate(ceed, op_apply); 89 90 // --Apply Sub-Operator for the volume 91 CeedCompositeOperatorAddSub(*op_apply, op_apply_vol); 92 93 // -- Create Sub-Operator for in/outflow BCs 94 if (phys->has_neumann || 1) { 95 // --- Setup 96 ierr = DMGetLabel(dm, "Face Sets", &domain_label); CHKERRQ(ierr); 97 //ierr = DMGetDimension(dm, &dim); CHKERRQ(ierr); 98 99 // --- Get number of quadrature points for the boundaries 100 CeedInt num_qpts_sur; 101 CeedBasisGetNumQuadraturePoints(ceed_data->basis_q_sur, &num_qpts_sur); 102 103 // --- Create Sub-Operator for inflow boundaries 104 for (CeedInt i=0; i < bc->num_inflow; i++) { 105 CeedVector q_data_sur; 106 CeedOperator op_setup_sur, op_apply_inflow; 107 CeedElemRestriction elem_restr_x_sur, elem_restr_q_sur, elem_restr_qd_i_sur; 108 109 // ---- CEED Restriction 110 ierr = GetRestrictionForDomain(ceed, dm, height, domain_label, bc->inflows[i], 111 Q_sur, q_data_size_sur, &elem_restr_q_sur, &elem_restr_x_sur, 112 &elem_restr_qd_i_sur); 113 CHKERRQ(ierr); 114 115 // ---- CEED Vector 116 PetscInt loc_num_elem_sur; 117 CeedElemRestrictionGetNumElements(elem_restr_q_sur, &loc_num_elem_sur); 118 CeedVectorCreate(ceed, q_data_size_sur*loc_num_elem_sur*num_qpts_sur, 119 &q_data_sur); 120 121 // ---- CEED Operator 122 // ----- CEED Operator for Setup (geometric factors) 123 CeedOperatorCreate(ceed, ceed_data->qf_setup_sur, NULL, NULL, &op_setup_sur); 124 CeedOperatorSetField(op_setup_sur, "dx", elem_restr_x_sur, 125 ceed_data->basis_x_sur, CEED_VECTOR_ACTIVE); 126 CeedOperatorSetField(op_setup_sur, "weight", CEED_ELEMRESTRICTION_NONE, 127 ceed_data->basis_x_sur, CEED_VECTOR_NONE); 128 CeedOperatorSetField(op_setup_sur, "surface qdata", elem_restr_qd_i_sur, 129 CEED_BASIS_COLLOCATED, CEED_VECTOR_ACTIVE); 130 131 // ----- CEED Operator for Physics 132 CeedOperatorCreate(ceed, ceed_data->qf_apply_inflow, NULL, NULL, 133 &op_apply_inflow); 134 CeedOperatorSetField(op_apply_inflow, "q", elem_restr_q_sur, 135 ceed_data->basis_q_sur, CEED_VECTOR_ACTIVE); 136 CeedOperatorSetField(op_apply_inflow, "surface qdata", elem_restr_qd_i_sur, 137 CEED_BASIS_COLLOCATED, q_data_sur); 138 CeedOperatorSetField(op_apply_inflow, "x", elem_restr_x_sur, 139 ceed_data->basis_x_sur, ceed_data->x_coord); 140 CeedOperatorSetField(op_apply_inflow, "v", elem_restr_q_sur, 141 ceed_data->basis_q_sur, CEED_VECTOR_ACTIVE); 142 143 // ----- Apply CEED operator for Setup 144 CeedOperatorApply(op_setup_sur, ceed_data->x_coord, q_data_sur, 145 CEED_REQUEST_IMMEDIATE); 146 147 // ----- Apply Sub-Operator for Physics 148 CeedCompositeOperatorAddSub(*op_apply, op_apply_inflow); 149 150 // ----- Cleanup 151 CeedVectorDestroy(&q_data_sur); 152 CeedElemRestrictionDestroy(&elem_restr_q_sur); 153 CeedElemRestrictionDestroy(&elem_restr_x_sur); 154 CeedElemRestrictionDestroy(&elem_restr_qd_i_sur); 155 CeedOperatorDestroy(&op_setup_sur); 156 CeedOperatorDestroy(&op_apply_inflow); 157 } 158 159 // --- Create Sub-Operator for outflow boundaries 160 for (CeedInt i=0; i < bc->num_outflow; i++) { 161 CeedVector q_data_sur; 162 CeedOperator op_setup_sur, op_apply_outflow; 163 CeedElemRestriction elem_restr_x_sur, elem_restr_q_sur, elem_restr_qd_i_sur; 164 165 // ---- CEED Restriction 166 ierr = GetRestrictionForDomain(ceed, dm, height, domain_label, bc->outflows[i], 167 Q_sur, q_data_size_sur, &elem_restr_q_sur, &elem_restr_x_sur, 168 &elem_restr_qd_i_sur); 169 CHKERRQ(ierr); 170 171 // ---- CEED Vector 172 PetscInt loc_num_elem_sur; 173 CeedElemRestrictionGetNumElements(elem_restr_q_sur, &loc_num_elem_sur); 174 CeedVectorCreate(ceed, q_data_size_sur*loc_num_elem_sur*num_qpts_sur, 175 &q_data_sur); 176 177 // ---- CEED Operator 178 // ----- CEED Operator for Setup (geometric factors) 179 CeedOperatorCreate(ceed, ceed_data->qf_setup_sur, NULL, NULL, &op_setup_sur); 180 CeedOperatorSetField(op_setup_sur, "dx", elem_restr_x_sur, 181 ceed_data->basis_x_sur, CEED_VECTOR_ACTIVE); 182 CeedOperatorSetField(op_setup_sur, "weight", CEED_ELEMRESTRICTION_NONE, 183 ceed_data->basis_x_sur, CEED_VECTOR_NONE); 184 CeedOperatorSetField(op_setup_sur, "surface qdata", elem_restr_qd_i_sur, 185 CEED_BASIS_COLLOCATED, CEED_VECTOR_ACTIVE); 186 187 // ----- CEED Operator for Physics 188 CeedOperatorCreate(ceed, ceed_data->qf_apply_outflow, NULL, NULL, 189 &op_apply_outflow); 190 CeedOperatorSetField(op_apply_outflow, "q", elem_restr_q_sur, 191 ceed_data->basis_q_sur, CEED_VECTOR_ACTIVE); 192 CeedOperatorSetField(op_apply_outflow, "surface qdata", elem_restr_qd_i_sur, 193 CEED_BASIS_COLLOCATED, q_data_sur); 194 CeedOperatorSetField(op_apply_outflow, "x", elem_restr_x_sur, 195 ceed_data->basis_x_sur, ceed_data->x_coord); 196 CeedOperatorSetField(op_apply_outflow, "v", elem_restr_q_sur, 197 ceed_data->basis_q_sur, CEED_VECTOR_ACTIVE); 198 199 // ----- Apply CEED operator for Setup 200 CeedOperatorApply(op_setup_sur, ceed_data->x_coord, q_data_sur, 201 CEED_REQUEST_IMMEDIATE); 202 203 // ----- Apply Sub-Operator for Physics 204 CeedCompositeOperatorAddSub(*op_apply, op_apply_outflow); 205 206 // ----- Cleanup 207 CeedVectorDestroy(&q_data_sur); 208 CeedElemRestrictionDestroy(&elem_restr_q_sur); 209 CeedElemRestrictionDestroy(&elem_restr_x_sur); 210 CeedElemRestrictionDestroy(&elem_restr_qd_i_sur); 211 CeedOperatorDestroy(&op_setup_sur); 212 CeedOperatorDestroy(&op_apply_outflow); 213 } 214 } 215 216 // ----- Get Context Labels for Operator 217 CeedOperatorContextGetFieldLabel(*op_apply, "solution time", 218 &phys->solution_time_label); 219 CeedOperatorContextGetFieldLabel(*op_apply, "timestep size", 220 &phys->timestep_size_label); 221 222 PetscFunctionReturn(0); 223 } 224 225 PetscErrorCode SetupLibceed(Ceed ceed, CeedData ceed_data, DM dm, User user, 226 AppCtx app_ctx, ProblemData *problem, SimpleBC bc) { 227 PetscErrorCode ierr; 228 PetscFunctionBeginUser; 229 230 // ***************************************************************************** 231 // Set up CEED objects for the interior domain (volume) 232 // ***************************************************************************** 233 const PetscInt num_comp_q = 5; 234 const CeedInt dim = problem->dim, 235 num_comp_x = problem->dim, 236 q_data_size_vol = problem->q_data_size_vol, 237 P = app_ctx->degree + 1, 238 Q = P + app_ctx->q_extra; 239 240 // ----------------------------------------------------------------------------- 241 // CEED Bases 242 // ----------------------------------------------------------------------------- 243 CeedBasisCreateTensorH1Lagrange(ceed, dim, num_comp_q, P, Q, CEED_GAUSS, 244 &ceed_data->basis_q); 245 CeedBasisCreateTensorH1Lagrange(ceed, dim, num_comp_x, 2, Q, CEED_GAUSS, 246 &ceed_data->basis_x); 247 CeedBasisCreateTensorH1Lagrange(ceed, dim, num_comp_x, 2, P, 248 CEED_GAUSS_LOBATTO, &ceed_data->basis_xc); 249 250 // ----------------------------------------------------------------------------- 251 // CEED Restrictions 252 // ----------------------------------------------------------------------------- 253 // -- Create restriction 254 ierr = GetRestrictionForDomain(ceed, dm, 0, 0, 0, Q, q_data_size_vol, 255 &ceed_data->elem_restr_q, &ceed_data->elem_restr_x, 256 &ceed_data->elem_restr_qd_i); CHKERRQ(ierr); 257 // -- Create E vectors 258 CeedElemRestrictionCreateVector(ceed_data->elem_restr_q, &user->q_ceed, NULL); 259 CeedElemRestrictionCreateVector(ceed_data->elem_restr_q, &user->q_dot_ceed, 260 NULL); 261 CeedElemRestrictionCreateVector(ceed_data->elem_restr_q, &user->g_ceed, NULL); 262 263 // ----------------------------------------------------------------------------- 264 // CEED QFunctions 265 // ----------------------------------------------------------------------------- 266 // -- Create QFunction for quadrature data 267 CeedQFunctionCreateInterior(ceed, 1, problem->setup_vol.qfunction, 268 problem->setup_vol.qfunction_loc, 269 &ceed_data->qf_setup_vol); 270 if (problem->setup_vol.qfunction_context) { 271 CeedQFunctionSetContext(ceed_data->qf_setup_vol, 272 problem->setup_vol.qfunction_context); 273 CeedQFunctionContextDestroy(&problem->setup_vol.qfunction_context); 274 } 275 CeedQFunctionAddInput(ceed_data->qf_setup_vol, "dx", num_comp_x*dim, 276 CEED_EVAL_GRAD); 277 CeedQFunctionAddInput(ceed_data->qf_setup_vol, "weight", 1, CEED_EVAL_WEIGHT); 278 CeedQFunctionAddOutput(ceed_data->qf_setup_vol, "qdata", q_data_size_vol, 279 CEED_EVAL_NONE); 280 281 // -- Create QFunction for ICs 282 CeedQFunctionCreateInterior(ceed, 1, problem->ics.qfunction, 283 problem->ics.qfunction_loc, 284 &ceed_data->qf_ics); 285 CeedQFunctionSetContext(ceed_data->qf_ics, problem->ics.qfunction_context); 286 CeedQFunctionContextDestroy(&problem->ics.qfunction_context); 287 CeedQFunctionAddInput(ceed_data->qf_ics, "x", num_comp_x, CEED_EVAL_INTERP); 288 CeedQFunctionAddOutput(ceed_data->qf_ics, "q0", num_comp_q, CEED_EVAL_NONE); 289 290 // -- Create QFunction for RHS 291 if (problem->apply_vol_rhs.qfunction) { 292 CeedQFunctionCreateInterior(ceed, 1, problem->apply_vol_rhs.qfunction, 293 problem->apply_vol_rhs.qfunction_loc, &ceed_data->qf_rhs_vol); 294 CeedQFunctionSetContext(ceed_data->qf_rhs_vol, 295 problem->apply_vol_rhs.qfunction_context); 296 CeedQFunctionContextDestroy(&problem->apply_vol_rhs.qfunction_context); 297 CeedQFunctionAddInput(ceed_data->qf_rhs_vol, "q", num_comp_q, CEED_EVAL_INTERP); 298 CeedQFunctionAddInput(ceed_data->qf_rhs_vol, "dq", num_comp_q*dim, 299 CEED_EVAL_GRAD); 300 CeedQFunctionAddInput(ceed_data->qf_rhs_vol, "qdata", q_data_size_vol, 301 CEED_EVAL_NONE); 302 CeedQFunctionAddInput(ceed_data->qf_rhs_vol, "x", num_comp_x, CEED_EVAL_INTERP); 303 CeedQFunctionAddOutput(ceed_data->qf_rhs_vol, "v", num_comp_q, 304 CEED_EVAL_INTERP); 305 CeedQFunctionAddOutput(ceed_data->qf_rhs_vol, "dv", num_comp_q*dim, 306 CEED_EVAL_GRAD); 307 } 308 309 // -- Create QFunction for IFunction 310 if (problem->apply_vol_ifunction.qfunction) { 311 CeedQFunctionCreateInterior(ceed, 1, problem->apply_vol_ifunction.qfunction, 312 problem->apply_vol_ifunction.qfunction_loc, &ceed_data->qf_ifunction_vol); 313 CeedQFunctionSetContext(ceed_data->qf_ifunction_vol, 314 problem->apply_vol_ifunction.qfunction_context); 315 CeedQFunctionContextDestroy(&problem->apply_vol_ifunction.qfunction_context); 316 CeedQFunctionAddInput(ceed_data->qf_ifunction_vol, "q", num_comp_q, 317 CEED_EVAL_INTERP); 318 CeedQFunctionAddInput(ceed_data->qf_ifunction_vol, "dq", num_comp_q*dim, 319 CEED_EVAL_GRAD); 320 CeedQFunctionAddInput(ceed_data->qf_ifunction_vol, "q dot", num_comp_q, 321 CEED_EVAL_INTERP); 322 CeedQFunctionAddInput(ceed_data->qf_ifunction_vol, "qdata", q_data_size_vol, 323 CEED_EVAL_NONE); 324 CeedQFunctionAddInput(ceed_data->qf_ifunction_vol, "x", num_comp_x, 325 CEED_EVAL_INTERP); 326 CeedQFunctionAddOutput(ceed_data->qf_ifunction_vol, "v", num_comp_q, 327 CEED_EVAL_INTERP); 328 CeedQFunctionAddOutput(ceed_data->qf_ifunction_vol, "dv", num_comp_q*dim, 329 CEED_EVAL_GRAD); 330 } 331 332 // --------------------------------------------------------------------------- 333 // Element coordinates 334 // --------------------------------------------------------------------------- 335 // -- Create CEED vector 336 CeedElemRestrictionCreateVector(ceed_data->elem_restr_x, &ceed_data->x_coord, 337 NULL); 338 339 // -- Copy PETSc vector in CEED vector 340 Vec X_loc; 341 const PetscScalar *X_loc_array; 342 ierr = DMGetCoordinatesLocal(dm, &X_loc); CHKERRQ(ierr); 343 ierr = VecScale(X_loc, problem->dm_scale); CHKERRQ(ierr); 344 ierr = VecGetArrayRead(X_loc, &X_loc_array); CHKERRQ(ierr); 345 CeedVectorSetArray(ceed_data->x_coord, CEED_MEM_HOST, CEED_COPY_VALUES, 346 (PetscScalar *)X_loc_array); 347 ierr = VecRestoreArrayRead(X_loc, &X_loc_array); CHKERRQ(ierr); 348 349 // ----------------------------------------------------------------------------- 350 // CEED vectors 351 // ----------------------------------------------------------------------------- 352 // -- Create CEED vector for geometric data 353 CeedInt num_qpts_vol; 354 PetscInt loc_num_elem_vol; 355 CeedBasisGetNumQuadraturePoints(ceed_data->basis_q, &num_qpts_vol); 356 CeedElemRestrictionGetNumElements(ceed_data->elem_restr_q, &loc_num_elem_vol); 357 CeedVectorCreate(ceed, q_data_size_vol*loc_num_elem_vol*num_qpts_vol, 358 &ceed_data->q_data); 359 360 // ----------------------------------------------------------------------------- 361 // CEED Operators 362 // ----------------------------------------------------------------------------- 363 // -- Create CEED operator for quadrature data 364 CeedOperatorCreate(ceed, ceed_data->qf_setup_vol, NULL, NULL, 365 &ceed_data->op_setup_vol); 366 CeedOperatorSetField(ceed_data->op_setup_vol, "dx", ceed_data->elem_restr_x, 367 ceed_data->basis_x, CEED_VECTOR_ACTIVE); 368 CeedOperatorSetField(ceed_data->op_setup_vol, "weight", 369 CEED_ELEMRESTRICTION_NONE, ceed_data->basis_x, CEED_VECTOR_NONE); 370 CeedOperatorSetField(ceed_data->op_setup_vol, "qdata", 371 ceed_data->elem_restr_qd_i, CEED_BASIS_COLLOCATED, CEED_VECTOR_ACTIVE); 372 373 // -- Create CEED operator for ICs 374 CeedOperatorCreate(ceed, ceed_data->qf_ics, NULL, NULL, &ceed_data->op_ics); 375 CeedOperatorSetField(ceed_data->op_ics, "x", ceed_data->elem_restr_x, 376 ceed_data->basis_xc, CEED_VECTOR_ACTIVE); 377 CeedOperatorSetField(ceed_data->op_ics, "q0", ceed_data->elem_restr_q, 378 CEED_BASIS_COLLOCATED, CEED_VECTOR_ACTIVE); 379 CeedOperatorContextGetFieldLabel(ceed_data->op_ics, "evaluation time", 380 &user->phys->ics_time_label); 381 382 // Create CEED operator for RHS 383 if (ceed_data->qf_rhs_vol) { 384 CeedOperator op; 385 CeedOperatorCreate(ceed, ceed_data->qf_rhs_vol, NULL, NULL, &op); 386 CeedOperatorSetField(op, "q", ceed_data->elem_restr_q, ceed_data->basis_q, 387 CEED_VECTOR_ACTIVE); 388 CeedOperatorSetField(op, "dq", ceed_data->elem_restr_q, ceed_data->basis_q, 389 CEED_VECTOR_ACTIVE); 390 CeedOperatorSetField(op, "qdata", ceed_data->elem_restr_qd_i, 391 CEED_BASIS_COLLOCATED, ceed_data->q_data); 392 CeedOperatorSetField(op, "x", ceed_data->elem_restr_x, ceed_data->basis_x, 393 ceed_data->x_coord); 394 CeedOperatorSetField(op, "v", ceed_data->elem_restr_q, ceed_data->basis_q, 395 CEED_VECTOR_ACTIVE); 396 CeedOperatorSetField(op, "dv", ceed_data->elem_restr_q, ceed_data->basis_q, 397 CEED_VECTOR_ACTIVE); 398 user->op_rhs_vol = op; 399 } 400 401 // -- CEED operator for IFunction 402 if (ceed_data->qf_ifunction_vol) { 403 CeedOperator op; 404 CeedOperatorCreate(ceed, ceed_data->qf_ifunction_vol, NULL, NULL, &op); 405 CeedOperatorSetField(op, "q", ceed_data->elem_restr_q, ceed_data->basis_q, 406 CEED_VECTOR_ACTIVE); 407 CeedOperatorSetField(op, "dq", ceed_data->elem_restr_q, ceed_data->basis_q, 408 CEED_VECTOR_ACTIVE); 409 CeedOperatorSetField(op, "q dot", ceed_data->elem_restr_q, ceed_data->basis_q, 410 user->q_dot_ceed); 411 CeedOperatorSetField(op, "qdata", ceed_data->elem_restr_qd_i, 412 CEED_BASIS_COLLOCATED, ceed_data->q_data); 413 CeedOperatorSetField(op, "x", ceed_data->elem_restr_x, ceed_data->basis_x, 414 ceed_data->x_coord); 415 CeedOperatorSetField(op, "v", ceed_data->elem_restr_q, ceed_data->basis_q, 416 CEED_VECTOR_ACTIVE); 417 CeedOperatorSetField(op, "dv", ceed_data->elem_restr_q, ceed_data->basis_q, 418 CEED_VECTOR_ACTIVE); 419 user->op_ifunction_vol = op; 420 } 421 422 // ***************************************************************************** 423 // Set up CEED objects for the exterior domain (surface) 424 // ***************************************************************************** 425 CeedInt height = 1, 426 dim_sur = dim - height, 427 P_sur = app_ctx->degree + 1, 428 Q_sur = P_sur + app_ctx->q_extra; 429 const CeedInt q_data_size_sur = problem->q_data_size_sur; 430 431 // ----------------------------------------------------------------------------- 432 // CEED Bases 433 // ----------------------------------------------------------------------------- 434 CeedBasisCreateTensorH1Lagrange(ceed, dim_sur, num_comp_q, P_sur, Q_sur, 435 CEED_GAUSS, &ceed_data->basis_q_sur); 436 CeedBasisCreateTensorH1Lagrange(ceed, dim_sur, num_comp_x, 2, Q_sur, CEED_GAUSS, 437 &ceed_data->basis_x_sur); 438 439 // ----------------------------------------------------------------------------- 440 // CEED QFunctions 441 // ----------------------------------------------------------------------------- 442 // -- Create QFunction for quadrature data 443 CeedQFunctionCreateInterior(ceed, 1, problem->setup_sur.qfunction, 444 problem->setup_sur.qfunction_loc, 445 &ceed_data->qf_setup_sur); 446 if (problem->setup_sur.qfunction_context) { 447 CeedQFunctionSetContext(ceed_data->qf_setup_sur, 448 problem->setup_sur.qfunction_context); 449 CeedQFunctionContextDestroy(&problem->setup_sur.qfunction_context); 450 } 451 CeedQFunctionAddInput(ceed_data->qf_setup_sur, "dx", num_comp_x*dim_sur, 452 CEED_EVAL_GRAD); 453 CeedQFunctionAddInput(ceed_data->qf_setup_sur, "weight", 1, CEED_EVAL_WEIGHT); 454 CeedQFunctionAddOutput(ceed_data->qf_setup_sur, "surface qdata", 455 q_data_size_sur, CEED_EVAL_NONE); 456 457 // -- Creat QFunction for inflow boundaries 458 if (problem->apply_inflow.qfunction) { 459 CeedQFunctionCreateInterior(ceed, 1, problem->apply_inflow.qfunction, 460 problem->apply_inflow.qfunction_loc, &ceed_data->qf_apply_inflow); 461 CeedQFunctionSetContext(ceed_data->qf_apply_inflow, 462 problem->apply_inflow.qfunction_context); 463 CeedQFunctionContextDestroy(&problem->apply_inflow.qfunction_context); 464 CeedQFunctionAddInput(ceed_data->qf_apply_inflow, "q", num_comp_q, 465 CEED_EVAL_INTERP); 466 CeedQFunctionAddInput(ceed_data->qf_apply_inflow, "surface qdata", 467 q_data_size_sur, CEED_EVAL_NONE); 468 CeedQFunctionAddInput(ceed_data->qf_apply_inflow, "x", num_comp_x, 469 CEED_EVAL_INTERP); 470 CeedQFunctionAddOutput(ceed_data->qf_apply_inflow, "v", num_comp_q, 471 CEED_EVAL_INTERP); 472 } 473 474 // -- Creat QFunction for outflow boundaries 475 if (problem->apply_outflow.qfunction) { 476 CeedQFunctionCreateInterior(ceed, 1, problem->apply_outflow.qfunction, 477 problem->apply_outflow.qfunction_loc, &ceed_data->qf_apply_outflow); 478 CeedQFunctionSetContext(ceed_data->qf_apply_outflow, 479 problem->apply_outflow.qfunction_context); 480 CeedQFunctionContextDestroy(&problem->apply_outflow.qfunction_context); 481 CeedQFunctionAddInput(ceed_data->qf_apply_outflow, "q", num_comp_q, 482 CEED_EVAL_INTERP); 483 CeedQFunctionAddInput(ceed_data->qf_apply_outflow, "surface qdata", 484 q_data_size_sur, CEED_EVAL_NONE); 485 CeedQFunctionAddInput(ceed_data->qf_apply_outflow, "x", num_comp_x, 486 CEED_EVAL_INTERP); 487 CeedQFunctionAddOutput(ceed_data->qf_apply_outflow, "v", num_comp_q, 488 CEED_EVAL_INTERP); 489 } 490 491 // ***************************************************************************** 492 // CEED Operator Apply 493 // ***************************************************************************** 494 // -- Apply CEED Operator for the geometric data 495 CeedOperatorApply(ceed_data->op_setup_vol, ceed_data->x_coord, 496 ceed_data->q_data, CEED_REQUEST_IMMEDIATE); 497 498 // -- Create and apply CEED Composite Operator for the entire domain 499 if (!user->phys->implicit) { // RHS 500 ierr = CreateOperatorForDomain(ceed, dm, bc, ceed_data, user->phys, 501 user->op_rhs_vol, height, P_sur, Q_sur, 502 q_data_size_sur, &user->op_rhs); CHKERRQ(ierr); 503 } else { // IFunction 504 ierr = CreateOperatorForDomain(ceed, dm, bc, ceed_data, user->phys, 505 user->op_ifunction_vol, height, P_sur, Q_sur, 506 q_data_size_sur, &user->op_ifunction); CHKERRQ(ierr); 507 } 508 509 PetscFunctionReturn(0); 510 } 511