Proposed improved API?

[tzanio]

Below is a mockup of some proposed "improvements" in the current API. (This is all up for debate of course.)

Some of the rationale for these changes:

• The Create + Set + Complete workflow is self-documenting and flexible if we decide to introduce additional options in the future (e.g. a rectangular operator).
• CeedBasis separates dof and qpt descriptions, allows Lagrange elements e.g. in Chebyshev Lobatto points
• The G and B operators are not independent, so bundling them in a space makes sense
• CeedOperator is an abstraction for just operator action (as opposed to action and assembly). The D matrix (qdata) is stored internally in it.

@jedbrown, what do you think? Are you open to any of these changes?

{
   // FE Basis = Operator B, B^t for a set of Eval modes
   CeedBasis basis;
   CeedBasisCreate(ceed, &basis);
   CeedBasisSetElement(basis, CEED_HEX); // implies 3D
   // Nodal scalar Lagrange basis, H1 or L2
   CeedBasisSetType(basis, CEED_LAGRANGE_BASIS, 3, CEED_GAUSS_LOBATTO);
   CeedBasisSetQuadrature(basis, 7, CEED_GAUSS); // order 7 = 4 points in 1D
   ierr = CeedBasisComplete(basis); // constructs B1d, etc., calls backend

   // Element restriction = Operator G, G^t
   CeedRestriction restriction;
   CeedRestrictionCreate(ceed, &restriction);
   CeedRestrictionSetElemSize(restriction, 64);
   CeedRestrictionSetNumElems(restriction, 8);
   CeedRestrictionSetNumDofs(restriction, 343);
   CeedInt *element_indices;
   // Set element_indices from mesh topology
   CeedRestrictionSetIndexType(restriction, element_indices,
                               CEED_MEM_HOST, CEED_USE_POINTER); // calls backend
   ierr = CeedRestrictionComplete(restriction); // error handling

   // Isoparametric case, otherwise define mesh-specific restriction & basis
   CeedSpace mesh_space;
   CeedSpaceCreate(ceed, &mesh_space);
   CeedSpaceSetRestriction(mesh_space, &restriction);
   CeedSpaceSetBasis(mesh_space, &basis);
   CeedSpaceSetNumComponents(mesh_space, 3);
   ierr = CeedSpaceComplete(mesh_space); // can call backend
   CeedVector mesh_nodes;
   // Set mesh nodes...

   // FE Space = Operators G, B, ...
   CeedSpace space;
   CeedSpaceCreate(ceed, &space);
   CeedSpaceSetRestriction(space, &restriction);
   CeedSpaceSetBasis(space, &basis);
   ierr = CeedSpaceComplete(space); // can call backend

   // FE Operator = Combines the mesh, the space, and kernels for assembly and
   // action of Operator D
   CeedOperator poisson;
   CeedOperatorCreate(ceed, &poisson);
   CeedOperatorSetMesh(poisson, mesh_space, mesh_nodes);
   CeedOperatorSetSpace(poisson, space);

   // Default kernels (set of Q-functions provided by libCEED)
   CeedOperatorSetKernels(poisson, "poisson");

   // Alternative: Custom kernels based on Q-functions
   CeedQFunction qf_poisson3d, qf_buildcoeffs;
   CeedQFunctionCreateInterior(ceed, 8, 1, 10*sizeof(CeedScalar), CEED_EVAL_GRAD,
                               CEED_EVAL_GRAD, f_poisson3d, "ex1.c:f_poisson3d",
                               &qf_poisson3d);
   CeedQFunctionCreateInterior(ceed, 1, 3, 10*sizeof(CeedScalar),
                               CEED_EVAL_INTERP | CEED_EVAL_GRAD, CEED_EVAL_NONE,
                               f_buildcoeffs,
                               "ex1.c:f_buildcoeffs", &qf_buildcoeffs);
   CeedOperatorSetAssemblyKernel(poisson, qf_my_poisson3d_assemble);
   CeedOperatorSetActionKernel(poisson, qf_my_poisson3d_apply);

   ierr = CeedOperatorComplete(poisson); // can call backend

   ierr = CeedOperatorSetup(poisson, CEED_REQUEST_IMMEDIATE); // calls backend

   ierr = CeedOperatorApply(poisson, u, r, CEED_REQUEST_IMMEDIATE); // calls backend
}

