step-33 with DG?

[Giorgio Giangaspero]

Hello everyone,

I would like to solve step-33 using DG only, instead of mixed cG and DG discretization.

Attached you will find my modified version of step-33, which of course diverges quickly, thus this post.

Here are the main modifications to the original code:

1) we use DG elements and a generic mapping, i.e.

#include <deal.II/fe/mapping.h>

#include <deal.II/fe/fe_dgq.h>

and in the constructor of the ConservationLaw class

fe (FE_DGQ<dim>(degree), EulerEquations<dim>::n_components)

2) the function setup_system uses the function make_flux_sparsity_pattern (as in step-12) and reads

template <int dim>

  void ConservationLaw<dim>::setup_system ()

  {

    DynamicSparsityPattern dsp(dof_handler.n_dofs());

    DoFTools::make_flux_sparsity_pattern (dof_handler, dsp);

    sparsity_pattern.copy_from(dsp);

    system_matrix.reinit (sparsity_pattern);

  }

3) after line 1602 of the original file, I have naively added (I am new to DG) the following lines which should take into account 

   the fact that in DG the solution is discontinuous also across faces shared by cells at the same refinement level.

               else

                {

                  // This is 'normal' face, that is a face between two cells

                  // at the same refinement level. This term is needed for

                  // pure DG.

                  const typename DoFHandler<dim>::cell_iterator

                  neighbor = cell->neighbor(face_no);

                  Assert(neighbor->level() == cell->level(),

                         ExcInternalError());

                  const unsigned int neighbor_face_no=

                    cell->neighbor_of_neighbor(face_no);

                  neighbor->get_dof_indices (dof_indices_neighbor);

                  fe_v_face.reinit (cell, face_no);

                  fe_v_face_neighbor.reinit (neighbor, neighbor_face_no);

                  assemble_face_term (face_no, fe_v_face,

                                      fe_v_face_neighbor,

                                      dof_indices,

                                      dof_indices_neighbor,

                                      false,

                                      numbers::invalid_unsigned_int,

                                      cell->face(face_no)->diameter());

                }

4) I have commented out the diffusion stabilization term proportional to h^(eta) which, as far as I know, is not needed for DG.

To keep things as easy as possible I have not changed anything else. Also, the input.prm file is the same.

I am aware that for DG we could use the MeshWorker class but I prefer the current implementation because I find it clearer and I would like to retain the possibility to switch to the mixed cG/DG implementation easily. 

So the big questions: what am I missing? Is the discretization correct? if so, which solver should I use?

Thanks a lot!

best

giorgio

 

[Praveen C]

Hi

If you have discontinuous solutions, you will need some stabilization even for DG (for degree > 0). The initial condition in step-33 is discontinuous. If you dont want to use stabilization, first try some problem with a smooth solution. Just put a smooth initial density distribution.

Another thing to try is degree=0 which is a finite volume scheme. No stabilization needed here and use theta=0 which is a forward euler scheme. But the code is setup to solve an implicit system, so I dont know if this will work.

I have a pure DG code starting from step-33 and using MeshWorker here

https://github.com/cpraveen/dflo

It is based on Runge-Kutta scheme and limiters. Code for the theta scheme as in step-33 is also there but it may or may not work.

Best

praveen

--
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see https://groups.google.com/d/forum/dealii?hl=en
---
You received this message because you are subscribed to the Google Groups "deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email to dealii+un...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[Giorgio Giangaspero]

Hi praveen,

thank you very much for you answer.

Indeed if I use elements with degree 0 the code runs smoothly. And if I set an initial smooth density the code works for higher-order elements as well. I am using the implicit time integration with theta=0.5. The only thing I have changed is the numerical flux function, I am now using the Roe flux, the implementation of which I copied and pasted from your code.

Thanks again,

best

giorgio