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 46 CeedVectorCreate(ceed, x_dim * dim, &x_corners); 47 { 48 CeedScalar x_array[x_dim * dim]; 49 50 for (CeedInt d = 0; d < dim; d++) { 51 for (CeedInt i = 0; i < x_dim; i++) x_array[x_dim * d + i] = (i % CeedIntPow(2, dim - d)) / CeedIntPow(2, dim - d - 1) ? 1 : -1; 52 } 53 CeedVectorSetArray(x_corners, CEED_MEM_HOST, CEED_COPY_VALUES, x_array); 54 } 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, CEED_GAUSS_LOBATTO, &basis_x); 66 CeedBasisApply(basis_x, 1, CEED_NOTRANSPOSE, CEED_EVAL_INTERP, x_corners, x_from); 67 CeedBasisDestroy(&basis_x); 68 CeedBasisCreateTensorH1Lagrange(ceed, dim, dim, 2, p_to, CEED_GAUSS_LOBATTO, &basis_x); 69 CeedBasisApply(basis_x, 1, CEED_NOTRANSPOSE, CEED_EVAL_INTERP, x_corners, x_to); 70 CeedBasisDestroy(&basis_x); 71 72 // Create U and projection bases 73 CeedBasisCreateTensorH1Lagrange(ceed, dim, 1, p_from, q, CEED_GAUSS, &basis_from); 74 CeedBasisCreateTensorH1Lagrange(ceed, dim, 1, p_to, q, CEED_GAUSS, &basis_to); 75 CeedBasisCreateProjection(basis_from, basis_to, &basis_project); 76 77 // Setup coarse solution 78 { 79 const CeedScalar *x_array; 80 CeedScalar u_array[p_from_dim]; 81 82 CeedVectorGetArrayRead(x_from, CEED_MEM_HOST, &x_array); 83 for (CeedInt i = 0; i < p_from_dim; i++) { 84 CeedScalar coord[dim]; 85 for (CeedInt d = 0; d < dim; d++) coord[d] = x_array[p_from_dim * d + i]; 86 u_array[i] = Eval(dim, coord); 87 } 88 CeedVectorRestoreArrayRead(x_from, &x_array); 89 CeedVectorSetArray(u_from, CEED_MEM_HOST, CEED_COPY_VALUES, u_array); 90 } 91 92 // Project to fine basis 93 CeedBasisApply(basis_project, 1, CEED_NOTRANSPOSE, CEED_EVAL_INTERP, u_from, u_to); 94 95 // Check solution 96 CeedScalar tol = GetTolerance(CEED_SCALAR_TYPE, dim); 97 { 98 const CeedScalar *x_array, *u_array; 99 100 CeedVectorGetArrayRead(x_to, CEED_MEM_HOST, &x_array); 101 CeedVectorGetArrayRead(u_to, CEED_MEM_HOST, &u_array); 102 for (CeedInt i = 0; i < p_to_dim; i++) { 103 CeedScalar coord[dim]; 104 for (CeedInt d = 0; d < dim; d++) coord[d] = x_array[d * p_to_dim + i]; 105 const CeedScalar u = Eval(dim, coord); 106 if (fabs(u - u_array[i]) > tol) printf("[%" CeedInt_FMT ", %" CeedInt_FMT "] %f != %f\n", dim, i, u_array[i], u); 107 } 108 CeedVectorRestoreArrayRead(x_to, &x_array); 109 CeedVectorRestoreArrayRead(u_to, &u_array); 110 } 111 112 // Project and take gradient 113 CeedBasisApply(basis_project, 1, CEED_NOTRANSPOSE, CEED_EVAL_GRAD, u_from, du_to); 114 115 // Check solution 116 { 117 const CeedScalar *x_array, *du_array; 118 119 CeedVectorGetArrayRead(x_to, CEED_MEM_HOST, &x_array); 120 CeedVectorGetArrayRead(du_to, CEED_MEM_HOST, &du_array); 121 for (CeedInt i = 0; i < p_to_dim; i++) { 122 CeedScalar coord[dim]; 123 for (CeedInt d = 0; d < dim; d++) coord[d] = x_array[p_to_dim * d + i]; 124 for (CeedInt d = 0; d < dim; d++) { 125 const CeedScalar du = EvalGrad(d, coord); 126 if (fabs(du - du_array[p_to_dim * (dim - 1 - d) + i]) > tol) { 127 // LCOV_EXCL_START 128 printf("[%" CeedInt_FMT ", %" CeedInt_FMT ", %" CeedInt_FMT "] %f != %f\n", dim, i, d, du_array[p_to_dim * (dim - 1 - d) + i], du); 129 // LCOV_EXCL_STOP 130 } 131 } 132 } 133 CeedVectorRestoreArrayRead(x_to, &x_array); 134 CeedVectorRestoreArrayRead(du_to, &du_array); 135 } 136 137 CeedVectorDestroy(&x_corners); 138 CeedVectorDestroy(&x_from); 139 CeedVectorDestroy(&x_to); 140 CeedVectorDestroy(&u_from); 141 CeedVectorDestroy(&u_to); 142 CeedVectorDestroy(&du_to); 143 CeedBasisDestroy(&basis_from); 144 CeedBasisDestroy(&basis_to); 145 CeedBasisDestroy(&basis_project); 146 } 147 CeedDestroy(&ceed); 148 return 0; 149 } 150