tests: add operator assembly tests for non-trivial oriented restrictions

Author: zatkins-dev

tests/t574-operator.c

@@ -0,0 +1,307 @@
+/// @file
+/// Test full assembly of multi-basis asymmetric mass-like operator with oriented and curl-oriented element restrictions
+/// \test Test full assembly of multi-basis asymmetric mass-like operator with oriented and curl-oriented element restrictions
+#include "t571-operator.h"
+
+#include <ceed.h>
+#include <math.h>
+#include <stdio.h>
+#include <stdlib.h>
+
+int main(int argc, char **argv) {
+  Ceed                ceed;
+  CeedElemRestriction elem_restriction_x, elem_restriction_q_data;
+  CeedElemRestriction oriented_elem_restriction_u, curl_oriented_elem_restriction_u;
+  CeedElemRestriction oriented_elem_restriction_p, curl_oriented_elem_restriction_p;
+  CeedBasis           basis_x, basis_u, basis_p;
+  CeedQFunction       qf_setup, qf_mass;
+  CeedOperator        op_setup, op_mass_oriented, op_mass_curl_oriented;
+  CeedVector          q_data, x, u, v;
+  CeedInt             p_u = 3, p_p = 2, q = 4, dim = 2;
+  CeedInt             n_x = 1, n_y = 1;
+  CeedInt             num_elem   = n_x * n_y;
+  CeedInt             num_comp_u = 1, num_comp_p = 1;  // Must be 1
+  CeedInt             num_dofs_u = (n_x * (p_u - 1) + 1) * (n_y * (p_u - 1) + 1), num_dofs_p = (n_x * (p_p - 1) + 1) * (n_y * (p_p - 1) + 1),
+          num_qpts = num_elem * q * q;
+  CeedInt    ind_x[num_elem * p_u * p_u], ind_p[num_elem * p_p * p_p];
+  CeedInt    l_size = num_comp_u * num_dofs_u + num_comp_p * num_dofs_p;
+  bool       orients_u[num_elem * p_u * p_u], orients_p[num_elem * p_p * p_p];
+  CeedInt8   curl_orients_u[3 * num_elem * p_u * p_u], curl_orients_p[3 * num_elem * p_p * p_p];
+  CeedScalar assembled_values_oriented[l_size * l_size];
+  CeedScalar assembled_values_curl_oriented[l_size * l_size];
+  CeedScalar assembled_true[l_size * l_size];
+
+  CeedInit(argv[1], &ceed);
+
+  // Vectors
+  CeedVectorCreate(ceed, dim * num_dofs_u, &x);
+  {
+    CeedScalar x_array[dim * num_dofs_u];
+
+    for (CeedInt i = 0; i < n_x * (p_u - 1) + 1; i++) {
+      for (CeedInt j = 0; j < n_y * (p_u - 1) + 1; j++) {
+        x_array[i + j * (n_x * 2 + 1) + 0 * num_dofs_u] = (CeedScalar)i / (n_x * (p_u - 1));
+        x_array[i + j * (n_x * 2 + 1) + 1 * num_dofs_u] = (CeedScalar)j / (n_y * (p_u - 1));
+      }
+    }
+    CeedVectorSetArray(x, CEED_MEM_HOST, CEED_COPY_VALUES, x_array);
+  }
+  CeedVectorCreate(ceed, l_size, &u);
+  CeedVectorCreate(ceed, l_size, &v);
+  CeedVectorCreate(ceed, num_qpts, &q_data);
+
+  // Restrictions
+  for (CeedInt i = 0; i < num_elem; i++) {
+    CeedInt col, row, offset;
+
+    col    = i % n_x;
+    row    = i / n_x;
+    offset = col * (p_u - 1) + row * (n_x * (p_u - 1) + 1) * (p_u - 1);
+    for (CeedInt j = 0; j < p_u; j++) {
+      for (CeedInt k = 0; k < p_u; k++) {
+        bool normal = !((j % 2 > 0 && k % 2 == 0) || (j % 2 == 0 && k % 2 > 0));
+
+        ind_x[p_u * (p_u * i + k) + j]                    = offset + k * p_u + j;
+        orients_u[p_u * (p_u * i + k) + j]                = !normal;
+        curl_orients_u[3 * (p_u * (p_u * i + k) + j) + 0] = j % 2 > 0 && k % 2 == 0 ? -1 : 0;
+        curl_orients_u[3 * (p_u * (p_u * i + k) + j) + 1] = normal ? 1 : 0;
+        curl_orients_u[3 * (p_u * (p_u * i + k) + j) + 2] = j % 2 == 0 && k % 2 > 0 ? -1 : 0;
+      }
+    }
+  }
+  const CeedInt offset_p = num_comp_u * num_dofs_u;
+  for (CeedInt i = 0; i < num_elem; i++) {
+    CeedInt col, row, offset;
+
+    col    = i % n_x;
+    row    = i / n_x;
+    offset = col * (p_p - 1) + row * (n_x * (p_p - 1) + 1) * (p_p - 1);
+    for (CeedInt j = 0; j < p_p; j++) {
+      for (CeedInt k = 0; k < p_p; k++) {
+        bool normal = !((j % 2 > 0 && k % 2 == 0) || (j % 2 == 0 && k % 2 > 0));
+
+        ind_p[p_p * (p_p * i + k) + j]                    = offset + k * p_p + j + offset_p;
+        orients_p[p_p * (p_p * i + k) + j]                = !normal;
+        curl_orients_p[3 * (p_p * (p_p * i + k) + j) + 0] = j % 2 == 0 && k % 2 > 0 ? -1 : 0;
+        curl_orients_p[3 * (p_p * (p_p * i + k) + j) + 1] = normal ? 1 : 0;
+        curl_orients_p[3 * (p_p * (p_p * i + k) + j) + 2] = j % 2 > 0 && k % 2 == 0 ? -1 : 0;
+      }
+    }
+  }
+  CeedElemRestrictionCreate(ceed, num_elem, p_u * p_u, dim, num_dofs_u, dim * num_dofs_u, CEED_MEM_HOST, CEED_USE_POINTER, ind_x,
+                            &elem_restriction_x);
+  CeedElemRestrictionCreateOriented(ceed, num_elem, p_u * p_u, num_comp_u, num_dofs_u, l_size, CEED_MEM_HOST, CEED_USE_POINTER, ind_x, orients_u,
+                                    &oriented_elem_restriction_u);
+  CeedElemRestrictionCreateCurlOriented(ceed, num_elem, p_u * p_u, num_comp_u, num_dofs_u, l_size, CEED_MEM_HOST, CEED_USE_POINTER, ind_x,
+                                        curl_orients_u, &curl_oriented_elem_restriction_u);
+
+  CeedElemRestrictionCreateOriented(ceed, num_elem, p_p * p_p, num_comp_p, num_dofs_p, l_size, CEED_MEM_HOST, CEED_USE_POINTER, ind_x, orients_p,
+                                    &oriented_elem_restriction_p);
+  CeedElemRestrictionCreateCurlOriented(ceed, num_elem, p_p * p_p, num_comp_p, num_dofs_p, l_size, CEED_MEM_HOST, CEED_USE_POINTER, ind_x,
+                                        curl_orients_p, &curl_oriented_elem_restriction_p);
+
+  CeedInt strides_q_data[3] = {1, q * q * num_elem, q * q};
+  CeedElemRestrictionCreateStrided(ceed, num_elem, q * q, 1, num_qpts * 1, strides_q_data, &elem_restriction_q_data);
+
+  // Bases
+  CeedBasisCreateTensorH1Lagrange(ceed, dim, dim, p_u, q, CEED_GAUSS, &basis_x);
+  CeedBasisCreateTensorH1Lagrange(ceed, dim, num_comp_u, p_u, q, CEED_GAUSS, &basis_u);
+  CeedBasisCreateTensorH1Lagrange(ceed, dim, num_comp_p, p_p, q, CEED_GAUSS, &basis_p);
+
+  // QFunctions
+  CeedQFunctionCreateInterior(ceed, 1, setup, setup_loc, &qf_setup);
+  CeedQFunctionAddInput(qf_setup, "weight", 1, CEED_EVAL_WEIGHT);
+  CeedQFunctionAddInput(qf_setup, "dx", dim * dim, CEED_EVAL_GRAD);
+  CeedQFunctionAddOutput(qf_setup, "qdata", 1, CEED_EVAL_NONE);
+
+  CeedQFunctionCreateInterior(ceed, 1, multi_basis, multi_basis_loc, &qf_mass);
+  CeedQFunctionAddInput(qf_mass, "qdata", 1, CEED_EVAL_NONE);
+  CeedQFunctionAddInput(qf_mass, "u", num_comp_u, CEED_EVAL_INTERP);
+  CeedQFunctionAddInput(qf_mass, "p", num_comp_p, CEED_EVAL_INTERP);
+  CeedQFunctionAddOutput(qf_mass, "v", num_comp_u, CEED_EVAL_INTERP);
+  CeedQFunctionAddOutput(qf_mass, "q", num_comp_p, CEED_EVAL_INTERP);
+
+  Context_t            context = {num_comp_u, num_comp_p};
+  CeedQFunctionContext ctx;
+
+  CeedQFunctionContextCreate(ceed, &ctx);
+  CeedQFunctionContextSetData(ctx, CEED_MEM_HOST, CEED_COPY_VALUES, sizeof(Context_t), &context);
+  CeedQFunctionSetContext(qf_mass, ctx);
+  CeedQFunctionContextDestroy(&ctx);
+
+  // Operators
+  CeedOperatorCreate(ceed, qf_setup, CEED_QFUNCTION_NONE, CEED_QFUNCTION_NONE, &op_setup);
+  CeedOperatorSetField(op_setup, "weight", CEED_ELEMRESTRICTION_NONE, basis_x, CEED_VECTOR_NONE);
+  CeedOperatorSetField(op_setup, "dx", elem_restriction_x, basis_x, CEED_VECTOR_ACTIVE);
+  CeedOperatorSetField(op_setup, "qdata", elem_restriction_q_data, CEED_BASIS_NONE, CEED_VECTOR_ACTIVE);
+
+  CeedOperatorCreate(ceed, qf_mass, CEED_QFUNCTION_NONE, CEED_QFUNCTION_NONE, &op_mass_oriented);
+  CeedOperatorSetField(op_mass_oriented, "qdata", elem_restriction_q_data, CEED_BASIS_NONE, q_data);
+  CeedOperatorSetField(op_mass_oriented, "u", oriented_elem_restriction_u, basis_u, CEED_VECTOR_ACTIVE);
+  CeedOperatorSetField(op_mass_oriented, "p", oriented_elem_restriction_p, basis_p, CEED_VECTOR_ACTIVE);
+  CeedOperatorSetField(op_mass_oriented, "v", oriented_elem_restriction_u, basis_u, CEED_VECTOR_ACTIVE);
+  CeedOperatorSetField(op_mass_oriented, "q", oriented_elem_restriction_p, basis_p, CEED_VECTOR_ACTIVE);
+
+  CeedOperatorCreate(ceed, qf_mass, CEED_QFUNCTION_NONE, CEED_QFUNCTION_NONE, &op_mass_curl_oriented);
+  CeedOperatorSetField(op_mass_curl_oriented, "qdata", elem_restriction_q_data, CEED_BASIS_NONE, q_data);
+  CeedOperatorSetField(op_mass_curl_oriented, "u", curl_oriented_elem_restriction_u, basis_u, CEED_VECTOR_ACTIVE);
+  CeedOperatorSetField(op_mass_curl_oriented, "p", curl_oriented_elem_restriction_p, basis_p, CEED_VECTOR_ACTIVE);
+  CeedOperatorSetField(op_mass_curl_oriented, "v", curl_oriented_elem_restriction_u, basis_u, CEED_VECTOR_ACTIVE);
+  CeedOperatorSetField(op_mass_curl_oriented, "q", curl_oriented_elem_restriction_p, basis_p, CEED_VECTOR_ACTIVE);
+
+  // Apply Setup Operator
+  CeedOperatorApply(op_setup, x, q_data, CEED_REQUEST_IMMEDIATE);
+
+  // Fully assemble operators
+  CeedSize   num_entries_oriented, num_entries_curl_oriented;
+  CeedInt   *rows_oriented, *rows_curl_oriented;
+  CeedInt   *cols_oriented, *cols_curl_oriented;
+  CeedVector assembled_oriented, assembled_curl_oriented;
+
+  for (CeedInt k = 0; k < l_size * l_size; ++k) {
+    assembled_values_oriented[k]      = 0.0;
+    assembled_values_curl_oriented[k] = 0.0;
+  }
+  CeedOperatorLinearAssembleSymbolic(op_mass_oriented, &num_entries_oriented, &rows_oriented, &cols_oriented);
+  CeedOperatorLinearAssembleSymbolic(op_mass_curl_oriented, &num_entries_curl_oriented, &rows_curl_oriented, &cols_curl_oriented);
+  CeedVectorCreate(ceed, num_entries_oriented, &assembled_oriented);
+  CeedVectorCreate(ceed, num_entries_curl_oriented, &assembled_curl_oriented);
+  CeedOperatorLinearAssemble(op_mass_oriented, assembled_oriented);
+  CeedOperatorLinearAssemble(op_mass_curl_oriented, assembled_curl_oriented);
+  {
+    const CeedScalar *assembled_array;
+
+    CeedVectorGetArrayRead(assembled_oriented, CEED_MEM_HOST, &assembled_array);
+    for (CeedInt k = 0; k < num_entries_oriented; ++k) {
+      assembled_values_oriented[rows_oriented[k] * l_size + cols_oriented[k]] += assembled_array[k];
+    }
+    CeedVectorRestoreArrayRead(assembled_oriented, &assembled_array);
+  }
+  {
+    const CeedScalar *assembled_array;
+
+    CeedVectorGetArrayRead(assembled_curl_oriented, CEED_MEM_HOST, &assembled_array);
+    for (CeedInt k = 0; k < num_entries_curl_oriented; ++k) {
+      assembled_values_curl_oriented[rows_curl_oriented[k] * l_size + cols_curl_oriented[k]] += assembled_array[k];
+    }
+    CeedVectorRestoreArrayRead(assembled_curl_oriented, &assembled_array);
+  }
+
+  // Manually assemble oriented operator
+  CeedInt old_index = -1;
+
+  CeedVectorSetValue(u, 0.0);
+  for (CeedInt j = 0; j < l_size; j++) {
+    CeedScalar       *u_array;
+    const CeedScalar *v_array;
+
+    // Set input
+    CeedVectorGetArray(u, CEED_MEM_HOST, &u_array);
+    u_array[j] = 1.0;
+    if (j > 0) u_array[old_index] = 0.0;
+    old_index = j;
+    CeedVectorRestoreArray(u, &u_array);
+
+    // Compute effect of DoF ind
+    CeedOperatorApply(op_mass_oriented, u, v, CEED_REQUEST_IMMEDIATE);
+
+    CeedVectorGetArrayRead(v, CEED_MEM_HOST, &v_array);
+    for (CeedInt i = 0; i < l_size; i++) assembled_true[i * l_size + j] = v_array[i];
+    CeedVectorRestoreArrayRead(v, &v_array);
+  }
+
+  // Check output
+  for (CeedInt i = 0; i < l_size; i++) {
+    for (CeedInt j = 0; j < l_size; j++) {
+      const CeedScalar assembled_value      = assembled_values_oriented[i * l_size + j];
+      const CeedScalar assembled_true_value = assembled_true[i * l_size + j];
+      const CeedScalar error                = fabs(assembled_value - assembled_true_value) / (fmax(fabs(assembled_true_value), 1.0));
+
+      if (!(error < 100. * CEED_EPSILON)) {
+        // LCOV_EXCL_START
+        const CeedInt node_out = (i < num_comp_u * num_dofs_u) ? i % num_dofs_u : (i - num_comp_u * num_dofs_u) % num_dofs_p;
+        const CeedInt comp_out = (i < num_comp_u * num_dofs_u) ? i / num_dofs_u : (i - num_comp_u * num_dofs_u) / num_dofs_p;
+        const CeedInt node_in  = (j < num_comp_u * num_dofs_u) ? j % num_dofs_u : (j - num_comp_u * num_dofs_u) % num_dofs_p;
+        const CeedInt comp_in  = (j < num_comp_u * num_dofs_u) ? j / num_dofs_u : (j - num_comp_u * num_dofs_u) / num_dofs_p;
+
+        printf("[(%s, %" CeedInt_FMT ", %" CeedInt_FMT "), (%s, %" CeedInt_FMT ", %" CeedInt_FMT
+               ")] Error in oriented assembly: %g != %g (error: %g)\n",
+               (i < num_comp_u * num_dofs_u ? "u" : "p"), node_out, comp_out, (j < num_comp_u * num_dofs_u ? "u" : "p"), node_in, comp_in,
+               assembled_value, assembled_true_value, error);
+        // LCOV_EXCL_STOP
+      }
+    }
+  }
+
+  // Manually assemble curl-oriented operator
+  old_index = -1;
+
+  CeedVectorSetValue(u, 0.0);
+  for (CeedInt j = 0; j < l_size; j++) {
+    CeedScalar       *u_array;
+    const CeedScalar *v_array;
+
+    // Set input
+    CeedVectorGetArray(u, CEED_MEM_HOST, &u_array);
+    u_array[j] = 1.0;
+    if (j > 0) u_array[old_index] = 0.0;
+    old_index = j;
+    CeedVectorRestoreArray(u, &u_array);
+
+    // Compute effect of DoF ind
+    CeedOperatorApply(op_mass_curl_oriented, u, v, CEED_REQUEST_IMMEDIATE);
+
+    CeedVectorGetArrayRead(v, CEED_MEM_HOST, &v_array);
+    for (CeedInt i = 0; i < l_size; i++) assembled_true[i * l_size + j] = v_array[i];
+    CeedVectorRestoreArrayRead(v, &v_array);
+  }
+  for (CeedInt i = 0; i < l_size; i++) {
+    for (CeedInt j = 0; j < l_size; j++) {
+      const CeedScalar assembled_value      = assembled_values_curl_oriented[i * l_size + j];
+      const CeedScalar assembled_true_value = assembled_true[i * l_size + j];
+      const CeedScalar error                = fabs(assembled_value - assembled_true_value) / (fmax(fabs(assembled_true_value), 1.0));
+
+      if (!(error < 100. * CEED_EPSILON)) {
+        // LCOV_EXCL_START
+        const CeedInt node_out = (i < num_comp_u * num_dofs_u) ? i % num_dofs_u : (i - num_comp_u * num_dofs_u) % num_dofs_p;
+        const CeedInt comp_out = (i < num_comp_u * num_dofs_u) ? i / num_dofs_u : (i - num_comp_u * num_dofs_u) / num_dofs_p;
+        const CeedInt node_in  = (j < num_comp_u * num_dofs_u) ? j % num_dofs_u : (j - num_comp_u * num_dofs_u) % num_dofs_p;
+        const CeedInt comp_in  = (j < num_comp_u * num_dofs_u) ? j / num_dofs_u : (j - num_comp_u * num_dofs_u) / num_dofs_p;
+
+        printf("[(%s, %" CeedInt_FMT ", %" CeedInt_FMT "), (%s, %" CeedInt_FMT ", %" CeedInt_FMT
+               ")] Error in curl-oriented assembly: %g != %g (error: %g)\n",
+               (i < num_comp_u * num_dofs_u ? "u" : "p"), node_out, comp_out, (j < num_comp_u * num_dofs_u ? "u" : "p"), node_in, comp_in,
+               assembled_value, assembled_true_value, error);
+        // LCOV_EXCL_STOP
+      }
+    }
+  }
+
+  // Cleanup
+  free(rows_oriented);
+  free(cols_oriented);
+  free(rows_curl_oriented);
+  free(cols_curl_oriented);
+  CeedVectorDestroy(&x);
+  CeedVectorDestroy(&q_data);
+  CeedVectorDestroy(&u);
+  CeedVectorDestroy(&v);
+  CeedVectorDestroy(&assembled_oriented);
+  CeedVectorDestroy(&assembled_curl_oriented);
+  CeedElemRestrictionDestroy(&oriented_elem_restriction_u);
+  CeedElemRestrictionDestroy(&curl_oriented_elem_restriction_u);
+  CeedElemRestrictionDestroy(&oriented_elem_restriction_p);
+  CeedElemRestrictionDestroy(&curl_oriented_elem_restriction_p);
+  CeedElemRestrictionDestroy(&elem_restriction_x);
+  CeedElemRestrictionDestroy(&elem_restriction_q_data);
+  CeedBasisDestroy(&basis_u);
+  CeedBasisDestroy(&basis_p);
+  CeedBasisDestroy(&basis_x);
+  CeedQFunctionDestroy(&qf_setup);
+  CeedQFunctionDestroy(&qf_mass);
+  CeedOperatorDestroy(&op_setup);
+  CeedOperatorDestroy(&op_mass_oriented);
+  CeedOperatorDestroy(&op_mass_curl_oriented);
+  CeedDestroy(&ceed);
+  return 0;
+}