Mesh refinement with periodic boundary conditions

[Andrew Davis]

HI all, I'm trying to implement a time-dependent convection equation using DG elements with periodic boundary conditions on and adaptive mesh. However, when I try to adapt the mesh using periodic using this code:   // estimate the error in each cell   Vector<float> estimatedErrorPerCell(triangulation.n_active_cells());   //KellyErrorEstimator<dim>::estimate(dofHandler, QGauss<dim-1>(feSystem.degree + 1), {}, solution.block(0), estimatedErrorPerCell);   KellyErrorEstimator<dim>::estimate(dofHandler, QGauss<dim-1>(feSystem.degree + 1), {}, solution, estimatedErrorPerCell);   // refine cells with the largest estimated error (80 per cent of the error) and coarsens those cells with the smallest error (10 per cent of the error)   GridRefinement::refine_and_coarsen_fixed_fraction(triangulation, estimatedErrorPerCell, 0.8, 0.1);   // to prevent the decrease in time step size limit the maximal refinement depth of the meshlimit the maximal refinement depth of the mesh   if( triangulation.n_levels()>maxGridLevel ) {     for( auto& cell : triangulation.active_cell_iterators_on_level(maxGridLevel) ) { cell->clear_refine_flag(); }   }   // store previous and current solution on the coarser mesh   std::vector<BlockVector<double> > tempSolution(2);   tempSolution[0] = oldSolution; tempSolution[1] = solution;   // transfer the solution vectors from the old mesh to the new one   SolutionTransfer<dim, BlockVector<double> > transfer(dofHandler);   // prepare for the refinement---we have to prepare the transfer as well   triangulation.prepare_coarsening_and_refinement(); // for some reason, this doesn't tell the cells across a periodic boundary that they need to be refined   transfer.prepare_for_coarsening_and_refinement(tempSolution);   // do the refinement and reset the DoFs   triangulation.execute_coarsening_and_refinement();   SetupSystem();   std::vector<BlockVector<double> > temp(2);   ReinitializeVector(temp[0]); ReinitializeVector(temp[1]);   transfer.interpolate(tempSolution, temp);   oldSolution = temp[0]; solution = temp[1];   constraints.distribute(oldSolution);   constraints.distribute(solution); I get this error: -------------------------------------------------------- An error occurred in line <992> of file </home/add8536/Software/dealii/source/grid/tria.cc> in function     void dealii::{anonymous}::update_periodic_face_map_recursively(const typename dealii::Triangulation<dim, spacedim>::cell_iterator&, const typename dealii::Triangulation<dim, spacedim>::cell_iterator&, unsigned int, unsigned int, const std::bitset<3>&, std::map<std::pair<typename dealii::Triangulation<dim, spacedim>::cell_iterator, unsigned int>, std::pair<std::pair<typename dealii::Triangulation<dim, spacedim>::cell_iterator, unsigned int>, std::bitset<3> > >&) [with int dim = 2; int spacedim = 2; typename dealii::Triangulation<dim, spacedim>::cell_iterator = dealii::TriaIterator<dealii::CellAccessor<2, 2> >] The violated condition was:     dim == 1 || std::abs(cell_1->level() - cell_2->level()) < 2 Additional information:     This exception -- which is used in many places in the library -- usually indicates that some condition which the author of the code thought must be satisfied at a certain point in an algorithm, is not fulfilled. An example would be that the first part of an algorithm sorts elements of an array in ascending order, and a second part of the algorithm later encounters an element that is not larger than the previous one. There is usually not very much you can do if you encounter such an exception since it indicates an error in deal.II, not in your own program. Try to come up with the smallest possible program that still demonstrates the error and contact the deal.II mailing lists with it to obtain help. For some reason the refinement does not seem to know that the boundaries are periodic and that it should refine the cells on the other side of the domain. I tried changing the refinement step so that instead of just calling triangulation.prepare_coarsening_and_refinement();, I created my own function that mimics the distributed triangulation (calling the distributed version of enforce_mesh_balance_over_periodic_boundaries) that does the following: while( true ) {     triangulation.prepare_coarsening_and_refinement();     triangulation.update_periodic_face_map();     bool changed = enforce_mesh_balance_over_periodic_boundaries(triangulation);     if( !changed ) { break; }   } This fixed the error but I'm getting the strange behavior---it looks like the cell face fluxes are no longer being computed correctly Does anyone have any idea about what the issue could be? I know I left out a lot of code, but hopefully I included the important bits. Thanks, Andy ReplyForward

