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