1 // Copyright (c) 2017-2018, Lawrence Livermore National Security, LLC. 2 // Produced at the Lawrence Livermore National Laboratory. LLNL-CODE-734707. 3 // All Rights 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 #include <ceed/ceed.h> 18 #include <ceed/backend.h> 19 #include <math.h> 20 #include <stdbool.h> 21 #include <stddef.h> 22 #include <stdint.h> 23 #include "ceed-ref.h" 24 25 //------------------------------------------------------------------------------ 26 // Setup Input/Output Fields 27 //------------------------------------------------------------------------------ 28 static int CeedOperatorSetupFields_Ref(CeedQFunction qf, CeedOperator op, 29 bool inOrOut, 30 CeedVector *full_evecs, CeedVector *e_vecs, 31 CeedVector *q_vecs, CeedInt starte, 32 CeedInt num_fields, CeedInt Q) { 33 CeedInt dim, ierr, size, P; 34 Ceed ceed; 35 ierr = CeedOperatorGetCeed(op, &ceed); CeedChkBackend(ierr); 36 CeedBasis basis; 37 CeedElemRestriction elem_restr; 38 CeedOperatorField *op_fields; 39 CeedQFunctionField *qf_fields; 40 if (inOrOut) { 41 ierr = CeedOperatorGetFields(op, NULL, NULL, NULL, &op_fields); 42 CeedChkBackend(ierr); 43 ierr = CeedQFunctionGetFields(qf, NULL, NULL, NULL, &qf_fields); 44 CeedChkBackend(ierr); 45 } else { 46 ierr = CeedOperatorGetFields(op, NULL, &op_fields, NULL, NULL); 47 CeedChkBackend(ierr); 48 ierr = CeedQFunctionGetFields(qf, NULL, &qf_fields, NULL, NULL); 49 CeedChkBackend(ierr); 50 } 51 52 // Loop over fields 53 for (CeedInt i=0; i<num_fields; i++) { 54 CeedEvalMode eval_mode; 55 ierr = CeedQFunctionFieldGetEvalMode(qf_fields[i], &eval_mode); 56 CeedChkBackend(ierr); 57 58 if (eval_mode != CEED_EVAL_WEIGHT) { 59 ierr = CeedOperatorFieldGetElemRestriction(op_fields[i], &elem_restr); 60 CeedChkBackend(ierr); 61 ierr = CeedElemRestrictionCreateVector(elem_restr, NULL, 62 &full_evecs[i+starte]); 63 CeedChkBackend(ierr); 64 } 65 66 switch(eval_mode) { 67 case CEED_EVAL_NONE: 68 ierr = CeedQFunctionFieldGetSize(qf_fields[i], &size); CeedChkBackend(ierr); 69 ierr = CeedVectorCreate(ceed, Q*size, &q_vecs[i]); CeedChkBackend(ierr); 70 break; 71 case CEED_EVAL_INTERP: 72 ierr = CeedQFunctionFieldGetSize(qf_fields[i], &size); CeedChkBackend(ierr); 73 ierr = CeedElemRestrictionGetElementSize(elem_restr, &P); 74 CeedChkBackend(ierr); 75 ierr = CeedVectorCreate(ceed, P*size, &e_vecs[i]); CeedChkBackend(ierr); 76 ierr = CeedVectorCreate(ceed, Q*size, &q_vecs[i]); CeedChkBackend(ierr); 77 break; 78 case CEED_EVAL_GRAD: 79 ierr = CeedOperatorFieldGetBasis(op_fields[i], &basis); CeedChkBackend(ierr); 80 ierr = CeedQFunctionFieldGetSize(qf_fields[i], &size); CeedChkBackend(ierr); 81 ierr = CeedBasisGetDimension(basis, &dim); CeedChkBackend(ierr); 82 ierr = CeedElemRestrictionGetElementSize(elem_restr, &P); 83 CeedChkBackend(ierr); 84 ierr = CeedVectorCreate(ceed, P*size/dim, &e_vecs[i]); CeedChkBackend(ierr); 85 ierr = CeedVectorCreate(ceed, Q*size, &q_vecs[i]); CeedChkBackend(ierr); 86 break; 87 case CEED_EVAL_WEIGHT: // Only on input fields 88 ierr = CeedOperatorFieldGetBasis(op_fields[i], &basis); CeedChkBackend(ierr); 89 ierr = CeedVectorCreate(ceed, Q, &q_vecs[i]); CeedChkBackend(ierr); 90 ierr = CeedBasisApply(basis, 1, CEED_NOTRANSPOSE, CEED_EVAL_WEIGHT, 91 CEED_VECTOR_NONE, q_vecs[i]); CeedChkBackend(ierr); 92 break; 93 case CEED_EVAL_DIV: 94 break; // Not implemented 95 case CEED_EVAL_CURL: 96 break; // Not implemented 97 } 98 } 99 return CEED_ERROR_SUCCESS; 100 } 101 102 //------------------------------------------------------------------------------ 103 // Setup Operator 104 //------------------------------------------------------------------------------/* 105 static int CeedOperatorSetup_Ref(CeedOperator op) { 106 int ierr; 107 bool setup_done; 108 ierr = CeedOperatorIsSetupDone(op, &setup_done); CeedChkBackend(ierr); 109 if (setup_done) return CEED_ERROR_SUCCESS; 110 Ceed ceed; 111 ierr = CeedOperatorGetCeed(op, &ceed); CeedChkBackend(ierr); 112 CeedOperator_Ref *impl; 113 ierr = CeedOperatorGetData(op, &impl); CeedChkBackend(ierr); 114 CeedQFunction qf; 115 ierr = CeedOperatorGetQFunction(op, &qf); CeedChkBackend(ierr); 116 CeedInt Q, num_input_fields, num_output_fields; 117 ierr = CeedOperatorGetNumQuadraturePoints(op, &Q); CeedChkBackend(ierr); 118 ierr = CeedQFunctionIsIdentity(qf, &impl->is_identity_qf); CeedChkBackend(ierr); 119 CeedOperatorField *op_input_fields, *op_output_fields; 120 ierr = CeedOperatorGetFields(op, &num_input_fields, &op_input_fields, 121 &num_output_fields, &op_output_fields); 122 CeedChkBackend(ierr); 123 CeedQFunctionField *qf_input_fields, *qf_output_fields; 124 ierr = CeedQFunctionGetFields(qf, NULL, &qf_input_fields, NULL, 125 &qf_output_fields); 126 CeedChkBackend(ierr); 127 128 // Allocate 129 ierr = CeedCalloc(num_input_fields + num_output_fields, &impl->e_vecs); 130 CeedChkBackend(ierr); 131 ierr = CeedCalloc(num_input_fields + num_output_fields, &impl->e_data); 132 CeedChkBackend(ierr); 133 134 ierr = CeedCalloc(16, &impl->input_state); CeedChkBackend(ierr); 135 ierr = CeedCalloc(16, &impl->e_vecs_in); CeedChkBackend(ierr); 136 ierr = CeedCalloc(16, &impl->e_vecs_out); CeedChkBackend(ierr); 137 ierr = CeedCalloc(16, &impl->q_vecs_in); CeedChkBackend(ierr); 138 ierr = CeedCalloc(16, &impl->q_vecs_out); CeedChkBackend(ierr); 139 140 impl->num_e_vecs_in = num_input_fields; 141 impl->num_e_vecs_out = num_output_fields; 142 143 // Set up infield and outfield e_vecs and q_vecs 144 // Infields 145 ierr = CeedOperatorSetupFields_Ref(qf, op, 0, impl->e_vecs, 146 impl->e_vecs_in, impl->q_vecs_in, 0, 147 num_input_fields, Q); 148 CeedChkBackend(ierr); 149 // Outfields 150 ierr = CeedOperatorSetupFields_Ref(qf, op, 1, impl->e_vecs, 151 impl->e_vecs_out, impl->q_vecs_out, 152 num_input_fields, num_output_fields, Q); 153 CeedChkBackend(ierr); 154 155 // Identity QFunctions 156 if (impl->is_identity_qf) { 157 CeedEvalMode in_mode, out_mode; 158 CeedQFunctionField *in_fields, *out_fields; 159 ierr = CeedQFunctionGetFields(qf, NULL, &in_fields, NULL, &out_fields); 160 CeedChkBackend(ierr); 161 ierr = CeedQFunctionFieldGetEvalMode(in_fields[0], &in_mode); 162 CeedChkBackend(ierr); 163 ierr = CeedQFunctionFieldGetEvalMode(out_fields[0], &out_mode); 164 CeedChkBackend(ierr); 165 166 if (in_mode == CEED_EVAL_NONE && out_mode == CEED_EVAL_NONE) { 167 impl->is_identity_restr_op = true; 168 } else { 169 ierr = CeedVectorDestroy(&impl->q_vecs_out[0]); CeedChkBackend(ierr); 170 impl->q_vecs_out[0] = impl->q_vecs_in[0]; 171 ierr = CeedVectorAddReference(impl->q_vecs_in[0]); CeedChkBackend(ierr); 172 } 173 } 174 175 ierr = CeedOperatorSetSetupDone(op); CeedChkBackend(ierr); 176 177 return CEED_ERROR_SUCCESS; 178 } 179 180 //------------------------------------------------------------------------------ 181 // Setup Operator Inputs 182 //------------------------------------------------------------------------------ 183 static inline int CeedOperatorSetupInputs_Ref(CeedInt num_input_fields, 184 CeedQFunctionField *qf_input_fields, CeedOperatorField *op_input_fields, 185 CeedVector in_vec, const bool skip_active, CeedOperator_Ref *impl, 186 CeedRequest *request) { 187 CeedInt ierr; 188 CeedEvalMode eval_mode; 189 CeedVector vec; 190 CeedElemRestriction elem_restr; 191 uint64_t state; 192 193 for (CeedInt i=0; i<num_input_fields; i++) { 194 // Get input vector 195 ierr = CeedOperatorFieldGetVector(op_input_fields[i], &vec); 196 CeedChkBackend(ierr); 197 if (vec == CEED_VECTOR_ACTIVE) { 198 if (skip_active) 199 continue; 200 else 201 vec = in_vec; 202 } 203 204 ierr = CeedQFunctionFieldGetEvalMode(qf_input_fields[i], &eval_mode); 205 CeedChkBackend(ierr); 206 // Restrict and Evec 207 if (eval_mode == CEED_EVAL_WEIGHT) { // Skip 208 } else { 209 // Restrict 210 ierr = CeedVectorGetState(vec, &state); CeedChkBackend(ierr); 211 // Skip restriction if input is unchanged 212 if (state != impl->input_state[i] || vec == in_vec) { 213 ierr = CeedOperatorFieldGetElemRestriction(op_input_fields[i], &elem_restr); 214 CeedChkBackend(ierr); 215 ierr = CeedElemRestrictionApply(elem_restr, CEED_NOTRANSPOSE, vec, 216 impl->e_vecs[i], request); CeedChkBackend(ierr); 217 impl->input_state[i] = state; 218 } 219 // Get evec 220 ierr = CeedVectorGetArrayRead(impl->e_vecs[i], CEED_MEM_HOST, 221 (const CeedScalar **) &impl->e_data[i]); 222 CeedChkBackend(ierr); 223 } 224 } 225 return CEED_ERROR_SUCCESS; 226 } 227 228 //------------------------------------------------------------------------------ 229 // Input Basis Action 230 //------------------------------------------------------------------------------ 231 static inline int CeedOperatorInputBasis_Ref(CeedInt e, CeedInt Q, 232 CeedQFunctionField *qf_input_fields, CeedOperatorField *op_input_fields, 233 CeedInt num_input_fields, const bool skip_active, CeedOperator_Ref *impl) { 234 CeedInt ierr; 235 CeedInt dim, elem_size, size; 236 CeedElemRestriction elem_restr; 237 CeedEvalMode eval_mode; 238 CeedBasis basis; 239 240 for (CeedInt i=0; i<num_input_fields; i++) { 241 // Skip active input 242 if (skip_active) { 243 CeedVector vec; 244 ierr = CeedOperatorFieldGetVector(op_input_fields[i], &vec); 245 CeedChkBackend(ierr); 246 if (vec == CEED_VECTOR_ACTIVE) 247 continue; 248 } 249 // Get elem_size, eval_mode, size 250 ierr = CeedOperatorFieldGetElemRestriction(op_input_fields[i], &elem_restr); 251 CeedChkBackend(ierr); 252 ierr = CeedElemRestrictionGetElementSize(elem_restr, &elem_size); 253 CeedChkBackend(ierr); 254 ierr = CeedQFunctionFieldGetEvalMode(qf_input_fields[i], &eval_mode); 255 CeedChkBackend(ierr); 256 ierr = CeedQFunctionFieldGetSize(qf_input_fields[i], &size); 257 CeedChkBackend(ierr); 258 // Basis action 259 switch(eval_mode) { 260 case CEED_EVAL_NONE: 261 ierr = CeedVectorSetArray(impl->q_vecs_in[i], CEED_MEM_HOST, 262 CEED_USE_POINTER, &impl->e_data[i][e*Q*size]); 263 CeedChkBackend(ierr); 264 break; 265 case CEED_EVAL_INTERP: 266 ierr = CeedOperatorFieldGetBasis(op_input_fields[i], &basis); 267 CeedChkBackend(ierr); 268 ierr = CeedVectorSetArray(impl->e_vecs_in[i], CEED_MEM_HOST, 269 CEED_USE_POINTER, &impl->e_data[i][e*elem_size*size]); 270 CeedChkBackend(ierr); 271 ierr = CeedBasisApply(basis, 1, CEED_NOTRANSPOSE, CEED_EVAL_INTERP, 272 impl->e_vecs_in[i], impl->q_vecs_in[i]); CeedChkBackend(ierr); 273 break; 274 case CEED_EVAL_GRAD: 275 ierr = CeedOperatorFieldGetBasis(op_input_fields[i], &basis); 276 CeedChkBackend(ierr); 277 ierr = CeedBasisGetDimension(basis, &dim); CeedChkBackend(ierr); 278 ierr = CeedVectorSetArray(impl->e_vecs_in[i], CEED_MEM_HOST, 279 CEED_USE_POINTER, &impl->e_data[i][e*elem_size*size/dim]); 280 CeedChkBackend(ierr); 281 ierr = CeedBasisApply(basis, 1, CEED_NOTRANSPOSE, 282 CEED_EVAL_GRAD, impl->e_vecs_in[i], 283 impl->q_vecs_in[i]); CeedChkBackend(ierr); 284 break; 285 case CEED_EVAL_WEIGHT: 286 break; // No action 287 // LCOV_EXCL_START 288 case CEED_EVAL_DIV: 289 case CEED_EVAL_CURL: { 290 ierr = CeedOperatorFieldGetBasis(op_input_fields[i], &basis); 291 CeedChkBackend(ierr); 292 Ceed ceed; 293 ierr = CeedBasisGetCeed(basis, &ceed); CeedChkBackend(ierr); 294 return CeedError(ceed, CEED_ERROR_BACKEND, 295 "Ceed evaluation mode not implemented"); 296 // LCOV_EXCL_STOP 297 } 298 } 299 } 300 return CEED_ERROR_SUCCESS; 301 } 302 303 //------------------------------------------------------------------------------ 304 // Output Basis Action 305 //------------------------------------------------------------------------------ 306 static inline int CeedOperatorOutputBasis_Ref(CeedInt e, CeedInt Q, 307 CeedQFunctionField *qf_output_fields, CeedOperatorField *op_output_fields, 308 CeedInt num_input_fields, CeedInt num_output_fields, CeedOperator op, 309 CeedOperator_Ref *impl) { 310 CeedInt ierr; 311 CeedInt dim, elem_size, size; 312 CeedElemRestriction elem_restr; 313 CeedEvalMode eval_mode; 314 CeedBasis basis; 315 316 for (CeedInt i=0; i<num_output_fields; i++) { 317 // Get elem_size, eval_mode, size 318 ierr = CeedOperatorFieldGetElemRestriction(op_output_fields[i], &elem_restr); 319 CeedChkBackend(ierr); 320 ierr = CeedElemRestrictionGetElementSize(elem_restr, &elem_size); 321 CeedChkBackend(ierr); 322 ierr = CeedQFunctionFieldGetEvalMode(qf_output_fields[i], &eval_mode); 323 CeedChkBackend(ierr); 324 ierr = CeedQFunctionFieldGetSize(qf_output_fields[i], &size); 325 CeedChkBackend(ierr); 326 // Basis action 327 switch(eval_mode) { 328 case CEED_EVAL_NONE: 329 break; // No action 330 case CEED_EVAL_INTERP: 331 ierr = CeedOperatorFieldGetBasis(op_output_fields[i], &basis); 332 CeedChkBackend(ierr); 333 ierr = CeedVectorSetArray(impl->e_vecs_out[i], CEED_MEM_HOST, 334 CEED_USE_POINTER, 335 &impl->e_data[i + num_input_fields][e*elem_size*size]); 336 CeedChkBackend(ierr); 337 ierr = CeedBasisApply(basis, 1, CEED_TRANSPOSE, 338 CEED_EVAL_INTERP, impl->q_vecs_out[i], 339 impl->e_vecs_out[i]); CeedChkBackend(ierr); 340 break; 341 case CEED_EVAL_GRAD: 342 ierr = CeedOperatorFieldGetBasis(op_output_fields[i], &basis); 343 CeedChkBackend(ierr); 344 ierr = CeedBasisGetDimension(basis, &dim); CeedChkBackend(ierr); 345 ierr = CeedVectorSetArray(impl->e_vecs_out[i], CEED_MEM_HOST, 346 CEED_USE_POINTER, 347 &impl->e_data[i + num_input_fields][e*elem_size*size/dim]); 348 CeedChkBackend(ierr); 349 ierr = CeedBasisApply(basis, 1, CEED_TRANSPOSE, 350 CEED_EVAL_GRAD, impl->q_vecs_out[i], 351 impl->e_vecs_out[i]); CeedChkBackend(ierr); 352 break; 353 // LCOV_EXCL_START 354 case CEED_EVAL_WEIGHT: { 355 Ceed ceed; 356 ierr = CeedOperatorGetCeed(op, &ceed); CeedChkBackend(ierr); 357 return CeedError(ceed, CEED_ERROR_BACKEND, 358 "CEED_EVAL_WEIGHT cannot be an output " 359 "evaluation mode"); 360 } 361 case CEED_EVAL_DIV: 362 case CEED_EVAL_CURL: { 363 Ceed ceed; 364 ierr = CeedOperatorGetCeed(op, &ceed); CeedChkBackend(ierr); 365 return CeedError(ceed, CEED_ERROR_BACKEND, 366 "Ceed evaluation mode not implemented"); 367 // LCOV_EXCL_STOP 368 } 369 } 370 } 371 return CEED_ERROR_SUCCESS; 372 } 373 374 //------------------------------------------------------------------------------ 375 // Restore Input Vectors 376 //------------------------------------------------------------------------------ 377 static inline int CeedOperatorRestoreInputs_Ref(CeedInt num_input_fields, 378 CeedQFunctionField *qf_input_fields, CeedOperatorField *op_input_fields, 379 const bool skip_active, CeedOperator_Ref *impl) { 380 CeedInt ierr; 381 CeedEvalMode eval_mode; 382 383 for (CeedInt i=0; i<num_input_fields; i++) { 384 // Skip active inputs 385 if (skip_active) { 386 CeedVector vec; 387 ierr = CeedOperatorFieldGetVector(op_input_fields[i], &vec); 388 CeedChkBackend(ierr); 389 if (vec == CEED_VECTOR_ACTIVE) 390 continue; 391 } 392 // Restore input 393 ierr = CeedQFunctionFieldGetEvalMode(qf_input_fields[i], &eval_mode); 394 CeedChkBackend(ierr); 395 if (eval_mode == CEED_EVAL_WEIGHT) { // Skip 396 } else { 397 ierr = CeedVectorRestoreArrayRead(impl->e_vecs[i], 398 (const CeedScalar **) &impl->e_data[i]); 399 CeedChkBackend(ierr); 400 } 401 } 402 return CEED_ERROR_SUCCESS; 403 } 404 405 //------------------------------------------------------------------------------ 406 // Operator Apply 407 //------------------------------------------------------------------------------ 408 static int CeedOperatorApplyAdd_Ref(CeedOperator op, CeedVector in_vec, 409 CeedVector out_vec, CeedRequest *request) { 410 int ierr; 411 CeedOperator_Ref *impl; 412 ierr = CeedOperatorGetData(op, &impl); CeedChkBackend(ierr); 413 CeedQFunction qf; 414 ierr = CeedOperatorGetQFunction(op, &qf); CeedChkBackend(ierr); 415 CeedInt Q, num_elem, num_input_fields, num_output_fields, size; 416 ierr = CeedOperatorGetNumQuadraturePoints(op, &Q); CeedChkBackend(ierr); 417 ierr = CeedOperatorGetNumElements(op, &num_elem); CeedChkBackend(ierr); 418 CeedOperatorField *op_input_fields, *op_output_fields; 419 ierr = CeedOperatorGetFields(op, &num_input_fields, &op_input_fields, 420 &num_output_fields, &op_output_fields); 421 CeedChkBackend(ierr); 422 CeedQFunctionField *qf_input_fields, *qf_output_fields; 423 ierr = CeedQFunctionGetFields(qf, NULL, &qf_input_fields, NULL, 424 &qf_output_fields); 425 CeedChkBackend(ierr); 426 CeedEvalMode eval_mode; 427 CeedVector vec; 428 CeedElemRestriction elem_restr; 429 430 // Setup 431 ierr = CeedOperatorSetup_Ref(op); CeedChkBackend(ierr); 432 433 // Restriction only operator 434 if (impl->is_identity_restr_op) { 435 ierr = CeedOperatorFieldGetElemRestriction(op_input_fields[0], &elem_restr); 436 CeedChkBackend(ierr); 437 ierr = CeedElemRestrictionApply(elem_restr, CEED_NOTRANSPOSE, in_vec, 438 impl->e_vecs[0], request); CeedChkBackend(ierr); 439 ierr = CeedOperatorFieldGetElemRestriction(op_output_fields[0], &elem_restr); 440 CeedChkBackend(ierr); 441 ierr = CeedElemRestrictionApply(elem_restr, CEED_TRANSPOSE, impl->e_vecs[0], 442 out_vec, request); CeedChkBackend(ierr); 443 return CEED_ERROR_SUCCESS; 444 } 445 446 // Input Evecs and Restriction 447 ierr = CeedOperatorSetupInputs_Ref(num_input_fields, qf_input_fields, 448 op_input_fields, in_vec, false, impl, 449 request); CeedChkBackend(ierr); 450 451 // Output Evecs 452 for (CeedInt i=0; i<num_output_fields; i++) { 453 ierr = CeedVectorGetArray(impl->e_vecs[i+impl->num_e_vecs_in], CEED_MEM_HOST, 454 &impl->e_data[i + num_input_fields]); CeedChkBackend(ierr); 455 } 456 457 // Loop through elements 458 for (CeedInt e=0; e<num_elem; e++) { 459 // Output pointers 460 for (CeedInt i=0; i<num_output_fields; i++) { 461 ierr = CeedQFunctionFieldGetEvalMode(qf_output_fields[i], &eval_mode); 462 CeedChkBackend(ierr); 463 if (eval_mode == CEED_EVAL_NONE) { 464 ierr = CeedQFunctionFieldGetSize(qf_output_fields[i], &size); 465 CeedChkBackend(ierr); 466 ierr = CeedVectorSetArray(impl->q_vecs_out[i], CEED_MEM_HOST, 467 CEED_USE_POINTER, 468 &impl->e_data[i + num_input_fields][e*Q*size]); 469 CeedChkBackend(ierr); 470 } 471 } 472 473 // Input basis apply 474 ierr = CeedOperatorInputBasis_Ref(e, Q, qf_input_fields, op_input_fields, 475 num_input_fields, false, impl); 476 CeedChkBackend(ierr); 477 478 // Q function 479 if (!impl->is_identity_qf) { 480 ierr = CeedQFunctionApply(qf, Q, impl->q_vecs_in, impl->q_vecs_out); 481 CeedChkBackend(ierr); 482 } 483 484 // Output basis apply 485 ierr = CeedOperatorOutputBasis_Ref(e, Q, qf_output_fields, op_output_fields, 486 num_input_fields, num_output_fields, op, impl); 487 CeedChkBackend(ierr); 488 } 489 490 // Output restriction 491 for (CeedInt i=0; i<num_output_fields; i++) { 492 // Restore Evec 493 ierr = CeedVectorRestoreArray(impl->e_vecs[i+impl->num_e_vecs_in], 494 &impl->e_data[i + num_input_fields]); 495 CeedChkBackend(ierr); 496 // Get output vector 497 ierr = CeedOperatorFieldGetVector(op_output_fields[i], &vec); 498 CeedChkBackend(ierr); 499 // Active 500 if (vec == CEED_VECTOR_ACTIVE) 501 vec = out_vec; 502 // Restrict 503 ierr = CeedOperatorFieldGetElemRestriction(op_output_fields[i], &elem_restr); 504 CeedChkBackend(ierr); 505 ierr = CeedElemRestrictionApply(elem_restr, CEED_TRANSPOSE, 506 impl->e_vecs[i+impl->num_e_vecs_in], vec, request); 507 CeedChkBackend(ierr); 508 } 509 510 // Restore input arrays 511 ierr = CeedOperatorRestoreInputs_Ref(num_input_fields, qf_input_fields, 512 op_input_fields, false, impl); 513 CeedChkBackend(ierr); 514 515 return CEED_ERROR_SUCCESS; 516 } 517 518 //------------------------------------------------------------------------------ 519 // Core code for assembling linear QFunction 520 //------------------------------------------------------------------------------ 521 static inline int CeedOperatorLinearAssembleQFunctionCore_Ref(CeedOperator op, 522 bool build_objects, CeedVector *assembled, CeedElemRestriction *rstr, 523 CeedRequest *request) { 524 int ierr; 525 CeedOperator_Ref *impl; 526 ierr = CeedOperatorGetData(op, &impl); CeedChkBackend(ierr); 527 CeedQFunction qf; 528 ierr = CeedOperatorGetQFunction(op, &qf); CeedChkBackend(ierr); 529 CeedInt Q, num_elem, num_input_fields, num_output_fields, size; 530 ierr = CeedOperatorGetNumQuadraturePoints(op, &Q); CeedChkBackend(ierr); 531 ierr = CeedOperatorGetNumElements(op, &num_elem); CeedChkBackend(ierr); 532 CeedOperatorField *op_input_fields, *op_output_fields; 533 ierr = CeedOperatorGetFields(op, &num_input_fields, &op_input_fields, 534 &num_output_fields, &op_output_fields); 535 CeedChkBackend(ierr); 536 CeedQFunctionField *qf_input_fields, *qf_output_fields; 537 ierr = CeedQFunctionGetFields(qf, NULL, &qf_input_fields, NULL, 538 &qf_output_fields); 539 CeedChkBackend(ierr); 540 CeedVector vec; 541 CeedInt num_active_in = impl->qf_num_active_in, 542 num_active_out = impl->qf_num_active_out; 543 CeedVector *active_in = impl->qf_active_in; 544 CeedScalar *a, *tmp; 545 Ceed ceed, ceed_parent; 546 ierr = CeedOperatorGetCeed(op, &ceed); CeedChkBackend(ierr); 547 ierr = CeedGetOperatorFallbackParentCeed(ceed, &ceed_parent); 548 CeedChkBackend(ierr); 549 ceed_parent = ceed_parent ? ceed_parent : ceed; 550 551 // Setup 552 ierr = CeedOperatorSetup_Ref(op); CeedChkBackend(ierr); 553 554 // Check for identity 555 if (impl->is_identity_qf) 556 // LCOV_EXCL_START 557 return CeedError(ceed, CEED_ERROR_BACKEND, 558 "Assembling identity QFunctions not supported"); 559 // LCOV_EXCL_STOP 560 561 // Input Evecs and Restriction 562 ierr = CeedOperatorSetupInputs_Ref(num_input_fields, qf_input_fields, 563 op_input_fields, NULL, true, impl, request); 564 CeedChkBackend(ierr); 565 566 // Count number of active input fields 567 if (!num_active_in) { 568 for (CeedInt i=0; i<num_input_fields; i++) { 569 // Get input vector 570 ierr = CeedOperatorFieldGetVector(op_input_fields[i], &vec); 571 CeedChkBackend(ierr); 572 // Check if active input 573 if (vec == CEED_VECTOR_ACTIVE) { 574 ierr = CeedQFunctionFieldGetSize(qf_input_fields[i], &size); 575 CeedChkBackend(ierr); 576 ierr = CeedVectorSetValue(impl->q_vecs_in[i], 0.0); CeedChkBackend(ierr); 577 ierr = CeedVectorGetArray(impl->q_vecs_in[i], CEED_MEM_HOST, &tmp); 578 CeedChkBackend(ierr); 579 ierr = CeedRealloc(num_active_in + size, &active_in); CeedChkBackend(ierr); 580 for (CeedInt field=0; field<size; field++) { 581 ierr = CeedVectorCreate(ceed, Q, &active_in[num_active_in+field]); 582 CeedChkBackend(ierr); 583 ierr = CeedVectorSetArray(active_in[num_active_in+field], CEED_MEM_HOST, 584 CEED_USE_POINTER, &tmp[field*Q]); 585 CeedChkBackend(ierr); 586 } 587 num_active_in += size; 588 ierr = CeedVectorRestoreArray(impl->q_vecs_in[i], &tmp); CeedChkBackend(ierr); 589 } 590 } 591 impl->qf_num_active_in = num_active_in; 592 impl->qf_active_in = active_in; 593 } 594 595 // Count number of active output fields 596 if (!num_active_out) { 597 for (CeedInt i=0; i<num_output_fields; i++) { 598 // Get output vector 599 ierr = CeedOperatorFieldGetVector(op_output_fields[i], &vec); 600 CeedChkBackend(ierr); 601 // Check if active output 602 if (vec == CEED_VECTOR_ACTIVE) { 603 ierr = CeedQFunctionFieldGetSize(qf_output_fields[i], &size); 604 CeedChkBackend(ierr); 605 num_active_out += size; 606 } 607 } 608 impl->qf_num_active_out = num_active_out; 609 } 610 611 // Check sizes 612 if (!num_active_in || !num_active_out) 613 // LCOV_EXCL_START 614 return CeedError(ceed, CEED_ERROR_BACKEND, 615 "Cannot assemble QFunction without active inputs " 616 "and outputs"); 617 // LCOV_EXCL_STOP 618 619 // Build objects if needed 620 if (build_objects) { 621 // Create output restriction 622 CeedInt strides[3] = {1, Q, num_active_in*num_active_out*Q}; /* *NOPAD* */ 623 ierr = CeedElemRestrictionCreateStrided(ceed_parent, num_elem, Q, 624 num_active_in*num_active_out, 625 num_active_in*num_active_out*num_elem*Q, 626 strides, rstr); CeedChkBackend(ierr); 627 // Create assembled vector 628 ierr = CeedVectorCreate(ceed_parent, num_elem*Q*num_active_in*num_active_out, 629 assembled); CeedChkBackend(ierr); 630 } 631 // Clear output vector 632 ierr = CeedVectorSetValue(*assembled, 0.0); CeedChkBackend(ierr); 633 ierr = CeedVectorGetArray(*assembled, CEED_MEM_HOST, &a); CeedChkBackend(ierr); 634 635 // Loop through elements 636 for (CeedInt e=0; e<num_elem; e++) { 637 // Input basis apply 638 ierr = CeedOperatorInputBasis_Ref(e, Q, qf_input_fields, op_input_fields, 639 num_input_fields, true, impl); 640 CeedChkBackend(ierr); 641 642 // Assemble QFunction 643 for (CeedInt in=0; in<num_active_in; in++) { 644 // Set Inputs 645 ierr = CeedVectorSetValue(active_in[in], 1.0); CeedChkBackend(ierr); 646 if (num_active_in > 1) { 647 ierr = CeedVectorSetValue(active_in[(in+num_active_in-1)%num_active_in], 648 0.0); CeedChkBackend(ierr); 649 } 650 // Set Outputs 651 for (CeedInt out=0; out<num_output_fields; out++) { 652 // Get output vector 653 ierr = CeedOperatorFieldGetVector(op_output_fields[out], &vec); 654 CeedChkBackend(ierr); 655 // Check if active output 656 if (vec == CEED_VECTOR_ACTIVE) { 657 CeedVectorSetArray(impl->q_vecs_out[out], CEED_MEM_HOST, 658 CEED_USE_POINTER, a); CeedChkBackend(ierr); 659 ierr = CeedQFunctionFieldGetSize(qf_output_fields[out], &size); 660 CeedChkBackend(ierr); 661 a += size*Q; // Advance the pointer by the size of the output 662 } 663 } 664 // Apply QFunction 665 ierr = CeedQFunctionApply(qf, Q, impl->q_vecs_in, impl->q_vecs_out); 666 CeedChkBackend(ierr); 667 } 668 } 669 670 // Un-set output Qvecs to prevent accidental overwrite of Assembled 671 for (CeedInt out=0; out<num_output_fields; out++) { 672 // Get output vector 673 ierr = CeedOperatorFieldGetVector(op_output_fields[out], &vec); 674 CeedChkBackend(ierr); 675 // Check if active output 676 if (vec == CEED_VECTOR_ACTIVE) { 677 CeedVectorTakeArray(impl->q_vecs_out[out], CEED_MEM_HOST, NULL); 678 CeedChkBackend(ierr); 679 } 680 } 681 682 // Restore input arrays 683 ierr = CeedOperatorRestoreInputs_Ref(num_input_fields, qf_input_fields, 684 op_input_fields, true, impl); 685 CeedChkBackend(ierr); 686 687 // Restore output 688 ierr = CeedVectorRestoreArray(*assembled, &a); CeedChkBackend(ierr); 689 690 return CEED_ERROR_SUCCESS; 691 } 692 693 //------------------------------------------------------------------------------ 694 // Assemble Linear QFunction 695 //------------------------------------------------------------------------------ 696 static int CeedOperatorLinearAssembleQFunction_Ref(CeedOperator op, 697 CeedVector *assembled, CeedElemRestriction *rstr, CeedRequest *request) { 698 return CeedOperatorLinearAssembleQFunctionCore_Ref(op, true, assembled, rstr, 699 request); 700 } 701 702 //------------------------------------------------------------------------------ 703 // Update Assembled Linear QFunction 704 //------------------------------------------------------------------------------ 705 static int CeedOperatorLinearAssembleQFunctionUpdate_Ref(CeedOperator op, 706 CeedVector assembled, CeedElemRestriction rstr, CeedRequest *request) { 707 return CeedOperatorLinearAssembleQFunctionCore_Ref(op, false, &assembled, 708 &rstr, request); 709 } 710 711 //------------------------------------------------------------------------------ 712 // Operator Destroy 713 //------------------------------------------------------------------------------ 714 static int CeedOperatorDestroy_Ref(CeedOperator op) { 715 int ierr; 716 CeedOperator_Ref *impl; 717 ierr = CeedOperatorGetData(op, &impl); CeedChkBackend(ierr); 718 719 for (CeedInt i=0; i<impl->num_e_vecs_in+impl->num_e_vecs_out; i++) { 720 ierr = CeedVectorDestroy(&impl->e_vecs[i]); CeedChkBackend(ierr); 721 } 722 ierr = CeedFree(&impl->e_vecs); CeedChkBackend(ierr); 723 ierr = CeedFree(&impl->e_data); CeedChkBackend(ierr); 724 ierr = CeedFree(&impl->input_state); CeedChkBackend(ierr); 725 726 for (CeedInt i=0; i<impl->num_e_vecs_in; i++) { 727 ierr = CeedVectorDestroy(&impl->e_vecs_in[i]); CeedChkBackend(ierr); 728 ierr = CeedVectorDestroy(&impl->q_vecs_in[i]); CeedChkBackend(ierr); 729 } 730 ierr = CeedFree(&impl->e_vecs_in); CeedChkBackend(ierr); 731 ierr = CeedFree(&impl->q_vecs_in); CeedChkBackend(ierr); 732 733 for (CeedInt i=0; i<impl->num_e_vecs_out; i++) { 734 ierr = CeedVectorDestroy(&impl->e_vecs_out[i]); CeedChkBackend(ierr); 735 ierr = CeedVectorDestroy(&impl->q_vecs_out[i]); CeedChkBackend(ierr); 736 } 737 ierr = CeedFree(&impl->e_vecs_out); CeedChkBackend(ierr); 738 ierr = CeedFree(&impl->q_vecs_out); CeedChkBackend(ierr); 739 740 // QFunction assembly 741 for (CeedInt i=0; i<impl->qf_num_active_in; i++) { 742 ierr = CeedVectorDestroy(&impl->qf_active_in[i]); CeedChkBackend(ierr); 743 } 744 ierr = CeedFree(&impl->qf_active_in); CeedChkBackend(ierr); 745 746 ierr = CeedFree(&impl); CeedChkBackend(ierr); 747 return CEED_ERROR_SUCCESS; 748 } 749 750 //------------------------------------------------------------------------------ 751 // Operator Create 752 //------------------------------------------------------------------------------ 753 int CeedOperatorCreate_Ref(CeedOperator op) { 754 int ierr; 755 Ceed ceed; 756 ierr = CeedOperatorGetCeed(op, &ceed); CeedChkBackend(ierr); 757 CeedOperator_Ref *impl; 758 759 ierr = CeedCalloc(1, &impl); CeedChkBackend(ierr); 760 ierr = CeedOperatorSetData(op, impl); CeedChkBackend(ierr); 761 762 ierr = CeedSetBackendFunction(ceed, "Operator", op, "LinearAssembleQFunction", 763 CeedOperatorLinearAssembleQFunction_Ref); 764 CeedChkBackend(ierr); 765 ierr = CeedSetBackendFunction(ceed, "Operator", op, 766 "LinearAssembleQFunctionUpdate", 767 CeedOperatorLinearAssembleQFunctionUpdate_Ref); 768 CeedChkBackend(ierr); 769 ierr = CeedSetBackendFunction(ceed, "Operator", op, "ApplyAdd", 770 CeedOperatorApplyAdd_Ref); CeedChkBackend(ierr); 771 ierr = CeedSetBackendFunction(ceed, "Operator", op, "Destroy", 772 CeedOperatorDestroy_Ref); CeedChkBackend(ierr); 773 return CEED_ERROR_SUCCESS; 774 } 775