Distribute constraints when using Workstream class

[Jaekwang Kim]

Hi all.... 

I am practicing the use of Workstream class to paralleize my assemble_system function. 

As a step stone, I am trying to solve simple possion equation with Workstream class.

To learn how to use Workstream class, I am referring step-9 tutorial.

However, I become little bit confused in distributing my constraints matrix to system matrix..

For example, my original set_up function and assemble system function is as follow...

template <int dim>

void Step4<dim>::setup_system ()

{

  dof_handler.distribute_dofs (fe);

 

  std::cout << "   Number of degrees of freedom: "

            << dof_handler.n_dofs()

            << std::endl;

 

    

  DynamicSparsityPattern dsp(dof_handler.n_dofs());

  DoFTools::make_sparsity_pattern (dof_handler, dsp);

  sparsity_pattern.copy_from(dsp);

    

  system_matrix.reinit (sparsity_pattern);

  solution.reinit (dof_handler.n_dofs());

  system_rhs.reinit (dof_handler.n_dofs());

    

    

  constraints.clear ();

    

  DoFTools::make_hanging_node_constraints (dof_handler,

                                             constraints);

    

  VectorTools::interpolate_boundary_values (dof_handler,0,

                                              BoundaryValues<dim>(),

                                              constraints);

    

  constraints.close ();

 

}

 

template <int dim>

void Step4<dim>::assemble_system ()

{

  QGauss<dim>  quadrature_formula(2);

 

  const RightHandSide<dim> right_hand_side;

 

  FEValues<dim> fe_values (fe, quadrature_formula,

                           update_values   | update_gradients |

                           update_quadrature_points | update_JxW_values);

 

  const unsigned int   dofs_per_cell = fe.dofs_per_cell;

  const unsigned int   n_q_points    = quadrature_formula.size();

 

  FullMatrix<double>   cell_matrix (dofs_per_cell, dofs_per_cell);

  Vector<double>       cell_rhs (dofs_per_cell);

 

  std::vector<types::global_dof_index> local_dof_indices (dofs_per_cell);

 

  typename DoFHandler<dim>::active_cell_iterator

  cell = dof_handler.begin_active(),

  endc = dof_handler.end();

 

  for (; cell!=endc; ++cell)

    {

      fe_values.reinit (cell);

      cell_matrix = 0;

      cell_rhs = 0;

 

      for (unsigned int q_index=0; q_index<n_q_points; ++q_index)

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

          {

            for (unsigned int j=0; j<dofs_per_cell; ++j)

              cell_matrix(i,j) += (fe_values.shape_grad (i, q_index) *

                                   fe_values.shape_grad (j, q_index) *

                                   fe_values.JxW (q_index));

 

            cell_rhs(i) += (fe_values.shape_value (i, q_index) *

                            right_hand_side.value (fe_values.quadrature_point (q_index)) *

                            fe_values.JxW (q_index));

          }

 

      cell->get_dof_indices (local_dof_indices);

  

     constraints.distribute_local_to_global (cell_matrix,cell_rhs,local_dof_indices,system_matrix,system_rhs);

    }

    

  

}

Maintaining same set_up function, where my constraints matrix is being generated, 

I tried to make another assemble_system function, that is being operated with Workstream class, like as follow

template <int dim>

void Step4<dim>::assemble_system_new ()

{

    WorkStream::run(dof_handler.begin_active(),

                    dof_handler.end(),

                    *this,

                    &Step4::local_assemble_system,

                    &Step4::copy_local_to_global,

                    AssemblyScratchData(fe),

                    AssemblyCopyData());

    

}

 

template <int dim>

Step4<dim>::AssemblyScratchData::

AssemblyScratchData (const FiniteElement<dim> &fe)

:

fe_values (fe,

           QGauss<dim>(2),

           update_values   | update_gradients |

           update_quadrature_points | update_JxW_values)

{}

 

 

template <int dim>

Step4<dim>::AssemblyScratchData::

AssemblyScratchData (const AssemblyScratchData &scratch_data)

:

fe_values (scratch_data.fe_values.get_fe(),

           scratch_data.fe_values.get_quadrature(),

           update_values   | update_gradients |

           update_quadrature_points | update_JxW_values)

