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 Basis functionality 19 20import os 21import math 22import libceed 23import numpy as np 24import buildmats as bm 25import check 26 27TOL = np.finfo(float).eps * 256 28 29# ------------------------------------------------------------------------------- 30# Utilities 31# ------------------------------------------------------------------------------- 32 33 34def eval(dim, x): 35 result, center = 1, 0.1 36 for d in range(dim): 37 result *= math.tanh(x[d] - center) 38 center += 0.1 39 return result 40 41 42def feval(x1, x2): 43 return x1 * x1 + x2 * x2 + x1 * x2 + 1 44 45 46def dfeval(x1, x2): 47 return 2 * x1 + x2 48 49# ------------------------------------------------------------------------------- 50# Test creation and distruction of a H1Lagrange basis 51# ------------------------------------------------------------------------------- 52 53 54def test_300(ceed_resource, capsys): 55 ceed = libceed.Ceed(ceed_resource) 56 57 b = ceed.BasisTensorH1Lagrange(1, 1, 4, 4, libceed.GAUSS_LOBATTO) 58 print(b) 59 del b 60 61 b = ceed.BasisTensorH1Lagrange(1, 1, 4, 4, libceed.GAUSS) 62 print(b) 63 del b 64 65 stdout, stderr, ref_stdout = check.output(capsys) 66 assert not stderr 67 assert stdout == ref_stdout 68 69# ------------------------------------------------------------------------------- 70# Test QR factorization 71# ------------------------------------------------------------------------------- 72 73 74def test_301(ceed_resource, capsys): 75 ceed = libceed.Ceed(ceed_resource) 76 77 qr = np.array([1, -1, 4, 1, 4, -2, 1, 4, 2, 1, -1, 0], dtype="float64") 78 tau = np.empty(3, dtype="float64") 79 80 qr, tau = ceed.qr_factorization(qr, tau, 4, 3) 81 82 for i in range(len(qr)): 83 if qr[i] <= TOL and qr[i] >= -TOL: 84 qr[i] = 0 85 print("%12.8f" % qr[i]) 86 87 for i in range(len(tau)): 88 if tau[i] <= TOL and tau[i] >= -TOL: 89 tau[i] = 0 90 print("%12.8f" % tau[i]) 91 92 stdout, stderr, ref_stdout = check.output(capsys) 93 assert not stderr 94 assert stdout == ref_stdout 95 96# ------------------------------------------------------------------------------- 97# Test Symmetric Schur Decomposition 98# ------------------------------------------------------------------------------- 99 100 101def test_304(ceed_resource, capsys): 102 ceed = libceed.Ceed(ceed_resource) 103 104 A = np.array([0.19996678, 0.0745459, -0.07448852, 0.0332866, 105 0.0745459, 1., 0.16666509, -0.07448852, 106 -0.07448852, 0.16666509, 1., 0.0745459, 107 0.0332866, -0.07448852, 0.0745459, 0.19996678], dtype="float64") 108 109 lam = ceed.symmetric_schur_decomposition(A, 4) 110 111 print("Q: ") 112 for i in range(4): 113 for j in range(4): 114 if A[j + 4 * i] <= TOL and A[j + 4 * i] >= -TOL: 115 A[j + 4 * i] = 0 116 print("%12.8f" % A[j + 4 * i]) 117 118 print("lambda: ") 119 for i in range(4): 120 if lam[i] <= TOL and lam[i] >= -TOL: 121 lam[i] = 0 122 print("%12.8f" % lam[i]) 123 124 stdout, stderr, ref_stdout = check.output(capsys) 125 assert not stderr 126 assert stdout == ref_stdout 127 128# ------------------------------------------------------------------------------- 129# Test Simultaneous Diagonalization 130# ------------------------------------------------------------------------------- 131 132 133def test_305(ceed_resource, capsys): 134 ceed = libceed.Ceed(ceed_resource) 135 136 M = np.array([0.19996678, 0.0745459, -0.07448852, 0.0332866, 137 0.0745459, 1., 0.16666509, -0.07448852, 138 -0.07448852, 0.16666509, 1., 0.0745459, 139 0.0332866, -0.07448852, 0.0745459, 0.19996678], dtype="float64") 140 K = np.array([3.03344425, -3.41501767, 0.49824435, -0.11667092, 141 -3.41501767, 5.83354662, -2.9167733, 0.49824435, 142 0.49824435, -2.9167733, 5.83354662, -3.41501767, 143 -0.11667092, 0.49824435, -3.41501767, 3.03344425], dtype="float64") 144 145 x, lam = ceed.simultaneous_diagonalization(K, M, 4) 146 147 print("x: ") 148 for i in range(4): 149 for j in range(4): 150 if x[j + 4 * i] <= TOL and x[j + 4 * i] >= -TOL: 151 x[j + 4 * i] = 0 152 print("%12.8f" % x[j + 4 * i]) 153 154 print("lambda: ") 155 for i in range(4): 156 if lam[i] <= TOL and lam[i] >= -TOL: 157 lam[i] = 0 158 print("%12.8f" % lam[i]) 159 160 stdout, stderr, ref_stdout = check.output(capsys) 161 assert not stderr 162 assert stdout == ref_stdout 163 164# ------------------------------------------------------------------------------- 165# Test GetNumNodes and GetNumQuadraturePoints for basis 166# ------------------------------------------------------------------------------- 167 168 169def test_306(ceed_resource): 170 ceed = libceed.Ceed(ceed_resource) 171 172 b = ceed.BasisTensorH1Lagrange(3, 1, 4, 5, libceed.GAUSS_LOBATTO) 173 174 p = b.get_num_nodes() 175 q = b.get_num_quadrature_points() 176 177 assert p == 64 178 assert q == 125 179 180# ------------------------------------------------------------------------------- 181# Test interpolation in multiple dimensions 182# ------------------------------------------------------------------------------- 183 184 185def test_313(ceed_resource): 186 ceed = libceed.Ceed(ceed_resource) 187 188 for dim in range(1, 4): 189 Q = 10 190 Qdim = Q**dim 191 Xdim = 2**dim 192 x = np.empty(Xdim * dim, dtype="float64") 193 uq = np.empty(Qdim, dtype="float64") 194 195 for d in range(dim): 196 for i in range(Xdim): 197 x[d * Xdim + i] = 1 if (i % 198 (2**(dim - d))) // (2**(dim - d - 1)) else -1 199 200 X = ceed.Vector(Xdim * dim) 201 X.set_array(x, cmode=libceed.USE_POINTER) 202 Xq = ceed.Vector(Qdim * dim) 203 Xq.set_value(0) 204 U = ceed.Vector(Qdim) 205 U.set_value(0) 206 Uq = ceed.Vector(Qdim) 207 208 bxl = ceed.BasisTensorH1Lagrange(dim, dim, 2, Q, libceed.GAUSS_LOBATTO) 209 bul = ceed.BasisTensorH1Lagrange(dim, 1, Q, Q, libceed.GAUSS_LOBATTO) 210 211 bxl.apply(1, libceed.EVAL_INTERP, X, Xq) 212 213 with Xq.array_read() as xq: 214 for i in range(Qdim): 215 xx = np.empty(dim, dtype="float64") 216 for d in range(dim): 217 xx[d] = xq[d * Qdim + i] 218 uq[i] = eval(dim, xx) 219 220 Uq.set_array(uq, cmode=libceed.USE_POINTER) 221 222 # This operation is the identity because the quadrature is collocated 223 bul.T.apply(1, libceed.EVAL_INTERP, Uq, U) 224 225 bxg = ceed.BasisTensorH1Lagrange(dim, dim, 2, Q, libceed.GAUSS) 226 bug = ceed.BasisTensorH1Lagrange(dim, 1, Q, Q, libceed.GAUSS) 227 228 bxg.apply(1, libceed.EVAL_INTERP, X, Xq) 229 bug.apply(1, libceed.EVAL_INTERP, U, Uq) 230 231 with Xq.array_read() as xq, Uq.array_read() as u: 232 for i in range(Qdim): 233 xx = np.empty(dim, dtype="float64") 234 for d in range(dim): 235 xx[d] = xq[d * Qdim + i] 236 fx = eval(dim, xx) 237 assert math.fabs(u[i] - fx) < 1E-4 238 239# ------------------------------------------------------------------------------- 240# Test grad in multiple dimensions 241# ------------------------------------------------------------------------------- 242 243 244def test_314(ceed_resource): 245 ceed = libceed.Ceed(ceed_resource) 246 247 for dim in range(1, 4): 248 P, Q = 8, 10 249 Pdim = P**dim 250 Qdim = Q**dim 251 Xdim = 2**dim 252 sum1 = sum2 = 0 253 x = np.empty(Xdim * dim, dtype="float64") 254 u = np.empty(Pdim, dtype="float64") 255 256 for d in range(dim): 257 for i in range(Xdim): 258 x[d * Xdim + i] = 1 if (i % 259 (2**(dim - d))) // (2**(dim - d - 1)) else -1 260 261 X = ceed.Vector(Xdim * dim) 262 X.set_array(x, cmode=libceed.USE_POINTER) 263 Xq = ceed.Vector(Pdim * dim) 264 Xq.set_value(0) 265 U = ceed.Vector(Pdim) 266 Uq = ceed.Vector(Qdim * dim) 267 Uq.set_value(0) 268 Ones = ceed.Vector(Qdim * dim) 269 Ones.set_value(1) 270 Gtposeones = ceed.Vector(Pdim) 271 Gtposeones.set_value(0) 272 273 # Get function values at quadrature points 274 bxl = ceed.BasisTensorH1Lagrange(dim, dim, 2, P, libceed.GAUSS_LOBATTO) 275 bxl.apply(1, libceed.EVAL_INTERP, X, Xq) 276 277 with Xq.array_read() as xq: 278 for i in range(Pdim): 279 xx = np.empty(dim, dtype="float64") 280 for d in range(dim): 281 xx[d] = xq[d * Pdim + i] 282 u[i] = eval(dim, xx) 283 284 U.set_array(u, cmode=libceed.USE_POINTER) 285 286 # Calculate G u at quadrature points, G' * 1 at dofs 287 bug = ceed.BasisTensorH1Lagrange(dim, 1, P, Q, libceed.GAUSS) 288 bug.apply(1, libceed.EVAL_GRAD, U, Uq) 289 bug.T.apply(1, libceed.EVAL_GRAD, Ones, Gtposeones) 290 291 # Check if 1' * G * u = u' * (G' * 1) 292 with Gtposeones.array_read() as gtposeones, Uq.array_read() as uq: 293 for i in range(Pdim): 294 sum1 += gtposeones[i] * u[i] 295 for i in range(dim * Qdim): 296 sum2 += uq[i] 297 298 assert math.fabs(sum1 - sum2) < 1E-10 299 300# ------------------------------------------------------------------------------- 301# Test creation and destruction of a 2D Simplex non-tensor H1 basis 302# ------------------------------------------------------------------------------- 303 304 305def test_320(ceed_resource): 306 ceed = libceed.Ceed(ceed_resource) 307 308 P, Q, dim = 6, 4, 2 309 310 in_array = np.empty(P, dtype="float64") 311 qref = np.empty(dim * Q, dtype="float64") 312 qweight = np.empty(Q, dtype="float64") 313 314 interp, grad = bm.buildmats(qref, qweight) 315 316 b = ceed.BasisH1(libceed.TRIANGLE, 1, P, Q, interp, grad, qref, qweight) 317 318 print(b) 319 del b 320 321# ------------------------------------------------------------------------------- 322# Test integration with a 2D Simplex non-tensor H1 basis 323# ------------------------------------------------------------------------------- 324 325 326def test_322(ceed_resource): 327 ceed = libceed.Ceed(ceed_resource) 328 329 P, Q, dim = 6, 4, 2 330 331 xr = np.array([0., 0.5, 1., 0., 0.5, 0., 0., 0., 332 0., 0.5, 0.5, 1.], dtype="float64") 333 in_array = np.empty(P, dtype="float64") 334 qref = np.empty(dim * Q, dtype="float64") 335 qweight = np.empty(Q, dtype="float64") 336 337 interp, grad = bm.buildmats(qref, qweight) 338 339 b = ceed.BasisH1(libceed.TRIANGLE, 1, P, Q, interp, grad, qref, qweight) 340 341 # Interpolate function to quadrature points 342 for i in range(P): 343 in_array[i] = feval(xr[0 * P + i], xr[1 * P + i]) 344 345 in_vec = ceed.Vector(P) 346 in_vec.set_array(in_array, cmode=libceed.USE_POINTER) 347 out_vec = ceed.Vector(Q) 348 out_vec.set_value(0) 349 weights_vec = ceed.Vector(Q) 350 weights_vec.set_value(0) 351 352 b.apply(1, libceed.EVAL_INTERP, in_vec, out_vec) 353 b.apply(1, libceed.EVAL_WEIGHT, libceed.VECTOR_NONE, weights_vec) 354 355 # Check values at quadrature points 356 with out_vec.array_read() as out_array, weights_vec.array_read() as weights_array: 357 sum = 0 358 for i in range(Q): 359 sum += out_array[i] * weights_array[i] 360 assert math.fabs(sum - 17. / 24.) < 1E-10 361 362# ------------------------------------------------------------------------------- 363# Test grad with a 2D Simplex non-tensor H1 basis 364# ------------------------------------------------------------------------------- 365 366 367def test_323(ceed_resource): 368 ceed = libceed.Ceed(ceed_resource) 369 370 P, Q, dim = 6, 4, 2 371 372 xq = np.array([0.2, 0.6, 1. / 3., 0.2, 0.2, 0.2, 373 1. / 3., 0.6], dtype="float64") 374 xr = np.array([0., 0.5, 1., 0., 0.5, 0., 0., 0., 375 0., 0.5, 0.5, 1.], dtype="float64") 376 in_array = np.empty(P, dtype="float64") 377 qref = np.empty(dim * Q, dtype="float64") 378 qweight = np.empty(Q, dtype="float64") 379 380 interp, grad = bm.buildmats(qref, qweight) 381 382 b = ceed.BasisH1(libceed.TRIANGLE, 1, P, Q, interp, grad, qref, qweight) 383 384 # Interpolate function to quadrature points 385 for i in range(P): 386 in_array[i] = feval(xr[0 * P + i], xr[1 * P + i]) 387 388 in_vec = ceed.Vector(P) 389 in_vec.set_array(in_array, cmode=libceed.USE_POINTER) 390 out_vec = ceed.Vector(Q * dim) 391 out_vec.set_value(0) 392 393 b.apply(1, libceed.EVAL_GRAD, in_vec, out_vec) 394 395 # Check values at quadrature points 396 with out_vec.array_read() as out_array: 397 for i in range(Q): 398 value = dfeval(xq[0 * Q + i], xq[1 * Q + i]) 399 assert math.fabs(out_array[0 * Q + i] - value) < 1E-10 400 401 value = dfeval(xq[1 * Q + i], xq[0 * Q + i]) 402 assert math.fabs(out_array[1 * Q + i] - value) < 1E-10 403 404# ------------------------------------------------------------------------------- 405