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