Discrete Galekin with periodic boundary condition

[Jaekwang Kim]

I am trying to solve advection equation using DG method with Meshworker. 

Yet, I have a problem in implementing periodic boundary condition. 

For example, I consider a 1D advection equation with periodic boundary condition....

with source term f=-sin(x), I assume that I should get cos(x) as a steady state solution, and I have confirmed this set up for the toy problem using Finite Difference.  

I derived my weak form as follow...

I am assembling equation (11) every time step,

Mass matrix for n+1 step solution 

 local_matrix(i,j) += fe_v.shape_value(i,point) *

                      fe_v.shape_value(j,point) *

                      JxW[point];

and right hand side is 

(cell integration) 

local_vector(i) += (fe_v.shape_value(i,point)

                    *(sol_phi[point] + dt * source_values[point])

                    +

                    dt*(advection_values[point]

                    *fe_v.shape_grad (i,point))

                    *sol_phi[point]

                    ) *JxW[point];

(Face integration) 

     for (unsigned int i=0; i<fe_v.dofs_per_cell; ++i)

            {

                if(advection_dot_normal >0) //outflux

                {    local_vector(i) -= dt * advection_dot_normal *

                                        (fe_v.shape_value(i,point)-fe_v_neighbor.shape_value(i,point)) *

                                        my_sol_phi[point] *

                                        JxW[point];

                }

                else //influx

                {

                    local_vector(i) -= dt * advection_dot_normal *

                                       (fe_v.shape_value(i,point)-fe_v_neighbor.shape_value(i,point)) *

                                       neighbor_sol_phi[point] *

                                       JxW[point];

                }

            }

(Boundary Integration)  <-think this is may be the problem.... 

                for (unsigned int i=0; i<fe_v.dofs_per_cell; ++i)

                {                       //sign changed here

                    if(advection_dot_normal >0) //outflux

                    { local_vector(i)-= dt*advection_dot_normal *

                                        fe_v.shape_value(i,point) *

                                        sol_phi[point] *

                                        JxW[point];

                    }

                    else

                    {

                        //influx this is problematic...it should be connected to last end cell (but..how?)

                      local_vector(i)-= dt*advection_dot_normal *

                                        fe_v.shape_value(i,point) *

                                        sol_phi[point] * it should be connected to last end cell, but the boundary integrate function cannot call information from reference cell...  (but..how?)

                                        JxW[point];

                    }

                }

As I am solving this x in [0,4pi] I should see two period of cosine wave after enough time step, but it is not..

Is there anyone who have the idea to use Meshworker for 

1D advection equation with periodic boundary condition?

Thanks in advance 

Jaekwang 

[Praveen C]

You have to do some extra work to use periodic bc with MeshWorker. I have a Euler DG code here which can do that

https://github.com/cpraveen/dflo/tree/master/src_mpi

The periodic bc is setup in these files (not my code, authors are mentioned in the files)

DealiiExtensions.cc

DealiiExtensions.h

Then see how these functions are used in claw.cc

This code is somewhat old now and I have not run it for quite some time. So it may or may not compile/run with latest deal.II, you have to try it out.

I have a simpler DG code for linear advection without MeshWorker here

https://bitbucket.org/cpraveen/deal_ii/src/master/dg/1d_scalar_legendre/

Best

praveen

On 08-Nov-2018, at 8:31 AM, Jaekwang Kim <jaekw...@gmail.com> wrote:

I am trying to solve advection equation using DG method with Meshworker. 

Yet, I have a problem in implementing periodic boundary condition. 

For example, I consider a 1D advection equation with periodic boundary condition....

with source term f=-sin(x), I assume that I should get cos(x) as a steady state solution, and I have confirmed this set up for the toy problem using Finite Difference.  

I derived my weak form as follow...

<Screen Shot 2018-11-07 at 6.15.31 PM.png>

<Screen Shot 2018-11-07 at 6.29.26 PM.png>

Is there anyone who have the idea to use Meshworker for 

1D advection equation with periodic boundary condition?

Thanks in advance 

Jaekwang 

--
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.
<CMakeLists.txt><dg_advection_periodic.cc><Screen Shot 2018-11-07 at 6.15.31 PM.png><Screen Shot 2018-11-07 at 6.29.26 PM.png>

[Daniel Arndt]

Jaekwang,

In general, using periodic boundary conditions with MeshWorker should be pretty easy.

All you would have to do is calling Triangulation.add_periodicity()

with a parameter constructed by GridTools::collect_periodic_faces().

Doing so will let MeshWorker treat all periodic faces as inner faces.

The 1D case is less well tested though. If you encounter additional difficulties feel free to ask.

Best,

Daniel

[Jaekwang Kim]

Thanks for the reply! 

I will take a look

Jaekwang 

[Jaekwang Kim]

Thanks, for the reply..

Hi am using following function to construct periodicity when I set up dof 

DoFTools::make_hanging_node_constraints (dof_handler, constraints);

DoFTools::make_periodicity_constraints(dof_handler,6,5 //Boundary ID 6 is outlet and 5 is inlet 

                                       ,0,

                                        constraints

                                        );

             

Question 1. 

I assume that this does same thing as Triangulation.add_periodicity () + Gridtools::collect_periodic_faces(). 

Isn't it? 

Question 2. 

What I don't know now is when I call 'integrate_boundary_term' function in the mesh worker 

how can I assess neighbor cell information?

For example, when I assemble for left end cell - the influx should be described by right-end cell at previous time solution solution at previous time.

To summary, I don't know how to call neighbor (described by periodic) cell solution within the 'integrate_boundary term' function

Thanks !! 

Jaekwang 

   template <int dim>

    void AdvectionProblem<dim>::integrate_boundary_term (DoFInfo &dinfo,

                                                         CellInfo &info)

    { ....

        FullMatrix<double> &local_matrix = dinfo.matrix(0).matrix;

        Vector<double> &local_vector = dinfo.vector(0).block(0);

        

        const FEValuesBase<dim> &fe_v = info.fe_values();

        const std::vector<double> &JxW = fe_v.get_JxW_values ();

        const std::vector<Tensor<1,dim> > &normals = fe_v.get_all_normal_vectors ();

        

        const AdvectionField<dim> advection_field;

        std::vector<Tensor<1,dim> > advection_values (fe_v.n_quadrature_points);

        advection_field.value_list (fe_v.get_quadrature_points(),advection_values);

....

[Daniel Arndt]

Jaekwang,

DoFTools::make_hanging_node_constraints (dof_handler, constraints);

DoFTools::make_periodicity_constraints(dof_handler,6,5 //Boundary ID 6 is outlet and 5 is inlet 

                                       ,0,

                                        constraints

                                        );

             

Question 1. 

I assume that this does same thing as Triangulation.add_periodicity () + Gridtools::collect_periodic_faces(). 

Isn't it? 

No, in general it isn't. This identifies degrees of freedom across periodic boundaries by writing the corresponding constraints in a ConstraintMatrix (or AffineConstraints) object.

However, for discontinuous approaches you normally don't want to enforce continuity across periodic boundaries strongly, but in a weak sense.

Question 2. 

What I don't know now is when I call 'integrate_boundary_term' function in the mesh worker 

how can I assess neighbor cell information?

For example, when I assemble for left end cell - the influx should be described by right-end cell at previous time solution solution at previous time.

To summary, I don't know how to call neighbor (described by periodic) cell solution within the 'integrate_boundary term' function

If you call add_periodicity() on a Triangulation object, the periodic faces are treated as internal faces in MeshWorker.

This means that you will not access them in a "integrate_boundary_term" function but in a "integrate_face_term" function.

Apart from that, you can always ask a cell iterator for its periodic neighbor on a given face via CellAccessor::periodic_neighbor().

Best,

Daniel

[Jaekwang Kim]

Thanks! Now it becomes clear to me. 

By the way, when I 'add_periodicity', I have compiler error on the red line of the code. 

Can you give me a clue to resolve this error? 

  //original argument

        std::vector<dealii::GridTools::

        PeriodicFacePair<typename DoFHandler<dim>::cell_iterator>> pairs;

    

        GridTools::collect_periodic_faces(dof_handler,

                                          /*b_id1*/ 6,

                                          /*b_id2*/ 5,

                                          /*direction*/ 0,

                                          pairs);

        

        triangulation.add_periodicity(pairs);

Error messages says that 

no viable conversion from

      'vector<PeriodicFacePair<TriaIterator<dealii::internal::DoFHandler::Iterators<dealii::DoFHandler<2,

      2>, false>::CellAccessor>>>' to 'const

      vector<PeriodicFacePair<TriaIterator<CellAccessor<2, 2>>>>'

        triangulation.add_periodicity(pairs);

                                      ^~~~~

also says that 

      candidate constructor not viable: no known conversion from

      'std::vector<dealii::GridTools::PeriodicFacePair<typename DoFHandler<2>::cell_iterator>

      >' (aka 'vector<PeriodicFacePair<TriaIterator<DoFCellAccessor<DoFHandler<2, 2>, false> >

      > >') to 'const

      std::__1::vector<dealii::GridTools::PeriodicFacePair<dealii::TriaIterator<dealii::CellAccessor<2,

      2> > >,

      std::__1::allocator<dealii::GridTools::PeriodicFacePair<dealii::TriaIterator<dealii::CellAccessor<2,

      2> > > > > &' for 1st argument

    vector(const vector& __x);

    ^

[Jaekwang Kim]

It seems that I don't need to add_periodicity to the triangulation unless it is parallel distributed triangulation? 

[Daniel Arndt]

Jaekwang,

It seems that I don't need to add_periodicity to the triangulation unless it is parallel distributed triangulation? 

I would be heavily surprised if you get correct results without calling Triangulation::add_periodicity() if you are using discontinuous elements and a serial Triangulation.

Otherwise, you are not providing any kind of information to these objects that you want to treat specific boundaries as periodic ones.

 

By the way, when I 'add_periodicity', I have compiler error on the red line of the code. 

Can you give me a clue to resolve this error? 

  //original argument

        std::vector<dealii::GridTools::

        PeriodicFacePair<typename DoFHandler<dim>::cell_iterator>> pairs;

    

        GridTools::collect_periodic_faces(dof_handler,

                                          /*b_id1*/ 6,

                                          /*b_id2*/ 5,

                                          /*direction*/ 0,

                                          pairs);

        

        triangulation.add_periodicity(pairs);

You need to call GridTools::collect_periodic_faces with a parameter of type

std::vector<dealii::GridTools::PeriodicFacePair<typename Triangulation<dim>::cell_iterator>> pairs;

if you want to use the resulting object for Triangulation::add_periodicity().

 

Best,

Daniel