MeshWorker vs Kelly (refined cells)?

[Denis Davydov]

Dear developers,

I was playing around with MeshWorker, in particular I tried to use it in step-16 instead of Kelly estimator, following ideas of Step-39.

Interesting part is that i was not able to make them identical... The local integrator is as follows

  template <int dim>

  void EstimatorIntegrator<dim>::face(MeshWorker::DoFInfo<dim> &dinfo1,

                                      MeshWorker::DoFInfo<dim> &dinfo2,

                                      typename MeshWorker::IntegrationInfo<dim> &info1,

                                      typename MeshWorker::IntegrationInfo<dim> &info2) const

  {

    const FEValuesBase<dim> &fe = info1.fe_values();

    const std::vector<Tensor<1,dim> > &Duh1 = info1.gradients[0][0];

    const std::vector<Tensor<1,dim> > &Duh2 = info2.gradients[0][0];

    const double h_K1 = dinfo1.cell->measure();

    const double h_K2 = dinfo2.cell->measure();

    double integral = 0.;

    for (unsigned k=0; k<fe.n_quadrature_points; ++k)

      {

        const double diff = fe.normal_vector(k) * Duh1[k] - fe.normal_vector(k) * Duh2[k];

        integral += diff * diff * fe.JxW(k);

      }

    dinfo1.value(0) = std::sqrt((h_K1/24.)*integral);

    dinfo2.value(0) = std::sqrt((h_K2/24.)*integral);

  }

The documentation of MeshWorker classes is a bit limited. I could find out how one treats integrals on faces when there are two elements with different refinement level. This is likely a reason for discrepancy, perhaps someone who knows the MeshWorker classes better can explain this or better yet update the documentation.

As for the Kelly estimator class, this case is considered in integrate_over_irregular_face() using FESubFaceValues() and this case is handled from the perspective of a coarser cell. Corresponding parts of the face contribution to the error are stored with each child face.

Btw, Step-39 should be updated as volumetric part of the estimator uses dinfo.cell->diameter(), whereas jump terms are scaled with dinfo1.face->measure(); 

Regards,

Denis.

[Wolfgang Bangerth]

On 04/04/2016 10:34 AM, Denis Davydov wrote:
>
> I was playing around with MeshWorker, in particular I tried to use it in
> step-16 instead of Kelly estimator, following ideas of Step-39.
> Interesting part is that i was not able to make them identical...

Are they just a bit different, or completely different? Does it yield the same
results on uniformly refined meshes?

The thing is that, collectively, we have forgotten how MeshWorker works. As
you (and others) found out, the documentation by itself is inadequate, and
nobody who's active on the mailing list knows how it all works :-(

Best
Wolfgang


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

[Denis Davydov]

Hi Wolfgang,

On Monday, April 4, 2016 at 6:09:25 PM UTC+2, bangerth wrote:

On 04/04/2016 10:34 AM, Denis Davydov wrote:
>
> I was playing around with MeshWorker, in particular I tried to use it in
> step-16 instead of Kelly estimator, following ideas of Step-39.
> Interesting part is that i was not able to make them identical...

Are they just a bit different, or completely different? Does it yield the same

they are completely different! 

 

results on uniformly refined meshes?

