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