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 jac_data_size_vol = num_comp_q + 3, 238 P = app_ctx->degree + 1, 239 Q = P + app_ctx->q_extra; 240 CeedElemRestriction elem_restr_jd_i; 241 CeedVector jac_data; 242 243 // ----------------------------------------------------------------------------- 244 // CEED Bases 245 // ----------------------------------------------------------------------------- 246 CeedBasisCreateTensorH1Lagrange(ceed, dim, num_comp_q, P, Q, CEED_GAUSS, 247 &ceed_data->basis_q); 248 CeedBasisCreateTensorH1Lagrange(ceed, dim, num_comp_x, 2, Q, CEED_GAUSS, 249 &ceed_data->basis_x); 250 CeedBasisCreateTensorH1Lagrange(ceed, dim, num_comp_x, 2, P, 251 CEED_GAUSS_LOBATTO, &ceed_data->basis_xc); 252 253 // ----------------------------------------------------------------------------- 254 // CEED Restrictions 255 // ----------------------------------------------------------------------------- 256 // -- Create restriction 257 ierr = GetRestrictionForDomain(ceed, dm, 0, 0, 0, Q, q_data_size_vol, 258 &ceed_data->elem_restr_q, &ceed_data->elem_restr_x, 259 &ceed_data->elem_restr_qd_i); CHKERRQ(ierr); 260 261 ierr = GetRestrictionForDomain(ceed, dm, 0, 0, 0, Q, jac_data_size_vol, 262 NULL, NULL, 263 &elem_restr_jd_i); CHKERRQ(ierr); 264 // -- Create E vectors 265 CeedElemRestrictionCreateVector(ceed_data->elem_restr_q, &user->q_ceed, NULL); 266 CeedElemRestrictionCreateVector(ceed_data->elem_restr_q, &user->q_dot_ceed, 267 NULL); 268 CeedElemRestrictionCreateVector(ceed_data->elem_restr_q, &user->g_ceed, NULL); 269 270 // ----------------------------------------------------------------------------- 271 // CEED QFunctions 272 // ----------------------------------------------------------------------------- 273 // -- Create QFunction for quadrature data 274 CeedQFunctionCreateInterior(ceed, 1, problem->setup_vol.qfunction, 275 problem->setup_vol.qfunction_loc, 276 &ceed_data->qf_setup_vol); 277 if (problem->setup_vol.qfunction_context) { 278 CeedQFunctionSetContext(ceed_data->qf_setup_vol, 279 problem->setup_vol.qfunction_context); 280 CeedQFunctionContextDestroy(&problem->setup_vol.qfunction_context); 281 } 282 CeedQFunctionAddInput(ceed_data->qf_setup_vol, "dx", num_comp_x*dim, 283 CEED_EVAL_GRAD); 284 CeedQFunctionAddInput(ceed_data->qf_setup_vol, "weight", 1, CEED_EVAL_WEIGHT); 285 CeedQFunctionAddOutput(ceed_data->qf_setup_vol, "qdata", q_data_size_vol, 286 CEED_EVAL_NONE); 287 288 // -- Create QFunction for ICs 289 CeedQFunctionCreateInterior(ceed, 1, problem->ics.qfunction, 290 problem->ics.qfunction_loc, 291 &ceed_data->qf_ics); 292 CeedQFunctionSetContext(ceed_data->qf_ics, problem->ics.qfunction_context); 293 CeedQFunctionContextDestroy(&problem->ics.qfunction_context); 294 CeedQFunctionAddInput(ceed_data->qf_ics, "x", num_comp_x, CEED_EVAL_INTERP); 295 CeedQFunctionAddOutput(ceed_data->qf_ics, "q0", num_comp_q, CEED_EVAL_NONE); 296 297 // -- Create QFunction for RHS 298 if (problem->apply_vol_rhs.qfunction) { 299 CeedQFunctionCreateInterior(ceed, 1, problem->apply_vol_rhs.qfunction, 300 problem->apply_vol_rhs.qfunction_loc, &ceed_data->qf_rhs_vol); 301 CeedQFunctionSetContext(ceed_data->qf_rhs_vol, 302 problem->apply_vol_rhs.qfunction_context); 303 CeedQFunctionContextDestroy(&problem->apply_vol_rhs.qfunction_context); 304 CeedQFunctionAddInput(ceed_data->qf_rhs_vol, "q", num_comp_q, CEED_EVAL_INTERP); 305 CeedQFunctionAddInput(ceed_data->qf_rhs_vol, "Grad_q", num_comp_q*dim, 306 CEED_EVAL_GRAD); 307 CeedQFunctionAddInput(ceed_data->qf_rhs_vol, "qdata", q_data_size_vol, 308 CEED_EVAL_NONE); 309 CeedQFunctionAddInput(ceed_data->qf_rhs_vol, "x", num_comp_x, CEED_EVAL_INTERP); 310 CeedQFunctionAddOutput(ceed_data->qf_rhs_vol, "v", num_comp_q, 311 CEED_EVAL_INTERP); 312 CeedQFunctionAddOutput(ceed_data->qf_rhs_vol, "Grad_v", num_comp_q*dim, 313 CEED_EVAL_GRAD); 314 } 315 316 // -- Create QFunction for IFunction 317 if (problem->apply_vol_ifunction.qfunction) { 318 CeedQFunctionCreateInterior(ceed, 1, problem->apply_vol_ifunction.qfunction, 319 problem->apply_vol_ifunction.qfunction_loc, &ceed_data->qf_ifunction_vol); 320 CeedQFunctionSetContext(ceed_data->qf_ifunction_vol, 321 problem->apply_vol_ifunction.qfunction_context); 322 CeedQFunctionContextDestroy(&problem->apply_vol_ifunction.qfunction_context); 323 CeedQFunctionAddInput(ceed_data->qf_ifunction_vol, "q", num_comp_q, 324 CEED_EVAL_INTERP); 325 CeedQFunctionAddInput(ceed_data->qf_ifunction_vol, "Grad_q", num_comp_q*dim, 326 CEED_EVAL_GRAD); 327 CeedQFunctionAddInput(ceed_data->qf_ifunction_vol, "q dot", num_comp_q, 328 CEED_EVAL_INTERP); 329 CeedQFunctionAddInput(ceed_data->qf_ifunction_vol, "qdata", q_data_size_vol, 330 CEED_EVAL_NONE); 331 CeedQFunctionAddInput(ceed_data->qf_ifunction_vol, "x", num_comp_x, 332 CEED_EVAL_INTERP); 333 CeedQFunctionAddOutput(ceed_data->qf_ifunction_vol, "v", num_comp_q, 334 CEED_EVAL_INTERP); 335 CeedQFunctionAddOutput(ceed_data->qf_ifunction_vol, "Grad_v", num_comp_q*dim, 336 CEED_EVAL_GRAD); 337 CeedQFunctionAddOutput(ceed_data->qf_ifunction_vol, "jac_data", 338 jac_data_size_vol, CEED_EVAL_NONE); 339 } 340 341 // --------------------------------------------------------------------------- 342 // Element coordinates 343 // --------------------------------------------------------------------------- 344 // -- Create CEED vector 345 CeedElemRestrictionCreateVector(ceed_data->elem_restr_x, &ceed_data->x_coord, 346 NULL); 347 348 // -- Copy PETSc vector in CEED vector 349 Vec X_loc; 350 const PetscScalar *X_loc_array; 351 ierr = DMGetCoordinatesLocal(dm, &X_loc); CHKERRQ(ierr); 352 ierr = VecScale(X_loc, problem->dm_scale); CHKERRQ(ierr); 353 ierr = VecGetArrayRead(X_loc, &X_loc_array); CHKERRQ(ierr); 354 CeedVectorSetArray(ceed_data->x_coord, CEED_MEM_HOST, CEED_COPY_VALUES, 355 (PetscScalar *)X_loc_array); 356 ierr = VecRestoreArrayRead(X_loc, &X_loc_array); CHKERRQ(ierr); 357 358 // ----------------------------------------------------------------------------- 359 // CEED vectors 360 // ----------------------------------------------------------------------------- 361 // -- Create CEED vector for geometric data 362 CeedInt num_qpts_vol; 363 PetscInt loc_num_elem_vol; 364 CeedBasisGetNumQuadraturePoints(ceed_data->basis_q, &num_qpts_vol); 365 CeedElemRestrictionGetNumElements(ceed_data->elem_restr_q, &loc_num_elem_vol); 366 CeedVectorCreate(ceed, q_data_size_vol*loc_num_elem_vol*num_qpts_vol, 367 &ceed_data->q_data); 368 369 CeedElemRestrictionCreateVector(elem_restr_jd_i, &jac_data, NULL); 370 // ----------------------------------------------------------------------------- 371 // CEED Operators 372 // ----------------------------------------------------------------------------- 373 // -- Create CEED operator for quadrature data 374 CeedOperatorCreate(ceed, ceed_data->qf_setup_vol, NULL, NULL, 375 &ceed_data->op_setup_vol); 376 CeedOperatorSetField(ceed_data->op_setup_vol, "dx", ceed_data->elem_restr_x, 377 ceed_data->basis_x, CEED_VECTOR_ACTIVE); 378 CeedOperatorSetField(ceed_data->op_setup_vol, "weight", 379 CEED_ELEMRESTRICTION_NONE, ceed_data->basis_x, CEED_VECTOR_NONE); 380 CeedOperatorSetField(ceed_data->op_setup_vol, "qdata", 381 ceed_data->elem_restr_qd_i, CEED_BASIS_COLLOCATED, CEED_VECTOR_ACTIVE); 382 383 // -- Create CEED operator for ICs 384 CeedOperatorCreate(ceed, ceed_data->qf_ics, NULL, NULL, &ceed_data->op_ics); 385 CeedOperatorSetField(ceed_data->op_ics, "x", ceed_data->elem_restr_x, 386 ceed_data->basis_xc, CEED_VECTOR_ACTIVE); 387 CeedOperatorSetField(ceed_data->op_ics, "q0", ceed_data->elem_restr_q, 388 CEED_BASIS_COLLOCATED, CEED_VECTOR_ACTIVE); 389 CeedOperatorContextGetFieldLabel(ceed_data->op_ics, "evaluation time", 390 &user->phys->ics_time_label); 391 392 // Create CEED operator for RHS 393 if (ceed_data->qf_rhs_vol) { 394 CeedOperator op; 395 CeedOperatorCreate(ceed, ceed_data->qf_rhs_vol, NULL, NULL, &op); 396 CeedOperatorSetField(op, "q", ceed_data->elem_restr_q, ceed_data->basis_q, 397 CEED_VECTOR_ACTIVE); 398 CeedOperatorSetField(op, "Grad_q", ceed_data->elem_restr_q, ceed_data->basis_q, 399 CEED_VECTOR_ACTIVE); 400 CeedOperatorSetField(op, "qdata", ceed_data->elem_restr_qd_i, 401 CEED_BASIS_COLLOCATED, ceed_data->q_data); 402 CeedOperatorSetField(op, "x", ceed_data->elem_restr_x, ceed_data->basis_x, 403 ceed_data->x_coord); 404 CeedOperatorSetField(op, "v", ceed_data->elem_restr_q, ceed_data->basis_q, 405 CEED_VECTOR_ACTIVE); 406 CeedOperatorSetField(op, "Grad_v", ceed_data->elem_restr_q, ceed_data->basis_q, 407 CEED_VECTOR_ACTIVE); 408 user->op_rhs_vol = op; 409 } 410 411 // -- CEED operator for IFunction 412 if (ceed_data->qf_ifunction_vol) { 413 CeedOperator op; 414 CeedOperatorCreate(ceed, ceed_data->qf_ifunction_vol, NULL, NULL, &op); 415 CeedOperatorSetField(op, "q", ceed_data->elem_restr_q, ceed_data->basis_q, 416 CEED_VECTOR_ACTIVE); 417 CeedOperatorSetField(op, "Grad_q", ceed_data->elem_restr_q, ceed_data->basis_q, 418 CEED_VECTOR_ACTIVE); 419 CeedOperatorSetField(op, "q dot", ceed_data->elem_restr_q, ceed_data->basis_q, 420 user->q_dot_ceed); 421 CeedOperatorSetField(op, "qdata", ceed_data->elem_restr_qd_i, 422 CEED_BASIS_COLLOCATED, ceed_data->q_data); 423 CeedOperatorSetField(op, "x", ceed_data->elem_restr_x, ceed_data->basis_x, 424 ceed_data->x_coord); 425 CeedOperatorSetField(op, "v", ceed_data->elem_restr_q, ceed_data->basis_q, 426 CEED_VECTOR_ACTIVE); 427 CeedOperatorSetField(op, "Grad_v", ceed_data->elem_restr_q, ceed_data->basis_q, 428 CEED_VECTOR_ACTIVE); 429 CeedOperatorSetField(op, "jac_data", elem_restr_jd_i, 430 CEED_BASIS_COLLOCATED, jac_data); 431 432 user->op_ifunction_vol = op; 433 } 434 435 // ***************************************************************************** 436 // Set up CEED objects for the exterior domain (surface) 437 // ***************************************************************************** 438 CeedInt height = 1, 439 dim_sur = dim - height, 440 P_sur = app_ctx->degree + 1, 441 Q_sur = P_sur + app_ctx->q_extra; 442 const CeedInt q_data_size_sur = problem->q_data_size_sur; 443 444 // ----------------------------------------------------------------------------- 445 // CEED Bases 446 // ----------------------------------------------------------------------------- 447 CeedBasisCreateTensorH1Lagrange(ceed, dim_sur, num_comp_q, P_sur, Q_sur, 448 CEED_GAUSS, &ceed_data->basis_q_sur); 449 CeedBasisCreateTensorH1Lagrange(ceed, dim_sur, num_comp_x, 2, Q_sur, CEED_GAUSS, 450 &ceed_data->basis_x_sur); 451 452 // ----------------------------------------------------------------------------- 453 // CEED QFunctions 454 // ----------------------------------------------------------------------------- 455 // -- Create QFunction for quadrature data 456 CeedQFunctionCreateInterior(ceed, 1, problem->setup_sur.qfunction, 457 problem->setup_sur.qfunction_loc, 458 &ceed_data->qf_setup_sur); 459 if (problem->setup_sur.qfunction_context) { 460 CeedQFunctionSetContext(ceed_data->qf_setup_sur, 461 problem->setup_sur.qfunction_context); 462 CeedQFunctionContextDestroy(&problem->setup_sur.qfunction_context); 463 } 464 CeedQFunctionAddInput(ceed_data->qf_setup_sur, "dx", num_comp_x*dim_sur, 465 CEED_EVAL_GRAD); 466 CeedQFunctionAddInput(ceed_data->qf_setup_sur, "weight", 1, CEED_EVAL_WEIGHT); 467 CeedQFunctionAddOutput(ceed_data->qf_setup_sur, "surface qdata", 468 q_data_size_sur, CEED_EVAL_NONE); 469 470 // -- Creat QFunction for inflow boundaries 471 if (problem->apply_inflow.qfunction) { 472 CeedQFunctionCreateInterior(ceed, 1, problem->apply_inflow.qfunction, 473 problem->apply_inflow.qfunction_loc, &ceed_data->qf_apply_inflow); 474 CeedQFunctionSetContext(ceed_data->qf_apply_inflow, 475 problem->apply_inflow.qfunction_context); 476 CeedQFunctionContextDestroy(&problem->apply_inflow.qfunction_context); 477 CeedQFunctionAddInput(ceed_data->qf_apply_inflow, "q", num_comp_q, 478 CEED_EVAL_INTERP); 479 CeedQFunctionAddInput(ceed_data->qf_apply_inflow, "surface qdata", 480 q_data_size_sur, CEED_EVAL_NONE); 481 CeedQFunctionAddInput(ceed_data->qf_apply_inflow, "x", num_comp_x, 482 CEED_EVAL_INTERP); 483 CeedQFunctionAddOutput(ceed_data->qf_apply_inflow, "v", num_comp_q, 484 CEED_EVAL_INTERP); 485 } 486 487 // -- Creat QFunction for outflow boundaries 488 if (problem->apply_outflow.qfunction) { 489 CeedQFunctionCreateInterior(ceed, 1, problem->apply_outflow.qfunction, 490 problem->apply_outflow.qfunction_loc, &ceed_data->qf_apply_outflow); 491 CeedQFunctionSetContext(ceed_data->qf_apply_outflow, 492 problem->apply_outflow.qfunction_context); 493 CeedQFunctionContextDestroy(&problem->apply_outflow.qfunction_context); 494 CeedQFunctionAddInput(ceed_data->qf_apply_outflow, "q", num_comp_q, 495 CEED_EVAL_INTERP); 496 CeedQFunctionAddInput(ceed_data->qf_apply_outflow, "surface qdata", 497 q_data_size_sur, CEED_EVAL_NONE); 498 CeedQFunctionAddInput(ceed_data->qf_apply_outflow, "x", num_comp_x, 499 CEED_EVAL_INTERP); 500 CeedQFunctionAddOutput(ceed_data->qf_apply_outflow, "v", num_comp_q, 501 CEED_EVAL_INTERP); 502 } 503 504 // ***************************************************************************** 505 // CEED Operator Apply 506 // ***************************************************************************** 507 // -- Apply CEED Operator for the geometric data 508 CeedOperatorApply(ceed_data->op_setup_vol, ceed_data->x_coord, 509 ceed_data->q_data, CEED_REQUEST_IMMEDIATE); 510 511 // -- Create and apply CEED Composite Operator for the entire domain 512 if (!user->phys->implicit) { // RHS 513 ierr = CreateOperatorForDomain(ceed, dm, bc, ceed_data, user->phys, 514 user->op_rhs_vol, height, P_sur, Q_sur, 515 q_data_size_sur, &user->op_rhs); CHKERRQ(ierr); 516 } else { // IFunction 517 ierr = CreateOperatorForDomain(ceed, dm, bc, ceed_data, user->phys, 518 user->op_ifunction_vol, height, P_sur, Q_sur, 519 q_data_size_sur, &user->op_ifunction); CHKERRQ(ierr); 520 } 521 CeedElemRestrictionDestroy(&elem_restr_jd_i); 522 CeedVectorDestroy(&jac_data); 523 PetscFunctionReturn(0); 524 } 525