Use libCEED for element-wise error estimate integrals

[sebastiangrimberg]

Resolves #213. Full error estimate calculation now takes place on GPU without host-device transfer.

Performance comparison for estimation time (the loop over mesh elements) on a particularly large example problem (on CPU):

• main: 1209.6 s
• This PR: 6.8 s

[hughcars]

LGTM. Some minor suggestions on function signatures in integrator.cpp, pushed to hughcars/estimator-libceed-dev if of interest. There's some libceed function/files I don't think needed.

[hughcars]

Each branch of these input if is the same, so I think the branch can be removed. The only changing on input seems to occur on the name above. The same thing can be done with the true/false argument on the other method, with a double dispatch, and dropping the input bool from AddOperatorActiveFields:

void AddOperatorActiveInputFields(unsigned int ops, bool input, Ceed ceed,
                             CeedElemRestriction restr, CeedBasis basis, CeedOperator op,
                             std::string name = "u", CeedVector v = CEED_VECTOR_ACTIVE)
{
 AddOperatorActiveFields(ops, ceed, restr, basis, op, name, v);
}
void AddOperatorActiveOutputFields(unsigned int ops, bool input, Ceed ceed,
                             CeedElemRestriction restr, CeedBasis basis, CeedOperator op,
                             std::string name = "v", CeedVector v = CEED_VECTOR_ACTIVE)
{
 AddOperatorActiveFields(ops, ceed, restr, basis, op, name, v);
}

the call sites then drop the true/false similarly.

[sebastiangrimberg]

Adopted these in 70aabf6.

[hughcars]

Every call for these I can see has a compile time known input value, so this could be split to avoid the nested if statements. Then the call sites are more explicit about what's going on, without needing to check the signature for what the true/false is doing.

void AddQFunctionActiveInputs(unsigned int ops, Ceed ceed, CeedBasis basis,
                              CeedQFunction qf, std::string name = "u")
{
  // Add inputs or outputs with evaluation modes for the active vector of a QFunction.
  CeedInt num_comp;
  PalaceCeedCall(ceed, CeedBasisGetNumComponents(basis, &num_comp));
  if (ops & EvalMode::None)
  {
    PalaceCeedCall(ceed, CeedQFunctionAddInput(qf, name.c_str(), num_comp, CEED_EVAL_NONE));
  }
  if (ops & EvalMode::Interp)
  {
    CeedInt q_comp;
    PalaceCeedCall(ceed,
                   CeedBasisGetNumQuadratureComponents(basis, CEED_EVAL_INTERP, &q_comp));
    PalaceCeedCall(
        ceed, CeedQFunctionAddInput(qf, name.c_str(), num_comp * q_comp, CEED_EVAL_INTERP));
  }
  if (ops & EvalMode::Grad)
  {
    CeedInt q_comp;
    PalaceCeedCall(ceed,
                   CeedBasisGetNumQuadratureComponents(basis, CEED_EVAL_GRAD, &q_comp));
    PalaceCeedCall(ceed, CeedQFunctionAddInput(qf, (std::string("grad_") + name).c_str(),
                                               num_comp * q_comp, CEED_EVAL_GRAD));
  }
  if (ops & EvalMode::Div)
  {
    CeedInt q_comp;
    PalaceCeedCall(ceed,
                   CeedBasisGetNumQuadratureComponents(basis, CEED_EVAL_DIV, &q_comp));
    PalaceCeedCall(ceed, CeedQFunctionAddInput(qf, (std::string("div_") + name).c_str(),
                                               num_comp * q_comp, CEED_EVAL_DIV));
  }
  if (ops & EvalMode::Curl)
  {
    CeedInt q_comp;
    PalaceCeedCall(ceed,
                   CeedBasisGetNumQuadratureComponents(basis, CEED_EVAL_CURL, &q_comp));
    PalaceCeedCall(ceed, CeedQFunctionAddInput(qf, (std::string("curl_") + name).c_str(),
                                               num_comp * q_comp, CEED_EVAL_CURL));
  }
}

void AddQFunctionActiveOutputs(unsigned int ops, Ceed ceed, CeedBasis basis,
                               CeedQFunction qf, std::string name = "v")
{
  // Add inputs or outputs with evaluation modes for the active vector of a QFunction.
  CeedInt num_comp;
  PalaceCeedCall(ceed, CeedBasisGetNumComponents(basis, &num_comp));
  if (ops & EvalMode::None)
  {
    PalaceCeedCall(ceed,
                   CeedQFunctionAddOutput(qf, name.c_str(), num_comp, CEED_EVAL_NONE));
  }
  if (ops & EvalMode::Interp)
  {
    CeedInt q_comp;
    PalaceCeedCall(ceed,
                   CeedBasisGetNumQuadratureComponents(basis, CEED_EVAL_INTERP, &q_comp));
    PalaceCeedCall(ceed, CeedQFunctionAddOutput(qf, name.c_str(), num_comp * q_comp,
                                                CEED_EVAL_INTERP));
  }
  if (ops & EvalMode::Grad)
  {
    CeedInt q_comp;
    PalaceCeedCall(ceed,
                   CeedBasisGetNumQuadratureComponents(basis, CEED_EVAL_GRAD, &q_comp));
    PalaceCeedCall(ceed, CeedQFunctionAddOutput(qf, (std::string("grad_") + name).c_str(),
                                                num_comp * q_comp, CEED_EVAL_GRAD));
  }
  if (ops & EvalMode::Div)
  {
    CeedInt q_comp;
    PalaceCeedCall(ceed,
                   CeedBasisGetNumQuadratureComponents(basis, CEED_EVAL_DIV, &q_comp));
    PalaceCeedCall(ceed, CeedQFunctionAddOutput(qf, (std::string("div_") + name).c_str(),
                                                num_comp * q_comp, CEED_EVAL_DIV));
  }
  if (ops & EvalMode::Curl)
  {
    CeedInt q_comp;
    PalaceCeedCall(ceed,
                   CeedBasisGetNumQuadratureComponents(basis, CEED_EVAL_CURL, &q_comp));
    PalaceCeedCall(ceed, CeedQFunctionAddOutput(qf, (std::string("curl_") + name).c_str(),
                                                num_comp * q_comp, CEED_EVAL_CURL));
  }
}

[sebastiangrimberg]

Adopted these in 70aabf6.

[hughcars]

Can coeff[6] be hoisted too? Provided the MultBAx33 calls are done before each load.

[sebastiangrimberg]

By "hoisted" do you mean put into the outer scope? I would think any compiler would not stack allocate two coeff[6] objects here since only one is in scope at a time. Is that not true? I hadn't thought about it but certainly moving it out could reduce local stack/register usage if that's not true.

[hughcars]

I think a sane compiler would notice, you're right with the scope. I was being dumb and thinking of a heap variable. Should be fine as is.

[hughcars]

For these two functions in here, given the coefficients are allowed to be discontinuous, I think the domain spaces are both L2d rather than hdiv and hcurl, for CurlFlux and GradFlux respectively. Judging by how these functions are all written however, I don't believe this matters? I.e. the discontinuity is coming from the coefficient which is evaluated after the basis at the quadrature points. Am I understanding that correctly?

[sebastiangrimberg]

The naming is addressing not the space of integration but just the space of the input functions for the interpolation/differentiation operations. That is, interpolate into physical space (requires the right Piola transformation), and only then multiply by the coefficient. Does that make sense? So the naming addresses the u space, not the space of C * u.

[hughcars]

I believe these functions are holdovers from when the grad flux estimation was h1d yes? Not sure they're still needed given they'd only be applicable if there was only a constant material property coefficient. I don't think specializing the estimation class on that assumption seems worth it.

[sebastiangrimberg]

They aren't needed but I figured they might be worth keeping in case down the line this type of element-wise integration is desired. There are a few other QFunctions that are also currently unused but left just for helpfulness later on.

[hughcars]

Passive here means that neither is modified by the operations correct? so there's no possibility of race conditions/can be easily vectorized.

[sebastiangrimberg]

Exactly, inputs are always read-only anyway. We have an operator d = f(u, v) which isn't nicely representable with a single active input vector. So both are just labeled as "passive" inputs.