Implement interior penalty discontinuous Galerkin method for temperature field.

[spco]

Hi guys,

This is my implementation of IPDG for the temperature field, as discussed a couple of months ago. There are several TODOs around, which may need some tweaking, but I'd rather show you what I have for now so you can feedback.

A brief overview:
parameters.use_discontinuous_temperature_discretization is the control variable (defaults to false).
When parameters.use_discontinuous_temperature_discretization==true:

(1) Set temperature to use FE_DGQ instead of FE_Q in setup_fes().
(2) DoFTools::make_flux_sparsity_pattern is called in setup_system_preconditioner() and setup_system_matrix().
(3) Dirichlet-type boundary constraints are skipped in compute_current_constraints() - we impose boundary conditions weakly
(4) Face terms are calculated via local_assemble_advection_face_terms() at the end of local_assemble_advection_terms().
(5) Face-term local matrices are added into system_matrix.
(6) compute_viscosity() returns zero (no stabilisation desired)

Changes I’ve had to make: compute_material_model_input_values() now takes an FEValuesBase object as argument to allow passing of FEFaceValues and FESubfaceValues

Scratch::AdvectionSystem now has a number of members to hold data from the faces, and its constructor requires a quadrature object for the faces.

CopyData::AdvectionSystem has vectors of local matrices to hold the internal-external, ext-int and ext-ext pairings (each element of the vector holds a matrix associated with each face or subface). A new member local_neighbor_dof_indices holds vectors of dof indices for each face (or subface) of the cell.

Any feedback gratefully received, particularly re whether these changes to AdvectionSystem are acceptable. I've tried to keep to the current style as far as I can - there's probably a lot of indenting tweaks needed, but I can do that later. I haven't done anything to attempt to optimize behaviour, e.g. minimising the number of local int_ext etc matrices by pre-calculating face/subface counts or various other parts as I don't know what returns we'd get.

I haven't included any DG parameter files - should these be included, eg convection-box-dg.prm?

Sam Cox

[gassmoeller]

Hi Sam,
great stuff. Thanks a lot for your contribution. I can take a closer look at your code later this week, but let me start by giving you a few things you can look into immediately:

• We usually handle all the indentation stuff with a program called astyle, preferably something around version 2.04. Many linux distributions have it as a package, or you can compile it yourself, it is pretty small and self-contained. If you have it installed, indenting all of ASPECT's source code and include files is as simple as typing ./doc/indent from the main directory of ASPECT. This will avoid the nasty red X next to your pull request.
• Providing some example models is nice, but what would be even better would to add some tests that prove your functionality works and acts as a reference for us to make sure nothing breaks it in the future. You can add tests by adding some (small model) prm files in the tests folder and adding a folder with the same name that contains some screen output and maybe a statistics file. This will include these models in our automated test suite. Of course the best option would be to have both: examples and tests 😄.
• A large and interesting feature like this calls for a section in the manual ;-), we can discuss later where and how we best include that, but if you have already an idea, feel free to go ahead. The manual would also be the place to explain the examples from the point above.
• A general question: From what I have seen so far you only implement the DG method for the temperature field. We strive to keep the handling of the temperature and compositional field as similar as possible to avoid a lot of if-loops in the advection system and to not introduce unexpected differences between the two. How much work would it be to enable DG for all advection fields? If very little: Would you be willing to include that as well? If medium to a lot: Could you leave a TODO or some explanation what would be necessary?

[spco]

Hi Rene, thanks for your feedback! I've corrected the formatting above.

On the question of DG for composition fields, are you suggesting all advection fields should use the same discretization, or that we should allow the user to pick and choose which fields are (dis)continuous?

[tjhei]

How do you handle assembly of face terms for faces on interfaces between processors (especially in combination with hanging nodes)?

[spco]

The layer of ghost cells provides the information on neighbouring cells owned by other processors. The issues that arise from the face-selection logic used in step-30 regarding faces on interfaces between processors are handled by modifying the face-selection logic for the case of equally-refined cells to rely on cell->id() rather than cell->level() and cell->index(). Have I missed another issue?

[bangerth]

On 11/04/2015 03:47 AM, Sam Cox wrote:

The layer of ghost cells provides the information on neighbouring cells owned
by other processors. The issues that arise from the face-selection logic used
in step-30 regarding faces on interfaces between processors are handled by
modifying the face-selection logic for the case of equally-refined cells to
rely on cell->id() rather than cell->level() and cell->index().

