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