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