1 /// @file 2 /// Test gradient transpose ApplyAdd in multiple dimensions at arbitrary points 3 /// \test Test gradient transpose ApplyAdd in multiple dimensions at arbitrary points 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 = tanh(x[0] + 0.1); 10 if (dim > 1) result += atan(x[1] + 0.2); 11 if (dim > 2) result += exp(-(x[2] + 0.3) * (x[2] + 0.3)); 12 return result; 13 } 14 15 static CeedScalar GetTolerance(CeedScalarType scalar_type, int dim) { 16 CeedScalar tol; 17 if (scalar_type == CEED_SCALAR_FP32) { 18 if (dim == 3) tol = 0.005; 19 else tol = 1.e-4; 20 } else { 21 tol = 1.e-11; 22 } 23 return tol; 24 } 25 26 int main(int argc, char **argv) { 27 Ceed ceed; 28 29 CeedInit(argv[1], &ceed); 30 31 for (CeedInt dim = 1; dim <= 3; dim++) { 32 CeedVector x, x_nodes, x_points, u, u_points, v, ones; 33 CeedBasis basis_x, basis_u; 34 const CeedInt p = 9, q = 9, num_points = 4, x_dim = CeedIntPow(2, dim), p_dim = CeedIntPow(p, dim); 35 CeedScalar sum_1 = 0, sum_2 = 0; 36 37 CeedVectorCreate(ceed, x_dim * dim, &x); 38 CeedVectorCreate(ceed, p_dim * dim, &x_nodes); 39 CeedVectorCreate(ceed, num_points * dim, &x_points); 40 CeedVectorCreate(ceed, p_dim, &u); 41 CeedVectorCreate(ceed, num_points * dim, &u_points); 42 CeedVectorCreate(ceed, p_dim, &v); 43 CeedVectorCreate(ceed, num_points * dim, &ones); 44 45 CeedVectorSetValue(ones, 1); 46 CeedVectorSetValue(v, 0); 47 48 // Get nodal coordinates 49 CeedBasisCreateTensorH1Lagrange(ceed, dim, dim, 2, p, CEED_GAUSS_LOBATTO, &basis_x); 50 { 51 CeedScalar x_array[x_dim * dim]; 52 53 for (CeedInt d = 0; d < dim; d++) { 54 for (CeedInt i = 0; i < x_dim; i++) x_array[d * x_dim + i] = (i % CeedIntPow(2, d + 1)) / CeedIntPow(2, d) ? 1 : -1; 55 } 56 CeedVectorSetArray(x, CEED_MEM_HOST, CEED_COPY_VALUES, x_array); 57 } 58 CeedBasisApply(basis_x, 1, CEED_NOTRANSPOSE, CEED_EVAL_INTERP, x, x_nodes); 59 60 // Set values of u at nodes 61 { 62 const CeedScalar *x_array; 63 CeedScalar u_array[p_dim]; 64 65 CeedVectorGetArrayRead(x_nodes, CEED_MEM_HOST, &x_array); 66 for (CeedInt i = 0; i < p_dim; i++) { 67 CeedScalar coord[dim]; 68 69 for (CeedInt d = 0; d < dim; d++) coord[d] = x_array[d * p_dim + i]; 70 u_array[i] = Eval(dim, coord); 71 } 72 CeedVectorRestoreArrayRead(x_nodes, &x_array); 73 CeedVectorSetArray(u, CEED_MEM_HOST, CEED_COPY_VALUES, (CeedScalar *)&u_array); 74 } 75 76 // Interpolate to arbitrary points 77 CeedBasisCreateTensorH1Lagrange(ceed, dim, 1, p, q, CEED_GAUSS, &basis_u); 78 { 79 CeedScalar x_array[12] = {-0.33, -0.65, 0.16, 0.99, -0.65, 0.16, 0.99, -0.33, 0.16, 0.99, -0.33, -0.65}; 80 81 CeedVectorSetArray(x_points, CEED_MEM_HOST, CEED_COPY_VALUES, x_array); 82 } 83 84 // Calculate G u at arbitrary points, G' * 1 at dofs 85 CeedBasisApplyAtPoints(basis_u, 1, &num_points, CEED_NOTRANSPOSE, CEED_EVAL_GRAD, x_points, u, u_points); 86 CeedBasisApplyAtPoints(basis_u, 1, &num_points, CEED_TRANSPOSE, CEED_EVAL_GRAD, x_points, ones, v); 87 // Double it 88 CeedBasisApplyAddAtPoints(basis_u, 1, &num_points, CEED_TRANSPOSE, CEED_EVAL_GRAD, x_points, ones, v); 89 { 90 const CeedScalar *u_array, *v_array, *u_points_array; 91 92 CeedVectorGetArrayRead(u, CEED_MEM_HOST, &u_array); 93 CeedVectorGetArrayRead(v, CEED_MEM_HOST, &v_array); 94 CeedVectorGetArrayRead(u_points, CEED_MEM_HOST, &u_points_array); 95 for (CeedInt i = 0; i < p_dim; i++) sum_1 += v_array[i] * u_array[i]; 96 for (CeedInt i = 0; i < num_points * dim; i++) sum_2 += u_points_array[i]; 97 CeedVectorRestoreArrayRead(u, &u_array); 98 CeedVectorRestoreArrayRead(v, &v_array); 99 CeedVectorRestoreArrayRead(u_points, &u_points_array); 100 } 101 { 102 CeedScalarType scalar_type; 103 104 CeedGetScalarType(&scalar_type); 105 106 CeedScalar tol = GetTolerance(scalar_type, dim); 107 108 if (fabs(sum_1 - 2.0 * sum_2) > tol) printf("[%" CeedInt_FMT "] %f != %f\n", dim, sum_1, 2.0 * sum_2); 109 } 110 111 CeedVectorDestroy(&x); 112 CeedVectorDestroy(&x_nodes); 113 CeedVectorDestroy(&x_points); 114 CeedVectorDestroy(&u); 115 CeedVectorDestroy(&u_points); 116 CeedVectorDestroy(&ones); 117 CeedVectorDestroy(&v); 118 CeedBasisDestroy(&basis_x); 119 CeedBasisDestroy(&basis_u); 120 } 121 CeedDestroy(&ceed); 122 return 0; 123 } 124