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