Implement periodic boundary condition on distributed mesh

[Zhicheng Wang]

Hi,
I'm doing large eddy simulation of a fully developed channel flow. I encountered a problem of imposing periodic boundary condition on faces which are distributed to different processors. I know how to enforce periodic condition by adding dofs constraints to A ConstraintMatrix like step-45. However, the problem is the global dofs are distributed to several cores, and usually the two sides of a periodic boundary are distributed to different cores, so how to setup the CONSTRAINTs between dofs located at different cores?

Thanks

Zhicheng

[Markus Bürg]

Hello Zhicheng,

there is also the function DoFTools::make_periodicity_constraints.
According to the documentation that one should be working for
distributed grids, too.

Best,
Markus

[Timo Heister]

Hi Markus,

> there is also the function DoFTools::make_periodicity_constraints. According
> to the documentation that one should be working for distributed grids, too.

where did you find the comment that it should work for distributed
grids? The code recursively goes over the faces and subfaces and I can
not see how this can work for a distributed Triangulation.

It is a difficult problem to handle.


--
Timo Heister
http://www.math.tamu.edu/~heister/

[Markus Bürg]

Hello Timo,

> where did you find the comment that it should work for distributed
> grids? The code recursively goes over the faces and subfaces and I can
> not see how this can work for a distributed Triangulation.

There is no comment statiting that it does not work. Thus, I supposed it
should.

Best,
Markus

[Timo Heister]

> There is no comment statiting that it does not work. Thus, I supposed it
> should.

Now it says that it won't work.

@Zhicheng:
I fear you have to figure out how to implement this. Right now I do
not see a simpler version than:
1. find all dof in own cells that correspond to the periodicity and
store them with their global position.
2. dublicate this information using MPI to each machine.
3. resolve the constraints using this information with the owned cells
on each machine by looking up in that datastructure (somewhat similar
to the implementation of make_periodicity_constraints()).

[Zhicheng Wang]

thank you all of you,

   i'll try Timo's method.

   as i in some fem codes , a pair of faces that are to be imposed periodic boundary condition on, are treated as two internal neighbors faces, then the mesh partitioning tool (Metis) could set the correspond element to be ghost cell. 

[Timo Heister]

> as i in some fem codes , a pair of faces that are to be imposed periodic
> boundary condition on, are treated as two internal neighbors faces, then the
> mesh partitioning tool (Metis) could set the correspond element to be ghost
> cell.

Metis? Are you using a parallel::distributed::Triangulation or a
regular Triangulation? Metis is only used in the latter case and there
the deal.II routine for periodic conditions should just work out of
the box.

[Zhicheng Wang]

No, I use parallel::distributed::Triangulation
I saw other fem code set faces of periodic as internal boundary, and use Metis partition the domain.

发自我的 iPhone

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

[Wolfgang Bangerth]

> No, I use parallel::distributed::Triangulation

Using this class together with METIS does not work.

> I saw other fem code set faces of periodic as internal boundary, and use Metis partition the domain.

But the FEM code has to be prepared to deal with this. deal.II can't
set periodic faces as internal. I'm afraid what you're trying to
do is not as easy in deal.II as you think. Timo already outlined
how to do what you want with parallel::distributed::Triangulation.

Best
W.

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

[Zhicheng Wang]

Hello, Tim,

     could you please give a fragment of the code or pseudo-code about periodic boundary ? I still have no idea how to implement it on a computer with shared memory. 

     

     BTW, does deal.ii@googlegroups have an another archive, since it is really hard to login google apps in China.

zhicheng

--
Zhicheng Wang
Ph.D. Candidate
RM 402, No.11 Beisihuan West Road
Institute of Engineering Thermophysics, Chinese Academy of Sciences
Beijing, 100190, P.R. China
Tel: (86) 10-82543092
E-mail: wangzhi...@gmail.com; wangzh...@mail.etp.ac.cn

[Timo Heister]

> could you please give a fragment of the code or pseudo-code about
> periodic boundary ? I still have no idea how to implement it on a computer
> with shared memory.

If you are not using a parallel::distributed::Triangulation, you can
use the method described in step-45:
http://www.dealii.org/developer/doxygen/deal.II/step_45.html