no :( see below.

 

The thing is that, collectively, we have forgotten how MeshWorker works. As

that's a pity... If that is the case, perhaps we should rework the examples so that

we show how to use classes for which we still remember what they do exactly 

and which are well documented?

 

you (and others) found out, the documentation by itself is inadequate, and
nobody who's active on the mailing list knows how it all works :-(

oh well... Perhaps Martin or Timo can comment on the best practices for assembly of geometric multigrids?

Looking at the more recent Step-37, i guess i should do assembly as usual, by hand.

I suppose one can use Worstream class. I will eventually move to MPI, however now I am doing a quick prototype. 

But MPI+TBB sounds good for assembly of MG matrices. What do you think?

p.s. Here's the output for h-adaptive refinement with Kelly and MeshWorker (first cell is homogeneous, even there the L2 norm of the error is different by a factor of two at least):

Kelly:

-------

DEAL::Cycle 0

DEAL::   Number of active cells:       20

DEAL::   Number of degrees of freedom: 25 (by level: 8, 25)

DEAL:cg::Starting value 0.510691

DEAL:cg::Convergence step 6 value 4.59193e-14

DEAL::Cycle 1

DEAL::error: 1.02041

DEAL::   Number of active cells:       41

DEAL::   Number of degrees of freedom: 52 (by level: 8, 25, 41)

DEAL:cg::Starting value 0.455356

DEAL:cg::Convergence step 8 value 3.09682e-13

DEAL::Cycle 2

DEAL::error: 0.778599

DEAL::   Number of active cells:       80

DEAL::   Number of degrees of freedom: 100 (by level: 8, 25, 61, 52)

DEAL:cg::Starting value 0.394469

DEAL:cg::Convergence step 9 value 1.96993e-13

DEAL::Cycle 3

DEAL::error: 0.668332

DEAL::   Number of active cells:       161

DEAL::   Number of degrees of freedom: 190 (by level: 8, 25, 77, 160)

DEAL:cg::Starting value 0.322156

DEAL:cg::Convergence step 9 value 2.94418e-13

MeshWorker:

-----------

DEAL::Cycle 0

DEAL::   Number of active cells:       20

DEAL::   Number of degrees of freedom: 25 (by level: 8, 25)

DEAL:cg::Starting value 0.510691

DEAL:cg::Convergence step 6 value 4.59193e-14

DEAL::Cycle 1

DEAL::error: 0.438336

DEAL::   Number of active cells:       41

DEAL::   Number of degrees of freedom: 52 (by level: 8, 25, 41)

DEAL:cg::Starting value 0.455356

DEAL:cg::Convergence step 8 value 3.09682e-13

DEAL::Cycle 2

DEAL::error: 0.294554

DEAL::   Number of active cells:       80

DEAL::   Number of degrees of freedom: 100 (by level: 8, 25, 57, 55)

DEAL:cg::Starting value 0.398498

DEAL:cg::Convergence step 9 value 2.80944e-13

DEAL::Cycle 3

DEAL::error: 0.223672

DEAL::   Number of active cells:       155

DEAL::   Number of degrees of freedom: 181 (by level: 8, 25, 77, 147)

DEAL:cg::Starting value 0.322165

DEAL:cg::Convergence step 9 value 5.09767e-13

[Denis Davydov]

On Monday, April 4, 2016 at 6:09:25 PM UTC+2, bangerth wrote:

The thing is that, collectively, we have forgotten how MeshWorker works. As
you (and others) found out, the documentation by itself is inadequate, and
nobody who's active on the mailing list knows how it all works :-(

It seems that step-16 is the only tutorial which does h-adaptive geometric MG and is currently bound to MeshWorker.

I will try to re-write mg-assembly loops explicitly and hopefully will make it work without a MeshWorker. 

But if anyone here already did that, please let me know.

Cheers,

Denis.

[Martin Kronbichler]

Dear Denis,

The unfinished step-50 should contain the integrals from step-16 but without MeshWorker and does GMG. It is not DG, but neither is step-16.

Best,
Martin=

[Denis Davydov]

Dear Martin,

Thanks for your reply, that will certainly save me a lot of time. 

I do not need DG, I am playing with h-daptive GMG on Laplace problem for now.

I will have a closer loot at step-50…

Have a good day,

Denis.

--
You received this message because you are subscribed to a topic in the Google Groups "deal.II developers" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/dealii-developers/66vk8nk7SKA/unsubscribe.
To unsubscribe from this group and all its topics, send an email to dealii-develop...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[Wolfgang Bangerth]

On 04/05/2016 03:35 AM, Denis Davydov wrote:
>
> It seems that step-16 is the only tutorial which does h-adaptive geometric MG
> and is currently bound to MeshWorker.
> I will try to re-write mg-assembly loops explicitly and hopefully will make it
> work without a MeshWorker.
> But if anyone here already did that, please let me know.

step-16 didn't use to use MeshWorker but was re-written with it a few years
ago. If you go back in history for this file, you may discover a version with
explicit assembly loops.

Best
W.

[Wolfgang Bangerth]

On 04/04/2016 01:47 PM, Denis Davydov wrote:
> Hi Wolfgang,
>
> On Monday, April 4, 2016 at 6:09:25 PM UTC+2, bangerth wrote:
>
> On 04/04/2016 10:34 AM, Denis Davydov wrote:
> >
> > I was playing around with MeshWorker, in particular I tried to use it in
> > step-16 instead of Kelly estimator, following ideas of Step-39.
> > Interesting part is that i was not able to make them identical...
>
> Are they just a bit different, or completely different? Does it yield the
> same
>
>
> they are completely different!

That's too bad. I don't think I can be of any help. For example, in your code here

const std::vector<Tensor<1,dim> > &Duh1 = info1.gradients[0][0];
const std::vector<Tensor<1,dim> > &Duh2 = info2.gradients[0][0];

try this with a linear solution, for which you know what the gradients should
be. Then make sure that the jump in gradients in this function is really
computed to be zero. I don't think I have any better suggestions than just
single stepping through the function and debugging what you have and what you get.

> that's a pity... If that is the case, perhaps we should rework the examples so
> that
> we show how to use classes for which we still remember what they do exactly
> and which are well documented?

I'll let others make this call.

> p.s. Here's the output for h-adaptive refinement with Kelly
> and MeshWorker (first cell is homogeneous, even there the L2 norm of the error
> is different by a factor of two at least):

Start with a mesh that has only one cell and single step through your
integration function. That's the best I can offer :-(

Best
W.

[Denis Davydov]

Thanks for your suggestions, Wolfgang.

I already made it work without MeshWorker following Step-50.

I also ported most of the assembly to Workstream so i will stick with that solution for now

as it is more transparent and will eventually play nicely with MPI as shown in Step-32.

On Tuesday, April 5, 2016 at 1:57:07 PM UTC+2, bangerth wrote:

On 04/04/2016 01:47 PM, Denis Davydov wrote:
> Hi Wolfgang,
>
> On Monday, April 4, 2016 at 6:09:25 PM UTC+2, bangerth wrote:
>
> 
> that's a pity... If that is the case, perhaps we should rework the examples so
> that
> we show how to use classes for which we still remember what they do exactly
> and which are well documented?

I'll let others make this call.

Fair enough.

Cheers,

Denis.