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