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 // Get Basis Emode Pointer 684 //------------------------------------------------------------------------------ 685 static inline void CeedOperatorGetBasisPointer_Ref(const CeedScalar **basis_ptr, 686 CeedEvalMode eval_mode, const CeedScalar *identity, const CeedScalar *interp, 687 const CeedScalar *grad) { 688 switch (eval_mode) { 689 case CEED_EVAL_NONE: 690 *basis_ptr = identity; 691 break; 692 case CEED_EVAL_INTERP: 693 *basis_ptr = interp; 694 break; 695 case CEED_EVAL_GRAD: 696 *basis_ptr = grad; 697 break; 698 case CEED_EVAL_WEIGHT: 699 case CEED_EVAL_DIV: 700 case CEED_EVAL_CURL: 701 break; // Caught by QF Assembly 702 } 703 } 704 705 //------------------------------------------------------------------------------ 706 // Create point block restriction 707 //------------------------------------------------------------------------------ 708 static int CreatePBRestriction_Ref(CeedElemRestriction rstr, 709 CeedElemRestriction *pb_rstr) { 710 int ierr; 711 Ceed ceed; 712 ierr = CeedElemRestrictionGetCeed(rstr, &ceed); CeedChkBackend(ierr); 713 const CeedInt *offsets; 714 ierr = CeedElemRestrictionGetOffsets(rstr, CEED_MEM_HOST, &offsets); 715 CeedChkBackend(ierr); 716 717 // Expand offsets 718 CeedInt num_elem, num_comp, elem_size, comp_stride, max = 1, *pbOffsets; 719 ierr = CeedElemRestrictionGetNumElements(rstr, &num_elem); CeedChkBackend(ierr); 720 ierr = CeedElemRestrictionGetNumComponents(rstr, &num_comp); 721 CeedChkBackend(ierr); 722 ierr = CeedElemRestrictionGetElementSize(rstr, &elem_size); 723 CeedChkBackend(ierr); 724 ierr = CeedElemRestrictionGetCompStride(rstr, &comp_stride); 725 CeedChkBackend(ierr); 726 CeedInt shift = num_comp; 727 if (comp_stride != 1) 728 shift *= num_comp; 729 ierr = CeedCalloc(num_elem*elem_size, &pbOffsets); CeedChkBackend(ierr); 730 for (CeedInt i = 0; i < num_elem*elem_size; i++) { 731 pbOffsets[i] = offsets[i]*shift; 732 if (pbOffsets[i] > max) 733 max = pbOffsets[i]; 734 } 735 736 // Create new restriction 737 ierr = CeedElemRestrictionCreate(ceed, num_elem, elem_size, num_comp*num_comp, 738 1, 739 max + num_comp*num_comp, CEED_MEM_HOST, 740 CEED_OWN_POINTER, pbOffsets, pb_rstr); 741 CeedChkBackend(ierr); 742 743 // Cleanup 744 ierr = CeedElemRestrictionRestoreOffsets(rstr, &offsets); CeedChkBackend(ierr); 745 746 return CEED_ERROR_SUCCESS; 747 } 748 749 //------------------------------------------------------------------------------ 750 // Assemble diagonal common code 751 //------------------------------------------------------------------------------ 752 static inline int CeedOperatorAssembleAddDiagonalCore_Ref(CeedOperator op, 753 CeedVector assembled, CeedRequest *request, const bool is_pointblock) { 754 int ierr; 755 Ceed ceed; 756 ierr = CeedOperatorGetCeed(op, &ceed); CeedChkBackend(ierr); 757 758 // Assemble QFunction 759 CeedQFunction qf; 760 ierr = CeedOperatorGetQFunction(op, &qf); CeedChkBackend(ierr); 761 CeedInt num_input_fields, num_output_fields; 762 ierr= CeedQFunctionGetNumArgs(qf, &num_input_fields, &num_output_fields); 763 CeedChkBackend(ierr); 764 CeedVector assembled_qf; 765 CeedElemRestriction rstr; 766 ierr = CeedOperatorLinearAssembleQFunction(op, &assembled_qf, &rstr, request); 767 CeedChkBackend(ierr); 768 ierr = CeedElemRestrictionDestroy(&rstr); CeedChkBackend(ierr); 769 CeedScalar max_norm = 0; 770 ierr = CeedVectorNorm(assembled_qf, CEED_NORM_MAX, &max_norm); 771 CeedChkBackend(ierr); 772 773 // Determine active input basis 774 CeedOperatorField *op_fields; 775 CeedQFunctionField *qf_fields; 776 ierr = CeedOperatorGetFields(op, &op_fields, NULL); CeedChkBackend(ierr); 777 ierr = CeedQFunctionGetFields(qf, &qf_fields, NULL); CeedChkBackend(ierr); 778 CeedInt num_eval_mode_in = 0, num_comp, dim = 1; 779 CeedEvalMode *eval_mode_in = NULL; 780 CeedBasis basis_in = NULL; 781 CeedElemRestriction rstr_in = NULL; 782 for (CeedInt i=0; i<num_input_fields; i++) { 783 CeedVector vec; 784 ierr = CeedOperatorFieldGetVector(op_fields[i], &vec); CeedChkBackend(ierr); 785 if (vec == CEED_VECTOR_ACTIVE) { 786 CeedElemRestriction rstr; 787 ierr = CeedOperatorFieldGetBasis(op_fields[i], &basis_in); CeedChkBackend(ierr); 788 ierr = CeedBasisGetNumComponents(basis_in, &num_comp); CeedChkBackend(ierr); 789 ierr = CeedBasisGetDimension(basis_in, &dim); CeedChkBackend(ierr); 790 ierr = CeedOperatorFieldGetElemRestriction(op_fields[i], &rstr); 791 CeedChkBackend(ierr); 792 if (rstr_in && rstr_in != rstr) 793 // LCOV_EXCL_START 794 return CeedError(ceed, CEED_ERROR_BACKEND, 795 "Multi-field non-composite operator diagonal assembly not supported"); 796 // LCOV_EXCL_STOP 797 rstr_in = rstr; 798 CeedEvalMode eval_mode; 799 ierr = CeedQFunctionFieldGetEvalMode(qf_fields[i], &eval_mode); 800 CeedChkBackend(ierr); 801 switch (eval_mode) { 802 case CEED_EVAL_NONE: 803 case CEED_EVAL_INTERP: 804 ierr = CeedRealloc(num_eval_mode_in + 1, &eval_mode_in); CeedChkBackend(ierr); 805 eval_mode_in[num_eval_mode_in] = eval_mode; 806 num_eval_mode_in += 1; 807 break; 808 case CEED_EVAL_GRAD: 809 ierr = CeedRealloc(num_eval_mode_in + dim, &eval_mode_in); CeedChkBackend(ierr); 810 for (CeedInt d=0; d<dim; d++) 811 eval_mode_in[num_eval_mode_in+d] = eval_mode; 812 num_eval_mode_in += dim; 813 break; 814 case CEED_EVAL_WEIGHT: 815 case CEED_EVAL_DIV: 816 case CEED_EVAL_CURL: 817 break; // Caught by QF Assembly 818 } 819 } 820 } 821 822 // Determine active output basis 823 ierr = CeedOperatorGetFields(op, NULL, &op_fields); CeedChkBackend(ierr); 824 ierr = CeedQFunctionGetFields(qf, NULL, &qf_fields); CeedChkBackend(ierr); 825 CeedInt num_eval_mode_out = 0; 826 CeedEvalMode *eval_mode_out = NULL; 827 CeedBasis basis_out = NULL; 828 CeedElemRestriction rstr_out = NULL; 829 for (CeedInt i=0; i<num_output_fields; i++) { 830 CeedVector vec; 831 ierr = CeedOperatorFieldGetVector(op_fields[i], &vec); CeedChkBackend(ierr); 832 if (vec == CEED_VECTOR_ACTIVE) { 833 CeedElemRestriction rstr; 834 ierr = CeedOperatorFieldGetBasis(op_fields[i], &basis_out); 835 CeedChkBackend(ierr); 836 ierr = CeedOperatorFieldGetElemRestriction(op_fields[i], &rstr); 837 CeedChkBackend(ierr); 838 if (rstr_out && rstr_out != rstr) 839 // LCOV_EXCL_START 840 return CeedError(ceed, CEED_ERROR_BACKEND, 841 "Multi-field non-composite operator diagonal assembly not supported"); 842 // LCOV_EXCL_STOP 843 rstr_out = rstr; 844 CeedEvalMode eval_mode; 845 ierr = CeedQFunctionFieldGetEvalMode(qf_fields[i], &eval_mode); 846 CeedChkBackend(ierr); 847 switch (eval_mode) { 848 case CEED_EVAL_NONE: 849 case CEED_EVAL_INTERP: 850 ierr = CeedRealloc(num_eval_mode_out + 1, &eval_mode_out); CeedChkBackend(ierr); 851 eval_mode_out[num_eval_mode_out] = eval_mode; 852 num_eval_mode_out += 1; 853 break; 854 case CEED_EVAL_GRAD: 855 ierr = CeedRealloc(num_eval_mode_out + dim, &eval_mode_out); 856 CeedChkBackend(ierr); 857 for (CeedInt d=0; d<dim; d++) 858 eval_mode_out[num_eval_mode_out+d] = eval_mode; 859 num_eval_mode_out += dim; 860 break; 861 case CEED_EVAL_WEIGHT: 862 case CEED_EVAL_DIV: 863 case CEED_EVAL_CURL: 864 break; // Caught by QF Assembly 865 } 866 } 867 } 868 869 // Assemble point-block diagonal restriction, if needed 870 CeedElemRestriction diag_rstr = rstr_out; 871 if (is_pointblock) { 872 ierr = CreatePBRestriction_Ref(rstr_out, &diag_rstr); CeedChkBackend(ierr); 873 } 874 875 // Create diagonal vector 876 CeedVector elem_diag; 877 ierr = CeedElemRestrictionCreateVector(diag_rstr, NULL, &elem_diag); 878 CeedChkBackend(ierr); 879 880 // Assemble element operator diagonals 881 CeedScalar *elem_diag_array, *assembled_qf_array; 882 ierr = CeedVectorSetValue(elem_diag, 0.0); CeedChkBackend(ierr); 883 ierr = CeedVectorGetArray(elem_diag, CEED_MEM_HOST, &elem_diag_array); 884 CeedChkBackend(ierr); 885 ierr = CeedVectorGetArray(assembled_qf, CEED_MEM_HOST, &assembled_qf_array); 886 CeedChkBackend(ierr); 887 CeedInt num_elem, num_nodes, num_qpts; 888 ierr = CeedElemRestrictionGetNumElements(diag_rstr, &num_elem); 889 CeedChkBackend(ierr); 890 ierr = CeedBasisGetNumNodes(basis_in, &num_nodes); CeedChkBackend(ierr); 891 ierr = CeedBasisGetNumQuadraturePoints(basis_in, &num_qpts); 892 CeedChkBackend(ierr); 893 // Basis matrices 894 const CeedScalar *interp_in, *interp_out, *grad_in, *grad_out; 895 CeedScalar *identity = NULL; 896 bool evalNone = false; 897 for (CeedInt i=0; i<num_eval_mode_in; i++) 898 evalNone = evalNone || (eval_mode_in[i] == CEED_EVAL_NONE); 899 for (CeedInt i=0; i<num_eval_mode_out; i++) 900 evalNone = evalNone || (eval_mode_out[i] == CEED_EVAL_NONE); 901 if (evalNone) { 902 ierr = CeedCalloc(num_qpts*num_nodes, &identity); CeedChkBackend(ierr); 903 for (CeedInt i=0; i<(num_nodes<num_qpts?num_nodes:num_qpts); i++) 904 identity[i*num_nodes+i] = 1.0; 905 } 906 ierr = CeedBasisGetInterp(basis_in, &interp_in); CeedChkBackend(ierr); 907 ierr = CeedBasisGetInterp(basis_out, &interp_out); CeedChkBackend(ierr); 908 ierr = CeedBasisGetGrad(basis_in, &grad_in); CeedChkBackend(ierr); 909 ierr = CeedBasisGetGrad(basis_out, &grad_out); CeedChkBackend(ierr); 910 // Compute the diagonal of B^T D B 911 // Each element 912 const CeedScalar qf_value_bound = max_norm*1e-12; 913 for (CeedInt e=0; e<num_elem; e++) { 914 CeedInt d_out = -1; 915 // Each basis eval mode pair 916 for (CeedInt e_out=0; e_out<num_eval_mode_out; e_out++) { 917 const CeedScalar *bt = NULL; 918 if (eval_mode_out[e_out] == CEED_EVAL_GRAD) 919 d_out += 1; 920 CeedOperatorGetBasisPointer_Ref(&bt, eval_mode_out[e_out], identity, interp_out, 921 &grad_out[d_out*num_qpts*num_nodes]); 922 CeedInt d_in = -1; 923 for (CeedInt e_in=0; e_in<num_eval_mode_in; e_in++) { 924 const CeedScalar *b = NULL; 925 if (eval_mode_in[e_in] == CEED_EVAL_GRAD) 926 d_in += 1; 927 CeedOperatorGetBasisPointer_Ref(&b, eval_mode_in[e_in], identity, interp_in, 928 &grad_in[d_in*num_qpts*num_nodes]); 929 // Each component 930 for (CeedInt c_out=0; c_out<num_comp; c_out++) 931 // Each qpoint/node pair 932 for (CeedInt q=0; q<num_qpts; q++) 933 if (is_pointblock) { 934 // Point Block Diagonal 935 for (CeedInt c_in=0; c_in<num_comp; c_in++) { 936 const CeedScalar qf_value = 937 assembled_qf_array[((((e*num_eval_mode_in+e_in)*num_comp+c_in)* 938 num_eval_mode_out+e_out)*num_comp+c_out)*num_qpts+q]; 939 if (fabs(qf_value) > qf_value_bound) 940 for (CeedInt n=0; n<num_nodes; n++) 941 elem_diag_array[((e*num_comp+c_out)*num_comp+c_in)*num_nodes+n] += 942 bt[q*num_nodes+n] * qf_value * b[q*num_nodes+n]; 943 } 944 } else { 945 // Diagonal Only 946 const CeedScalar qf_value = 947 assembled_qf_array[((((e*num_eval_mode_in+e_in)*num_comp+c_out)* 948 num_eval_mode_out+e_out)*num_comp+c_out)*num_qpts+q]; 949 if (fabs(qf_value) > qf_value_bound) 950 for (CeedInt n=0; n<num_nodes; n++) 951 elem_diag_array[(e*num_comp+c_out)*num_nodes+n] += 952 bt[q*num_nodes+n] * qf_value * b[q*num_nodes+n]; 953 } 954 } 955 } 956 } 957 ierr = CeedVectorRestoreArray(elem_diag, &elem_diag_array); 958 CeedChkBackend(ierr); 959 ierr = CeedVectorRestoreArray(assembled_qf, &assembled_qf_array); 960 CeedChkBackend(ierr); 961 962 // Assemble local operator diagonal 963 ierr = CeedElemRestrictionApply(diag_rstr, CEED_TRANSPOSE, elem_diag, 964 assembled, request); CeedChkBackend(ierr); 965 966 // Cleanup 967 if (is_pointblock) { 968 ierr = CeedElemRestrictionDestroy(&diag_rstr); CeedChkBackend(ierr); 969 } 970 ierr = CeedVectorDestroy(&assembled_qf); CeedChkBackend(ierr); 971 ierr = CeedVectorDestroy(&elem_diag); CeedChkBackend(ierr); 972 ierr = CeedFree(&eval_mode_in); CeedChkBackend(ierr); 973 ierr = CeedFree(&eval_mode_out); CeedChkBackend(ierr); 974 ierr = CeedFree(&identity); CeedChkBackend(ierr); 975 976 return CEED_ERROR_SUCCESS; 977 } 978 979 //------------------------------------------------------------------------------ 980 // Assemble composite diagonal common code 981 //------------------------------------------------------------------------------ 982 static inline int CeedOperatorLinearAssembleAddDiagonalCompositeCore_Ref( 983 CeedOperator op, CeedVector assembled, CeedRequest *request, 984 const bool is_pointblock) { 985 int ierr; 986 CeedInt num_sub; 987 CeedOperator *suboperators; 988 ierr = CeedOperatorGetNumSub(op, &num_sub); CeedChkBackend(ierr); 989 ierr = CeedOperatorGetSubList(op, &suboperators); CeedChkBackend(ierr); 990 for (CeedInt i = 0; i < num_sub; i++) { 991 ierr = CeedOperatorAssembleAddDiagonalCore_Ref(suboperators[i], assembled, 992 request, is_pointblock); CeedChkBackend(ierr); 993 } 994 return CEED_ERROR_SUCCESS; 995 } 996 997 //------------------------------------------------------------------------------ 998 // Assemble Linear Diagonal 999 //------------------------------------------------------------------------------ 1000 static int CeedOperatorLinearAssembleAddDiagonal_Ref(CeedOperator op, 1001 CeedVector assembled, CeedRequest *request) { 1002 int ierr; 1003 bool is_composite; 1004 ierr = CeedOperatorIsComposite(op, &is_composite); CeedChkBackend(ierr); 1005 if (is_composite) { 1006 return CeedOperatorLinearAssembleAddDiagonalCompositeCore_Ref(op, assembled, 1007 request, false); 1008 } else { 1009 return CeedOperatorAssembleAddDiagonalCore_Ref(op, assembled, request, false); 1010 } 1011 } 1012 1013 //------------------------------------------------------------------------------ 1014 // Assemble Linear Point Block Diagonal 1015 //------------------------------------------------------------------------------ 1016 static int CeedOperatorLinearAssembleAddPointBlockDiagonal_Ref(CeedOperator op, 1017 CeedVector assembled, CeedRequest *request) { 1018 int ierr; 1019 bool is_composite; 1020 ierr = CeedOperatorIsComposite(op, &is_composite); CeedChkBackend(ierr); 1021 if (is_composite) { 1022 return CeedOperatorLinearAssembleAddDiagonalCompositeCore_Ref(op, assembled, 1023 request, true); 1024 } else { 1025 return CeedOperatorAssembleAddDiagonalCore_Ref(op, assembled, request, true); 1026 } 1027 } 1028 1029 //------------------------------------------------------------------------------ 1030 // Create FDM Element Inverse 1031 //------------------------------------------------------------------------------ 1032 int CeedOperatorCreateFDMElementInverse_Ref(CeedOperator op, 1033 CeedOperator *fdm_inv, CeedRequest *request) { 1034 int ierr; 1035 Ceed ceed, ceed_parent; 1036 ierr = CeedOperatorGetCeed(op, &ceed); CeedChkBackend(ierr); 1037 ierr = CeedGetOperatorFallbackParentCeed(ceed, &ceed_parent); 1038 CeedChkBackend(ierr); 1039 ceed_parent = ceed_parent ? ceed_parent : ceed; 1040 CeedQFunction qf; 1041 ierr = CeedOperatorGetQFunction(op, &qf); CeedChkBackend(ierr); 1042 1043 // Determine active input basis 1044 bool interp = false, grad = false; 1045 CeedBasis basis = NULL; 1046 CeedElemRestriction rstr = NULL; 1047 CeedOperatorField *op_fields; 1048 CeedQFunctionField *qf_fields; 1049 ierr = CeedOperatorGetFields(op, &op_fields, NULL); CeedChkBackend(ierr); 1050 ierr = CeedQFunctionGetFields(qf, &qf_fields, NULL); CeedChkBackend(ierr); 1051 CeedInt num_input_fields; 1052 ierr = CeedQFunctionGetNumArgs(qf, &num_input_fields, NULL); 1053 CeedChkBackend(ierr); 1054 for (CeedInt i=0; i<num_input_fields; i++) { 1055 CeedVector vec; 1056 ierr = CeedOperatorFieldGetVector(op_fields[i], &vec); CeedChkBackend(ierr); 1057 if (vec == CEED_VECTOR_ACTIVE) { 1058 CeedEvalMode eval_mode; 1059 ierr = CeedQFunctionFieldGetEvalMode(qf_fields[i], &eval_mode); 1060 CeedChkBackend(ierr); 1061 interp = interp || eval_mode == CEED_EVAL_INTERP; 1062 grad = grad || eval_mode == CEED_EVAL_GRAD; 1063 ierr = CeedOperatorFieldGetBasis(op_fields[i], &basis); CeedChkBackend(ierr); 1064 ierr = CeedOperatorFieldGetElemRestriction(op_fields[i], &rstr); 1065 CeedChkBackend(ierr); 1066 } 1067 } 1068 if (!basis) 1069 // LCOV_EXCL_START 1070 return CeedError(ceed, CEED_ERROR_BACKEND, "No active field set"); 1071 // LCOV_EXCL_STOP 1072 CeedInt P_1d, Q_1d, elem_size, num_qpts, dim, num_comp = 1, num_elem = 1, 1073 l_size = 1; 1074 ierr = CeedBasisGetNumNodes1D(basis, &P_1d); CeedChkBackend(ierr); 1075 ierr = CeedBasisGetNumNodes(basis, &elem_size); CeedChkBackend(ierr); 1076 ierr = CeedBasisGetNumQuadraturePoints1D(basis, &Q_1d); CeedChkBackend(ierr); 1077 ierr = CeedBasisGetNumQuadraturePoints(basis, &num_qpts); CeedChkBackend(ierr); 1078 ierr = CeedBasisGetDimension(basis, &dim); CeedChkBackend(ierr); 1079 ierr = CeedBasisGetNumComponents(basis, &num_comp); CeedChkBackend(ierr); 1080 ierr = CeedElemRestrictionGetNumElements(rstr, &num_elem); CeedChkBackend(ierr); 1081 ierr = CeedElemRestrictionGetLVectorSize(rstr, &l_size); CeedChkBackend(ierr); 1082 1083 // Build and diagonalize 1D Mass and Laplacian 1084 bool tensor_basis; 1085 ierr = CeedBasisIsTensor(basis, &tensor_basis); CeedChkBackend(ierr); 1086 if (!tensor_basis) 1087 // LCOV_EXCL_START 1088 return CeedError(ceed, CEED_ERROR_BACKEND, 1089 "FDMElementInverse only supported for tensor " 1090 "bases"); 1091 // LCOV_EXCL_STOP 1092 CeedScalar *work, *mass, *laplace, *x, *x2, *lambda; 1093 ierr = CeedMalloc(Q_1d*P_1d, &work); CeedChkBackend(ierr); 1094 ierr = CeedMalloc(P_1d*P_1d, &mass); CeedChkBackend(ierr); 1095 ierr = CeedMalloc(P_1d*P_1d, &laplace); CeedChkBackend(ierr); 1096 ierr = CeedMalloc(P_1d*P_1d, &x); CeedChkBackend(ierr); 1097 ierr = CeedMalloc(P_1d*P_1d, &x2); CeedChkBackend(ierr); 1098 ierr = CeedMalloc(P_1d, &lambda); CeedChkBackend(ierr); 1099 // -- Mass 1100 const CeedScalar *interp_1d, *grad_1d, *q_weight_1d; 1101 ierr = CeedBasisGetInterp1D(basis, &interp_1d); CeedChkBackend(ierr); 1102 ierr = CeedBasisGetGrad1D(basis, &grad_1d); CeedChkBackend(ierr); 1103 ierr = CeedBasisGetQWeights(basis, &q_weight_1d); CeedChkBackend(ierr); 1104 for (CeedInt i=0; i<Q_1d; i++) 1105 for (CeedInt j=0; j<P_1d; j++) 1106 work[i+j*Q_1d] = interp_1d[i*P_1d+j]*q_weight_1d[i]; 1107 ierr = CeedMatrixMultiply(ceed, (const CeedScalar *)work, 1108 (const CeedScalar *)interp_1d, mass, P_1d, P_1d, Q_1d); 1109 CeedChkBackend(ierr); 1110 // -- Laplacian 1111 for (CeedInt i=0; i<Q_1d; i++) 1112 for (CeedInt j=0; j<P_1d; j++) 1113 work[i+j*Q_1d] = grad_1d[i*P_1d+j]*q_weight_1d[i]; 1114 ierr = CeedMatrixMultiply(ceed, (const CeedScalar *)work, 1115 (const CeedScalar *)grad_1d, laplace, P_1d, P_1d, Q_1d); 1116 CeedChkBackend(ierr); 1117 // -- Diagonalize 1118 ierr = CeedSimultaneousDiagonalization(ceed, laplace, mass, x, lambda, P_1d); 1119 CeedChkBackend(ierr); 1120 ierr = CeedFree(&work); CeedChkBackend(ierr); 1121 ierr = CeedFree(&mass); CeedChkBackend(ierr); 1122 ierr = CeedFree(&laplace); CeedChkBackend(ierr); 1123 for (CeedInt i=0; i<P_1d; i++) 1124 for (CeedInt j=0; j<P_1d; j++) 1125 x2[i+j*P_1d] = x[j+i*P_1d]; 1126 ierr = CeedFree(&x); CeedChkBackend(ierr); 1127 1128 // Assemble QFunction 1129 CeedVector assembled; 1130 CeedElemRestriction rstr_qf; 1131 ierr = CeedOperatorLinearAssembleQFunction(op, &assembled, &rstr_qf, 1132 request); CeedChkBackend(ierr); 1133 ierr = CeedElemRestrictionDestroy(&rstr_qf); CeedChkBackend(ierr); 1134 CeedScalar max_norm = 0; 1135 ierr = CeedVectorNorm(assembled, CEED_NORM_MAX, &max_norm); 1136 CeedChkBackend(ierr); 1137 1138 // Calculate element averages 1139 CeedInt num_fields = ((interp?1:0) + (grad?dim:0))*((interp?1:0) + 1140 (grad?dim:0)); 1141 CeedScalar *elem_avg; 1142 const CeedScalar *assembledarray, *q_weight_array; 1143 CeedVector q_weight; 1144 ierr = CeedVectorCreate(ceed_parent, num_qpts, &q_weight); CeedChkBackend(ierr); 1145 ierr = CeedBasisApply(basis, 1, CEED_NOTRANSPOSE, CEED_EVAL_WEIGHT, 1146 CEED_VECTOR_NONE, q_weight); CeedChkBackend(ierr); 1147 ierr = CeedVectorGetArrayRead(assembled, CEED_MEM_HOST, &assembledarray); 1148 CeedChkBackend(ierr); 1149 ierr = CeedVectorGetArrayRead(q_weight, CEED_MEM_HOST, &q_weight_array); 1150 CeedChkBackend(ierr); 1151 ierr = CeedCalloc(num_elem, &elem_avg); CeedChkBackend(ierr); 1152 for (CeedInt e=0; e<num_elem; e++) { 1153 CeedInt count = 0; 1154 for (CeedInt q=0; q<num_qpts; q++) 1155 for (CeedInt i=0; i<num_comp*num_comp*num_fields; i++) 1156 if (fabs(assembledarray[e*num_elem*num_qpts*num_comp*num_comp*num_fields + 1157 i*num_qpts + q]) > max_norm*1e-12) { 1158 elem_avg[e] += assembledarray[e*num_elem*num_qpts*num_comp*num_comp*num_fields + 1159 i*num_qpts + q] / q_weight_array[q]; 1160 count++; 1161 } 1162 if (count) 1163 elem_avg[e] /= count; 1164 } 1165 ierr = CeedVectorRestoreArrayRead(assembled, &assembledarray); 1166 CeedChkBackend(ierr); 1167 ierr = CeedVectorDestroy(&assembled); CeedChkBackend(ierr); 1168 ierr = CeedVectorRestoreArrayRead(q_weight, &q_weight_array); 1169 CeedChkBackend(ierr); 1170 ierr = CeedVectorDestroy(&q_weight); CeedChkBackend(ierr); 1171 1172 // Build FDM diagonal 1173 CeedVector q_data; 1174 CeedScalar *q_data_array; 1175 ierr = CeedVectorCreate(ceed_parent, num_elem*num_comp*l_size, &q_data); 1176 CeedChkBackend(ierr); 1177 ierr = CeedVectorSetArray(q_data, CEED_MEM_HOST, CEED_COPY_VALUES, NULL); 1178 CeedChkBackend(ierr); 1179 ierr = CeedVectorGetArray(q_data, CEED_MEM_HOST, &q_data_array); 1180 CeedChkBackend(ierr); 1181 for (CeedInt e=0; e<num_elem; e++) 1182 for (CeedInt c=0; c<num_comp; c++) 1183 for (CeedInt n=0; n<l_size; n++) { 1184 if (interp) 1185 q_data_array[(e*num_comp+c)*l_size+n] = 1; 1186 if (grad) 1187 for (CeedInt d=0; d<dim; d++) { 1188 CeedInt i = (n / CeedIntPow(P_1d, d)) % P_1d; 1189 q_data_array[(e*num_comp+c)*l_size+n] += lambda[i]; 1190 } 1191 q_data_array[(e*num_comp+c)*l_size+n] = 1 / (elem_avg[e] * 1192 q_data_array[(e*num_comp+c)*l_size+n]); 1193 } 1194 ierr = CeedFree(&elem_avg); CeedChkBackend(ierr); 1195 ierr = CeedVectorRestoreArray(q_data, &q_data_array); CeedChkBackend(ierr); 1196 1197 // Setup FDM operator 1198 // -- Basis 1199 CeedBasis fdm_basis; 1200 CeedScalar *grad_dummy, *q_ref_dummy, *q_weight_dummy; 1201 ierr = CeedCalloc(P_1d*P_1d, &grad_dummy); CeedChkBackend(ierr); 1202 ierr = CeedCalloc(P_1d, &q_ref_dummy); CeedChkBackend(ierr); 1203 ierr = CeedCalloc(P_1d, &q_weight_dummy); CeedChkBackend(ierr); 1204 ierr = CeedBasisCreateTensorH1(ceed_parent, dim, num_comp, P_1d, P_1d, x2, 1205 grad_dummy, q_ref_dummy, q_weight_dummy, 1206 &fdm_basis); CeedChkBackend(ierr); 1207 ierr = CeedFree(&grad_dummy); CeedChkBackend(ierr); 1208 ierr = CeedFree(&q_ref_dummy); CeedChkBackend(ierr); 1209 ierr = CeedFree(&q_weight_dummy); CeedChkBackend(ierr); 1210 ierr = CeedFree(&x2); CeedChkBackend(ierr); 1211 ierr = CeedFree(&lambda); CeedChkBackend(ierr); 1212 1213 // -- Restriction 1214 CeedElemRestriction rstr_i; 1215 CeedInt strides[3] = {1, l_size, l_size*num_comp}; 1216 ierr = CeedElemRestrictionCreateStrided(ceed_parent, num_elem, l_size, num_comp, 1217 l_size*num_elem*num_comp, strides, &rstr_i); 1218 CeedChkBackend(ierr); 1219 // -- QFunction 1220 CeedQFunction mass_qf; 1221 ierr = CeedQFunctionCreateInteriorByName(ceed_parent, "MassApply", &mass_qf); 1222 CeedChkBackend(ierr); 1223 // -- Operator 1224 ierr = CeedOperatorCreate(ceed_parent, mass_qf, NULL, NULL, fdm_inv); 1225 CeedChkBackend(ierr); 1226 CeedOperatorSetField(*fdm_inv, "u", rstr_i, fdm_basis, CEED_VECTOR_ACTIVE); 1227 CeedChkBackend(ierr); 1228 CeedOperatorSetField(*fdm_inv, "qdata", rstr_i, CEED_BASIS_COLLOCATED, q_data); 1229 CeedChkBackend(ierr); 1230 CeedOperatorSetField(*fdm_inv, "v", rstr_i, fdm_basis, CEED_VECTOR_ACTIVE); 1231 CeedChkBackend(ierr); 1232 1233 // Cleanup 1234 ierr = CeedVectorDestroy(&q_data); CeedChkBackend(ierr); 1235 ierr = CeedBasisDestroy(&fdm_basis); CeedChkBackend(ierr); 1236 ierr = CeedElemRestrictionDestroy(&rstr_i); CeedChkBackend(ierr); 1237 ierr = CeedQFunctionDestroy(&mass_qf); CeedChkBackend(ierr); 1238 1239 return CEED_ERROR_SUCCESS; 1240 } 1241 1242 //------------------------------------------------------------------------------ 1243 // Operator Destroy 1244 //------------------------------------------------------------------------------ 1245 static int CeedOperatorDestroy_Ref(CeedOperator op) { 1246 int ierr; 1247 CeedOperator_Ref *impl; 1248 ierr = CeedOperatorGetData(op, &impl); CeedChkBackend(ierr); 1249 1250 for (CeedInt i=0; i<impl->num_e_vecs_in+impl->num_e_vecs_out; i++) { 1251 ierr = CeedVectorDestroy(&impl->e_vecs[i]); CeedChkBackend(ierr); 1252 } 1253 ierr = CeedFree(&impl->e_vecs); CeedChkBackend(ierr); 1254 ierr = CeedFree(&impl->e_data); CeedChkBackend(ierr); 1255 ierr = CeedFree(&impl->input_state); CeedChkBackend(ierr); 1256 1257 for (CeedInt i=0; i<impl->num_e_vecs_in; i++) { 1258 ierr = CeedVectorDestroy(&impl->e_vecs_in[i]); CeedChkBackend(ierr); 1259 ierr = CeedVectorDestroy(&impl->q_vecs_in[i]); CeedChkBackend(ierr); 1260 } 1261 ierr = CeedFree(&impl->e_vecs_in); CeedChkBackend(ierr); 1262 ierr = CeedFree(&impl->q_vecs_in); CeedChkBackend(ierr); 1263 1264 for (CeedInt i=0; i<impl->num_e_vecs_out; i++) { 1265 ierr = CeedVectorDestroy(&impl->e_vecs_out[i]); CeedChkBackend(ierr); 1266 ierr = CeedVectorDestroy(&impl->q_vecs_out[i]); CeedChkBackend(ierr); 1267 } 1268 ierr = CeedFree(&impl->e_vecs_out); CeedChkBackend(ierr); 1269 ierr = CeedFree(&impl->q_vecs_out); CeedChkBackend(ierr); 1270 1271 ierr = CeedFree(&impl); CeedChkBackend(ierr); 1272 return CEED_ERROR_SUCCESS; 1273 } 1274 1275 //------------------------------------------------------------------------------ 1276 // Operator Create 1277 //------------------------------------------------------------------------------ 1278 int CeedOperatorCreate_Ref(CeedOperator op) { 1279 int ierr; 1280 Ceed ceed; 1281 ierr = CeedOperatorGetCeed(op, &ceed); CeedChkBackend(ierr); 1282 CeedOperator_Ref *impl; 1283 1284 ierr = CeedCalloc(1, &impl); CeedChkBackend(ierr); 1285 ierr = CeedOperatorSetData(op, impl); CeedChkBackend(ierr); 1286 1287 ierr = CeedSetBackendFunction(ceed, "Operator", op, "LinearAssembleQFunction", 1288 CeedOperatorLinearAssembleQFunction_Ref); 1289 CeedChkBackend(ierr); 1290 ierr = CeedSetBackendFunction(ceed, "Operator", op, "LinearAssembleAddDiagonal", 1291 CeedOperatorLinearAssembleAddDiagonal_Ref); 1292 CeedChkBackend(ierr); 1293 ierr = CeedSetBackendFunction(ceed, "Operator", op, 1294 "LinearAssembleAddPointBlockDiagonal", 1295 CeedOperatorLinearAssembleAddPointBlockDiagonal_Ref); 1296 CeedChkBackend(ierr); 1297 ierr = CeedSetBackendFunction(ceed, "Operator", op, "CreateFDMElementInverse", 1298 CeedOperatorCreateFDMElementInverse_Ref); 1299 CeedChkBackend(ierr); 1300 ierr = CeedSetBackendFunction(ceed, "Operator", op, "ApplyAdd", 1301 CeedOperatorApplyAdd_Ref); CeedChkBackend(ierr); 1302 ierr = CeedSetBackendFunction(ceed, "Operator", op, "Destroy", 1303 CeedOperatorDestroy_Ref); CeedChkBackend(ierr); 1304 return CEED_ERROR_SUCCESS; 1305 } 1306 1307 //------------------------------------------------------------------------------ 1308 // Composite Operator Create 1309 //------------------------------------------------------------------------------ 1310 int CeedCompositeOperatorCreate_Ref(CeedOperator op) { 1311 int ierr; 1312 Ceed ceed; 1313 ierr = CeedOperatorGetCeed(op, &ceed); CeedChkBackend(ierr); 1314 ierr = CeedSetBackendFunction(ceed, "Operator", op, "LinearAssembleAddDiagonal", 1315 CeedOperatorLinearAssembleAddDiagonal_Ref); 1316 CeedChkBackend(ierr); 1317 ierr = CeedSetBackendFunction(ceed, "Operator", op, 1318 "LinearAssembleAddPointBlockDiagonal", 1319 CeedOperatorLinearAssembleAddPointBlockDiagonal_Ref); 1320 CeedChkBackend(ierr); 1321 return CEED_ERROR_SUCCESS; 1322 } 1323 //------------------------------------------------------------------------------ 1324