{}

 

 

 

template <int dim>

void

Step4<dim>::

local_assemble_system (const typename DoFHandler<dim>::active_cell_iterator &cell,

                       AssemblyScratchData                                  &scratch_data,

                       AssemblyCopyData                                     &copy_data)

{

    const unsigned int dofs_per_cell   = fe.dofs_per_cell;

    const unsigned int n_q_points  = scratch_data.fe_values.get_quadrature().size();

    

    copy_data.cell_matrix.reinit (dofs_per_cell, dofs_per_cell);

    copy_data.cell_rhs.reinit (dofs_per_cell);

    copy_data.local_dof_indices.resize(dofs_per_cell);

    

    

    const RightHandSide<dim>  right_hand_side;

    std::vector<double>         rhs_values (n_q_points);

 

    scratch_data.fe_values.reinit (cell);

    right_hand_side.value_list (scratch_data.fe_values.get_quadrature_points(),

                                rhs_values);

    

    

    for (unsigned int q_index=0; q_index<n_q_points; ++q_index)

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

        {

            for (unsigned int j=0; j<dofs_per_cell; ++j)

                

            { copy_data.cell_matrix(i,j) += scratch_data.fe_values.shape_grad(i,q_index)*

                                                scratch_data.fe_values.shape_grad(j,q_index)*

                                            scratch_data.fe_values.JxW(q_index);

            }

            

                copy_data.cell_rhs(i) += scratch_data.fe_values.shape_value(i,q_index)*

                                        right_hand_side.value(scratch_data.fe_values.quadrature_point(q_index))*

                                        scratch_data.fe_values.JxW(q_index);

            

        }

 

    cell->get_dof_indices (copy_data.local_dof_indices);

}

 

 

template <int dim>

void

Step4<dim>::copy_local_to_global (const AssemblyCopyData &copy_data)

{

    for (unsigned int i=0; i<copy_data.local_dof_indices.size(); ++i)

    {

        for (unsigned int j=0; j<copy_data.local_dof_indices.size(); ++j)

            system_matrix.add (copy_data.local_dof_indices[i],

                               copy_data.local_dof_indices[j],

                               copy_data.cell_matrix(i,j));

        system_rhs(copy_data.local_dof_indices[i]) += copy_data.cell_rhs(i);

        

    }

    

     constraints.distribute_local_to_global (copy_data.cell_matrix,copy_data.cell_rhs,copy_data.local_dof_indices,system_matrix,system_rhs);

    // I am suspicious on this part...

}

However, those two different assemble_system function result different solution

While the result from original assemble function looks qualitatively correct, the result from new assemble function is weird. 

Is there anyone who can figure out what I was doing wrong...?

<original assemble function , qualitatively correct >

<assemble function modified, with Workstream Class, wrong

[Wolfgang Bangerth]

On 04/06/2017 02:55 PM, Jaekwang Kim wrote:
> template<intdim>
>
> void
>
> Step4<dim>::copy_local_to_global (constAssemblyCopyData &copy_data)
>
> {
>
> for(unsignedinti=0; i<copy_data.local_dof_indices.size(); ++i)
>
> {
>
> for(unsignedintj=0; j<copy_data.local_dof_indices.size(); ++j)

>
> system_matrix.add (copy_data.local_dof_indices[i],
>
> copy_data.local_dof_indices[j],
>
> copy_data.cell_matrix(i,j));
>
> system_rhs(copy_data.local_dof_indices[i]) += copy_data.cell_rhs(i);
>
>
>
> }
>
>
>

> *constraints.distribute_local_to_global
> (copy_data.cell_matrix,copy_data.cell_rhs,copy_data.local_dof_indices,system_matrix,system_rhs);*
>
> * // I am suspicious on this part...*
>
> }

You copy things twice into the global matrix and rhs, but only once taking
into accounts the constraints. Take a look at step-32, for example, to see how
this function is supposed to look.

Best
W.


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

[Jaekwang Kim]

I just understand what you meant ! 

Thanks !

2017년 4월 6일 목요일 오후 7시 24분 1초 UTC-5, Wolfgang Bangerth 님의 말: