Artifacts at the refinement level boundary for Step 48

[Stephen DeWitt]

Hello,

I've been working on a phase field code using the matrix free approach from Step 48. When I use adaptive meshing, I notice some artifacts at the boundary between levels of refinement and I'm trying to understand their origin.

For example, if I make the following modification to the "local_apply" function in Step 48, the solution should remain unchanged from time step to time step:

template <int dim, int fe_degree>
  void SineGordonOperation<dim, fe_degree>::
  local_apply (const MatrixFree<dim>                      &data,
               parallel::distributed::Vector<double>      &dst,
               const std::vector<parallel::distributed::Vector<double>*> &src,
               const std::pair<unsigned int,unsigned int> &cell_range) const
  {
    AssertDimension (src.size(), 2);
    FEEvaluation<dim,fe_degree> current (data), old (data);
    for (unsigned int cell=cell_range.first; cell<cell_range.second; ++cell)
      {
        current.reinit (cell);
        old.reinit (cell);


        current.read_dof_values (*src[0]);
        old.read_dof_values     (*src[1]);


        current.evaluate (true, true, false);
        old.evaluate (true, false, false);


        for (unsigned int q=0; q<current.n_q_points; ++q)
          {
            const VectorizedArray<double> current_value = current.get_value(q);
            const VectorizedArray<double> old_value     = old.get_value(q);


            // ORIGINAL
            // current.submit_value (2.*current_value - old_value -
            //                       delta_t_sqr * std::sin(current_value),q);
            // current.submit_gradient (- delta_t_sqr *
            //                          current.get_gradient(q), q);


             // CHANGES
             current.submit_value (current_value,q);
             current.submit_gradient (- delta_t_sqr *
                                      current.get_gradient(q)*0.0, q);
          }


        current.integrate (true,true);
        current.distribute_local_to_global (dst);
      }
  }

However, there is a non-trivial change in the solution (same plot with and without the mesh):

Has anyone noticed this before? In this example, the magnitude of the artifacts is low, but in some of my other tests, the magnitude has been on the order of the "real" change in the solution. When doing implicit time stepping (still in the matrix free framework), the artifacts don't show up, which suggests that they aren't just a consequence of having an adapted mesh. In the Step 48 description, it's noted that the inverse of the mass matrix isn't computed exactly. Could the artifacts be related to that?

For Step 48, the artifacts only appear from the initial time step to the first time step (in my application, I'm getting the artifacts every time step and I'm still in the midst of finding what the source of that difference is).

Thank you for any insight you can give.

Best,

Steve

Stephen DeWitt
Asst. Research Scientist
PRISMS Center
Dept. of Materials Science and Engineering
University of Michigan

[Martin Kronbichler]

Dear Stephen,

At hanging nodes, there is definitely going to be a larger error due to
the approximation of the diagonal mass matrix. I do not remember the
exact details but to get a diagonal mass matrix you need to assume an
interpolation in addition to the approximation leading to the mass
lumping. If I remember correctly my co-worker Katharina Kormann
describes this in her paper:
https://doi.org/10.4208/cicp.101214.021015a

So in general, I would assume that a change of 1e-8 as you have in your
picture is simply the approximation error that occurs when going from
the initially interpolated solution to the new mesh. I would assume that
it goes away if you compute an "approximate" L2 projection where you
compute a right hand side (v, u_init) for test functions v and the
initial field u_init from ExactSolution of step-48 rather than
VectorTools::interpolate.

Have you tried that?

If the change is larger in other cases as you describe, I would be more
worried. Have you tried to investigate more? Does the projection instead
of interpolation help in that case, too?

The implementation should be correct as far as I can tell - we have
verified it on a large series of applications. In my experience, the
most tricky thing to get right regarding constraints and boundary
conditions is for nonlinear problems or implicit solvers with
inhomogeneous boundary data, but that appears to not be the case here...

Best,
Martin

[Stephen DeWitt]

Dear Martin,

Thank you for your quick reply.

I tried as you suggested, and set the initial condition by computing a right hand side, and it lowered the difference between the initial condition and the first time step by about a factor of 5. So it helped, but didn't make a qualitative change.

While doing that test, I realized that I was wrong in my first post where I said that the artifacts are only present between the initial condition and the first time step -- I missed that the field is only outputted every 200 time steps. Changing it so that it outputs every time step, it turns out that the solution changes by up to 1e-8 between time steps for the first 40 time steps or so before settling down.

Assuming that this is just intrinsic error in the approximation of the mass matrix, given how you said that the implementation has been repeatedly verified, it must just be that I'm being too aggressive in choosing where to coarsen my mesh. 

Thank you for sending me that paper -- I'll read through it to better understand the approximations that go into the diagonal mass matrix.

Thanks again,

Steve