You can do it that way; or you could say that on faces between a locally owned
and a ghost cell, the processor with the lower subdomain_id/MPI rank is
responsible.

[spco]

@bangerth : the documentation stated id() takes O(level) time to compute - I guess subdomain_id() is known anyway - would it be worth changing to subdomain_id() then?

[bangerth]

Yes, that's likely more efficient. In the end you just need some kind of tie breaker; what exactly it is you do is not important as long as both processors agree on the procedure, so doing something cheap is better I suppose.

[tjhei]

I am asking because the logic is slightly more complicated with adaptive refinement, where you typically have to assemble coming from the finer cells. Using id() is certainly not a great idea but works in contrast to index(). It might be easier to use MeshWorker for the face assembly (or at least take a look at how the iteration is done based on the settings inside LoopControl).

[spco]

I think the current algorithm works, under a couple of assumptions:

• isotropic refinement (anisotropic would break it I'm sure, but is not currently implemented in ASPECT that I can find)
• that cell comes from a FilteredIterator<typename DoFHandler<dim>::active_cell_iterator> over IteratorFilters::LocallyOwnedCell() (as currently happens in assemble_advection_system()

A rough proof:

• loop over all active, locally-owned cells K
• loop over every face f of K
• if f is at domain boundary, then we’re happy (no hanging nodes or processor boundary)
• else
  • L is the neighbor() behind K and we have 6 cases:
  • (1) If neighbor L is coarser and locally-owned
    • we’re happy (due to mesh regularity, L must be active, and face f is not on a processor boundary).
  • (2) If neighbor L is coarser and not locally-owned
    • due to mesh regularity, L must be active, so compute from the other side on processor 2, where we will be in case (6)
  • (3) If neighbor L is the same size, active and locally-owned
    • we’re happy (face f is not on a processor boundary).
  • (4) If neighbor L is the same size, active but not-locally owned.
    • Then because L is adjacent to f, L is a ghost cell, and we can find its subdomain_id(). So compare these and compute from the side with lower subdomain_id (no hanging nodes).
  • (5) If neighbor L is the same size, inactive but all its children that touch f are locally owned
    • we’re happy (each subface is not on a processor boundary).
  • (6) L is the same size, inactive and some of its children that touch f are not locally owned.
    • For each cell A (child of L) lying adjacent to face f that is not locally-owned, A is a ghost cell so the subface between A and L is computable. All other children B lying adjacent to f are locally owned and the subface is not on a processor boundary, so treat as in (5)

Essentially, in any case where the neighbor (or neighbor_child as necessary) is not locally owned, it is a ghost cell so the face term is computable.

From what I can see, all the nasty cases being caught in LoopControl don't exist under the above assumptions, as it's been written for more general iterators and with more control over which faces are of interest.

[bangerth]

On 11/04/2015 10:56 AM, Timo Heister wrote:

I am asking because the logic is slightly more complicated with adaptive
refinement, where you typically have to assemble coming from the finer cells.

The tie breaker is only necessary between cells of the same level.

[spco]

So, my thoughts on what might need cleaning up in this code:

(1) Near-identical code in the equal-face and subface cases: Should this be cleaned up by a call to a function

void local_assemble_advection_individual_face(const AdvectionField     &advection_field,
                                              const typename DoFHandler<dim>::active_cell_iterator &cell,
                                              internal::Assembly::Scratch::AdvectionSystem<dim> &scratch,
                                              internal::Assembly::CopyData::AdvectionSystem<dim> &data,
                                              const FEFaceValuesBase<dim> face_finite_element_values)

(2) Removal of global_field_range, global_max_velocity, global_entropy_variation as input arguments to local_assemble_advection_face_terms() - these aren't used in the function.

(3) On the issue of composition fields - there shouldn’t be too much of an issue, although the only terms needed from the faces in the composition field (due to lack of diffusion) are the convection-dependent terms, so there will be a good number of (advection_field.is_temperature ? . : . ) calls. Also, allowing only some fields to be dG while others staying continuous would require some fix to the parameters, but I don’t think that would make a huge amount of sense - presumably, all advection fields are continuous or discontinuous, so we’d only want one parameter.

I am happy to implement the composition fields and will work on that in the next day or so when I get a chance.

[gassmoeller]

Hi Sam,
sorry for let you hanging so long. I do not know if people want to use some continuous and some discontinuous fields, I just know that I used fields for the weirdest different things, so I would like to keep the possibility. What about the following suggestion:
You could add another function to the struct AdvectionField in simulator.h that is called is_discontinuous() and decides whether any field (temperature or composition) is handled discontinuously. And we add a parameter to parameters.h (something like List of discontinuous fields). The parameter could be a multiple selection of all compositional fields (Patterns::MultipleSelection(pattern_of_names) Then is_discontinuous() could check if it is the temperature field and discontinuous (like you do right now) and if so return true, if not check if it is a compositional field, and the field is in this list and if so return true, and if none is true return false. This would keep the assembly more or less clean. Of course in assemble_face_terms you would still need to ask for is_temperature() to decide if you need particular terms or not.

[bangerth]

My initial inclination would have been to either make all fields continuous or all discontinuous. But @gassmoeller is right -- people have (ab)used the flexibility of compositional fields in ways we could never have imagined, and keeping things flexible would allow to continue this tradition. If it isn't too much work, it might be interesting to leave it up to the user whether a field is cG or dG.

[spco]

Hi guys,

sorry for leaving this hanging for such a long time - I had to go get my teeth into some theory for a couple of months!

I've fixed a couple of bugs in the original implementation, so that compositions can now be treated as continuous or discontinuous. However, I've hit an issue in the implementation of the flexible approach (allowing different fields to have cG or dG). The issue I have is that if we allow two composition fields, one cG and one dG, then this starts to break a number of assumptions in the code. For example, Introspection has a struct BaseElements, enumerating the FEs. If we have cG and dG compositional fields, they no longer share the same base element, so we need an unsigned int continuous_compositional_fields and unsigned int discontinuous_compositional_fields to replace the unsigned int compositional_fields. This is then going to break various other places: for example, in mesh_refinement/composition.cc, line 44 calls this->get_fe().base_element(this->introspection().base_elements.compositional_fields . This is going to need tweaking to check whether the compositional field is continuous or discontinuous.

Similarly, in introspection.cc, in setup_fes(), the two will need treating differently, and this will affect the structure of multiplicities and fes.

I'm asking for advice - is this the right way to go, or is there a better way?

[gassmoeller]

Thanks for coming back to this. You are right, the individual selection seems to cause more problems than I expected. Since the individual selection is a feature that is not needed right away, would it be much work to implement your new additions in a way that allows for flexible selection of discontinuous fields, but enforce either all fields being selected as discontinuous or none, for now?
That way we could leave the rest of the code untouched, but in case somebody wants to use mixed continuous/discontinuous elements he would only need to touch the other parts of the code, not your additions.
If that is also causing much trouble, then I would say let us just go with the all continuous / discontinuous fields for now and lets worry about the mixed case when somebody actually needs it.
Let us know if there are problems or questions.

[bangerth]

If that is also causing much trouble, then I would say let us just go with the all continuous / discontinuous fields for now and lets worry about the mixed case when somebody actually needs it.

I think that makes sense.

[spco]

Ok, I think I have a (nearly) finished product! I haven't yet added any test cases. To do. Any pointers on that front would be welcomed.

Let me know what needs tidying up - otherwise, I think this is about ready to go, from my POV.

I also haven't added anything to the manual yet - from a quick run through, I think it might be worth mentioning this in §2.9 Numerical methods: Accurate discretizations, as well as §5.1.2: Parameters for certain aspects of the numerical algorithm. The two parameters use_discontinuous_temperature_discretization and use_discontinuous_composition_discretization would need adding to §5.29. I don't know whether a more substantial section is either wanted or needed? There doesn't seem to be a parallel one for the continuous case, for example. It might be worth stating explicitly somewhere that, if temperature is discontinuous, then playing around with stabilization parameters will be pointless. Just a thought.

[gassmoeller]

/run-tests

[gassmoeller]

The overall structure looks good, it should also be compatible to the new structure of the advection assembly (as a follow-up to #718 for the Stokes system. @bangerth: This patch needs face terms similar to #722 but for the Advection equation).
I think you need to address a compile issue with deal.II 8.2 though, and I suspect it is because of the name change of boundary_indicator to boundary_id. For a possible fix look at assembly.cc:1237.

I do not feel competent enough to review the math of the patch, but I guess we need a benchmark and tests for that anyway.
Concerning the tests: It might be easiest to copy a few tests for the advection system and just make sure the discontinuous formulation gives (nearly) identical results to the continuous one for now. Tests that come to mind could be: diffusion.prm, time_dependent_temperature_bc_2.prm, composition_stabilization.prm.
A interesting benchmark case that is popular in the geosciences would be a variation of the test van_keken_smooth_tracer.prm with discontinuous compositional and temperature fields. It is not a perfect benchmark, but when comparing the fields at an end time of 2000 and with say global refinement 6, then it would be interesting to see how much the artificial diffusion differs to the continuous formulation.

I will leave a few minor comments in the code, but I think it might be worth merging this soon (when it compiles) to not block the restructuring of the assembly by @bangerth. Essentially so far nothing old has been changed, only the new functions are missing tests and documentation, so nobody will use it yet 😄

[gassmoeller]

I think this change will be merged into the developer version with #723 soon anyway, will be removed during the rebase before merging.

[bangerth]

The issue with the quadrature is in fact not that it's not possible to create a quadrature rule with no points, but that it's not possible to create a FEValuesBase object with such a quadrature formula:

--------------------------------------------------------
An error occurred in line <2164> of file </root/deal.II/dealii/source/fe/fe_values.cc> in function
    dealii::FEValuesBase<dim, spacedim>::FEValuesBase(unsigned int, unsigned int, dealii::UpdateFlags, const dealii::Mapping<dim, spacedim>&, const dealii::FiniteElement<dim, spacedim>&) [with int dim = 2; int spacedim = 2]
The violated condition was: 
    n_q_points > 0
The name and call sequence of the exception was:
    ExcMessage ("There is nothing useful you can do with an FEValues " "object when using a quadrature formula with zero " "quadrature points!")
Additional Information: 
There is nothing useful you can do with an FEValues object when using a quadrature formula with zero quadrature points!

I think that's actually quite reasonable.

What do you think of the following proposal: While the scratch objects have member variables for the FEValues object, they could have a std_cxx11::shared_ptr<FEFaceValues> member for the face integration. If the face quadrature formula passed down has no quadrature points, then that pointer would simply be initialized with NULL. Otherwise, you would create an object (via operator new) initialized with an appropriate FEFaceValues object. This way you're not just creating an invalid FEFaceValues object, you're in fact not creating one at all.

This should avoid the error above.

[tjhei]

may we require 8.4 once that is out?

I would be happy to require 8.4.0 after the next ASPECT release but I would prefer the next ASPECT release (soon!?) to be at 8.3 or 8.2.1.

I think using a shared_ptr for some of the objects that are unused is cleaner than constructing an invalid FEFaceValues.

[spco]

I've implemented the shared_ptr method. Let me know what you think! It all seems to work locally - just waiting for the tests to run here.

Re the errors for 8.2.0: It also didn't work in 8.4.0 either, I just hadn't tested it enough locally :)

[bangerth]

OK, great. Let me look through it one more time. In the meantime, would you mind squashing the 27 commits into one?

[bangerth]

Please break these lines (starting at line 368) differently so that an expression in either of the branches does not unnecessarily crosses lines. Maybe like

  field_is_discontinuous
  ?
  ...
  :
  0

?

[bangerth]

I had only two really minor issues left. Please address them, squash everything, and then this is ready to merge!

I really appreciate all of the work you've put into this, and your patience with both my slowness in getting this reviewed and with addressing all of the comments I've had. Many thanks, @spco !

[tjhei]

can you add a comment what these matrices are? What is the length of the vector?

[tjhei]

Sorry for only looking at it now. Looks good except my minor comments above.

[spco]

Thanks for your comments @bangerth and @tjhei . I've addressed all the issues apart from that of averaged material_model_outputs. On that front, there are face_material_model_outputs and neighbor_face_material_model_outputs. Most terms then use whichever is appropriate. Some instead (eg penalty term and density_c_P_and_latent_heat) use the maximum over the 2 cells, as that is more appropriate. Hopefully that makes sense? I will squash these once the tests finish.

[tjhei]

Yes, taking the maximum makes sense. 👍

[spco]

Right, all tests have finally passed!

[bangerth]

I've hit the "merge" button so we can get on with the other patches that were queued up after this.

@spco -- would you mind adding an entry to doc/modules/changes.h regarding the DG discretization that describes the functionality you implemented in one paragraph (as a follow-up patch)?

@tjhei -- what was your take on averaging? We can do this in a follow-up patch without having to hold this further up.

@jdannberg and @gassmoeller -- if you want, go ahead and also put the advection assemblers into the signals format.

[tjhei]

what was your take on averaging?

The correct thing to do would be to average, even if the assembly here takes the maximum of the values in a single quadrature point.

@jdannberg and @gassmoeller -- if you want, go ahead and also put the advection assemblers into the signals format.

Can you take a look at #751 before we continue with modifying the other assemblers?

[bangerth]

Some related discussion is in #808.