You might want to read this step and the implementation of
make_periodicity_constraints to get started.

> BTW, does deal.ii@googlegroups have an another archive, since it is
> really hard to login google apps in China.

No, we have no other archive at this time.

[Zhicheng Wang]

sorry, I still have a problem with applying periodic boundary constraint on distributed mesh.
bellow the email is an example that apply the constraints. I tested some cases with this piece of code, it works well on 2 dimensional, but when going to 3 dimensional cases,  it will get incorrect results, depends on the number of dofs: when the number of dofs is small, the result is right, when the number  of dofs is large, the results is wrong. 
It looks like the Constraint Matrix only knows 'locally_relevant_dofs',  if the dof I want to constraint is outside of this range,the ConstraintMatrix could not handle it.

Thanks in advance!    
*********************************************************************************************************************************************8

 template <int dim>
  void
  NavierStokesProjection<dim>::    setup_periodic_boundary_constraints ()
 {
        unsigned int myid = Utilities::MPI::this_mpi_process(mpi_communicator);
    const unsigned int n_processes = Utilities::System::
                                     get_n_mpi_processes(mpi_communicator);

    const QGauss<dim-1> quad_vel(deg+2);
    const QGauss<dim-1> quad_pres(deg+1);

    FEFaceValues<dim> feface_vel (fe_velocity,quad_vel,
                                              update_quadrature_points);
    FEFaceValues<dim> feface_pres (fe_pressure,quad_pres,
                                              update_quadrature_points);

    const unsigned int dpc_v = fe_velocity.dofs_per_cell;
    const unsigned int dpc_p = fe_pressure.dofs_per_cell;

        std::vector< double > global_dofs_positions_vel,global_dofs_positions_pres;
        std::vector< double > global_dofs_positions_vel_buff,global_dofs_positions_pres_buff;

        const unsigned int nfpc = GeometryInfo<dim>::faces_per_cell;
    const unsigned int nvpc = GeometryInfo<dim>::vertices_per_face;


    u_constraints.clear ();
    u_constraints.reinit(locally_relevant_dofs_velocity);
    DoFTools::make_hanging_node_constraints (dof_handler_velocity,
                       u_constraints);
//
    v_constraints.clear ();
    v_constraints.reinit(locally_relevant_dofs_velocity);
    DoFTools::make_hanging_node_constraints (dof_handler_velocity,
                       v_constraints);
//
    w_constraints.clear ();
    w_constraints.reinit(locally_relevant_dofs_velocity);
    DoFTools::make_hanging_node_constraints (dof_handler_velocity,
                       w_constraints);
//
//
    p_constraints.clear ();
    p_constraints.reinit(locally_relevant_dofs_pressure);
    DoFTools::make_hanging_node_constraints (dof_handler_pressure,
                       p_constraints);

     for(typename DoFHandler<dim>::active_cell_iterator
             cell = dof_handler_velocity.begin_active ();
             cell != dof_handler_velocity.end (); ++cell )
       if  (cell->is_locally_owned())
        for(unsigned int f=0; f<nfpc; ++f)
       {

          if( ( cell->at_boundary())
                 && (cell->face(f)->boundary_indicator() == 4 ) )
          {
                 typename DoFHandler<dim>::active_cell_iterator
                      pres_cell (&triangulation,
                                cell->level(),
                                cell->index(),
                                &dof_handler_pressure);

                 const typename DoFHandler<dim>::active_face_iterator
                        face_v = cell->face(f);
                 const typename DoFHandler<dim>::active_face_iterator
                        face_p = pres_cell->face(f);


                 const unsigned int dpfv
                    =face_v->get_dof_handler().get_fe().dofs_per_face;
                 const unsigned int dpfp
                    =face_p->get_dof_handler().get_fe().dofs_per_face;

                 std::vector<unsigned int> ldi_v (dpfv);
                 std::vector<unsigned int> ldi_p (dpfp);
                 face_v->get_dof_indices (ldi_v);
                 face_p->get_dof_indices (ldi_p);
                
                 feface_vel.reinit(cell,f);
                 feface_pres.reinit(pres_cell,f);

                 const std::vector< Point<dim> > points_v =
                     feface_vel.get_present_fe_values().get_quadrature_points();
                 const std::vector< Point<dim> > points_p =
                     feface_pres.get_present_fe_values().get_quadrature_points();


            for(unsigned int i=0; i<dpfv; ++i)
              {
                  global_dofs_positions_vel.push_back(ldi_v[i]);
                  global_dofs_positions_vel.push_back(points_v[i][0]);
                  global_dofs_positions_vel.push_back(points_v[i][1]);
              }

            for(unsigned int i=0; i<dpfp; ++i)
               {
                  global_dofs_positions_pres.push_back(ldi_p[i]);
                  global_dofs_positions_pres.push_back(points_p[i][0]);
                  global_dofs_positions_pres.push_back(points_p[i][1]);
               }
               
                              
           } //
         } //


                   MPI_Barrier(mpi_communicator);


                unsigned int len_v = global_dofs_positions_vel.size();
                const unsigned int max_len_v = Utilities::MPI::max(len_v,mpi_communicator);
                const unsigned int num_periodic_dofs_vel = Utilities::MPI::sum(len_v,mpi_communicator);

                unsigned int len_p = global_dofs_positions_pres.size();
                const unsigned int max_len_p = Utilities::MPI::max(len_p,mpi_communicator);
                const unsigned int num_periodic_dofs_pres = Utilities::MPI::sum(len_p,mpi_communicator);

                global_dofs_positions_vel_buff.resize(max_len_v*n_processes);
                global_dofs_positions_pres_buff.resize(max_len_p*n_processes);


                if(len_v<max_len_v)
                for(unsigned int k=0; k<(max_len_v-len_v); ++k)
                    global_dofs_positions_vel.push_back(-10e6);

                if(len_p<max_len_p)
                for(unsigned int k=0; k<(max_len_p-len_p); ++k)
                    global_dofs_positions_pres.push_back(-10e6);




//
//
                MPI_Allgather(&global_dofs_positions_vel[0], max_len_v, MPI_DOUBLE,
                          &global_dofs_positions_vel_buff[0], max_len_v, MPI_DOUBLE, mpi_communicator);


                MPI_Allgather(&global_dofs_positions_pres[0], max_len_p, MPI_DOUBLE,
                          &global_dofs_positions_pres_buff[0], max_len_p, MPI_DOUBLE, mpi_communicator);


                unsigned int buff_size_vel = global_dofs_positions_vel_buff.size();
                unsigned int buff_size_pres = global_dofs_positions_pres_buff.size();
               

                std::vector< double >velocity_dofs, pressure_dofs;

                for(unsigned int i=0; i<buff_size_vel; ++i)
                    if(global_dofs_positions_vel_buff[i]!=-10e6)
                        velocity_dofs.push_back(global_dofs_positions_vel_buff[i]);
                for(unsigned int i=0; i<buff_size_pres; ++i)
                    if(global_dofs_positions_pres_buff[i]!=-10e6)
                        pressure_dofs.push_back(global_dofs_positions_pres_buff[i]);

                unsigned int size_vel = velocity_dofs.size();
                unsigned int size_pres = pressure_dofs.size();


                     Assert(size_vel== num_periodic_dofs_vel,
                                ExcMessage("Wrong number of periodic dofs for velocity!"));
                     Assert(size_pres== num_periodic_dofs_pres,
                                ExcMessage("Wrong number of periodic dofs for pressure!"));

                     


     for(typename DoFHandler<dim>::active_cell_iterator
             cell = dof_handler_velocity.begin_active ();
             cell != dof_handler_velocity.end (); ++cell )
       if  (cell->is_locally_owned())
        for(unsigned int f=0; f<nfpc; ++f)
        {

          if( ( cell->at_boundary())
                 && (cell->face(f)->boundary_indicator() == 5 ) )
        {


             typename DoFHandler<dim>::active_cell_iterator
                      pres_cell (&triangulation,
                                cell->level(),
                                cell->index(),
                                &dof_handler_pressure);
                
                 const typename DoFHandler<dim>::active_face_iterator
                        face_v = cell->face(f);
                 const typename DoFHandler<dim>::active_face_iterator
                        face_p = pres_cell->face(f);

                 const unsigned int dpfv
                    =face_v->get_dof_handler().get_fe().dofs_per_face;
                 const unsigned int dpfp
                    =face_p->get_dof_handler().get_fe().dofs_per_face;

                 std::vector<unsigned int> ldi_v (dpfv),ldi_p(dpfp);
                 face_v->get_dof_indices (ldi_v);
                 face_p->get_dof_indices (ldi_p);

                 feface_vel.reinit(cell,f);
                 feface_pres.reinit(pres_cell,f);

                 const std::vector< Point<dim> > points_v =
                     feface_vel.get_present_fe_values().get_quadrature_points();
                 const std::vector< Point<dim> > points_p =
                     feface_pres.get_present_fe_values().get_quadrature_points();



                 unsigned int search_step = dim;

              for(unsigned int i=0; i<dpfv; ++i)
                for(unsigned int j=0; j<num_periodic_dofs_vel; j=j+search_step)
                   if( (std::fabs(points_v[i][0]-velocity_dofs[j+1])<1e-8) &&
                       (std::fabs(points_v[i][1]-velocity_dofs[j+2])<1e-8) )
//                    if(!u_constraints.is_inhomogeneously_constrained(ldi_v[i]))
                      {

                         u_constraints.add_line  (ldi_v[i]);
                         v_constraints.add_line  (ldi_v[i]);
                         w_constraints.add_line  (ldi_v[i]);

                         u_constraints.add_entry (ldi_v[i],
                                         (int)velocity_dofs[j],
                                     1.0);
                         v_constraints.add_entry (ldi_v[i],
                                         (int)velocity_dofs[j],
                                      1.0);
                         w_constraints.add_entry (ldi_v[i],
                                     (int)velocity_dofs[j],
                                      1.0);
                      }

              for(unsigned int i=0; i<dpfp; ++i)
                for(unsigned int j=0; j<num_periodic_dofs_pres; j=j+search_step)
                  if( (std::fabs(points_p[i][0]-pressure_dofs[j+1])<1e-8) &&
                      (std::fabs(points_p[i][1]-pressure_dofs[j+2])<1e-8) )
//                    if(!p_constraints.is_inhomogeneously_constrained(ldi_p[i]))
                       {
                         p_constraints.add_line  (ldi_p[i]);
                         p_constraints.add_entry (ldi_p[i],
                                         (int)pressure_dofs[j],
                                     1.0);
                       }

          } //

       }//


    u_constraints.close ();
    v_constraints.close ();
    w_constraints.close ();
    p_constraints.close ();
        

        }
 **************************************************************************************************************************************************************

