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