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