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