On Friday, September 28, 2012 2:31:29 AM UTC+8, Timo Heister wrote:

[Timo Heister]

> bellow the email is an example that apply the constraints. I tested some
> cases with this piece of code, it works well on 2 dimensional, but when
> going to 3 dimensional cases, it will get incorrect results, depends on the
> number of dofs: when the number of dofs is small, the result is right, when
> the number of dofs is large, the results is wrong.
> It looks like the Constraint Matrix only knows 'locally_relevant_dofs', if
> the dof I want to constraint is outside of this range,the ConstraintMatrix
> could not handle it.

Yes, you are initializing the ContraintMatrix using
u_constraints.reinit(locally_relevant_dofs_velocity);
This tells it to only consider lines that are locally relevant. You
can remove the parameter (so just call reinit()) to test if this is
your problem. But, this is done for a reason. A ConstraintMatrix
without an IndexSet will consume a very large amount of memory
(proportional to the total number of DoFs in your computation instead
of proportional to the local number of DoFs).

[Zhicheng Wang]

Hi, Timo, Thank you for your reply, but constraint.reinit() could not solve the problem. Is it possible 

the order of face->get_dofs_indices() is not agree with that of FeFaceValue.get_quadrature_points() ?

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

[Timo Heister]

On Thu, Jan 10, 2013 at 10:31 PM, Zhicheng Wang
<wangzhi...@gmail.com> wrote:
> Hi, Timo, Thank you for your reply, but constraint.reinit() could not solve
> the problem. Is it possible
> the order of face->get_dofs_indices() is not agree with that of
> FeFaceValue.get_quadrature_points() ?

unlikely. You should try to output the ConstraintMatrices for a
testproblem and check if they are consistent with each other
(information in there should be the same or only found in either one).