[Daniel Arndt]

Andrew,

Did you tell the triangulation that it should take periodic boundaries into account? My best bet (with the information given) would be that you didn't call parallel::distributed::Triangulation::add_periodicity()

as described in https://www.dealii.org/current/doxygen/deal.II/step_45.html.

Best,

Daniel

--
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.
To view this discussion on the web visit https://groups.google.com/d/msgid/dealii/fdcc9fe0-0a99-4db5-bf49-f447ed26aa4e%40googlegroups.com.

[Andrew Davis]

Yes, here is the code I used to tell the triangulation about the periodic boundaries: 

  // create a mesh on [0,1] with colorized boundaries
  GridGenerator::hyper_cube(triangulation, 0.0, 1.0, true);

  // collect the periodic boundary pairs
  std::vector<GridTools::PeriodicFacePair<typename parallel::distributed::Triangulation<dim>::cell_iterator> > periodicityVector;
  GridTools::collect_periodic_faces(triangulation, 0, 1, 0, periodicityVector);
  GridTools::collect_periodic_faces(triangulation, 2, 3, 1, periodicityVector);

  triangulation.add_periodicity(periodicityVector);

  // refine one less then maximum refinement level globally
  triangulation.refine_global(maxGridLevel-1);

Perhaps I'm doing it wrong?

On Thursday, February 27, 2020 at 10:47:09 PM UTC-5, Daniel Arndt wrote:

Andrew,

Did you tell the triangulation that it should take periodic boundaries into account? My best bet (with the information given) would be that you didn't call parallel::distributed::Triangulation::add_periodicity()

as described in https://www.dealii.org/current/doxygen/deal.II/step_45.html.

Best,

Daniel

To unsubscribe from this group and stop receiving emails from it, send an email to dea...@googlegroups.com.

[Andrew Davis]

I just say a slight problem and changed this line the declaration of periodicityVector to:

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

But I'm getting the same behavior.

[Daniel Arndt]

Andrew,

Still hard to say what the problem is. Which commit are you using?

We have quite a number of tests that make sure that adaptive grid refinement also works with periodic boundary conditions.

Can you provide a full minimal example that shows the problem you need your workaround for?

How are you assembling the linear system? Are you using MeshWorker, MatrixFree or do you assemble matrices on your own using FEFaceValues?

Is your code working in serial and for global refinement?

Best,

Daniel

To unsubscribe from this group and stop receiving emails from it, send an email to dealii+un...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/dealii/3c417b20-5136-4182-bd4d-c3d2b31c0007%40googlegroups.com.

[Andrew Davis]

The commit that I am using is: 489c9d198d4ff0654f7d1f69f09404a2613b7954

I will work on a minimal example, but in the mean time here are the functions I use to create the sparsity pattern and assemble the system. The code works perfectly if either I turn off mesh refinement but keep periodic BCs OR keep mesh refinement but instead prescribe in/outflow boundary conditions.  I'm solving the advection equation c_{t} + \nabla \cdot (u c) = f for c with prescribed velocity u and forcing f using Crank-Nicolson as my time integrator.

I'm using MeshWorker to assemble the system. 

The code is working in serial for global refinement.

Set up sparsity pattern:

template<unsigned int dim>
void SeaIceModel<dim>::SetupDoFs() {
  // set up the DoFs
  dofHandler.distribute_dofs(feSystem);

  // reset the constraints
  constraints.clear();
  

  DoFTools::make_hanging_node_constraints(dofHandler, constraints);

  std::vector<GridTools::PeriodicFacePair<typename DoFHandler<dim>::cell_iterator> > periodicityVector;
  GridTools::collect_periodic_faces(dofHandler, 0, 1, 0, periodicityVector);
  GridTools::collect_periodic_faces(dofHandler, 2, 3, 1, periodicityVector);
  DoFTools::make_periodicity_constraints<DoFHandler<dim> >(periodicityVector, constraints);

  constraints.close();

  // get the number of degrees of freedom for each block
  std::vector<unsigned int> ndofsPerBlock(nblocks);
  DoFTools::count_dofs_per_block(dofHandler, ndofsPerBlock);

  BlockDynamicSparsityPattern dsp(nblocks, nblocks);

  // initalize the sparsity pattern
  for( unsigned int i=0; i<nblocks; ++i ) {
    for( unsigned int j=0; j<nblocks; ++j ) {
      dsp.block(i, j).reinit(ndofsPerBlock[i], ndofsPerBlock[j]);
    }
  }
  dsp.collect_sizes();

  // generate the sparsity pattern
  DoFTools::make_flux_sparsity_pattern(dofHandler, dsp, constraints);
  sparsityPattern.copy_from(dsp);

}

