1# Copyright (c) 2017, Lawrence Livermore National Security, LLC. Produced at 2# the Lawrence Livermore National Laboratory. LLNL-CODE-734707. All Rights 3# reserved. See files LICENSE and NOTICE for details. 4# 5# This file is part of CEED, a collection of benchmarks, miniapps, software 6# libraries and APIs for efficient high-order finite element and spectral 7# element discretizations for exascale applications. For more information and 8# source code availability see http://github.com/ceed. 9# 10# The CEED research is supported by the Exascale Computing Project 17-SC-20-SC, 11# a collaborative effort of two U.S. Department of Energy organizations (Office 12# of Science and the National Nuclear Security Administration) responsible for 13# the planning and preparation of a capable exascale ecosystem, including 14# software, applications, hardware, advanced system engineering and early 15# testbed platforms, in support of the nation's exascale computing imperative. 16 17# @file 18# Test Ceed Operator functionality 19 20import os 21import libceed 22import numpy as np 23import check 24import buildmats as bm 25 26TOL = libceed.EPSILON * 256 27 28# ------------------------------------------------------------------------------- 29# Utility 30# ------------------------------------------------------------------------------- 31 32 33def load_qfs_so(): 34 from distutils.sysconfig import get_config_var 35 import ctypes 36 37 file_dir = os.path.dirname(os.path.abspath(__file__)) 38 qfs_so = os.path.join( 39 file_dir, 40 "libceed_qfunctions" + get_config_var("EXT_SUFFIX")) 41 42 # Load library 43 return ctypes.cdll.LoadLibrary(qfs_so) 44 45# ------------------------------------------------------------------------------- 46# Test creation, action, and destruction for mass matrix operator 47# ------------------------------------------------------------------------------- 48 49 50def test_500(ceed_resource): 51 ceed = libceed.Ceed(ceed_resource) 52 53 nelem = 15 54 p = 5 55 q = 8 56 nx = nelem + 1 57 nu = nelem * (p - 1) + 1 58 59 # Vectors 60 x = ceed.Vector(nx) 61 x_array = np.zeros(nx) 62 for i in range(nx): 63 x_array[i] = i / (nx - 1.0) 64 x.set_array(x_array, cmode=libceed.USE_POINTER) 65 66 qdata = ceed.Vector(nelem * q) 67 u = ceed.Vector(nu) 68 v = ceed.Vector(nu) 69 70 # Restrictions 71 indx = np.zeros(nx * 2, dtype="int32") 72 for i in range(nx): 73 indx[2 * i + 0] = i 74 indx[2 * i + 1] = i + 1 75 rx = ceed.ElemRestriction(nelem, 2, 1, 1, nx, indx, 76 cmode=libceed.USE_POINTER) 77 78 indu = np.zeros(nelem * p, dtype="int32") 79 for i in range(nelem): 80 for j in range(p): 81 indu[p * i + j] = i * (p - 1) + j 82 ru = ceed.ElemRestriction(nelem, p, 1, 1, nu, indu, 83 cmode=libceed.USE_POINTER) 84 strides = np.array([1, q, q], dtype="int32") 85 rui = ceed.StridedElemRestriction(nelem, q, 1, q * nelem, strides) 86 87 # Bases 88 bx = ceed.BasisTensorH1Lagrange(1, 1, 2, q, libceed.GAUSS) 89 bu = ceed.BasisTensorH1Lagrange(1, 1, p, q, libceed.GAUSS) 90 91 # QFunctions 92 file_dir = os.path.dirname(os.path.abspath(__file__)) 93 qfs = load_qfs_so() 94 95 qf_setup = ceed.QFunction(1, qfs.setup_mass, 96 os.path.join(file_dir, "test-qfunctions.h:setup_mass")) 97 qf_setup.add_input("weights", 1, libceed.EVAL_WEIGHT) 98 qf_setup.add_input("dx", 1, libceed.EVAL_GRAD) 99 qf_setup.add_output("rho", 1, libceed.EVAL_NONE) 100 101 qf_mass = ceed.QFunction(1, qfs.apply_mass, 102 os.path.join(file_dir, "test-qfunctions.h:apply_mass")) 103 qf_mass.add_input("rho", 1, libceed.EVAL_NONE) 104 qf_mass.add_input("u", 1, libceed.EVAL_INTERP) 105 qf_mass.add_output("v", 1, libceed.EVAL_INTERP) 106 107 # Operators 108 op_setup = ceed.Operator(qf_setup) 109 op_setup.set_field("weights", libceed.ELEMRESTRICTION_NONE, bx, 110 libceed.VECTOR_NONE) 111 op_setup.set_field("dx", rx, bx, libceed.VECTOR_ACTIVE) 112 op_setup.set_field("rho", rui, libceed.BASIS_COLLOCATED, 113 libceed.VECTOR_ACTIVE) 114 op_setup.check() 115 116 op_mass = ceed.Operator(qf_mass) 117 op_mass.set_field("rho", rui, libceed.BASIS_COLLOCATED, qdata) 118 op_mass.set_field("u", ru, bu, libceed.VECTOR_ACTIVE) 119 op_mass.set_field("v", ru, bu, libceed.VECTOR_ACTIVE) 120 op_mass.check() 121 122 # Setup 123 op_setup.apply(x, qdata) 124 125 # Apply mass matrix 126 u.set_value(0) 127 op_mass.apply(u, v) 128 129 # Check 130 with v.array_read() as v_array: 131 for i in range(q): 132 assert abs(v_array[i]) < TOL 133 134# ------------------------------------------------------------------------------- 135# Test creation, action, and destruction for mass matrix operator 136# ------------------------------------------------------------------------------- 137 138 139def test_501(ceed_resource): 140 ceed = libceed.Ceed(ceed_resource) 141 142 nelem = 15 143 p = 5 144 q = 8 145 nx = nelem + 1 146 nu = nelem * (p - 1) + 1 147 148 # Vectors 149 x = ceed.Vector(nx) 150 x_array = np.zeros(nx, dtype=ceed.scalar_type()) 151 for i in range(nx): 152 x_array[i] = i / (nx - 1.0) 153 x.set_array(x_array, cmode=libceed.USE_POINTER) 154 155 qdata = ceed.Vector(nelem * q) 156 u = ceed.Vector(nu) 157 v = ceed.Vector(nu) 158 159 # Restrictions 160 indx = np.zeros(nx * 2, dtype="int32") 161 for i in range(nx): 162 indx[2 * i + 0] = i 163 indx[2 * i + 1] = i + 1 164 rx = ceed.ElemRestriction(nelem, 2, 1, 1, nx, indx, 165 cmode=libceed.USE_POINTER) 166 167 indu = np.zeros(nelem * p, dtype="int32") 168 for i in range(nelem): 169 for j in range(p): 170 indu[p * i + j] = i * (p - 1) + j 171 ru = ceed.ElemRestriction(nelem, p, 1, 1, nu, indu, 172 cmode=libceed.USE_POINTER) 173 strides = np.array([1, q, q], dtype="int32") 174 rui = ceed.StridedElemRestriction(nelem, q, 1, q * nelem, strides) 175 176 # Bases 177 bx = ceed.BasisTensorH1Lagrange(1, 1, 2, q, libceed.GAUSS) 178 bu = ceed.BasisTensorH1Lagrange(1, 1, p, q, libceed.GAUSS) 179 180 # QFunctions 181 file_dir = os.path.dirname(os.path.abspath(__file__)) 182 qfs = load_qfs_so() 183 184 qf_setup = ceed.QFunction(1, qfs.setup_mass, 185 os.path.join(file_dir, "test-qfunctions.h:setup_mass")) 186 qf_setup.add_input("weights", 1, libceed.EVAL_WEIGHT) 187 qf_setup.add_input("dx", 1, libceed.EVAL_GRAD) 188 qf_setup.add_output("rho", 1, libceed.EVAL_NONE) 189 190 qf_mass = ceed.QFunction(1, qfs.apply_mass, 191 os.path.join(file_dir, "test-qfunctions.h:apply_mass")) 192 qf_mass.add_input("rho", 1, libceed.EVAL_NONE) 193 qf_mass.add_input("u", 1, libceed.EVAL_INTERP) 194 qf_mass.add_output("v", 1, libceed.EVAL_INTERP) 195 196 # Operators 197 op_setup = ceed.Operator(qf_setup) 198 op_setup.set_field("weights", libceed.ELEMRESTRICTION_NONE, bx, 199 libceed.VECTOR_NONE) 200 op_setup.set_field("dx", rx, bx, libceed.VECTOR_ACTIVE) 201 op_setup.set_field("rho", rui, libceed.BASIS_COLLOCATED, 202 libceed.VECTOR_ACTIVE) 203 204 op_mass = ceed.Operator(qf_mass) 205 op_mass.set_field("rho", rui, libceed.BASIS_COLLOCATED, qdata) 206 op_mass.set_field("u", ru, bu, libceed.VECTOR_ACTIVE) 207 op_mass.set_field("v", ru, bu, libceed.VECTOR_ACTIVE) 208 209 # Setup 210 op_setup.apply(x, qdata) 211 212 # Apply mass matrix 213 u.set_value(1.) 214 op_mass.apply(u, v) 215 216 # Check 217 with v.array_read() as v_array: 218 total = 0.0 219 for i in range(nu): 220 total = total + v_array[i] 221 assert abs(total - 1.0) < TOL 222 223# ------------------------------------------------------------------------------- 224# Test creation, action, and destruction for mass matrix operator with multiple 225# components 226# ------------------------------------------------------------------------------- 227 228 229def test_502(ceed_resource): 230 ceed = libceed.Ceed(ceed_resource) 231 232 nelem = 15 233 p = 5 234 q = 8 235 nx = nelem + 1 236 nu = nelem * (p - 1) + 1 237 238 # Vectors 239 x = ceed.Vector(nx) 240 x_array = np.zeros(nx, dtype=ceed.scalar_type()) 241 for i in range(nx): 242 x_array[i] = i / (nx - 1.0) 243 x.set_array(x_array, cmode=libceed.USE_POINTER) 244 245 qdata = ceed.Vector(nelem * q) 246 u = ceed.Vector(2 * nu) 247 v = ceed.Vector(2 * nu) 248 249 # Restrictions 250 indx = np.zeros(nx * 2, dtype="int32") 251 for i in range(nx): 252 indx[2 * i + 0] = i 253 indx[2 * i + 1] = i + 1 254 rx = ceed.ElemRestriction(nelem, 2, 1, 1, nx, indx, 255 cmode=libceed.USE_POINTER) 256 257 indu = np.zeros(nelem * p, dtype="int32") 258 for i in range(nelem): 259 for j in range(p): 260 indu[p * i + j] = 2 * (i * (p - 1) + j) 261 ru = ceed.ElemRestriction(nelem, p, 2, 1, 2 * nu, indu, 262 cmode=libceed.USE_POINTER) 263 strides = np.array([1, q, q], dtype="int32") 264 rui = ceed.StridedElemRestriction(nelem, q, 1, q * nelem, strides) 265 266 # Bases 267 bx = ceed.BasisTensorH1Lagrange(1, 1, 2, q, libceed.GAUSS) 268 bu = ceed.BasisTensorH1Lagrange(1, 2, p, q, libceed.GAUSS) 269 270 # QFunctions 271 file_dir = os.path.dirname(os.path.abspath(__file__)) 272 qfs = load_qfs_so() 273 274 qf_setup = ceed.QFunction(1, qfs.setup_mass, 275 os.path.join(file_dir, "test-qfunctions.h:setup_mass")) 276 qf_setup.add_input("weights", 1, libceed.EVAL_WEIGHT) 277 qf_setup.add_input("dx", 1, libceed.EVAL_GRAD) 278 qf_setup.add_output("rho", 1, libceed.EVAL_NONE) 279 280 qf_mass = ceed.QFunction(1, qfs.apply_mass_two, 281 os.path.join(file_dir, "test-qfunctions.h:apply_mass_two")) 282 qf_mass.add_input("rho", 1, libceed.EVAL_NONE) 283 qf_mass.add_input("u", 2, libceed.EVAL_INTERP) 284 qf_mass.add_output("v", 2, libceed.EVAL_INTERP) 285 286 # Operators 287 op_setup = ceed.Operator(qf_setup) 288 op_setup.set_field("weights", libceed.ELEMRESTRICTION_NONE, bx, 289 libceed.VECTOR_NONE) 290 op_setup.set_field("dx", rx, bx, libceed.VECTOR_ACTIVE) 291 op_setup.set_field("rho", rui, libceed.BASIS_COLLOCATED, 292 libceed.VECTOR_ACTIVE) 293 294 op_mass = ceed.Operator(qf_mass) 295 op_mass.set_field("rho", rui, libceed.BASIS_COLLOCATED, qdata) 296 op_mass.set_field("u", ru, bu, libceed.VECTOR_ACTIVE) 297 op_mass.set_field("v", ru, bu, libceed.VECTOR_ACTIVE) 298 299 # Setup 300 op_setup.apply(x, qdata) 301 302 # Apply mass matrix 303 with u.array() as u_array: 304 for i in range(nu): 305 u_array[2 * i] = 1. 306 u_array[2 * i + 1] = 2. 307 op_mass.apply(u, v) 308 309 # Check 310 with v.array_read() as v_array: 311 total_1 = 0.0 312 total_2 = 0.0 313 for i in range(nu): 314 total_1 = total_1 + v_array[2 * i] 315 total_2 = total_2 + v_array[2 * i + 1] 316 assert abs(total_1 - 1.0) < TOL 317 assert abs(total_2 - 2.0) < TOL 318 319# ------------------------------------------------------------------------------- 320# Test creation, action, and destruction for mass matrix operator with passive 321# inputs and outputs 322# ------------------------------------------------------------------------------- 323 324 325def test_503(ceed_resource): 326 ceed = libceed.Ceed(ceed_resource) 327 328 nelem = 15 329 p = 5 330 q = 8 331 nx = nelem + 1 332 nu = nelem * (p - 1) + 1 333 334 # Vectors 335 x = ceed.Vector(nx) 336 x_array = np.zeros(nx, dtype=ceed.scalar_type()) 337 for i in range(nx): 338 x_array[i] = i / (nx - 1.0) 339 x.set_array(x_array, cmode=libceed.USE_POINTER) 340 341 qdata = ceed.Vector(nelem * q) 342 u = ceed.Vector(nu) 343 v = ceed.Vector(nu) 344 345 # Restrictions 346 indx = np.zeros(nx * 2, dtype="int32") 347 for i in range(nx): 348 indx[2 * i + 0] = i 349 indx[2 * i + 1] = i + 1 350 rx = ceed.ElemRestriction(nelem, 2, 1, 1, nx, indx, 351 cmode=libceed.USE_POINTER) 352 353 indu = np.zeros(nelem * p, dtype="int32") 354 for i in range(nelem): 355 for j in range(p): 356 indu[p * i + j] = i * (p - 1) + j 357 ru = ceed.ElemRestriction(nelem, p, 1, 1, nu, indu, 358 cmode=libceed.USE_POINTER) 359 strides = np.array([1, q, q], dtype="int32") 360 rui = ceed.StridedElemRestriction(nelem, q, 1, q * nelem, strides) 361 362 # Bases 363 bx = ceed.BasisTensorH1Lagrange(1, 1, 2, q, libceed.GAUSS) 364 bu = ceed.BasisTensorH1Lagrange(1, 1, p, q, libceed.GAUSS) 365 366 # QFunctions 367 file_dir = os.path.dirname(os.path.abspath(__file__)) 368 qfs = load_qfs_so() 369 370 qf_setup = ceed.QFunction(1, qfs.setup_mass, 371 os.path.join(file_dir, "test-qfunctions.h:setup_mass")) 372 qf_setup.add_input("weights", 1, libceed.EVAL_WEIGHT) 373 qf_setup.add_input("dx", 1, libceed.EVAL_GRAD) 374 qf_setup.add_output("rho", 1, libceed.EVAL_NONE) 375 376 qf_mass = ceed.QFunction(1, qfs.apply_mass, 377 os.path.join(file_dir, "test-qfunctions.h:apply_mass")) 378 qf_mass.add_input("rho", 1, libceed.EVAL_NONE) 379 qf_mass.add_input("u", 1, libceed.EVAL_INTERP) 380 qf_mass.add_output("v", 1, libceed.EVAL_INTERP) 381 382 # Operators 383 op_setup = ceed.Operator(qf_setup) 384 op_setup.set_field("weights", libceed.ELEMRESTRICTION_NONE, bx, 385 libceed.VECTOR_NONE) 386 op_setup.set_field("dx", rx, bx, libceed.VECTOR_ACTIVE) 387 op_setup.set_field("rho", rui, libceed.BASIS_COLLOCATED, 388 libceed.VECTOR_ACTIVE) 389 390 op_mass = ceed.Operator(qf_mass) 391 op_mass.set_field("rho", rui, libceed.BASIS_COLLOCATED, qdata) 392 op_mass.set_field("u", ru, bu, u) 393 op_mass.set_field("v", ru, bu, v) 394 395 # Setup 396 op_setup.apply(x, qdata) 397 398 # Apply mass matrix 399 u.set_value(1) 400 op_mass.apply(libceed.VECTOR_NONE, libceed.VECTOR_NONE) 401 402 # Check 403 with v.array_read() as v_array: 404 total = 0.0 405 for i in range(nu): 406 total = total + v_array[i] 407 assert abs(total - 1.0) < TOL 408 409# ------------------------------------------------------------------------------- 410# Test viewing of mass matrix operator 411# ------------------------------------------------------------------------------- 412 413 414def test_504(ceed_resource, capsys): 415 ceed = libceed.Ceed(ceed_resource) 416 417 nelem = 15 418 p = 5 419 q = 8 420 nx = nelem + 1 421 nu = nelem * (p - 1) + 1 422 423 # Vectors 424 qdata = ceed.Vector(nelem * q) 425 426 # Restrictions 427 indx = np.zeros(nx * 2, dtype="int32") 428 for i in range(nx): 429 indx[2 * i + 0] = i 430 indx[2 * i + 1] = i + 1 431 rx = ceed.ElemRestriction(nelem, 2, 1, 1, nx, indx, 432 cmode=libceed.USE_POINTER) 433 434 indu = np.zeros(nelem * p, dtype="int32") 435 for i in range(nelem): 436 for j in range(p): 437 indu[p * i + j] = i * (p - 1) + j 438 ru = ceed.ElemRestriction(nelem, p, 1, 1, nu, indu, 439 cmode=libceed.USE_POINTER) 440 strides = np.array([1, q, q], dtype="int32") 441 rui = ceed.StridedElemRestriction(nelem, q, 1, q * nelem, strides) 442 443 # Bases 444 bx = ceed.BasisTensorH1Lagrange(1, 1, 2, q, libceed.GAUSS) 445 bu = ceed.BasisTensorH1Lagrange(1, 1, p, q, libceed.GAUSS) 446 447 # QFunctions 448 file_dir = os.path.dirname(os.path.abspath(__file__)) 449 qfs = load_qfs_so() 450 451 qf_setup = ceed.QFunction(1, qfs.setup_mass, 452 os.path.join(file_dir, "test-qfunctions.h:setup_mass")) 453 qf_setup.add_input("weights", 1, libceed.EVAL_WEIGHT) 454 qf_setup.add_input("dx", 1, libceed.EVAL_GRAD) 455 qf_setup.add_output("rho", 1, libceed.EVAL_NONE) 456 457 qf_mass = ceed.QFunction(1, qfs.apply_mass, 458 os.path.join(file_dir, "test-qfunctions.h:apply_mass")) 459 qf_mass.add_input("rho", 1, libceed.EVAL_NONE) 460 qf_mass.add_input("u", 1, libceed.EVAL_INTERP) 461 qf_mass.add_output("v", 1, libceed.EVAL_INTERP) 462 463 # Operators 464 op_setup = ceed.Operator(qf_setup) 465 op_setup.set_field("weights", libceed.ELEMRESTRICTION_NONE, bx, 466 libceed.VECTOR_NONE) 467 op_setup.set_field("dx", rx, bx, libceed.VECTOR_ACTIVE) 468 op_setup.set_field("rho", rui, libceed.BASIS_COLLOCATED, 469 libceed.VECTOR_ACTIVE) 470 471 op_mass = ceed.Operator(qf_mass) 472 op_mass.set_field("rho", rui, libceed.BASIS_COLLOCATED, qdata) 473 op_mass.set_field("u", ru, bu, libceed.VECTOR_ACTIVE) 474 op_mass.set_field("v", ru, bu, libceed.VECTOR_ACTIVE) 475 476 # View 477 print(op_setup) 478 print(op_mass) 479 480 stdout, stderr, ref_stdout = check.output(capsys) 481 assert not stderr 482 assert stdout == ref_stdout 483 484# ------------------------------------------------------------------------------- 485# Test CeedOperatorApplyAdd 486# ------------------------------------------------------------------------------- 487 488 489def test_505(ceed_resource): 490 ceed = libceed.Ceed(ceed_resource) 491 492 nelem = 15 493 p = 5 494 q = 8 495 nx = nelem + 1 496 nu = nelem * (p - 1) + 1 497 498 # Vectors 499 x = ceed.Vector(nx) 500 x_array = np.zeros(nx, dtype=ceed.scalar_type()) 501 for i in range(nx): 502 x_array[i] = i / (nx - 1.0) 503 x.set_array(x_array, cmode=libceed.USE_POINTER) 504 505 qdata = ceed.Vector(nelem * q) 506 u = ceed.Vector(nu) 507 v = ceed.Vector(nu) 508 509 # Restrictions 510 indx = np.zeros(nx * 2, dtype="int32") 511 for i in range(nx): 512 indx[2 * i + 0] = i 513 indx[2 * i + 1] = i + 1 514 rx = ceed.ElemRestriction(nelem, 2, 1, 1, nx, indx, 515 cmode=libceed.USE_POINTER) 516 517 indu = np.zeros(nelem * p, dtype="int32") 518 for i in range(nelem): 519 for j in range(p): 520 indu[p * i + j] = i * (p - 1) + j 521 ru = ceed.ElemRestriction(nelem, p, 1, 1, nu, indu, 522 cmode=libceed.USE_POINTER) 523 strides = np.array([1, q, q], dtype="int32") 524 rui = ceed.StridedElemRestriction(nelem, q, 1, q * nelem, strides) 525 526 # Bases 527 bx = ceed.BasisTensorH1Lagrange(1, 1, 2, q, libceed.GAUSS) 528 bu = ceed.BasisTensorH1Lagrange(1, 1, p, q, libceed.GAUSS) 529 530 # QFunctions 531 file_dir = os.path.dirname(os.path.abspath(__file__)) 532 qfs = load_qfs_so() 533 534 qf_setup = ceed.QFunction(1, qfs.setup_mass, 535 os.path.join(file_dir, "test-qfunctions.h:setup_mass")) 536 qf_setup.add_input("weights", 1, libceed.EVAL_WEIGHT) 537 qf_setup.add_input("dx", 1, libceed.EVAL_GRAD) 538 qf_setup.add_output("rho", 1, libceed.EVAL_NONE) 539 540 qf_mass = ceed.QFunction(1, qfs.apply_mass, 541 os.path.join(file_dir, "test-qfunctions.h:apply_mass")) 542 qf_mass.add_input("rho", 1, libceed.EVAL_NONE) 543 qf_mass.add_input("u", 1, libceed.EVAL_INTERP) 544 qf_mass.add_output("v", 1, libceed.EVAL_INTERP) 545 546 # Operators 547 op_setup = ceed.Operator(qf_setup) 548 op_setup.set_field("weights", libceed.ELEMRESTRICTION_NONE, bx, 549 libceed.VECTOR_NONE) 550 op_setup.set_field("dx", rx, bx, libceed.VECTOR_ACTIVE) 551 op_setup.set_field("rho", rui, libceed.BASIS_COLLOCATED, 552 libceed.VECTOR_ACTIVE) 553 554 op_mass = ceed.Operator(qf_mass) 555 op_mass.set_field("rho", rui, libceed.BASIS_COLLOCATED, qdata) 556 op_mass.set_field("u", ru, bu, libceed.VECTOR_ACTIVE) 557 op_mass.set_field("v", ru, bu, libceed.VECTOR_ACTIVE) 558 559 # Setup 560 op_setup.apply(x, qdata) 561 562 # Apply mass matrix with v = 0 563 u.set_value(1.) 564 v.set_value(0.) 565 op_mass.apply_add(u, v) 566 567 # Check 568 with v.array_read() as v_array: 569 total = 0.0 570 for i in range(nu): 571 total = total + v_array[i] 572 assert abs(total - 1.0) < TOL 573 574 # Apply mass matrix with v = 0 575 v.set_value(1.) 576 op_mass.apply_add(u, v) 577 578 # Check 579 with v.array_read() as v_array: 580 total = -nu 581 for i in range(nu): 582 total = total + v_array[i] 583 assert abs(total - 1.0) < 10. * TOL 584 585# ------------------------------------------------------------------------------- 586# Test creation, action, and destruction for mass matrix operator 587# ------------------------------------------------------------------------------- 588 589 590def test_510(ceed_resource): 591 ceed = libceed.Ceed(ceed_resource) 592 593 nelem = 12 594 dim = 2 595 p = 6 596 q = 4 597 nx, ny = 3, 2 598 ndofs = (nx * 2 + 1) * (ny * 2 + 1) 599 nqpts = nelem * q 600 601 # Vectors 602 x = ceed.Vector(dim * ndofs) 603 x_array = np.zeros(dim * ndofs) 604 for i in range(ndofs): 605 x_array[i] = (1. / (nx * 2)) * (i % (nx * 2 + 1)) 606 x_array[i + ndofs] = (1. / (ny * 2)) * (i / (nx * 2 + 1)) 607 x.set_array(x_array, cmode=libceed.USE_POINTER) 608 609 qdata = ceed.Vector(nqpts) 610 u = ceed.Vector(ndofs) 611 v = ceed.Vector(ndofs) 612 613 # Restrictions 614 indx = np.zeros(nelem * p, dtype="int32") 615 for i in range(nelem // 2): 616 col = i % nx 617 row = i // nx 618 offset = col * 2 + row * (nx * 2 + 1) * 2 619 620 indx[i * 2 * p + 0] = 2 + offset 621 indx[i * 2 * p + 1] = 9 + offset 622 indx[i * 2 * p + 2] = 16 + offset 623 indx[i * 2 * p + 3] = 1 + offset 624 indx[i * 2 * p + 4] = 8 + offset 625 indx[i * 2 * p + 5] = 0 + offset 626 627 indx[i * 2 * p + 6] = 14 + offset 628 indx[i * 2 * p + 7] = 7 + offset 629 indx[i * 2 * p + 8] = 0 + offset 630 indx[i * 2 * p + 9] = 15 + offset 631 indx[i * 2 * p + 10] = 8 + offset 632 indx[i * 2 * p + 11] = 16 + offset 633 634 rx = ceed.ElemRestriction(nelem, p, dim, ndofs, dim * ndofs, indx, 635 cmode=libceed.USE_POINTER) 636 637 ru = ceed.ElemRestriction(nelem, p, 1, 1, ndofs, indx, 638 cmode=libceed.USE_POINTER) 639 strides = np.array([1, q, q], dtype="int32") 640 rui = ceed.StridedElemRestriction(nelem, q, 1, nqpts, strides) 641 642 # Bases 643 qref = np.empty(dim * q, dtype=ceed.scalar_type()) 644 qweight = np.empty(q, dtype=ceed.scalar_type()) 645 interp, grad = bm.buildmats(qref, qweight, libceed.scalar_types[ 646 libceed.lib.CEED_SCALAR_TYPE]) 647 648 bx = ceed.BasisH1(libceed.TRIANGLE, dim, p, q, interp, grad, qref, qweight) 649 bu = ceed.BasisH1(libceed.TRIANGLE, 1, p, q, interp, grad, qref, qweight) 650 651 # QFunctions 652 file_dir = os.path.dirname(os.path.abspath(__file__)) 653 qfs = load_qfs_so() 654 655 qf_setup = ceed.QFunction(1, qfs.setup_mass_2d, 656 os.path.join(file_dir, "test-qfunctions.h:setup_mass_2d")) 657 qf_setup.add_input("weights", 1, libceed.EVAL_WEIGHT) 658 qf_setup.add_input("dx", dim * dim, libceed.EVAL_GRAD) 659 qf_setup.add_output("rho", 1, libceed.EVAL_NONE) 660 661 qf_mass = ceed.QFunction(1, qfs.apply_mass, 662 os.path.join(file_dir, "test-qfunctions.h:apply_mass")) 663 qf_mass.add_input("rho", 1, libceed.EVAL_NONE) 664 qf_mass.add_input("u", 1, libceed.EVAL_INTERP) 665 qf_mass.add_output("v", 1, libceed.EVAL_INTERP) 666 667 # Operators 668 op_setup = ceed.Operator(qf_setup) 669 op_setup.set_field("weights", libceed.ELEMRESTRICTION_NONE, bx, 670 libceed.VECTOR_NONE) 671 op_setup.set_field("dx", rx, bx, libceed.VECTOR_ACTIVE) 672 op_setup.set_field("rho", rui, libceed.BASIS_COLLOCATED, 673 libceed.VECTOR_ACTIVE) 674 675 op_mass = ceed.Operator(qf_mass) 676 op_mass.set_field("rho", rui, libceed.BASIS_COLLOCATED, qdata) 677 op_mass.set_field("u", ru, bu, libceed.VECTOR_ACTIVE) 678 op_mass.set_field("v", ru, bu, libceed.VECTOR_ACTIVE) 679 680 # Setup 681 op_setup.apply(x, qdata) 682 683 # Apply mass matrix 684 u.set_value(0.) 685 op_mass.apply(u, v) 686 687 # Check 688 with v.array_read() as v_array: 689 for i in range(ndofs): 690 assert abs(v_array[i]) < TOL 691 692# ------------------------------------------------------------------------------- 693# Test creation, action, and destruction for mass matrix operator 694# ------------------------------------------------------------------------------- 695 696 697def test_511(ceed_resource): 698 ceed = libceed.Ceed(ceed_resource) 699 700 nelem = 12 701 dim = 2 702 p = 6 703 q = 4 704 nx, ny = 3, 2 705 ndofs = (nx * 2 + 1) * (ny * 2 + 1) 706 nqpts = nelem * q 707 708 # Vectors 709 x = ceed.Vector(dim * ndofs) 710 x_array = np.zeros(dim * ndofs, dtype=ceed.scalar_type()) 711 for i in range(ndofs): 712 x_array[i] = (1. / (nx * 2)) * (i % (nx * 2 + 1)) 713 x_array[i + ndofs] = (1. / (ny * 2)) * (i / (nx * 2 + 1)) 714 x.set_array(x_array, cmode=libceed.USE_POINTER) 715 716 qdata = ceed.Vector(nqpts) 717 u = ceed.Vector(ndofs) 718 v = ceed.Vector(ndofs) 719 720 # Restrictions 721 indx = np.zeros(nelem * p, dtype="int32") 722 for i in range(nelem // 2): 723 col = i % nx 724 row = i // nx 725 offset = col * 2 + row * (nx * 2 + 1) * 2 726 727 indx[i * 2 * p + 0] = 2 + offset 728 indx[i * 2 * p + 1] = 9 + offset 729 indx[i * 2 * p + 2] = 16 + offset 730 indx[i * 2 * p + 3] = 1 + offset 731 indx[i * 2 * p + 4] = 8 + offset 732 indx[i * 2 * p + 5] = 0 + offset 733 734 indx[i * 2 * p + 6] = 14 + offset 735 indx[i * 2 * p + 7] = 7 + offset 736 indx[i * 2 * p + 8] = 0 + offset 737 indx[i * 2 * p + 9] = 15 + offset 738 indx[i * 2 * p + 10] = 8 + offset 739 indx[i * 2 * p + 11] = 16 + offset 740 741 rx = ceed.ElemRestriction(nelem, p, dim, ndofs, dim * ndofs, indx, 742 cmode=libceed.USE_POINTER) 743 744 ru = ceed.ElemRestriction(nelem, p, 1, 1, ndofs, indx, 745 cmode=libceed.USE_POINTER) 746 strides = np.array([1, q, q], dtype="int32") 747 rui = ceed.StridedElemRestriction(nelem, q, 1, nqpts, strides) 748 749 # Bases 750 qref = np.empty(dim * q, dtype=ceed.scalar_type()) 751 qweight = np.empty(q, dtype=ceed.scalar_type()) 752 interp, grad = bm.buildmats(qref, qweight, libceed.scalar_types[ 753 libceed.lib.CEED_SCALAR_TYPE]) 754 755 bx = ceed.BasisH1(libceed.TRIANGLE, dim, p, q, interp, grad, qref, qweight) 756 bu = ceed.BasisH1(libceed.TRIANGLE, 1, p, q, interp, grad, qref, qweight) 757 758 # QFunctions 759 file_dir = os.path.dirname(os.path.abspath(__file__)) 760 qfs = load_qfs_so() 761 762 qf_setup = ceed.QFunction(1, qfs.setup_mass_2d, 763 os.path.join(file_dir, "test-qfunctions.h:setup_mass_2d")) 764 qf_setup.add_input("weights", 1, libceed.EVAL_WEIGHT) 765 qf_setup.add_input("dx", dim * dim, libceed.EVAL_GRAD) 766 qf_setup.add_output("rho", 1, libceed.EVAL_NONE) 767 768 qf_mass = ceed.QFunction(1, qfs.apply_mass, 769 os.path.join(file_dir, "test-qfunctions.h:apply_mass")) 770 qf_mass.add_input("rho", 1, libceed.EVAL_NONE) 771 qf_mass.add_input("u", 1, libceed.EVAL_INTERP) 772 qf_mass.add_output("v", 1, libceed.EVAL_INTERP) 773 774 # Operators 775 op_setup = ceed.Operator(qf_setup) 776 op_setup.set_field("weights", libceed.ELEMRESTRICTION_NONE, bx, 777 libceed.VECTOR_NONE) 778 op_setup.set_field("dx", rx, bx, libceed.VECTOR_ACTIVE) 779 op_setup.set_field("rho", rui, libceed.BASIS_COLLOCATED, 780 libceed.VECTOR_ACTIVE) 781 782 op_mass = ceed.Operator(qf_mass) 783 op_mass.set_field("rho", rui, libceed.BASIS_COLLOCATED, qdata) 784 op_mass.set_field("u", ru, bu, libceed.VECTOR_ACTIVE) 785 op_mass.set_field("v", ru, bu, libceed.VECTOR_ACTIVE) 786 787 # Setup 788 op_setup.apply(x, qdata) 789 790 # Apply mass matrix 791 u.set_value(1.) 792 op_mass.apply(u, v) 793 794 # Check 795 with v.array_read() as v_array: 796 total = 0.0 797 for i in range(ndofs): 798 total = total + v_array[i] 799 assert abs(total - 1.0) < 10. * TOL 800 801# ------------------------------------------------------------------------------- 802# Test creation, action, and destruction for composite mass matrix operator 803# ------------------------------------------------------------------------------- 804 805 806def test_520(ceed_resource): 807 ceed = libceed.Ceed(ceed_resource) 808 809 nelem_tet, p_tet, q_tet = 6, 6, 4 810 nelem_hex, p_hex, q_hex = 6, 3, 4 811 nx, ny = 3, 3 812 dim = 2 813 nx_tet, ny_tet, nx_hex = 3, 1, 3 814 ndofs = (nx * 2 + 1) * (ny * 2 + 1) 815 nqpts_tet, nqpts_hex = nelem_tet * q_tet, nelem_hex * q_hex * q_hex 816 817 # Vectors 818 x = ceed.Vector(dim * ndofs) 819 x_array = np.zeros(dim * ndofs) 820 for i in range(ny * 2 + 1): 821 for j in range(nx * 2 + 1): 822 x_array[i + j * (ny * 2 + 1)] = i / (2 * ny) 823 x_array[i + j * (ny * 2 + 1) + ndofs] = j / (2 * nx) 824 x.set_array(x_array, cmode=libceed.USE_POINTER) 825 826 qdata_hex = ceed.Vector(nqpts_hex) 827 qdata_tet = ceed.Vector(nqpts_tet) 828 u = ceed.Vector(ndofs) 829 v = ceed.Vector(ndofs) 830 831 # ------------------------- Tet Elements ------------------------- 832 833 # Restrictions 834 indx_tet = np.zeros(nelem_tet * p_tet, dtype="int32") 835 for i in range(nelem_tet // 2): 836 col = i % nx 837 row = i // nx 838 offset = col * 2 + row * (nx * 2 + 1) * 2 839 840 indx_tet[i * 2 * p_tet + 0] = 2 + offset 841 indx_tet[i * 2 * p_tet + 1] = 9 + offset 842 indx_tet[i * 2 * p_tet + 2] = 16 + offset 843 indx_tet[i * 2 * p_tet + 3] = 1 + offset 844 indx_tet[i * 2 * p_tet + 4] = 8 + offset 845 indx_tet[i * 2 * p_tet + 5] = 0 + offset 846 847 indx_tet[i * 2 * p_tet + 6] = 14 + offset 848 indx_tet[i * 2 * p_tet + 7] = 7 + offset 849 indx_tet[i * 2 * p_tet + 8] = 0 + offset 850 indx_tet[i * 2 * p_tet + 9] = 15 + offset 851 indx_tet[i * 2 * p_tet + 10] = 8 + offset 852 indx_tet[i * 2 * p_tet + 11] = 16 + offset 853 854 rx_tet = ceed.ElemRestriction(nelem_tet, p_tet, dim, ndofs, dim * ndofs, 855 indx_tet, cmode=libceed.USE_POINTER) 856 857 ru_tet = ceed.ElemRestriction(nelem_tet, p_tet, 1, 1, ndofs, indx_tet, 858 cmode=libceed.USE_POINTER) 859 strides = np.array([1, q_tet, q_tet], dtype="int32") 860 rui_tet = ceed.StridedElemRestriction(nelem_tet, q_tet, 1, nqpts_tet, 861 strides) 862 863 # Bases 864 qref = np.empty(dim * q_tet, dtype=ceed.scalar_type()) 865 qweight = np.empty(q_tet, dtype=ceed.scalar_type()) 866 interp, grad = bm.buildmats(qref, qweight, libceed.scalar_types[ 867 libceed.lib.CEED_SCALAR_TYPE]) 868 869 bx_tet = ceed.BasisH1(libceed.TRIANGLE, dim, p_tet, q_hex, interp, grad, qref, 870 qweight) 871 bu_tet = ceed.BasisH1(libceed.TRIANGLE, 1, p_tet, q_hex, interp, grad, qref, 872 qweight) 873 874 # QFunctions 875 file_dir = os.path.dirname(os.path.abspath(__file__)) 876 qfs = load_qfs_so() 877 878 qf_setup_tet = ceed.QFunction(1, qfs.setup_mass_2d, 879 os.path.join(file_dir, "test-qfunctions.h:setup_mass_2d")) 880 qf_setup_tet.add_input("weights", 1, libceed.EVAL_WEIGHT) 881 qf_setup_tet.add_input("dx", dim * dim, libceed.EVAL_GRAD) 882 qf_setup_tet.add_output("rho", 1, libceed.EVAL_NONE) 883 884 qf_mass_tet = ceed.QFunction(1, qfs.apply_mass, 885 os.path.join(file_dir, "test-qfunctions.h:apply_mass")) 886 qf_mass_tet.add_input("rho", 1, libceed.EVAL_NONE) 887 qf_mass_tet.add_input("u", 1, libceed.EVAL_INTERP) 888 qf_mass_tet.add_output("v", 1, libceed.EVAL_INTERP) 889 890 # Operators 891 op_setup_tet = ceed.Operator(qf_setup_tet) 892 op_setup_tet.set_field("weights", libceed.ELEMRESTRICTION_NONE, bx_tet, 893 libceed.VECTOR_NONE) 894 op_setup_tet.set_field("dx", rx_tet, bx_tet, libceed.VECTOR_ACTIVE) 895 op_setup_tet.set_field("rho", rui_tet, libceed.BASIS_COLLOCATED, 896 qdata_tet) 897 898 op_mass_tet = ceed.Operator(qf_mass_tet) 899 op_mass_tet.set_field("rho", rui_tet, libceed.BASIS_COLLOCATED, qdata_tet) 900 op_mass_tet.set_field("u", ru_tet, bu_tet, libceed.VECTOR_ACTIVE) 901 op_mass_tet.set_field("v", ru_tet, bu_tet, libceed.VECTOR_ACTIVE) 902 903 # ------------------------- Hex Elements ------------------------- 904 905 # Restrictions 906 indx_hex = np.zeros(nelem_hex * p_hex * p_hex, dtype="int32") 907 for i in range(nelem_hex): 908 col = i % nx_hex 909 row = i // nx_hex 910 offset = (nx_tet * 2 + 1) * (ny_tet * 2) * (1 + row) + col * 2 911 912 for j in range(p_hex): 913 for k in range(p_hex): 914 indx_hex[p_hex * (p_hex * i + k) + j] = offset + \ 915 k * (nx_hex * 2 + 1) + j 916 917 rx_hex = ceed.ElemRestriction(nelem_hex, p_hex * p_hex, dim, ndofs, 918 dim * ndofs, indx_hex, cmode=libceed.USE_POINTER) 919 920 ru_hex = ceed.ElemRestriction(nelem_hex, p_hex * p_hex, 1, 1, ndofs, 921 indx_hex, cmode=libceed.USE_POINTER) 922 strides = np.array([1, q_hex * q_hex, q_hex * q_hex], dtype="int32") 923 rui_hex = ceed.StridedElemRestriction(nelem_hex, q_hex * q_hex, 1, 924 nqpts_hex, strides) 925 926 # Bases 927 bx_hex = ceed.BasisTensorH1Lagrange(dim, dim, p_hex, q_hex, libceed.GAUSS) 928 bu_hex = ceed.BasisTensorH1Lagrange(dim, 1, p_hex, q_hex, libceed.GAUSS) 929 930 # QFunctions 931 qf_setup_hex = ceed.QFunction(1, qfs.setup_mass_2d, 932 os.path.join(file_dir, "test-qfunctions.h:setup_mass_2d")) 933 qf_setup_hex.add_input("weights", 1, libceed.EVAL_WEIGHT) 934 qf_setup_hex.add_input("dx", dim * dim, libceed.EVAL_GRAD) 935 qf_setup_hex.add_output("rho", 1, libceed.EVAL_NONE) 936 937 qf_mass_hex = ceed.QFunction(1, qfs.apply_mass, 938 os.path.join(file_dir, "test-qfunctions.h:apply_mass")) 939 qf_mass_hex.add_input("rho", 1, libceed.EVAL_NONE) 940 qf_mass_hex.add_input("u", 1, libceed.EVAL_INTERP) 941 qf_mass_hex.add_output("v", 1, libceed.EVAL_INTERP) 942 943 # Operators 944 op_setup_hex = ceed.Operator(qf_setup_tet) 945 op_setup_hex.set_field("weights", libceed.ELEMRESTRICTION_NONE, bx_hex, 946 libceed.VECTOR_NONE) 947 op_setup_hex.set_field("dx", rx_hex, bx_hex, libceed.VECTOR_ACTIVE) 948 op_setup_hex.set_field("rho", rui_hex, libceed.BASIS_COLLOCATED, 949 qdata_hex) 950 951 op_mass_hex = ceed.Operator(qf_mass_hex) 952 op_mass_hex.set_field("rho", rui_hex, libceed.BASIS_COLLOCATED, qdata_hex) 953 op_mass_hex.set_field("u", ru_hex, bu_hex, libceed.VECTOR_ACTIVE) 954 op_mass_hex.set_field("v", ru_hex, bu_hex, libceed.VECTOR_ACTIVE) 955 956 # ------------------------- Composite Operators ------------------------- 957 958 # Setup 959 op_setup = ceed.CompositeOperator() 960 op_setup.add_sub(op_setup_tet) 961 op_setup.add_sub(op_setup_hex) 962 op_setup.apply(x, libceed.VECTOR_NONE) 963 964 # Apply mass matrix 965 op_mass = ceed.CompositeOperator() 966 op_mass.add_sub(op_mass_tet) 967 op_mass.add_sub(op_mass_hex) 968 969 u.set_value(0.) 970 op_mass.apply(u, v) 971 972 # Check 973 with v.array_read() as v_array: 974 for i in range(ndofs): 975 assert abs(v_array[i]) < TOL 976 977# ------------------------------------------------------------------------------- 978# Test creation, action, and destruction for composite mass matrix operator 979# ------------------------------------------------------------------------------- 980 981 982def test_521(ceed_resource): 983 ceed = libceed.Ceed(ceed_resource) 984 985 nelem_tet, p_tet, q_tet = 6, 6, 4 986 nelem_hex, p_hex, q_hex = 6, 3, 4 987 nx, ny = 3, 3 988 dim = 2 989 nx_tet, ny_tet, nx_hex = 3, 1, 3 990 ndofs = (nx * 2 + 1) * (ny * 2 + 1) 991 nqpts_tet, nqpts_hex = nelem_tet * q_tet, nelem_hex * q_hex * q_hex 992 993 # Vectors 994 x = ceed.Vector(dim * ndofs) 995 x_array = np.zeros(dim * ndofs, dtype=ceed.scalar_type()) 996 for i in range(ny * 2 + 1): 997 for j in range(nx * 2 + 1): 998 x_array[i + j * (ny * 2 + 1)] = i / (2 * ny) 999 x_array[i + j * (ny * 2 + 1) + ndofs] = j / (2 * nx) 1000 x.set_array(x_array, cmode=libceed.USE_POINTER) 1001 1002 qdata_hex = ceed.Vector(nqpts_hex) 1003 qdata_tet = ceed.Vector(nqpts_tet) 1004 u = ceed.Vector(ndofs) 1005 v = ceed.Vector(ndofs) 1006 1007 # ------------------------- Tet Elements ------------------------- 1008 1009 # Restrictions 1010 indx_tet = np.zeros(nelem_tet * p_tet, dtype="int32") 1011 for i in range(nelem_tet // 2): 1012 col = i % nx 1013 row = i // nx 1014 offset = col * 2 + row * (nx * 2 + 1) * 2 1015 1016 indx_tet[i * 2 * p_tet + 0] = 2 + offset 1017 indx_tet[i * 2 * p_tet + 1] = 9 + offset 1018 indx_tet[i * 2 * p_tet + 2] = 16 + offset 1019 indx_tet[i * 2 * p_tet + 3] = 1 + offset 1020 indx_tet[i * 2 * p_tet + 4] = 8 + offset 1021 indx_tet[i * 2 * p_tet + 5] = 0 + offset 1022 1023 indx_tet[i * 2 * p_tet + 6] = 14 + offset 1024 indx_tet[i * 2 * p_tet + 7] = 7 + offset 1025 indx_tet[i * 2 * p_tet + 8] = 0 + offset 1026 indx_tet[i * 2 * p_tet + 9] = 15 + offset 1027 indx_tet[i * 2 * p_tet + 10] = 8 + offset 1028 indx_tet[i * 2 * p_tet + 11] = 16 + offset 1029 1030 rx_tet = ceed.ElemRestriction(nelem_tet, p_tet, dim, ndofs, dim * ndofs, 1031 indx_tet, cmode=libceed.USE_POINTER) 1032 1033 ru_tet = ceed.ElemRestriction(nelem_tet, p_tet, 1, 1, ndofs, indx_tet, 1034 cmode=libceed.USE_POINTER) 1035 strides = np.array([1, q_tet, q_tet], dtype="int32") 1036 rui_tet = ceed.StridedElemRestriction(nelem_tet, q_tet, 1, nqpts_tet, 1037 strides) 1038 1039 # Bases 1040 qref = np.empty(dim * q_tet, dtype=ceed.scalar_type()) 1041 qweight = np.empty(q_tet, dtype=ceed.scalar_type()) 1042 interp, grad = bm.buildmats(qref, qweight, libceed.scalar_types[ 1043 libceed.lib.CEED_SCALAR_TYPE]) 1044 1045 bx_tet = ceed.BasisH1(libceed.TRIANGLE, dim, p_tet, q_hex, interp, grad, qref, 1046 qweight) 1047 bu_tet = ceed.BasisH1(libceed.TRIANGLE, 1, p_tet, q_hex, interp, grad, qref, 1048 qweight) 1049 1050 # QFunctions 1051 file_dir = os.path.dirname(os.path.abspath(__file__)) 1052 qfs = load_qfs_so() 1053 1054 qf_setup_tet = ceed.QFunction(1, qfs.setup_mass_2d, 1055 os.path.join(file_dir, "test-qfunctions.h:setup_mass_2d")) 1056 qf_setup_tet.add_input("weights", 1, libceed.EVAL_WEIGHT) 1057 qf_setup_tet.add_input("dx", dim * dim, libceed.EVAL_GRAD) 1058 qf_setup_tet.add_output("rho", 1, libceed.EVAL_NONE) 1059 1060 qf_mass_tet = ceed.QFunction(1, qfs.apply_mass, 1061 os.path.join(file_dir, "test-qfunctions.h:apply_mass")) 1062 qf_mass_tet.add_input("rho", 1, libceed.EVAL_NONE) 1063 qf_mass_tet.add_input("u", 1, libceed.EVAL_INTERP) 1064 qf_mass_tet.add_output("v", 1, libceed.EVAL_INTERP) 1065 1066 # Operators 1067 op_setup_tet = ceed.Operator(qf_setup_tet) 1068 op_setup_tet.set_field("weights", libceed.ELEMRESTRICTION_NONE, bx_tet, 1069 libceed.VECTOR_NONE) 1070 op_setup_tet.set_field("dx", rx_tet, bx_tet, libceed.VECTOR_ACTIVE) 1071 op_setup_tet.set_field("rho", rui_tet, libceed.BASIS_COLLOCATED, 1072 qdata_tet) 1073 1074 op_mass_tet = ceed.Operator(qf_mass_tet) 1075 op_mass_tet.set_field("rho", rui_tet, libceed.BASIS_COLLOCATED, qdata_tet) 1076 op_mass_tet.set_field("u", ru_tet, bu_tet, libceed.VECTOR_ACTIVE) 1077 op_mass_tet.set_field("v", ru_tet, bu_tet, libceed.VECTOR_ACTIVE) 1078 1079 # ------------------------- Hex Elements ------------------------- 1080 1081 # Restrictions 1082 indx_hex = np.zeros(nelem_hex * p_hex * p_hex, dtype="int32") 1083 for i in range(nelem_hex): 1084 col = i % nx_hex 1085 row = i // nx_hex 1086 offset = (nx_tet * 2 + 1) * (ny_tet * 2) * (1 + row) + col * 2 1087 1088 for j in range(p_hex): 1089 for k in range(p_hex): 1090 indx_hex[p_hex * (p_hex * i + k) + j] = offset + \ 1091 k * (nx_hex * 2 + 1) + j 1092 1093 rx_hex = ceed.ElemRestriction(nelem_hex, p_hex * p_hex, dim, ndofs, 1094 dim * ndofs, indx_hex, cmode=libceed.USE_POINTER) 1095 1096 ru_hex = ceed.ElemRestriction(nelem_hex, p_hex * p_hex, 1, 1, ndofs, 1097 indx_hex, cmode=libceed.USE_POINTER) 1098 strides = np.array([1, q_hex * q_hex, q_hex * q_hex], dtype="int32") 1099 rui_hex = ceed.StridedElemRestriction(nelem_hex, q_hex * q_hex, 1, 1100 nqpts_hex, strides) 1101 1102 # Bases 1103 bx_hex = ceed.BasisTensorH1Lagrange(dim, dim, p_hex, q_hex, libceed.GAUSS) 1104 bu_hex = ceed.BasisTensorH1Lagrange(dim, 1, p_hex, q_hex, libceed.GAUSS) 1105 1106 # QFunctions 1107 qf_setup_hex = ceed.QFunction(1, qfs.setup_mass_2d, 1108 os.path.join(file_dir, "test-qfunctions.h:setup_mass_2d")) 1109 qf_setup_hex.add_input("weights", 1, libceed.EVAL_WEIGHT) 1110 qf_setup_hex.add_input("dx", dim * dim, libceed.EVAL_GRAD) 1111 qf_setup_hex.add_output("rho", 1, libceed.EVAL_NONE) 1112 1113 qf_mass_hex = ceed.QFunction(1, qfs.apply_mass, 1114 os.path.join(file_dir, "test-qfunctions.h:apply_mass")) 1115 qf_mass_hex.add_input("rho", 1, libceed.EVAL_NONE) 1116 qf_mass_hex.add_input("u", 1, libceed.EVAL_INTERP) 1117 qf_mass_hex.add_output("v", 1, libceed.EVAL_INTERP) 1118 1119 # Operators 1120 op_setup_hex = ceed.Operator(qf_setup_tet) 1121 op_setup_hex.set_field("weights", libceed.ELEMRESTRICTION_NONE, bx_hex, 1122 libceed.VECTOR_NONE) 1123 op_setup_hex.set_field("dx", rx_hex, bx_hex, libceed.VECTOR_ACTIVE) 1124 op_setup_hex.set_field("rho", rui_hex, libceed.BASIS_COLLOCATED, 1125 qdata_hex) 1126 1127 op_mass_hex = ceed.Operator(qf_mass_hex) 1128 op_mass_hex.set_field("rho", rui_hex, libceed.BASIS_COLLOCATED, qdata_hex) 1129 op_mass_hex.set_field("u", ru_hex, bu_hex, libceed.VECTOR_ACTIVE) 1130 op_mass_hex.set_field("v", ru_hex, bu_hex, libceed.VECTOR_ACTIVE) 1131 1132 # ------------------------- Composite Operators ------------------------- 1133 1134 # Setup 1135 op_setup = ceed.CompositeOperator() 1136 op_setup.add_sub(op_setup_tet) 1137 op_setup.add_sub(op_setup_hex) 1138 op_setup.apply(x, libceed.VECTOR_NONE) 1139 1140 # Apply mass matrix 1141 op_mass = ceed.CompositeOperator() 1142 op_mass.add_sub(op_mass_tet) 1143 op_mass.add_sub(op_mass_hex) 1144 u.set_value(1.) 1145 op_mass.apply(u, v) 1146 1147 # Check 1148 with v.array_read() as v_array: 1149 total = 0.0 1150 for i in range(ndofs): 1151 total = total + v_array[i] 1152 assert abs(total - 1.0) < 10. * TOL 1153 1154# ------------------------------------------------------------------------------- 1155# Test viewing of composite mass matrix operator 1156# ------------------------------------------------------------------------------- 1157 1158 1159def test_523(ceed_resource, capsys): 1160 ceed = libceed.Ceed(ceed_resource) 1161 1162 nelem_tet, p_tet, q_tet = 6, 6, 4 1163 nelem_hex, p_hex, q_hex = 6, 3, 4 1164 nx, ny = 3, 3 1165 dim = 2 1166 nx_tet, ny_tet, nx_hex = 3, 1, 3 1167 ndofs = (nx * 2 + 1) * (ny * 2 + 1) 1168 nqpts_tet, nqpts_hex = nelem_tet * q_tet, nelem_hex * q_hex * q_hex 1169 1170 # Vectors 1171 qdata_hex = ceed.Vector(nqpts_hex) 1172 qdata_tet = ceed.Vector(nqpts_tet) 1173 1174 # ------------------------- Tet Elements ------------------------- 1175 1176 # Restrictions 1177 indx_tet = np.zeros(nelem_tet * p_tet, dtype="int32") 1178 for i in range(nelem_tet // 2): 1179 col = i % nx 1180 row = i // nx 1181 offset = col * 2 + row * (nx * 2 + 1) * 2 1182 1183 indx_tet[i * 2 * p_tet + 0] = 2 + offset 1184 indx_tet[i * 2 * p_tet + 1] = 9 + offset 1185 indx_tet[i * 2 * p_tet + 2] = 16 + offset 1186 indx_tet[i * 2 * p_tet + 3] = 1 + offset 1187 indx_tet[i * 2 * p_tet + 4] = 8 + offset 1188 indx_tet[i * 2 * p_tet + 5] = 0 + offset 1189 1190 indx_tet[i * 2 * p_tet + 6] = 14 + offset 1191 indx_tet[i * 2 * p_tet + 7] = 7 + offset 1192 indx_tet[i * 2 * p_tet + 8] = 0 + offset 1193 indx_tet[i * 2 * p_tet + 9] = 15 + offset 1194 indx_tet[i * 2 * p_tet + 10] = 8 + offset 1195 indx_tet[i * 2 * p_tet + 11] = 16 + offset 1196 1197 rx_tet = ceed.ElemRestriction(nelem_tet, p_tet, dim, ndofs, dim * ndofs, 1198 indx_tet, cmode=libceed.USE_POINTER) 1199 1200 ru_tet = ceed.ElemRestriction(nelem_tet, p_tet, 1, 1, ndofs, indx_tet, 1201 cmode=libceed.USE_POINTER) 1202 strides = np.array([1, q_tet, q_tet], dtype="int32") 1203 rui_tet = ceed.StridedElemRestriction(nelem_tet, q_tet, 1, nqpts_tet, 1204 strides) 1205 1206 # Bases 1207 qref = np.empty(dim * q_tet, dtype=ceed.scalar_type()) 1208 qweight = np.empty(q_tet, dtype=ceed.scalar_type()) 1209 interp, grad = bm.buildmats(qref, qweight, libceed.scalar_types[ 1210 libceed.lib.CEED_SCALAR_TYPE]) 1211 1212 bx_tet = ceed.BasisH1(libceed.TRIANGLE, dim, p_tet, q_hex, interp, grad, qref, 1213 qweight) 1214 bu_tet = ceed.BasisH1(libceed.TRIANGLE, 1, p_tet, q_hex, interp, grad, qref, 1215 qweight) 1216 1217 # QFunctions 1218 file_dir = os.path.dirname(os.path.abspath(__file__)) 1219 qfs = load_qfs_so() 1220 1221 qf_setup_tet = ceed.QFunction(1, qfs.setup_mass_2d, 1222 os.path.join(file_dir, "test-qfunctions.h:setup_mass_2d")) 1223 qf_setup_tet.add_input("weights", 1, libceed.EVAL_WEIGHT) 1224 qf_setup_tet.add_input("dx", dim * dim, libceed.EVAL_GRAD) 1225 qf_setup_tet.add_output("rho", 1, libceed.EVAL_NONE) 1226 1227 qf_mass_tet = ceed.QFunction(1, qfs.apply_mass, 1228 os.path.join(file_dir, "test-qfunctions.h:apply_mass")) 1229 qf_mass_tet.add_input("rho", 1, libceed.EVAL_NONE) 1230 qf_mass_tet.add_input("u", 1, libceed.EVAL_INTERP) 1231 qf_mass_tet.add_output("v", 1, libceed.EVAL_INTERP) 1232 1233 # Operators 1234 op_setup_tet = ceed.Operator(qf_setup_tet) 1235 op_setup_tet.set_field("weights", libceed.ELEMRESTRICTION_NONE, bx_tet, 1236 libceed.VECTOR_NONE) 1237 op_setup_tet.set_field("dx", rx_tet, bx_tet, libceed.VECTOR_ACTIVE) 1238 op_setup_tet.set_field("rho", rui_tet, libceed.BASIS_COLLOCATED, 1239 qdata_tet) 1240 1241 op_mass_tet = ceed.Operator(qf_mass_tet) 1242 op_mass_tet.set_field("rho", rui_tet, libceed.BASIS_COLLOCATED, qdata_tet) 1243 op_mass_tet.set_field("u", ru_tet, bu_tet, libceed.VECTOR_ACTIVE) 1244 op_mass_tet.set_field("v", ru_tet, bu_tet, libceed.VECTOR_ACTIVE) 1245 1246 # ------------------------- Hex Elements ------------------------- 1247 1248 # Restrictions 1249 indx_hex = np.zeros(nelem_hex * p_hex * p_hex, dtype="int32") 1250 for i in range(nelem_hex): 1251 col = i % nx_hex 1252 row = i // nx_hex 1253 offset = (nx_tet * 2 + 1) * (ny_tet * 2) * (1 + row) + col * 2 1254 1255 for j in range(p_hex): 1256 for k in range(p_hex): 1257 indx_hex[p_hex * (p_hex * i + k) + j] = offset + \ 1258 k * (nx_hex * 2 + 1) + j 1259 1260 rx_hex = ceed.ElemRestriction(nelem_hex, p_hex * p_hex, dim, ndofs, 1261 dim * ndofs, indx_hex, 1262 cmode=libceed.USE_POINTER) 1263 1264 ru_hex = ceed.ElemRestriction(nelem_hex, p_hex * p_hex, 1, 1, ndofs, 1265 indx_hex, cmode=libceed.USE_POINTER) 1266 strides = np.array([1, q_hex * q_hex, q_hex * q_hex], dtype="int32") 1267 rui_hex = ceed.StridedElemRestriction(nelem_hex, q_hex * q_hex, 1, 1268 nqpts_hex, strides) 1269 1270 # Bases 1271 bx_hex = ceed.BasisTensorH1Lagrange(dim, dim, p_hex, q_hex, libceed.GAUSS) 1272 bu_hex = ceed.BasisTensorH1Lagrange(dim, 1, p_hex, q_hex, libceed.GAUSS) 1273 1274 # QFunctions 1275 qf_setup_hex = ceed.QFunction(1, qfs.setup_mass_2d, 1276 os.path.join(file_dir, "test-qfunctions.h:setup_mass_2d")) 1277 qf_setup_hex.add_input("weights", 1, libceed.EVAL_WEIGHT) 1278 qf_setup_hex.add_input("dx", dim * dim, libceed.EVAL_GRAD) 1279 qf_setup_hex.add_output("rho", 1, libceed.EVAL_NONE) 1280 1281 qf_mass_hex = ceed.QFunction(1, qfs.apply_mass, 1282 os.path.join(file_dir, "test-qfunctions.h:apply_mass")) 1283 qf_mass_hex.add_input("rho", 1, libceed.EVAL_NONE) 1284 qf_mass_hex.add_input("u", 1, libceed.EVAL_INTERP) 1285 qf_mass_hex.add_output("v", 1, libceed.EVAL_INTERP) 1286 1287 # Operators 1288 op_setup_hex = ceed.Operator(qf_setup_tet) 1289 op_setup_hex.set_field("weights", libceed.ELEMRESTRICTION_NONE, bx_hex, 1290 libceed.VECTOR_NONE) 1291 op_setup_hex.set_field("dx", rx_hex, bx_hex, libceed.VECTOR_ACTIVE) 1292 op_setup_hex.set_field("rho", rui_hex, libceed.BASIS_COLLOCATED, 1293 qdata_hex) 1294 1295 op_mass_hex = ceed.Operator(qf_mass_hex) 1296 op_mass_hex.set_field("rho", rui_hex, libceed.BASIS_COLLOCATED, qdata_hex) 1297 op_mass_hex.set_field("u", ru_hex, bu_hex, libceed.VECTOR_ACTIVE) 1298 op_mass_hex.set_field("v", ru_hex, bu_hex, libceed.VECTOR_ACTIVE) 1299 1300 # ------------------------- Composite Operators ------------------------- 1301 1302 # Setup 1303 op_setup = ceed.CompositeOperator() 1304 op_setup.add_sub(op_setup_tet) 1305 op_setup.add_sub(op_setup_hex) 1306 1307 # Apply mass matrix 1308 op_mass = ceed.CompositeOperator() 1309 op_mass.add_sub(op_mass_tet) 1310 op_mass.add_sub(op_mass_hex) 1311 1312 # View 1313 print(op_setup) 1314 print(op_mass) 1315 1316 stdout, stderr, ref_stdout = check.output(capsys) 1317 assert not stderr 1318 assert stdout == ref_stdout 1319 1320# ------------------------------------------------------------------------------- 1321# CeedOperatorApplyAdd for composite operator 1322# ------------------------------------------------------------------------------- 1323 1324 1325def test_524(ceed_resource): 1326 ceed = libceed.Ceed(ceed_resource) 1327 1328 nelem_tet, p_tet, q_tet = 6, 6, 4 1329 nelem_hex, p_hex, q_hex = 6, 3, 4 1330 nx, ny = 3, 3 1331 dim = 2 1332 nx_tet, ny_tet, nx_hex = 3, 1, 3 1333 ndofs = (nx * 2 + 1) * (ny * 2 + 1) 1334 nqpts_tet, nqpts_hex = nelem_tet * q_tet, nelem_hex * q_hex * q_hex 1335 1336 # Vectors 1337 x = ceed.Vector(dim * ndofs) 1338 x_array = np.zeros(dim * ndofs, dtype=ceed.scalar_type()) 1339 for i in range(ny * 2 + 1): 1340 for j in range(nx * 2 + 1): 1341 x_array[i + j * (ny * 2 + 1)] = i / (2 * ny) 1342 x_array[i + j * (ny * 2 + 1) + ndofs] = j / (2 * nx) 1343 x.set_array(x_array, cmode=libceed.USE_POINTER) 1344 1345 qdata_hex = ceed.Vector(nqpts_hex) 1346 qdata_tet = ceed.Vector(nqpts_tet) 1347 u = ceed.Vector(ndofs) 1348 v = ceed.Vector(ndofs) 1349 1350 # ------------------------- Tet Elements ------------------------- 1351 1352 # Restrictions 1353 indx_tet = np.zeros(nelem_tet * p_tet, dtype="int32") 1354 for i in range(nelem_tet // 2): 1355 col = i % nx 1356 row = i // nx 1357 offset = col * 2 + row * (nx * 2 + 1) * 2 1358 1359 indx_tet[i * 2 * p_tet + 0] = 2 + offset 1360 indx_tet[i * 2 * p_tet + 1] = 9 + offset 1361 indx_tet[i * 2 * p_tet + 2] = 16 + offset 1362 indx_tet[i * 2 * p_tet + 3] = 1 + offset 1363 indx_tet[i * 2 * p_tet + 4] = 8 + offset 1364 indx_tet[i * 2 * p_tet + 5] = 0 + offset 1365 1366 indx_tet[i * 2 * p_tet + 6] = 14 + offset 1367 indx_tet[i * 2 * p_tet + 7] = 7 + offset 1368 indx_tet[i * 2 * p_tet + 8] = 0 + offset 1369 indx_tet[i * 2 * p_tet + 9] = 15 + offset 1370 indx_tet[i * 2 * p_tet + 10] = 8 + offset 1371 indx_tet[i * 2 * p_tet + 11] = 16 + offset 1372 1373 rx_tet = ceed.ElemRestriction(nelem_tet, p_tet, dim, ndofs, dim * ndofs, 1374 indx_tet, cmode=libceed.USE_POINTER) 1375 1376 ru_tet = ceed.ElemRestriction(nelem_tet, p_tet, 1, 1, ndofs, indx_tet, 1377 cmode=libceed.USE_POINTER) 1378 strides = np.array([1, q_tet, q_tet], dtype="int32") 1379 rui_tet = ceed.StridedElemRestriction(nelem_tet, q_tet, 1, nqpts_tet, 1380 strides) 1381 1382 # Bases 1383 qref = np.empty(dim * q_tet, dtype=ceed.scalar_type()) 1384 qweight = np.empty(q_tet, dtype=ceed.scalar_type()) 1385 interp, grad = bm.buildmats(qref, qweight, libceed.scalar_types[ 1386 libceed.lib.CEED_SCALAR_TYPE]) 1387 1388 bx_tet = ceed.BasisH1(libceed.TRIANGLE, dim, p_tet, q_hex, interp, grad, qref, 1389 qweight) 1390 bu_tet = ceed.BasisH1(libceed.TRIANGLE, 1, p_tet, q_hex, interp, grad, qref, 1391 qweight) 1392 1393 # QFunctions 1394 file_dir = os.path.dirname(os.path.abspath(__file__)) 1395 qfs = load_qfs_so() 1396 1397 qf_setup_tet = ceed.QFunction(1, qfs.setup_mass_2d, 1398 os.path.join(file_dir, "test-qfunctions.h:setup_mass_2d")) 1399 qf_setup_tet.add_input("weights", 1, libceed.EVAL_WEIGHT) 1400 qf_setup_tet.add_input("dx", dim * dim, libceed.EVAL_GRAD) 1401 qf_setup_tet.add_output("rho", 1, libceed.EVAL_NONE) 1402 1403 qf_mass_tet = ceed.QFunction(1, qfs.apply_mass, 1404 os.path.join(file_dir, "test-qfunctions.h:apply_mass")) 1405 qf_mass_tet.add_input("rho", 1, libceed.EVAL_NONE) 1406 qf_mass_tet.add_input("u", 1, libceed.EVAL_INTERP) 1407 qf_mass_tet.add_output("v", 1, libceed.EVAL_INTERP) 1408 1409 # Operators 1410 op_setup_tet = ceed.Operator(qf_setup_tet) 1411 op_setup_tet.set_field("weights", libceed.ELEMRESTRICTION_NONE, bx_tet, 1412 libceed.VECTOR_NONE) 1413 op_setup_tet.set_field("dx", rx_tet, bx_tet, libceed.VECTOR_ACTIVE) 1414 op_setup_tet.set_field("rho", rui_tet, libceed.BASIS_COLLOCATED, 1415 qdata_tet) 1416 1417 op_mass_tet = ceed.Operator(qf_mass_tet) 1418 op_mass_tet.set_field("rho", rui_tet, libceed.BASIS_COLLOCATED, qdata_tet) 1419 op_mass_tet.set_field("u", ru_tet, bu_tet, libceed.VECTOR_ACTIVE) 1420 op_mass_tet.set_field("v", ru_tet, bu_tet, libceed.VECTOR_ACTIVE) 1421 1422 # ------------------------- Hex Elements ------------------------- 1423 1424 # Restrictions 1425 indx_hex = np.zeros(nelem_hex * p_hex * p_hex, dtype="int32") 1426 for i in range(nelem_hex): 1427 col = i % nx_hex 1428 row = i // nx_hex 1429 offset = (nx_tet * 2 + 1) * (ny_tet * 2) * (1 + row) + col * 2 1430 1431 for j in range(p_hex): 1432 for k in range(p_hex): 1433 indx_hex[p_hex * (p_hex * i + k) + j] = offset + \ 1434 k * (nx_hex * 2 + 1) + j 1435 1436 rx_hex = ceed.ElemRestriction(nelem_hex, p_hex * p_hex, dim, ndofs, 1437 dim * ndofs, indx_hex, 1438 cmode=libceed.USE_POINTER) 1439 1440 ru_hex = ceed.ElemRestriction(nelem_hex, p_hex * p_hex, 1, 1, ndofs, 1441 indx_hex, cmode=libceed.USE_POINTER) 1442 strides = np.array([1, q_hex * q_hex, q_hex * q_hex], dtype="int32") 1443 rui_hex = ceed.StridedElemRestriction(nelem_hex, q_hex * q_hex, 1, 1444 nqpts_hex, strides) 1445 1446 # Bases 1447 bx_hex = ceed.BasisTensorH1Lagrange(dim, dim, p_hex, q_hex, libceed.GAUSS) 1448 bu_hex = ceed.BasisTensorH1Lagrange(dim, 1, p_hex, q_hex, libceed.GAUSS) 1449 1450 # QFunctions 1451 qf_setup_hex = ceed.QFunction(1, qfs.setup_mass_2d, 1452 os.path.join(file_dir, "test-qfunctions.h:setup_mass_2d")) 1453 qf_setup_hex.add_input("weights", 1, libceed.EVAL_WEIGHT) 1454 qf_setup_hex.add_input("dx", dim * dim, libceed.EVAL_GRAD) 1455 qf_setup_hex.add_output("rho", 1, libceed.EVAL_NONE) 1456 1457 qf_mass_hex = ceed.QFunction(1, qfs.apply_mass, 1458 os.path.join(file_dir, "test-qfunctions.h:apply_mass")) 1459 qf_mass_hex.add_input("rho", 1, libceed.EVAL_NONE) 1460 qf_mass_hex.add_input("u", 1, libceed.EVAL_INTERP) 1461 qf_mass_hex.add_output("v", 1, libceed.EVAL_INTERP) 1462 1463 # Operators 1464 op_setup_hex = ceed.Operator(qf_setup_tet) 1465 op_setup_hex.set_field("weights", libceed.ELEMRESTRICTION_NONE, bx_hex, 1466 libceed.VECTOR_NONE) 1467 op_setup_hex.set_field("dx", rx_hex, bx_hex, libceed.VECTOR_ACTIVE) 1468 op_setup_hex.set_field("rho", rui_hex, libceed.BASIS_COLLOCATED, 1469 qdata_hex) 1470 1471 op_mass_hex = ceed.Operator(qf_mass_hex) 1472 op_mass_hex.set_field("rho", rui_hex, libceed.BASIS_COLLOCATED, qdata_hex) 1473 op_mass_hex.set_field("u", ru_hex, bu_hex, libceed.VECTOR_ACTIVE) 1474 op_mass_hex.set_field("v", ru_hex, bu_hex, libceed.VECTOR_ACTIVE) 1475 1476 # ------------------------- Composite Operators ------------------------- 1477 1478 # Setup 1479 op_setup = ceed.CompositeOperator() 1480 op_setup.add_sub(op_setup_tet) 1481 op_setup.add_sub(op_setup_hex) 1482 op_setup.apply(x, libceed.VECTOR_NONE) 1483 1484 # Apply mass matrix 1485 op_mass = ceed.CompositeOperator() 1486 op_mass.add_sub(op_mass_tet) 1487 op_mass.add_sub(op_mass_hex) 1488 u.set_value(1.) 1489 op_mass.apply(u, v) 1490 1491 # Check 1492 with v.array_read() as v_array: 1493 total = 0.0 1494 for i in range(ndofs): 1495 total = total + v_array[i] 1496 assert abs(total - 1.0) < 10. * TOL 1497 1498 # ApplyAdd mass matrix 1499 v.set_value(1.) 1500 op_mass.apply_add(u, v) 1501 1502 # Check 1503 with v.array_read() as v_array: 1504 total = -ndofs 1505 for i in range(ndofs): 1506 total = total + v_array[i] 1507 assert abs(total - 1.0) < 10. * TOL 1508 1509# ------------------------------------------------------------------------------- 1510# Test assembly of mass matrix operator diagonal 1511# ------------------------------------------------------------------------------- 1512 1513 1514def test_533(ceed_resource): 1515 ceed = libceed.Ceed(ceed_resource) 1516 1517 nelem = 6 1518 p = 3 1519 q = 4 1520 dim = 2 1521 nx = 3 1522 ny = 2 1523 ndofs = (nx * 2 + 1) * (ny * 2 + 1) 1524 nqpts = nelem * q * q 1525 1526 # Vectors 1527 x = ceed.Vector(dim * ndofs) 1528 x_array = np.zeros(dim * ndofs) 1529 for i in range(nx * 2 + 1): 1530 for j in range(ny * 2 + 1): 1531 x_array[i + j * (nx * 2 + 1) + 0 * ndofs] = i / (2 * nx) 1532 x_array[i + j * (nx * 2 + 1) + 1 * ndofs] = j / (2 * ny) 1533 x.set_array(x_array, cmode=libceed.USE_POINTER) 1534 1535 qdata = ceed.Vector(nqpts) 1536 u = ceed.Vector(ndofs) 1537 v = ceed.Vector(ndofs) 1538 1539 # Restrictions 1540 indx = np.zeros(nelem * p * p, dtype="int32") 1541 for i in range(nelem): 1542 col = i % nx 1543 row = i // nx 1544 offset = col * (p - 1) + row * (nx * 2 + 1) * (p - 1) 1545 for j in range(p): 1546 for k in range(p): 1547 indx[p * (p * i + k) + j] = offset + k * (nx * 2 + 1) + j 1548 rx = ceed.ElemRestriction(nelem, p * p, dim, ndofs, dim * ndofs, 1549 indx, cmode=libceed.USE_POINTER) 1550 1551 ru = ceed.ElemRestriction(nelem, p * p, 1, 1, ndofs, indx, 1552 cmode=libceed.USE_POINTER) 1553 strides = np.array([1, q * q, q * q], dtype="int32") 1554 rui = ceed.StridedElemRestriction(nelem, q * q, 1, nqpts, strides) 1555 1556 # Bases 1557 bx = ceed.BasisTensorH1Lagrange(dim, dim, p, q, libceed.GAUSS) 1558 bu = ceed.BasisTensorH1Lagrange(dim, 1, p, q, libceed.GAUSS) 1559 1560 # QFunctions 1561 qf_setup = ceed.QFunctionByName("Mass2DBuild") 1562 1563# ------------------------------------------------------------------------------- 1564# Test creation, action, and destruction for mass matrix operator with 1565# multigrid level, tensor basis and interpolation basis generation 1566# ------------------------------------------------------------------------------- 1567 1568 1569def test_550(ceed_resource): 1570 ceed = libceed.Ceed(ceed_resource) 1571 1572 nelem = 15 1573 p_coarse = 3 1574 p_fine = 5 1575 q = 8 1576 ncomp = 2 1577 nx = nelem + 1 1578 nu_coarse = nelem * (p_coarse - 1) + 1 1579 nu_fine = nelem * (p_fine - 1) + 1 1580 1581 # Vectors 1582 x = ceed.Vector(nx) 1583 x_array = np.zeros(nx, dtype=ceed.scalar_type()) 1584 for i in range(nx): 1585 x_array[i] = i / (nx - 1.0) 1586 x.set_array(x_array, cmode=libceed.USE_POINTER) 1587 1588 qdata = ceed.Vector(nelem * q) 1589 u_coarse = ceed.Vector(ncomp * nu_coarse) 1590 v_coarse = ceed.Vector(ncomp * nu_coarse) 1591 u_fine = ceed.Vector(ncomp * nu_fine) 1592 v_fine = ceed.Vector(ncomp * nu_fine) 1593 1594 # Restrictions 1595 indx = np.zeros(nx * 2, dtype="int32") 1596 for i in range(nx): 1597 indx[2 * i + 0] = i 1598 indx[2 * i + 1] = i + 1 1599 rx = ceed.ElemRestriction(nelem, 2, 1, 1, nx, indx, 1600 cmode=libceed.USE_POINTER) 1601 1602 indu_coarse = np.zeros(nelem * p_coarse, dtype="int32") 1603 for i in range(nelem): 1604 for j in range(p_coarse): 1605 indu_coarse[p_coarse * i + j] = i * (p_coarse - 1) + j 1606 ru_coarse = ceed.ElemRestriction(nelem, p_coarse, ncomp, nu_coarse, 1607 ncomp * nu_coarse, indu_coarse, 1608 cmode=libceed.USE_POINTER) 1609 1610 indu_fine = np.zeros(nelem * p_fine, dtype="int32") 1611 for i in range(nelem): 1612 for j in range(p_fine): 1613 indu_fine[p_fine * i + j] = i * (p_fine - 1) + j 1614 ru_fine = ceed.ElemRestriction(nelem, p_fine, ncomp, nu_fine, 1615 ncomp * nu_fine, indu_fine, 1616 cmode=libceed.USE_POINTER) 1617 strides = np.array([1, q, q], dtype="int32") 1618 1619 rui = ceed.StridedElemRestriction(nelem, q, 1, q * nelem, strides) 1620 1621 # Bases 1622 bx = ceed.BasisTensorH1Lagrange(1, 1, 2, q, libceed.GAUSS) 1623 bu_coarse = ceed.BasisTensorH1Lagrange(1, ncomp, p_coarse, q, libceed.GAUSS) 1624 bu_fine = ceed.BasisTensorH1Lagrange(1, ncomp, p_fine, q, libceed.GAUSS) 1625 1626 # QFunctions 1627 file_dir = os.path.dirname(os.path.abspath(__file__)) 1628 qfs = load_qfs_so() 1629 1630 qf_setup = ceed.QFunction(1, qfs.setup_mass, 1631 os.path.join(file_dir, "test-qfunctions.h:setup_mass")) 1632 qf_setup.add_input("weights", 1, libceed.EVAL_WEIGHT) 1633 qf_setup.add_input("dx", 1, libceed.EVAL_GRAD) 1634 qf_setup.add_output("rho", 1, libceed.EVAL_NONE) 1635 1636 qf_mass = ceed.QFunction(1, qfs.apply_mass_two, 1637 os.path.join(file_dir, "test-qfunctions.h:apply_mass_two")) 1638 qf_mass.add_input("rho", 1, libceed.EVAL_NONE) 1639 qf_mass.add_input("u", ncomp, libceed.EVAL_INTERP) 1640 qf_mass.add_output("v", ncomp, libceed.EVAL_INTERP) 1641 1642 # Operators 1643 op_setup = ceed.Operator(qf_setup) 1644 op_setup.set_field("weights", libceed.ELEMRESTRICTION_NONE, bx, 1645 libceed.VECTOR_NONE) 1646 op_setup.set_field("dx", rx, bx, libceed.VECTOR_ACTIVE) 1647 op_setup.set_field("rho", rui, libceed.BASIS_COLLOCATED, 1648 libceed.VECTOR_ACTIVE) 1649 1650 op_mass_fine = ceed.Operator(qf_mass) 1651 op_mass_fine.set_field("rho", rui, libceed.BASIS_COLLOCATED, qdata) 1652 op_mass_fine.set_field("u", ru_fine, bu_fine, libceed.VECTOR_ACTIVE) 1653 op_mass_fine.set_field("v", ru_fine, bu_fine, libceed.VECTOR_ACTIVE) 1654 1655 # Setup 1656 op_setup.apply(x, qdata) 1657 1658 # Create multigrid level 1659 p_mult_fine = ceed.Vector(ncomp * nu_fine) 1660 p_mult_fine.set_value(1.0) 1661 [op_mass_coarse, op_prolong, op_restrict] = op_mass_fine.multigrid_create(p_mult_fine, 1662 ru_coarse, 1663 bu_coarse) 1664 1665 # Apply coarse mass matrix 1666 u_coarse.set_value(1.0) 1667 op_mass_coarse.apply(u_coarse, v_coarse) 1668 1669 # Check 1670 with v_coarse.array_read() as v_array: 1671 total = 0.0 1672 for i in range(nu_coarse * ncomp): 1673 total = total + v_array[i] 1674 assert abs(total - 2.0) < 10. * TOL 1675 1676 # Prolong coarse u 1677 op_prolong.apply(u_coarse, u_fine) 1678 1679 # Apply mass matrix 1680 op_mass_fine.apply(u_fine, v_fine) 1681 1682 # Check 1683 with v_fine.array_read() as v_array: 1684 total = 0.0 1685 for i in range(nu_fine * ncomp): 1686 total = total + v_array[i] 1687 assert abs(total - 2.0) < 10. * TOL 1688 1689 # Restrict state to coarse grid 1690 op_restrict.apply(v_fine, v_coarse) 1691 1692 # Check 1693 with v_coarse.array_read() as v_array: 1694 total = 0.0 1695 for i in range(nu_coarse * ncomp): 1696 total = total + v_array[i] 1697 assert abs(total - 2.0) < 10. * TOL 1698 1699# ------------------------------------------------------------------------------- 1700# Test creation, action, and destruction for mass matrix operator with 1701# multigrid level, tensor basis 1702# ------------------------------------------------------------------------------- 1703 1704 1705def test_552(ceed_resource): 1706 ceed = libceed.Ceed(ceed_resource) 1707 1708 nelem = 15 1709 p_coarse = 3 1710 p_fine = 5 1711 q = 8 1712 ncomp = 2 1713 nx = nelem + 1 1714 nu_coarse = nelem * (p_coarse - 1) + 1 1715 nu_fine = nelem * (p_fine - 1) + 1 1716 1717 # Vectors 1718 x = ceed.Vector(nx) 1719 x_array = np.zeros(nx, dtype=ceed.scalar_type()) 1720 for i in range(nx): 1721 x_array[i] = i / (nx - 1.0) 1722 x.set_array(x_array, cmode=libceed.USE_POINTER) 1723 1724 qdata = ceed.Vector(nelem * q) 1725 u_coarse = ceed.Vector(ncomp * nu_coarse) 1726 v_coarse = ceed.Vector(ncomp * nu_coarse) 1727 u_fine = ceed.Vector(ncomp * nu_fine) 1728 v_fine = ceed.Vector(ncomp * nu_fine) 1729 1730 # Restrictions 1731 indx = np.zeros(nx * 2, dtype="int32") 1732 for i in range(nx): 1733 indx[2 * i + 0] = i 1734 indx[2 * i + 1] = i + 1 1735 rx = ceed.ElemRestriction(nelem, 2, 1, 1, nx, indx, 1736 cmode=libceed.USE_POINTER) 1737 1738 indu_coarse = np.zeros(nelem * p_coarse, dtype="int32") 1739 for i in range(nelem): 1740 for j in range(p_coarse): 1741 indu_coarse[p_coarse * i + j] = i * (p_coarse - 1) + j 1742 ru_coarse = ceed.ElemRestriction(nelem, p_coarse, ncomp, nu_coarse, 1743 ncomp * nu_coarse, indu_coarse, 1744 cmode=libceed.USE_POINTER) 1745 1746 indu_fine = np.zeros(nelem * p_fine, dtype="int32") 1747 for i in range(nelem): 1748 for j in range(p_fine): 1749 indu_fine[p_fine * i + j] = i * (p_fine - 1) + j 1750 ru_fine = ceed.ElemRestriction(nelem, p_fine, ncomp, nu_fine, 1751 ncomp * nu_fine, indu_fine, 1752 cmode=libceed.USE_POINTER) 1753 strides = np.array([1, q, q], dtype="int32") 1754 1755 rui = ceed.StridedElemRestriction(nelem, q, 1, q * nelem, strides) 1756 1757 # Bases 1758 bx = ceed.BasisTensorH1Lagrange(1, 1, 2, q, libceed.GAUSS) 1759 bu_coarse = ceed.BasisTensorH1Lagrange(1, ncomp, p_coarse, q, libceed.GAUSS) 1760 bu_fine = ceed.BasisTensorH1Lagrange(1, ncomp, p_fine, q, libceed.GAUSS) 1761 1762 # QFunctions 1763 file_dir = os.path.dirname(os.path.abspath(__file__)) 1764 qfs = load_qfs_so() 1765 1766 qf_setup = ceed.QFunction(1, qfs.setup_mass, 1767 os.path.join(file_dir, "test-qfunctions.h:setup_mass")) 1768 qf_setup.add_input("weights", 1, libceed.EVAL_WEIGHT) 1769 qf_setup.add_input("dx", 1, libceed.EVAL_GRAD) 1770 qf_setup.add_output("rho", 1, libceed.EVAL_NONE) 1771 1772 qf_mass = ceed.QFunction(1, qfs.apply_mass_two, 1773 os.path.join(file_dir, "test-qfunctions.h:apply_mass_two")) 1774 qf_mass.add_input("rho", 1, libceed.EVAL_NONE) 1775 qf_mass.add_input("u", ncomp, libceed.EVAL_INTERP) 1776 qf_mass.add_output("v", ncomp, libceed.EVAL_INTERP) 1777 1778 # Operators 1779 op_setup = ceed.Operator(qf_setup) 1780 op_setup.set_field("weights", libceed.ELEMRESTRICTION_NONE, bx, 1781 libceed.VECTOR_NONE) 1782 op_setup.set_field("dx", rx, bx, libceed.VECTOR_ACTIVE) 1783 op_setup.set_field("rho", rui, libceed.BASIS_COLLOCATED, 1784 libceed.VECTOR_ACTIVE) 1785 1786 op_mass_fine = ceed.Operator(qf_mass) 1787 op_mass_fine.set_field("rho", rui, libceed.BASIS_COLLOCATED, qdata) 1788 op_mass_fine.set_field("u", ru_fine, bu_fine, libceed.VECTOR_ACTIVE) 1789 op_mass_fine.set_field("v", ru_fine, bu_fine, libceed.VECTOR_ACTIVE) 1790 1791 # Setup 1792 op_setup.apply(x, qdata) 1793 1794 # Create multigrid level 1795 p_mult_fine = ceed.Vector(ncomp * nu_fine) 1796 p_mult_fine.set_value(1.0) 1797 b_c_to_f = ceed.BasisTensorH1Lagrange( 1798 1, ncomp, p_coarse, p_fine, libceed.GAUSS_LOBATTO) 1799 interp_C_to_F = b_c_to_f.get_interp_1d() 1800 [op_mass_coarse, op_prolong, op_restrict] = op_mass_fine.multigrid_create_tensor_h1(p_mult_fine, 1801 ru_coarse, bu_coarse, interp_C_to_F) 1802 1803 # Apply coarse mass matrix 1804 u_coarse.set_value(1.0) 1805 op_mass_coarse.apply(u_coarse, v_coarse) 1806 1807 # Check 1808 with v_coarse.array_read() as v_array: 1809 total = 0.0 1810 for i in range(nu_coarse * ncomp): 1811 total = total + v_array[i] 1812 assert abs(total - 2.0) < TOL 1813 1814 # Prolong coarse u 1815 op_prolong.apply(u_coarse, u_fine) 1816 1817 # Apply mass matrix 1818 op_mass_fine.apply(u_fine, v_fine) 1819 1820 # Check 1821 with v_fine.array_read() as v_array: 1822 total = 0.0 1823 for i in range(nu_fine * ncomp): 1824 total = total + v_array[i] 1825 assert abs(total - 2.0) < TOL 1826 1827 # Restrict state to coarse grid 1828 op_restrict.apply(v_fine, v_coarse) 1829 1830 # Check 1831 with v_coarse.array_read() as v_array: 1832 total = 0.0 1833 for i in range(nu_coarse * ncomp): 1834 total = total + v_array[i] 1835 assert abs(total - 2.0) < TOL 1836 1837# ------------------------------------------------------------------------------- 1838# Test creation, action, and destruction for mass matrix operator with 1839# multigrid level, non-tensor basis 1840# ------------------------------------------------------------------------------- 1841 1842 1843def test_553(ceed_resource): 1844 ceed = libceed.Ceed(ceed_resource) 1845 1846 nelem = 15 1847 p_coarse = 3 1848 p_fine = 5 1849 q = 8 1850 ncomp = 1 1851 nx = nelem + 1 1852 nu_coarse = nelem * (p_coarse - 1) + 1 1853 nu_fine = nelem * (p_fine - 1) + 1 1854 1855 # Vectors 1856 x = ceed.Vector(nx) 1857 x_array = np.zeros(nx, dtype=ceed.scalar_type()) 1858 for i in range(nx): 1859 x_array[i] = i / (nx - 1.0) 1860 x.set_array(x_array, cmode=libceed.USE_POINTER) 1861 1862 qdata = ceed.Vector(nelem * q) 1863 u_coarse = ceed.Vector(ncomp * nu_coarse) 1864 v_coarse = ceed.Vector(ncomp * nu_coarse) 1865 u_fine = ceed.Vector(ncomp * nu_fine) 1866 v_fine = ceed.Vector(ncomp * nu_fine) 1867 1868 # Restrictions 1869 indx = np.zeros(nx * 2, dtype="int32") 1870 for i in range(nx): 1871 indx[2 * i + 0] = i 1872 indx[2 * i + 1] = i + 1 1873 rx = ceed.ElemRestriction(nelem, 2, 1, 1, nx, indx, 1874 cmode=libceed.USE_POINTER) 1875 1876 indu_coarse = np.zeros(nelem * p_coarse, dtype="int32") 1877 for i in range(nelem): 1878 for j in range(p_coarse): 1879 indu_coarse[p_coarse * i + j] = i * (p_coarse - 1) + j 1880 ru_coarse = ceed.ElemRestriction(nelem, p_coarse, ncomp, nu_coarse, 1881 ncomp * nu_coarse, indu_coarse, 1882 cmode=libceed.USE_POINTER) 1883 1884 indu_fine = np.zeros(nelem * p_fine, dtype="int32") 1885 for i in range(nelem): 1886 for j in range(p_fine): 1887 indu_fine[p_fine * i + j] = i * (p_fine - 1) + j 1888 ru_fine = ceed.ElemRestriction(nelem, p_fine, ncomp, nu_fine, 1889 ncomp * nu_fine, indu_fine, 1890 cmode=libceed.USE_POINTER) 1891 strides = np.array([1, q, q], dtype="int32") 1892 1893 rui = ceed.StridedElemRestriction(nelem, q, 1, q * nelem, strides) 1894 1895 # Bases 1896 bx = ceed.BasisTensorH1Lagrange(1, 1, 2, q, libceed.GAUSS) 1897 bu_coarse = ceed.BasisTensorH1Lagrange(1, ncomp, p_coarse, q, libceed.GAUSS) 1898 bu_fine = ceed.BasisTensorH1Lagrange(1, ncomp, p_fine, q, libceed.GAUSS) 1899 1900 # QFunctions 1901 file_dir = os.path.dirname(os.path.abspath(__file__)) 1902 qfs = load_qfs_so() 1903 1904 qf_setup = ceed.QFunctionByName("Mass1DBuild") 1905 qf_mass = ceed.QFunctionByName("MassApply") 1906 1907 # Operators 1908 op_setup = ceed.Operator(qf_setup) 1909 op_setup.set_field("weights", libceed.ELEMRESTRICTION_NONE, bx, 1910 libceed.VECTOR_NONE) 1911 op_setup.set_field("dx", rx, bx, libceed.VECTOR_ACTIVE) 1912 op_setup.set_field("qdata", rui, libceed.BASIS_COLLOCATED, 1913 libceed.VECTOR_ACTIVE) 1914 1915 op_mass_fine = ceed.Operator(qf_mass) 1916 op_mass_fine.set_field("qdata", rui, libceed.BASIS_COLLOCATED, qdata) 1917 op_mass_fine.set_field("u", ru_fine, bu_fine, libceed.VECTOR_ACTIVE) 1918 op_mass_fine.set_field("v", ru_fine, bu_fine, libceed.VECTOR_ACTIVE) 1919 1920 # Setup 1921 op_setup.apply(x, qdata) 1922 1923 # Create multigrid level 1924 p_mult_fine = ceed.Vector(ncomp * nu_fine) 1925 p_mult_fine.set_value(1.0) 1926 b_c_to_f = ceed.BasisTensorH1Lagrange( 1927 1, ncomp, p_coarse, p_fine, libceed.GAUSS_LOBATTO) 1928 interp_C_to_F = b_c_to_f.get_interp_1d() 1929 [op_mass_coarse, op_prolong, op_restrict] = op_mass_fine.multigrid_create_h1(p_mult_fine, 1930 ru_coarse, bu_coarse, interp_C_to_F) 1931 1932 # Apply coarse mass matrix 1933 u_coarse.set_value(1.0) 1934 op_mass_coarse.apply(u_coarse, v_coarse) 1935 1936 # Check 1937 with v_coarse.array_read() as v_array: 1938 total = 0.0 1939 for i in range(nu_coarse * ncomp): 1940 total = total + v_array[i] 1941 assert abs(total - 1.0) < TOL 1942 1943 # Prolong coarse u 1944 op_prolong.apply(u_coarse, u_fine) 1945 1946 # Apply mass matrix 1947 op_mass_fine.apply(u_fine, v_fine) 1948 1949 # Check 1950 with v_fine.array_read() as v_array: 1951 total = 0.0 1952 for i in range(nu_fine * ncomp): 1953 total = total + v_array[i] 1954 assert abs(total - 1.0) < TOL 1955 1956 # Restrict state to coarse grid 1957 op_restrict.apply(v_fine, v_coarse) 1958 1959 # Check 1960 with v_coarse.array_read() as v_array: 1961 total = 0.0 1962 for i in range(nu_coarse * ncomp): 1963 total = total + v_array[i] 1964 assert abs(total - 1.0) < TOL 1965 1966# ------------------------------------------------------------------------------- 1967