[jedbrown]

1. Split constructors into Create, Set, Complete (or Setup). We do this in PETSc and have "convenience" interfaces for some common configurations. I figured there was a good chance we would adopt this approach in Ceed. Changing it is easy and not disruptive (we can preserve convenience interfaces, perhaps with further name specialization). I don't see this as being in the critical path to a usable interface, but would be happy if someone else wants to make the transformation (and document it along the way).

2. I don't like CeedSpace as proposed because I want a single restriction to be able to restrict to many different element topologies. I envision supporting that by making a composite CeedBasis that is indexed by topology and order (for p-adaptivity), then having the CeedElemRestriction label the elements accordingly. I skipped this at the present stage in the interest of simplicity because previous discussions concluded that homogeneous elements was sufficient for now.

3. There are some ways in which I would like to use the same qdata with different operators. It can be something we set (so that it isn't in the CeedOperatorApply interface) but cannot be exclusive.

[tzanio]

1. We are on the same page -- thanks!

2. That will be nice (though none of our codes currently supports these cases). Wouldn't the single restriction you have in mind plus the composite basis be another example of a CeedSpace? The idea was to group the restriction and the basis, to emphasize that they are synchronized, which is still the case for more general elements/order. Can you clarify what reusability is lost?

3. Yes, I can also sharing (parts of) the qdata between different operators will be useful. I did not mean to preclude that in the above interface. After CeedOperatorAssemble one can grab poisson->qdata.

[jedbrown]

2. If it's just a non-exclusive pair, it's just more terminology for the user and more different objects to interact with. There are plenty of cases where different bases have the same restriction. Lagrange and spline bases on hexes is one. And restrictions change under refinement while the basis can remain fixed. I think having them separate isolates concerns in a natural way. How does bundling them together help the user?

3. You also need to be able to set the qdata without calling CeedOperatorAssemble. I envisioned keeping it as an explicit argument because it makes the interface less stateful.

[tzanio]

2. Okay, I'm flexible on this, we can maybe not require the space objects explicitly, but just adjust slightly the above API:

   
   // Describe the mesh space
   CeedOperatorSetMesh(poisson, basis, restriction, num_components, mesh_nodes); // can call backend

// Describe the solution space (input/output for the Operator). For rectangular operators // we can later add CeedOperatorSetDomainSpace and CeedOperatorSetRangeSpace CeedOperatorSetSpace(poisson, basis, restriction, num_components); // can call backend

Is that OK? @v-dobrev what do you think?

3. Can't this be handled with an optional `CeedOperatorSetQData`? This could be an extra factor to the qdata computed by the operator and we can even handle multiple calls that accumulate (multiply) the pointwise entries. After `CeedOperatorAssemble`, one can call `CeedOperatorGetQData` if they want.

Looking at the interface in ex1.c:

https://github.com/CEED/libCEED/blob/master/examples/ex1.c#L107-L110

my concerns are 1) that the Assembly and Action are both cast as "Operator", which may be confusing to new users 2) since it is external, it is not clear that the qdata is a part of the operator. The above design was meant to address these with minimal changes...

[jedbrown]

2. There is no reason coordinates need to be expressed using the same finite element space as the solution variables. Coordinates might be expressed analytically and are really only defined at quadrature points. Nodes are merely a dual space. I'm not wild about making an interface that doesn't cleanly support all of these variants. For that, we would need a concept of a discrete field decoupled from the finite element space, and it would need to handle various forms of periodic and symmetry conditions. I don't want to add it at this time.

3. CeedOperatorSetQData is possible but makes the interface more stateful. I have mixed feelings about it.

3b. I very frequently want to compute qdata as a side-effect of another computation. Caching information to evaluate a Jacobian or adjoint as a side-effect of computing a residual is an example, but also computing sub-systems in multiphysics and evaluating functionals of interest and other postprocessing uses. In the case of ensembles and adjoint time integration, the qdata needs to be a first-class entity and it will be different for every call to CeedOperatorApply.

[tzanio]

