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

[Wolfgang Bangerth]

Sam,

> 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.

Yes, that's almost certainly less work.

As a suggestion, work from the ASPECT github repository and try to keep 
your branch relatively up to date with regard to the main development 
line. This is going to make all discussions and modifications easier, 
and make merging easier when the code is ready.

If you want, create a pull request early on so that you can share (and 
we can see) your code. This allows for early feedback on design issues, 
and avoids you running into unnecessary hurdles later on.


> 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.

Yes, I think you should model this on the
   "Use locally conservative discretization"
flag that determines the element used for the pressure variable.


> (1) push_back an FE_DGQ object in setup_fes

Yes. Use the same code as for the pressure element.


> (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.

I think the two masks are simply there to indicate which vector 
components couple with which others. I suspect that, without knowing for 
sure, the arguments corresponding to the coupling matrix in the non-DG 
function:

http://dealii.org/developer/doxygen/deal.II/group__constraints.html#ga250c261f4dba6c77834b2fa8d3bf7f5c

We should probably update the documentation of 
make_flux_sparsity_pattern() :-(

My initial impulse would be to simply add the call to make_f_s_p() to 
the existing call to make_s_p() and see what happens. We can optimize 
this later.


> (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.

Yes, though they already have access to it. See how the code that 
assembles the Stokes matrix deals with Neumann (traction) type boundary 
conditions for the velocity, for example.


> (5) In assembly, change the bilinear form substantially, and add
> extra functions to compute the face term contributions to both the
> matrix and rhs.

Yes. My suggestion would be to look through 
local_assemble_advection_system() in assembly.cc. I would create a 
second function that is responsible only for the interface terms, and if 
use_discontinuous_temperature is set, call the other function from the 
end of local_assemble_advection_system(). It will make everyone's job 
easier during code review and merging if you manage to keep the code for 
interface assembly as separate from the existing code as possible.


> (6) Inclusion of a penalty parameter, either hard-coded or as an
> extra parameter from file.

Yes. If there's a good value for it, hardcode it for now. We can always 
make it into a parameter later on.


> 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.

All of the mesh refinement classes already have access to the material 
model (and all other model components) if you derive your class from 
SimulatorAccess. Take a look at the section in the manual that describes 
SimulatorAccess.


> 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?

None of this will be necessary if you use SimulatorAccess.


> Finally, would this be of interest to add to the code base?

Yes. It's a sizable enough functionality that we'll probably have to go 
through a couple or three rounds of review before merging, but you will 
see that while that may be a bit of a hassle, it also makes your code 
much better!

Best
   W.



-- 
------------------------------------------------------------------------
Wolfgang Bangerth               email:            bangerth at math.tamu.edu
                                 www: http://www.math.tamu.edu/~bangerth/