1 /// @file 2 /// Test projection interp and grad in multiple dimensions 3 /// \test Test projection interp and grad in multiple dimensions 4 #include <ceed.h> 5 #include <math.h> 6 7 static CeedScalar Eval(CeedInt dim, const CeedScalar x[]) { 8 CeedScalar result = (x[0] + 0.1) * (x[0] + 0.1); 9 if (dim > 1) result += (x[1] + 0.2) * (x[1] + 0.2); 10 if (dim > 2) result += -(x[2] + 0.3) * (x[2] + 0.3); 11 return result; 12 } 13 14 static CeedScalar EvalGrad(CeedInt dim, const CeedScalar x[]) { 15 switch (dim) { 16 case 0: 17 return 2 * x[0] + 0.2; 18 case 1: 19 return 2 * x[1] + 0.4; 20 default: 21 return -2 * x[2] - 0.6; 22 } 23 } 24 25 static CeedScalar GetTolerance(CeedScalarType scalar_type, int dim) { 26 CeedScalar tol; 27 if (scalar_type == CEED_SCALAR_FP32) { 28 if (dim == 3) tol = 1.e-4; 29 else tol = 1.e-5; 30 } else { 31 tol = 1.e-11; 32 } 33 return tol; 34 } 35 36 int main(int argc, char **argv) { 37 Ceed ceed; 38 39 CeedInit(argv[1], &ceed); 40 41 for (CeedInt dim = 1; dim <= 3; dim++) { 42 CeedVector X_corners, X_from, X_to, U_from, U_to, dU_to; 43 CeedBasis basis_x, basis_from, basis_to, basis_project; 44 CeedInt P_from = 5, P_to = 6, Q = 7, X_dim = CeedIntPow(2, dim), P_from_dim = CeedIntPow(P_from, dim), P_to_dim = CeedIntPow(P_to, dim); 45 CeedScalar x[X_dim * dim], u_from[P_from_dim]; 46 const CeedScalar *u_to, *du_to, *x_from, *x_to; 47 48 for (CeedInt d = 0; d < dim; d++) { 49 for (CeedInt i = 0; i < X_dim; i++) x[X_dim * d + i] = (i % CeedIntPow(2, dim - d)) / CeedIntPow(2, dim - d - 1) ? 1 : -1; 50 } 51 52 CeedVectorCreate(ceed, X_dim * dim, &X_corners); 53 CeedVectorSetArray(X_corners, CEED_MEM_HOST, CEED_USE_POINTER, x); 54 CeedVectorCreate(ceed, P_from_dim * dim, &X_from); 55 CeedVectorCreate(ceed, P_to_dim * dim, &X_to); 56 CeedVectorCreate(ceed, P_from_dim, &U_from); 57 CeedVectorSetValue(U_from, 0); 58 CeedVectorCreate(ceed, P_to_dim, &U_to); 59 CeedVectorSetValue(U_to, 0); 60 CeedVectorCreate(ceed, P_to_dim * dim, &dU_to); 61 CeedVectorSetValue(dU_to, 0); 62 63 // Get nodal coordinates 64 CeedBasisCreateTensorH1Lagrange(ceed, dim, dim, 2, P_from, CEED_GAUSS_LOBATTO, &basis_x); 65 CeedBasisApply(basis_x, 1, CEED_NOTRANSPOSE, CEED_EVAL_INTERP, X_corners, X_from); 66 CeedBasisDestroy(&basis_x); 67 CeedBasisCreateTensorH1Lagrange(ceed, dim, dim, 2, P_to, CEED_GAUSS_LOBATTO, &basis_x); 68 CeedBasisApply(basis_x, 1, CEED_NOTRANSPOSE, CEED_EVAL_INTERP, X_corners, X_to); 69 CeedBasisDestroy(&basis_x); 70 71 // Create U and projection bases 72 CeedBasisCreateTensorH1Lagrange(ceed, dim, 1, P_from, Q, CEED_GAUSS, &basis_from); 73 CeedBasisCreateTensorH1Lagrange(ceed, dim, 1, P_to, Q, CEED_GAUSS, &basis_to); 74 CeedBasisCreateProjection(basis_from, basis_to, &basis_project); 75 76 // Setup coarse solution 77 CeedVectorGetArrayRead(X_from, CEED_MEM_HOST, &x_from); 78 for (CeedInt i = 0; i < P_from_dim; i++) { 79 CeedScalar xx[dim]; 80 for (CeedInt d = 0; d < dim; d++) xx[d] = x_from[P_from_dim * d + i]; 81 u_from[i] = Eval(dim, xx); 82 } 83 CeedVectorRestoreArrayRead(X_from, &x_from); 84 CeedVectorSetArray(U_from, CEED_MEM_HOST, CEED_USE_POINTER, u_from); 85 86 // Project to fine basis 87 CeedBasisApply(basis_project, 1, CEED_NOTRANSPOSE, CEED_EVAL_INTERP, U_from, U_to); 88 89 // Check solution 90 CeedVectorGetArrayRead(U_to, CEED_MEM_HOST, &u_to); 91 CeedVectorGetArrayRead(X_to, CEED_MEM_HOST, &x_to); 92 CeedScalar tol = GetTolerance(CEED_SCALAR_TYPE, dim); 93 for (CeedInt i = 0; i < P_to_dim; i++) { 94 CeedScalar xx[dim]; 95 for (CeedInt d = 0; d < dim; d++) xx[d] = x_to[d * P_to_dim + i]; 96 const CeedScalar u = Eval(dim, xx); 97 if (fabs(u - u_to[i]) > tol) printf("[%" CeedInt_FMT ", %" CeedInt_FMT "] %f != %f\n", dim, i, u_to[i], u); 98 } 99 CeedVectorRestoreArrayRead(X_to, &x_to); 100 CeedVectorRestoreArrayRead(U_to, &u_to); 101 102 // Project and take gradient 103 CeedBasisApply(basis_project, 1, CEED_NOTRANSPOSE, CEED_EVAL_GRAD, U_from, dU_to); 104 105 // Check solution 106 CeedVectorGetArrayRead(dU_to, CEED_MEM_HOST, &du_to); 107 CeedVectorGetArrayRead(X_to, CEED_MEM_HOST, &x_to); 108 for (CeedInt i = 0; i < P_to_dim; i++) { 109 CeedScalar xx[dim]; 110 for (CeedInt d = 0; d < dim; d++) xx[d] = x_to[P_to_dim * d + i]; 111 for (CeedInt d = 0; d < dim; d++) { 112 const CeedScalar du = EvalGrad(d, xx); 113 if (fabs(du - du_to[P_to_dim * (dim - 1 - d) + i]) > tol) { 114 // LCOV_EXCL_START 115 printf("[%" CeedInt_FMT ", %" CeedInt_FMT ", %" CeedInt_FMT "] %f != %f\n", dim, i, d, du_to[P_to_dim * (dim - 1 - d) + i], du); 116 // LCOV_EXCL_STOP 117 } 118 } 119 } 120 CeedVectorRestoreArrayRead(X_to, &x_to); 121 CeedVectorRestoreArrayRead(dU_to, &du_to); 122 123 CeedVectorDestroy(&X_corners); 124 CeedVectorDestroy(&X_from); 125 CeedVectorDestroy(&X_to); 126 CeedVectorDestroy(&U_from); 127 CeedVectorDestroy(&U_to); 128 CeedVectorDestroy(&dU_to); 129 CeedBasisDestroy(&basis_from); 130 CeedBasisDestroy(&basis_to); 131 CeedBasisDestroy(&basis_project); 132 } 133 CeedDestroy(&ceed); 134 return 0; 135 } 136