2. To clarify -- this was just the isoparametric case, in general the coordinates will come with their own, mesh_basis and mesh_restriction. For this initial release, I think we should focus on the case where the mesh nodes are just another (vector) finite element function (in appropriate space). Handling geometry described by piecewise-smooth analytic function will be cool, but if libCEED cannot do e.g. fast Jacobian evaluations, then maybe that's out of scope?

[jedbrown]

So with my existing spec, they are just another field with no special treatment at all. We will need a reliable way to evaluate coordinates (with gradients) to quadrature points independent of the space that coordinates come from. That's the idea behind what happens in buildcoeffs. The library can provide a function for computing those coordinate transformations so that users don't have to reproduce it. I'd rather not include special handling for coordinates at this time and I'm not sure it will be justified later.

[v-dobrev]

Regarding CeedSpace: I think the CeedSpace expresses a simple concept - description of a discrete space on a mesh. The most common use is to define an L-space but it can definitely be used to define E- and Q-spaces as well. It gives a concrete meaning of a CeedVector as an L-, E-, or Q-vector. As a parameter of a CeedOperator, it tells the operator what type of vectors it operates on and what type of vector the output of the operator Apply function should be. It can be reused by multiple CeedOperators. I think that the meaning of CeedSpace will be much more transparent to users compared to CeedElemRestriction and the fact that different CeedSpaces can share the same CeedElemRestriction or CeedBasis does not diminish the usefulness of CeedSpace. As for the name, CeedSpace, I'm not attached to it; CeedDiscreteSpace may be another option.

[v-dobrev]

Regarding the Q-data: I understand that users may want to reuse/share Q-data (or portions of it) among multiple operators, so adding a type CeedQData that facilitates the sharing is perfectly fine with me. In many cases CeedQData will be simply a multi-component Q-vector, potentially using some specialized layout. However, it can be more general than that - in particular it may even be opaque to libCEED, i.e. it can be specific to a set of operators, sharable among them. On the other hand, there are many cases where users will not need to use shared CeedQData, so I also prefer to hide it inside the CeedOperator - invisible to beginner users but still available to more advanced users.

@jedbrown, I don't understand your objection: what you mean by "more stateful" interface?

[jedbrown]

So we create an Operator and several other objects (forward and adjoint solvers, multi-physics) get a reference to it. Now I'd rather that the interface to use it does not require that I modify the Operator and then need to remember to set it back to its prior configuration so that my usage is not observable to other clients.

Note that the qdata vector is currently opaque to Ceed and can have any layout within batches of quadrature points. I'd rather not impose more structure unless/until we have a very clear reason.

[v-dobrev]

I still don't understand: why do you need to modify the Operator and then set it back? Is this a case where you needs to have multiple Q-data (at the same moment), e.g. corresponding to different physical states and/or different values of t?

[jedbrown]

Yes, there are many scenarios where that is needed. I'd rather avoid this pattern:

CeedVector qdata_save;
CeedOperatorGetQData(op, &qdata_save);
CeedOperatorSetQData(op, qdata);
CeedOperatorApply(op, ...);
CeedOperatorSetQData(op, qdata_save);

[jeremylt]

Closed by PR #41 (I'm pretty sure)

[v-dobrev]

@jeremylt this is an old discussion with suggestions for API changes that #41 does not address.

[jedbrown]

I think too much of this discussion is obsolete given API changes that have occurred since. Can you create a new issue for any outstanding concerns that you would like to address?

[v-dobrev]

Most of the example above is still valid. The only thing that #41 changes is how the "assemble" and "apply" kernels (actually called "operators" in the master) are constructed - if constructed explicitly. The above example proposes (among other things) a way to access operators already defined by libCEED, e.g.

  CeedOperatorSetKernels(poisson, "poisson");

This type of functionality is something we still need to add.

[jedbrown]

See #37 ?

[v-dobrev]

That is just one of the things proposed in the example above.

[jedbrown]

This thread has too many non-orthogonal suggestions and the discussion is too biased by the most disruptive of those, which was addressed in #41. Please make new issues at a more modest granularity, as has been done for #37. (You have been largely silent in this project for the past few months. We welcome your comments, but this thread is too broad in scope and too stale for productive conversation.)