Periodic boundary condition by using ConstraintMatrix::condense()

[Chenchen Liu]

Hi all,

I want to do a time-dependent problem with periodic boundary condition. From the previous discussion https://groups.google.com/forum/#!searchin/dealii/apply$20constraintmatrix/dealii/kovM6_AGSf4/zo7c7NSvL1QJ, I understand that I need to use ConstraintMatrix::condense() after assembling to apply any constraints in my time dependent problem. 

As a practice, I want to modify step-45 as an equivalent single time step. As I mentioned before, I did not put any constraints before or during the assembling of system_matrix and system_rhs. After I finished assembling them in a way like step-4, I began to add periodic constraints (left and right boundary) and Dirichlet constraints (top and bottom boundary with indicator 1). The following is my code:

      constraints.clear ();

      make_periodicity_constraints ();

      std::map<types::global_dof_index,double> boundary_values;

      VectorTools::interpolate_boundary_values (dof_handler, 1,

                                                                            ZeroFunction<2> (),

                                                                            boundary_values);

      constraints.close ();

      constraints.condense (system_matrix, system_rhs);

      MatrixTools::apply_boundary_values (boundary_values,

                                                                   system_matrix,

                                                                   solution,

                                                                   system_rhs);

However, it gave me an error:

An error occurred in line <1723> of file <../include/deal.II/lac/sparse_matrix.h> in function

    void dealii::SparseMatrix<double>::add(const size_type, const size_type, const number) [number = double]

The violated condition was: 

    (index != SparsityPattern::invalid_entry) || (value == number()).

The name and call sequence of the exception was:

    ExcInvalidIndex(i, j)

Additional Information: 

The entry with index <389,1> does not exist.

I revised the code in the following,

      constraints.clear ();

      make_periodicity_constraints ();

      VectorTools::interpolate_boundary_values (dof_handler, 1,

                                                                            ZeroFunction<2> (),

                                                                            constraints);

      constraints.close ();

      constraints.condense (system_matrix, system_rhs);

which still gave me similar error.

Did I make any mistake when using constraints.condense() or should I also revise any other part of the code? Thanks!

Chenchen

[Daniel Arndt]

Chenchen,

from your description it's hard to tell what the problem is. It seems that there is something wrong in your make_periodicity_constraints.
Preferably, you just try to use the already implemented functions as explained in the recent version of step-45 [1].
The old version really restricts to Q1 elements.
In case this does not solve your problem, it would be helpful if you could provide some more information about the elements and the mesh you are using.

Bests,
Daniel

[1] https://www.dealii.org/developer/doxygen/deal.II/step_45.html

[Chenchen Liu]

Thanks Daniel.

I am pretty sure that the function "make_periodicity_constraints()" is correct, which is an existing function in old version of step-45 (https://www.dealii.org/8.3.0/doxygen/deal.II/step_45.html). I am dealing with time dependent problem, in which the mass matrix or system matrix remains the same for all time steps. In order to save computational time, I want to save the system_matrix at the beginning of the code, and call it at each time step. Then, at each time step, I want to apply periodic boundary condition and Dirichlet boundary condition to the system_matrix and system_rhs.

Now the problem is, when I use the following codes:

      constraints.clear ();

      make_periodicity_constraints ();

      constraints.close ();

if I put the above lines before I define "DoFTools::make_sparsity_pattern (dof_handler, dsp, constraints, true);", it works perfectly.

If I put the same codes later and don't consider "constraints" in "make_sparsity_pattern" (the code is listed in previous email), it does not work. It seems that if I want to apply periodic constraints, I have to put the constraints in the definition of sparsity pattern. In other words, in order to apply periodic boundary condition, at each time step, I have to re-assemble the system_matrix... However, this is what I want to avoid.

Best,

Chenchen

在 2015年11月19日星期四 UTC-5下午2:00:36，Daniel Arndt写道：

[Daniel Arndt]

Chenchen,

Now the problem is, when I use the following codes:

      constraints.clear ();

      make_periodicity_constraints ();

      constraints.close ();

if I put the above lines before I define "DoFTools::make_sparsity_pattern (dof_handler, dsp, constraints, true);", it works perfectly.

If I put the same codes later and don't consider "constraints" in "make_sparsity_pattern" (the code is listed in previous email), it does not work. It seems that if I want to apply periodic constraints, I have to put the constraints in the definition of sparsity pattern. In other words, in order to apply periodic boundary condition, at each time step, I have to re-assemble the system_matrix... However, this is what I want to avoid.

By default (without ConstraintMatrix) DoFTools::make_sparsity_pattern just creates entries for degrees of freedom that share at least one common cell. 

Therefore entries for DoFs coupling through the periodic boundary condition are not considered and you have to provide make_sparsity_pattern with the ConstraintMatrix.

As far as I understand in your problem neither the mesh nor the constraints change between time steps. Therefore ConstraintMatrix and SparsityPattern do not change.

You don't need to assemble system_matrix in each time step, but you can store the required unconstrained matrices and copy them in each time step. Have a look at step-23 for this.

Best,

Daniel

[Chenchen Liu]

Thanks Daniel. I see your point. So for sure, we have to make sparsity pattern with constraints. In the following, we can still modify the constraints in each time step without re-assembling system matrix.

Best,

Chenchen

在 2015年11月20日星期五 UTC-5上午5:23:06，Daniel Arndt写道：