Postprocessing through global vectors?

[kronbichler]

After reviewing #369 and looking through the code base for compressible + incompressible NS, I feel that we should discuss the design of the postprocessors. Right now, we create global vectors for every field we write into output. While this ok for DG elements, it feels less good for H(div) because we project onto the finite element space, rather than keeping the actual form of the quantity. Also, I do not like the idea of the big state that gets created from this setup. So the general question is whether we want to keep this in the long run. It feels we are duplicating a lot of effort from the deal.II class PostProcessor that can be directly attached to DataOut. Let me cite part of the comment https://github.com/exadg/exadg/pull/349#discussion_r1130660192:

This means that I would intuitively not create a global vector first and then write it somewhere, but use the DataPostprocessor facilities of deal.II for that, see https://www.dealii.org/developer/doxygen/deal.II/classDataPostprocessor.html Even though there might be some extensions (or performance optimizations) necessary in that module, the overall feeling is that we should not duplicate such functionality in ExaDG unless there are good reasons.

Let us discuss pros and cons on these two options here. Comparing the current design, I would classify it as follows:

• Pro: It is better performance than dealii::DataOut because it uses FEEvaluation rather than FEValues interfaces. Without checking closely, I think this could be fixed in deal.II.
• Con: We have a lot of infrastructure in the basic operator, which is not really nice. Let me give as example just the part in the compressible NS part and only in the operator, there is even more in the incompressible operator: https://github.com/exadg/exadg/blob/f8e533cad5fdd2fb1c6845657d5fda59228590e4/include/exadg/compressible_navier_stokes/spatial_discretization/operator.h#L127-L132 and https://github.com/exadg/exadg/blob/f8e533cad5fdd2fb1c6845657d5fda59228590e4/include/exadg/compressible_navier_stokes/spatial_discretization/operator.h#L136-L158 and https://github.com/exadg/exadg/blob/f8e533cad5fdd2fb1c6845657d5fda59228590e4/include/exadg/compressible_navier_stokes/spatial_discretization/operator.h#L216-L217 and https://github.com/exadg/exadg/blob/f8e533cad5fdd2fb1c6845657d5fda59228590e4/include/exadg/compressible_navier_stokes/spatial_discretization/operator.h#L225-L226 Of course, this code would not completely go away, but it could be reduced to local computation at the level of points and thus be more modular.
• The projection on the finite element space mentioned above is ad-hoc - one might get better quality sometimes due to superconvergence, but otherwise it can hide the actual form of a derived quantity.

So my question is:

• Are there uses where the global vectors are really necessary? It would be good to have a plan for all uses of the infrastructure, possibly beyond postprocessing/output, that require this interface?
• How much effort would we estimate for such a change? Of course only if we are in favor.

[kronbichler]

Let me add that I was more skeptical against dealii::PostProcessor some years ago (because part of the infrastructure is so incredibly slow for high orders), but I would like to have an open discussion if we should spend considerable effort here or if it is better spent elsewhere, and what the future costs would be.

[nfehn]

I understand this issue as related to classes that we call OutputGenerator in ExaDG. Maybe we should restrict the issue name to this type of postprocessing functionality (?), because postprocessing has a broader notion in ExaDG than OutputGenerator.

[nfehn]

Yes, you are right, computing a global vector can have many reasons, more than for writing vtu output. See e.g. the compressible Navier-Stokes postprocessing and the complex if-statements appearing there: https://github.com/exadg/exadg/blob/f8e533cad5fdd2fb1c6845657d5fda59228590e4/include/exadg/compressible_navier_stokes/postprocessor/postprocessor.cpp#L214-L218

[richardschu]

Let me summarize my pamphlet from #349 for postprocessing on codim 1 manifolds (boundaries or the interface in FSI):

1. the vectors might in fact be used: maybe for (i) timestepping (hemodynamic indicators are integrated in time for temporal averages), (ii) to project quantities (depending on the flow solver, the traction might be needed as a continuous quantitity, see the PPE solvers I used), (iii) in FSI coupling, working with smaller (interface-only) vectors, requires less work than using the volumetric vectors (I did not look at the implementation there!)

And in the present case, I think about output/postprocessing: 2) this happens only every ~1000 steps? So I think dumping a lot of time into that might not be the best investment? Does that really hold us back? I did in fact realize that output is expensive, but not overwhelming.

[nfehn]

Superconvergence is an important aspect. We have realized that the recently introduced normal-flux calculator (see #333) is not optimal in accuracy and that there seem to be better options with superconvergence property (@bergbauer mentioned the work https://onlinelibrary.wiley.com/doi/abs/10.1002/fld.1650070406, but I did not find time to look into it in detail). Wall-shear stress calculations seem to be similar in terms of computing "normal gradients", so this could be relevant for @richardschu as well (is something like this realized in #349?).

My general feeling is that postprocessing functionality is very much problem-dependent in my experience. So developing generic routines is something nice, but most probably very time consuming (?!). There would definitely be value in generic routines for data living on boundaries (with MPI-parallel data structures) and lifting operations to optimally project volumetric data to faces. I agree with @richardschu in the sense that I currently do not see a major bottleneck in global vectors used for postprocessing purposes in ExaDG. For past developments, I personally ranked modularity higher than optimal performance or memory efficiency. So the essential question for me would right now be whether the Postprocessor in deal.II can serve the requirements in terms of modularity and the different use cases that we have in ExaDG (edit: I probably use modularity here in a different meaning than @kronbichler above). Our hand-written functionality might have been some optimum in terms of "overall efficiency" (in the sense of reaching certain research goals) at least in the past, which maybe explains why the code is currently designed as it is.

[nfehn]

To keep the discussion focussed, let me raise the question whether the present issue is the right place for reduced dimensional output (e.g. data living on boundaries like heat fluxes or wall-shear stresses). What was your original intention @kronbichler regarding this aspect? We could potentially open another issue if reasonable.

[richardschu]

Wall-shear stress calculations seem to be similar in terms of computing "normal gradients", so this could be relevant for @richardschu as well (is something like this realized in #349?).

The time averaged wall shear stress is averaged in time, but not in space, so it is not integrated over some ID.

Mentioning the surface output etc. here was just to keep it in mind (and that specific question was raised). Assuming now, we are just talking about postprocessing and not using the output anymore, i.e., just storing the data to look at it in awe. In that sense, I do not think that output generation currently is i) generated that often (user dependent), and hence ii) does not contribute overwhelmingly to the computing time (up to 90% if one is being ridiculous with the output, but otherwise maybe 20%-40%).

[nfehn]

Maybe this was a misunderstanding or unclear description from my side. I did not mean integrating over the boundary, but how to obtain the flux across the boundary. Just taking the gradient of the FE solution might be sub-optimal and there might be techniques with super-convergence property. Integrating that quantity over the boundary would then be a second, optional step.

[richardschu]

Maybe this was a misunderstanding or unclear description from my side. I did not mean integrating over the boundary, but how to obtain the flux across the boundary. Just taking the gradient of the FE solution might be sub-optimal and there might be techniques with super-convergence property. Integrating that quantity over the boundary would then be a second, optional step.

yep that is true, also for my case actually. another reason for introducing a projection onto a surface. thanks for the heads up. lets stay focused on the current issue though, this goes a bit off-topic :)