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 // Assemble Linear QFunction 520 //------------------------------------------------------------------------------ 521 static int CeedOperatorLinearAssembleQFunction_Ref(CeedOperator op, 522 CeedVector *assembled, CeedElemRestriction *rstr, CeedRequest *request) { 523 int ierr; 524 CeedOperator_Ref *impl; 525 ierr = CeedOperatorGetData(op, &impl); CeedChkBackend(ierr); 526 CeedQFunction qf; 527 ierr = CeedOperatorGetQFunction(op, &qf); CeedChkBackend(ierr); 528 CeedInt Q, num_elem, num_input_fields, num_output_fields, size; 529 ierr = CeedOperatorGetNumQuadraturePoints(op, &Q); CeedChkBackend(ierr); 530 ierr = CeedOperatorGetNumElements(op, &num_elem); CeedChkBackend(ierr); 531 CeedOperatorField *op_input_fields, *op_output_fields; 532 ierr = CeedOperatorGetFields(op, &num_input_fields, &op_input_fields, 533 &num_output_fields, &op_output_fields); 534 CeedChkBackend(ierr); 535 CeedQFunctionField *qf_input_fields, *qf_output_fields; 536 ierr = CeedQFunctionGetFields(qf, NULL, &qf_input_fields, NULL, 537 &qf_output_fields); 538 CeedChkBackend(ierr); 539 CeedVector vec; 540 CeedInt num_active_in = 0, num_active_out = 0; 541 CeedVector *active_in = NULL; 542 CeedScalar *a, *tmp; 543 Ceed ceed, ceed_parent; 544 ierr = CeedOperatorGetCeed(op, &ceed); CeedChkBackend(ierr); 545 ierr = CeedGetOperatorFallbackParentCeed(ceed, &ceed_parent); 546 CeedChkBackend(ierr); 547 ceed_parent = ceed_parent ? ceed_parent : ceed; 548 549 // Setup 550 ierr = CeedOperatorSetup_Ref(op); CeedChkBackend(ierr); 551 552 // Check for identity 553 if (impl->is_identity_qf) 554 // LCOV_EXCL_START 555 return CeedError(ceed, CEED_ERROR_BACKEND, 556 "Assembling identity QFunctions not supported"); 557 // LCOV_EXCL_STOP 558 559 // Input Evecs and Restriction 560 ierr = CeedOperatorSetupInputs_Ref(num_input_fields, qf_input_fields, 561 op_input_fields, NULL, true, impl, request); 562 CeedChkBackend(ierr); 563 564 // Count number of active input fields 565 for (CeedInt i=0; i<num_input_fields; i++) { 566 // Get input vector 567 ierr = CeedOperatorFieldGetVector(op_input_fields[i], &vec); 568 CeedChkBackend(ierr); 569 // Check if active input 570 if (vec == CEED_VECTOR_ACTIVE) { 571 ierr = CeedQFunctionFieldGetSize(qf_input_fields[i], &size); 572 CeedChkBackend(ierr); 573 ierr = CeedVectorSetValue(impl->q_vecs_in[i], 0.0); CeedChkBackend(ierr); 574 ierr = CeedVectorGetArray(impl->q_vecs_in[i], CEED_MEM_HOST, &tmp); 575 CeedChkBackend(ierr); 576 ierr = CeedRealloc(num_active_in + size, &active_in); CeedChkBackend(ierr); 577 for (CeedInt field=0; field<size; field++) { 578 ierr = CeedVectorCreate(ceed, Q, &active_in[num_active_in+field]); 579 CeedChkBackend(ierr); 580 ierr = CeedVectorSetArray(active_in[num_active_in+field], CEED_MEM_HOST, 581 CEED_USE_POINTER, &tmp[field*Q]); 582 CeedChkBackend(ierr); 583 } 584 num_active_in += size; 585 ierr = CeedVectorRestoreArray(impl->q_vecs_in[i], &tmp); CeedChkBackend(ierr); 586 } 587 } 588 589 // Count number of active output fields 590 for (CeedInt i=0; i<num_output_fields; i++) { 591 // Get output vector 592 ierr = CeedOperatorFieldGetVector(op_output_fields[i], &vec); 593 CeedChkBackend(ierr); 594 // Check if active output 595 if (vec == CEED_VECTOR_ACTIVE) { 596 ierr = CeedQFunctionFieldGetSize(qf_output_fields[i], &size); 597 CeedChkBackend(ierr); 598 num_active_out += size; 599 } 600 } 601 602 // Check sizes 603 if (!num_active_in || !num_active_out) 604 // LCOV_EXCL_START 605 return CeedError(ceed, CEED_ERROR_BACKEND, 606 "Cannot assemble QFunction without active inputs " 607 "and outputs"); 608 // LCOV_EXCL_STOP 609 610 // Create output restriction 611 CeedInt strides[3] = {1, Q, num_active_in*num_active_out*Q}; /* *NOPAD* */ 612 ierr = CeedElemRestrictionCreateStrided(ceed_parent, num_elem, Q, 613 num_active_in*num_active_out, 614 num_active_in*num_active_out*num_elem*Q, 615 strides, rstr); CeedChkBackend(ierr); 616 // Create assembled vector 617 ierr = CeedVectorCreate(ceed_parent, num_elem*Q*num_active_in*num_active_out, 618 assembled); CeedChkBackend(ierr); 619 ierr = CeedVectorSetValue(*assembled, 0.0); CeedChkBackend(ierr); 620 ierr = CeedVectorGetArray(*assembled, CEED_MEM_HOST, &a); CeedChkBackend(ierr); 621 622 // Loop through elements 623 for (CeedInt e=0; e<num_elem; e++) { 624 // Input basis apply 625 ierr = CeedOperatorInputBasis_Ref(e, Q, qf_input_fields, op_input_fields, 626 num_input_fields, true, impl); 627 CeedChkBackend(ierr); 628 629 // Assemble QFunction 630 for (CeedInt in=0; in<num_active_in; in++) { 631 // Set Inputs 632 ierr = CeedVectorSetValue(active_in[in], 1.0); CeedChkBackend(ierr); 633 if (num_active_in > 1) { 634 ierr = CeedVectorSetValue(active_in[(in+num_active_in-1)%num_active_in], 635 0.0); CeedChkBackend(ierr); 636 } 637 // Set Outputs 638 for (CeedInt out=0; out<num_output_fields; out++) { 639 // Get output vector 640 ierr = CeedOperatorFieldGetVector(op_output_fields[out], &vec); 641 CeedChkBackend(ierr); 642 // Check if active output 643 if (vec == CEED_VECTOR_ACTIVE) { 644 CeedVectorSetArray(impl->q_vecs_out[out], CEED_MEM_HOST, 645 CEED_USE_POINTER, a); CeedChkBackend(ierr); 646 ierr = CeedQFunctionFieldGetSize(qf_output_fields[out], &size); 647 CeedChkBackend(ierr); 648 a += size*Q; // Advance the pointer by the size of the output 649 } 650 } 651 // Apply QFunction 652 ierr = CeedQFunctionApply(qf, Q, impl->q_vecs_in, impl->q_vecs_out); 653 CeedChkBackend(ierr); 654 } 655 } 656 657 // Un-set output Qvecs to prevent accidental overwrite of Assembled 658 for (CeedInt out=0; out<num_output_fields; out++) { 659 // Get output vector 660 ierr = CeedOperatorFieldGetVector(op_output_fields[out], &vec); 661 CeedChkBackend(ierr); 662 // Check if active output 663 if (vec == CEED_VECTOR_ACTIVE) { 664 CeedVectorTakeArray(impl->q_vecs_out[out], CEED_MEM_HOST, NULL); 665 CeedChkBackend(ierr); 666 } 667 } 668 669 // Restore input arrays 670 ierr = CeedOperatorRestoreInputs_Ref(num_input_fields, qf_input_fields, 671 op_input_fields, true, impl); 672 CeedChkBackend(ierr); 673 674 // Restore output 675 ierr = CeedVectorRestoreArray(*assembled, &a); CeedChkBackend(ierr); 676 677 // Cleanup 678 for (CeedInt i=0; i<num_active_in; i++) { 679 ierr = CeedVectorDestroy(&active_in[i]); CeedChkBackend(ierr); 680 } 681 ierr = CeedFree(&active_in); CeedChkBackend(ierr); 682 683 return CEED_ERROR_SUCCESS; 684 } 685 686 //------------------------------------------------------------------------------ 687 // Operator Destroy 688 //------------------------------------------------------------------------------ 689 static int CeedOperatorDestroy_Ref(CeedOperator op) { 690 int ierr; 691 CeedOperator_Ref *impl; 692 ierr = CeedOperatorGetData(op, &impl); CeedChkBackend(ierr); 693 694 for (CeedInt i=0; i<impl->num_e_vecs_in+impl->num_e_vecs_out; i++) { 695 ierr = CeedVectorDestroy(&impl->e_vecs[i]); CeedChkBackend(ierr); 696 } 697 ierr = CeedFree(&impl->e_vecs); CeedChkBackend(ierr); 698 ierr = CeedFree(&impl->e_data); CeedChkBackend(ierr); 699 ierr = CeedFree(&impl->input_state); CeedChkBackend(ierr); 700 701 for (CeedInt i=0; i<impl->num_e_vecs_in; i++) { 702 ierr = CeedVectorDestroy(&impl->e_vecs_in[i]); CeedChkBackend(ierr); 703 ierr = CeedVectorDestroy(&impl->q_vecs_in[i]); CeedChkBackend(ierr); 704 } 705 ierr = CeedFree(&impl->e_vecs_in); CeedChkBackend(ierr); 706 ierr = CeedFree(&impl->q_vecs_in); CeedChkBackend(ierr); 707 708 for (CeedInt i=0; i<impl->num_e_vecs_out; i++) { 709 ierr = CeedVectorDestroy(&impl->e_vecs_out[i]); CeedChkBackend(ierr); 710 ierr = CeedVectorDestroy(&impl->q_vecs_out[i]); CeedChkBackend(ierr); 711 } 712 ierr = CeedFree(&impl->e_vecs_out); CeedChkBackend(ierr); 713 ierr = CeedFree(&impl->q_vecs_out); CeedChkBackend(ierr); 714 715 ierr = CeedFree(&impl); CeedChkBackend(ierr); 716 return CEED_ERROR_SUCCESS; 717 } 718 719 //------------------------------------------------------------------------------ 720 // Operator Create 721 //------------------------------------------------------------------------------ 722 int CeedOperatorCreate_Ref(CeedOperator op) { 723 int ierr; 724 Ceed ceed; 725 ierr = CeedOperatorGetCeed(op, &ceed); CeedChkBackend(ierr); 726 CeedOperator_Ref *impl; 727 728 ierr = CeedCalloc(1, &impl); CeedChkBackend(ierr); 729 ierr = CeedOperatorSetData(op, impl); CeedChkBackend(ierr); 730 731 ierr = CeedSetBackendFunction(ceed, "Operator", op, "LinearAssembleQFunction", 732 CeedOperatorLinearAssembleQFunction_Ref); 733 CeedChkBackend(ierr); 734 ierr = CeedSetBackendFunction(ceed, "Operator", op, "ApplyAdd", 735 CeedOperatorApplyAdd_Ref); CeedChkBackend(ierr); 736 ierr = CeedSetBackendFunction(ceed, "Operator", op, "Destroy", 737 CeedOperatorDestroy_Ref); CeedChkBackend(ierr); 738 return CEED_ERROR_SUCCESS; 739 } 740