Assemble the system:

template<unsigned int dim>
void SeaIceModel<dim>::AssembleSystem(double const prevTime, double const currTime) {
  // reset the system matrix and rhs to zero
  systemMatrix = 0.0;
  systemRHS = 0.0;

  MeshWorker::IntegrationInfoBox<dim> infoBox;

  const unsigned int nGaussPoints = feSystem.degree + 1;
  infoBox.initialize_gauss_quadrature(nGaussPoints, nGaussPoints, nGaussPoints);

  infoBox.initialize_update_flags();
  infoBox.add_update_flags_cell(update_quadrature_points | update_values | update_gradients | update_JxW_values);
  infoBox.add_update_flags(update_normal_vectors | update_quadrature_points | update_values | update_JxW_values, false, true, true, true);

  infoBox.initialize(feSystem, mapping);

  MeshWorker::DoFInfo<dim> dofInfo(dofHandler);

  MeshWorker::Assembler::SystemSimple<BlockSparseMatrix<double>, BlockVector<double> >
  assembler;
  assembler.initialize(systemMatrix, systemRHS);

  // objects to integrate over cells and boundaries---these objects have overloaded operator() functions for the MeshWorker
  SystemCellIntegration<dim> cellIntegration(eqnBlocks, extractors, solutionIndexMap, oldSolution, solution, std::tuple<double, double, double>(prevTime, currTime, timeStep), nfields);
  SystemBoundaryIntegration<dim> boundaryIntegration(eqnBlocks, extractors, solutionIndexMap, oldSolution, solution, std::tuple<double, double, double>(prevTime, currTime, timeStep), nfields);
  SystemFaceIntegration<dim> faceIntegration(eqnBlocks, extractors, solutionIndexMap, oldSolution, solution, std::tuple<double, double, double>(prevTime, currTime, timeStep), nfields);

  MeshWorker::loop<dim, dim, MeshWorker::DoFInfo<dim>, MeshWorker::IntegrationInfoBox<dim> > (dofHandler.begin_active(), dofHandler.end(), dofInfo, infoBox,
    cellIntegration, boundaryIntegration, faceIntegration, assembler);
}

On Friday, February 28, 2020 at 11:19:57 AM UTC-5, Daniel Arndt wrote:

Andrew,

Still hard to say what the problem is. Which commit are you using?

We have quite a number of tests that make sure that adaptive grid refinement also works with periodic boundary conditions.

Can you provide a full minimal example that shows the problem you need your workaround for?

How are you assembling the linear system? Are you using MeshWorker, MatrixFree or do you assemble matrices on your own using FEFaceValues?

Is your code working in serial and for global refinement?

Best,

Daniel

To view this discussion on the web visit https://groups.google.com/d/msgid/dealii/3c417b20-5136-4182-bd4d-c3d2b31c0007%40googlegroups.com.

[Andrew Davis]

Here is a slight modification to example step-12 that fails.  This is the same sample (steady state) as step-12. The only difference is that I made the system periodic in the x direction by adding these lines:

            GridGenerator::hyper_cube(triangulation, 0.0, 1.0, true);

            // collect the periodic boundary pairs

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

            GridTools::collect_periodic_faces(triangulation, 0, 1, 0, periodicityVector);

            triangulation.add_periodicity(periodicityVector);

            triangulation.refine_global(3);

Attached is the full code.

[David]

Hi Andrew,

I'm facing exactly the same issue, which means either mesh refinement without periodic boundaries works fine or the other way around, but not both at the same time.

The code uses DG methods with matrix-free techniques.

It seems, that the refinement destroys the periodicity somehow. The fluxes are incorrect calculated after the refinement and (after a while) I get an exception referring to an inconsistent mesh refinement across the boundary as you already mentioned above.

do you managed to solve this issue?

Best regards
David

[Andrew Davis]

Hello,

I have not managed to solve this issue. I have gotten as far as diagnosing that the error occurs when the refinement needs to happen "across" the period boundary. For some reason, this does not happen and then when you go to solve the problem you end up with an inconsistent mesh. For the toy problems I'm using to develop my code I'm just using a globally refined mesh where this isn't a problem.  This will not be scalable so I will need to tackle this issue at some point, but I haven't gotten around to it yet.

Has anyone else managed to find a solution?

-Andy