[aspect-devel] Advice/discussion on implementing discontinuous temperature field in ASPECT

[Cox, Samuel P.]

Hi Wolfgang and Timo (and others),

I'm looking into adding a few features to ASPECT, and would appreciate your thoughts on my plans before I dive in. 

I have previously converted deal.II’s step-32 to use discontinuous (FE_DGQ) elements for the temperature field, using interior-penalty dG for the advection equation, and have developed a residual-based error estimator for the problem.

My plan is to re-implement the IP-DG and error estimator in ASPECT, rather than attempt to re-implement the functionality of ASPECT in my own code. As such, I'd like to ask for your advice on the feasibility of this, as well as the desirability of adding this to the ASPECT codebase.

In a little detail:

Discontinuous elements for temperature would, roughly, require some alteration to simulator/core.cc, simulator/assembly.cc and simulator/introspection.cc plus parameters.cc: The simplest method I can think of would be to have a boolean parameter use_discontinuous_temperature from the input file, and use this to handle cases where we would need a different setup. The following are my best shot at where changes would need to be made - I've probably missed a few:

(1) push_back an FE_DGQ object in setup_fes
(2) call DoFTools::make_flux_sparsity_pattern in Simulator::setup_system_matrix. This would probably need an improved method than simply calling make_flux_sparsity_method, as this would put in flux terms for all the variables, not just temperature. My dG code from step-32 has kept separate DoFHandlers and matrices for the Stokes and temperature parts, so I haven't worked through all the intricacies of doing this for just part of a system - please correct me if I'm wrong! It seems make_flux_sparsity_pattern can probably handle this through the int_mask and flux_mask arguments.
(3) A similar solution for skipping the temperature part of make_hanging_node_constraints in setup_dofs.
(4) Weak imposition of temperature boundary conditions means skipping the interpolate_boundary_values in compute_current_constraints, and instead passing the boundary temperature object to the assembly routines.
(5) In assembly, change the bilinear form substantially, and add extra functions to compute the face term contributions to both the matrix and rhs.
(6) Inclusion of a penalty parameter, either hard-coded or as an extra parameter from file.
(7) Skipping of artificial viscosity calculations, if it's worth it.


The error estimator is based on the framework of KellyErrorEstimator, which has already been easily used in MeshRefinement plugins, so this can be standalone. The only difference that I think may be problematic is that this requires access to the parameters of the equations, eg density, viscosity. In my current code these are passed in as CoefficientFunctions (a class I've written that looks like Function but allows input of temperature or other values at the given points).
I'm not sure how best to pass this information in from ASPECT - would I want to pass in the MaterialModel, and evaluate the functions at the necessary quadrature points, or would I want to build a full-size vector of, say, density for the owned part of the domain, and pass this in along with my solution vector?


Finally, would this be of interest to add to the code base? In particular, the changes to the files in simulator/ would otherwise render my code incompatible with simply upgrading to any new versions of ASPECT, unless they become part of the main code. I'm happy to invest some time in this, but would like to see how the land lies before diving in! Any thoughts on this, or corrections to my understanding of what would be involved, are most welcome!

Sam Cox
University of Leicester