Getting old solution values at quadrature points using MeshWorker::loop

[Andrew Davis]

I'm trying to implement a time-dependent solver that assembles the system using the MeshWorker::loop tool. However, I cannot figure out how get the values of the solution at the previous timestep at each quadrature point.

Currently, I store the old solution in a class called "AdvectionProblem" such that---I've left out some code in the interest of brevity.

template<unsigned int dim>

class AdvectionProblem {

public:

...

void AssembleSystem()

private:

...

Vector<double> old_solution;

DoFHandler<dim> dofHandler;

};

Since we need the old_solution values to assemble the system, I create a CellIntegrator class that MeshWorker::loop can use to integrate over each cell (I also create similar objects for faces/boundaries):

template<unsigned int dim>
class IntegrateCell {
public:
  IntegrateCell(dealii::Vector<double> const& old_solution);

  virtual ~IntegrateCell() = default;

  void operator()(dealii::MeshWorker::DoFInfo<dim>& dinfo, dealii::MeshWorker::IntegrationInfo<dim>& info) const;

private:

  /// Store a constant reference to the old solution
  const dealii::Vector<double>& old_solution;
};

template<unsigned int dim>
IntegrateCell<dim>::IntegrateCell(Vector<double> const& old_solution) : old_solution(old_solution) {}

template<unsigned int dim>
void IntegrateCell<dim>::operator()(MeshWorker::DoFInfo<dim>& dinfo, MeshWorker::IntegrationInfo<dim>& info) const {
  // get some useful objects
  const FEValuesBase<dim>& feValues = info.fe_values();
  FullMatrix<double>& localMatrix = dinfo.matrix(0).matrix;
  Vector<double>& localVector = dinfo.vector(0).block(0);
  const std::vector<double>& JxW = feValues.get_JxW_values();

  // get the values of the old solution
  std::vector<double> old_solution_values(feValues.n_quadrature_points);
  feValues.get_function_values(old_solution, old_solution_values);

  // these are non-zero, IngtegrateCell is correctly getting the old_solution
  std::cout << old_solution << std::endl;

  // loop over the quadrature points and compute the advection vector at the current point
  for( unsigned int point=0; point<feValues.n_quadrature_points; ++point ) {

    // these are always all zero

    std::cout << "value: " << old_solution_values[point] << std::endl;
  }
}

I then try to assemble the system using the MeshWorker::loop tool:

template<unsigned int dim>
void AdvectionProblem<dim>::AssembleSystem() {
  // an object that knows about data structures and local integration
  MeshWorker::IntegrationInfoBox<dim> infoBox;

  // initialize the quadrature formula
  const unsigned int nGauss = dofHandler.get_fe().degree + 1;
  infoBox.initialize_gauss_quadrature(nGauss, nGauss, nGauss);

  // these are the values we need to integrate our system
  infoBox.initialize_update_flags();
  const UpdateFlags updateFlags = update_quadrature_points | update_values | update_gradients | update_JxW_values;
  infoBox.add_update_flags(updateFlags, true, true, true, true);

  // initialize the finite element values
  infoBox.initialize(fe, mapping);

  // create an object that forwards local data to the assembler
  MeshWorker::DoFInfo<dim> dofInfo(dofHandler);

  // create teh assembler---tell it where to put local data
  MeshWorker::Assembler::SystemSimple<SparseMatrix<double>, Vector<double> > assembler;
  assembler.initialize(systemMatrix, rhs);

  IntegrateCell<dim> cellIntegration(old_solution);

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

For some reason the feValues.get_function_values(old_solution, old_solution_values); call in the CellIntegrator function always sets the old_solution_values to zero. Does anyone see what I'm doing wrong?

Thanks!

[Andrew Davis]

I'm not sure if this is helpful but as a diagnostic, replacing MeshWorker::loop with this loop that I wrote myself:

const QGauss<dim>  quadrature_formula(fe.degree + 1);
  const unsigned int n_q_points = quadrature_formula.size();
  FEValues<dim> fe_values(fe, quadrature_formula, update_values | update_JxW_values);
  std::vector<double> solution_values(n_q_points);
  for( const auto& cell : dofHandler.active_cell_iterators() ) {
    fe_values.reinit(cell);
    fe_values.get_function_values(old_solution, solution_values);
    for( unsigned int i=0; i<n_q_points; ++i ) {

      // these values are non-zero, as expected

      std::cout << "self loop value: " << solution_values[i] << std::endl;
    }
  }

Correctly prints the values that I expect.

[Wolfgang Bangerth]

On 1/29/20 1:14 PM, Andrew Davis wrote:
> *
> *
> *For some reason the feValues.get_function_values(old_solution,

> old_solution_values); call in the CellIntegrator function always sets the

> old_solution_values to zero. Does anyone see what I'm doing wrong?*

So the output of this function is wrong. Is the input correct? My guess would
be that if you looked at 'old_solution', you'd find that it's a zero vector.
Why that is so I don't know, but it's something you should be able to trace
back to its source.

Best
W.

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

[Andrew Davis]

I thought the same thing---that 'old_solution' would be a zero vector---but when I add the line `std::cout << old_solution << std::endl;` in the function IntegrateCell<dim>::operator() (just after the `get_function_values' call in the original post) prints a non-zero vector that I expect.

[Wolfgang Bangerth]

On 1/30/20 8:45 AM, Andrew Davis wrote:
> I thought the same thing---that 'old_solution' would be a zero vector---but
> when I add the line `std::cout << old_solution << std::endl;` in the function
> IntegrateCell<dim>::operator() (just after the `get_function_values' call in
> the original post) prints a non-zero vector that I expect.

That does not sound credible to me. Can you construct a small example that
does exactly that? I.e., print out, for example, the l2-norm of old_solution
right before the call to fe_values.get_function_values(), and the l2-norm of
old_solution_values right after.

[Andrew Davis]

Sure, attached is a *.cpp file that should run. The cout statement on line 87 prints non-zero values for `old_solution` but the cout statement online 92 always prints zero.

[Andrew Davis]

As another diagnostic---I added cout statements in the function `get_function_values` in fe_values.cc.  It looks like the values of the input variable `fe_function` are correct but the `dof_values` are always set to zero.  Specifically:

  // get function values of dofs on this cell
  Vector<Number> dof_values(dofs_per_cell);
  present_cell->get_interpolated_dof_values(fe_function, dof_values);
  for( unsigned int i=0; i<dof_values.size(); ++i ) {
    std::cout << "dof_val: " << dof_values[i] << std::endl;
  }

always prints zero.

[Andrew Davis]

I'm continuing to dig into this problem. For some reason the `get_function_values` is calling the function

for (VEC : VECTOR_TYPES)
  {
    template <int dim, int spacedim>
    void
    FEValuesBase<dim, spacedim>::TriaCellIterator::get_interpolated_dof_values(
      const VEC &, Vector<VEC::value_type> &) const
    { //
      Assert(false, ExcMessage(message_string));
    \}
  }

from the file fe_values.impl.2.inst.in. The Assert(false) seems strange to me---why doesn't it crash? (I'm in Release mode so that could answer that). However, should it be calling a different function?

Thanks!

[Bruno Turcksin]

On Thursday, January 30, 2020 at 2:28:55 PM UTC-5, Andrew Davis wrote:

from the file fe_values.impl.2.inst.in. The Assert(false) seems strange to me---why doesn't it crash? (I'm in Release mode so that could answer that). However, should it be calling a different function?

You should always debug your code in Debug mode. It's very possible that there is an assert that will catch a problem earlier in the code and so the code after that assert makes no sense because you are not supposed to get there.

Best,

Bruno

[Andrew Davis]

Fair point. I reinstalled in Debug mode and now I'm failing at that Assert.  Here is the error message (using the same code I attached previously):

An error occurred in line <2891> of file <dealii/source/fe/fe_values.cc> in function
    dealii::types::global_dof_index dealii::FEValuesBase<dim, spacedim>::TriaCellIterator::n_dofs_for_dof_handler() const [with int dim = 2; int spacedim = 2; dealii::types::global_dof_index = unsigned int]
The violated condition was:
    false
Additional information:
    You have previously called the FEValues::reinit function with a
cell iterator of type Triangulation<dim,spacedim>::cell_iterator. However,
when you do this, you cannot call some functions in the FEValues
class, such as the get_function_values/gradients/hessians/third_derivatives
functions. If you need these functions, then you need to call
FEValues::reinit with an iterator type that allows to extract
degrees of freedom, such as DoFHandler<dim,spacedim>::cell_iterator.

I'm going to dig into this further but does anyone know if this is a bug or user error? If it is user error does anyone know how to fix this?

Cheers,

Andy

[Andrew Davis]

I found a fix for this by editing the deal.ii source code. I am very new to deal.ii so it is extremely likely that there is a way to accomplish what I want without touching the deal.ii source or there might be a better way of editing the source. So, here is the fix that worked for me. Let me know if there is a better solution.

I first implemented a new reinit function for IntegrationInfo:

    /**
     * Reinitialize internal data structures for use on a cell.
     */
    template <typename number, class ITERATOR>
    void
    reinit(ITERATOR cell, const DoFInfo<dim, spacedim, number> &info);

  template <int dim, int spacedim>

  template <typename number, class ITERATOR>
  inline void
  IntegrationInfo<dim, spacedim>::reinit(ITERATOR cell,
    const DoFInfo<dim, spacedim, number> &info)
  {
    for (unsigned int i = 0; i < fevalv.size(); ++i)
      {
        FEValuesBase<dim, spacedim> &febase = *fevalv[i];
        if (info.sub_number != numbers::invalid_unsigned_int)
          {
            // This is a subface
            FESubfaceValues<dim, spacedim> &fe =
              dynamic_cast<FESubfaceValues<dim, spacedim> &>(febase);
            fe.reinit(cell, info.face_number, info.sub_number);
          }
        else if (info.face_number != numbers::invalid_unsigned_int)
          {
            // This is a face
            FEFaceValues<dim, spacedim> &fe =
              dynamic_cast<FEFaceValues<dim, spacedim> &>(febase);
            fe.reinit(cell, info.face_number);
          }
        else
          {
            // This is a cell
            FEValues<dim, spacedim> &fe =
              dynamic_cast<FEValues<dim, spacedim> &>(febase);
            fe.reinit(cell);
          }
      }

    const bool split_fevalues = info.block_info != nullptr;
    if (!global_data->empty())
      fill_local_data(info, split_fevalues);
  }

And then I changed the MeshWorker::loop function on lines 231-232 of deal.ii/meshworker/loop.h to

if (integrate_cell)
      info.cell.reinit(cell, dof_info.cell);

Cheers,

Andy

[Wolfgang Bangerth]

On 1/30/20 1:39 PM, Andrew Davis wrote:
> Fair point. I reinstalled in Debug mode and now I'm failing at that Assert.

So learning opportunity: You could have saved yourself a lot of searching of
where the problem actually is by always first running your program in debug
mode :-)

> Here is the error message (using the same code I attached previously):
>
> An error occurred in line <2891> of file <dealii/source/fe/fe_values.cc> in
> function
>     dealii::types::global_dof_index dealii::FEValuesBase<dim,
> spacedim>::TriaCellIterator::n_dofs_for_dof_handler() const [with int dim = 2;
> int spacedim = 2; dealii::types::global_dof_index = unsigned int]
> The violated condition was:
>     false
> Additional information:
>     You have previously called the FEValues::reinit function with a
> cell iterator of type Triangulation<dim,spacedim>::cell_iterator. However,
> when you do this, you cannot call some functions in the FEValues
> class, such as the get_function_values/gradients/hessians/third_derivatives
> functions. If you need these functions, then you need to call
> FEValues::reinit with an iterator type that allows to extract
> degrees of freedom, such as DoFHandler<dim,spacedim>::cell_iterator.
>
> I'm going to dig into this further but does anyone know if this is a bug or
> user error? If it is user error does anyone know how to fix this?

Well, at least it explains what the issue is: You can't evaluate a finite
element field if you don't have a DoFHandler. The question is why the
FEValues::reinit() function was called with only a
Triangulation::active_cell_iterator and not a
DoFHandler::active_cell_iterator. I have to admit that, as with so many things
that relate to MeshWorker, I don't know and nobody might actually know. The
place you should find by going up the call stack in a debugger is where
FEValues::reinit() is called, and what the iterator object there is.

[Wolfgang Bangerth]

On 1/30/20 2:48 PM, Andrew Davis wrote:
> I found a fix for this by editing the deal.ii source code. I am very new to
> deal.ii so it is extremely likely that there is a way to accomplish what I
> want without touching the deal.ii source or there might be a better way of
> editing the source. So, here is the fix that worked for me. Let me know if
> there is a better solution.

Interesting solution. I don't know if anyone really still understands the
MeshWorker code (as has been mentioned many times on the mailing list) -- I
think we are all confused by the differences between DoFInfo, DoFInfoBox,
IntegrationInfo, IntegrationInfoBox, or generally the term "Box" in this
combination to begin with :-) All of this predates the era when we did patch
reviews...

I think this could work. Have you looked whether the existing function is
called anywhere else? Why not just change the existing function to take the
"real" cell object as an additional argument? If there are (were) only a
handful of places where that function is called, maybe it's worth changing
this everywhere.

[Andrew Davis]

I have not looked deeply into whether this function is called in a lot of places. Attached is an implementation of a time-dependent advection equation with an adaptive mesh. It seems like it is working given the changes that I made but it does not work given the current master branch.

I can clean up my "fix" and submit a pull request, however, I figure I'd through out the example code here before putting to much work into that. I'm not to deal.ii and, therefore, there is a good chance that someone will come up with a better solution (or find a bug in my code).

[Wolfgang Bangerth]

On 2/3/20 12:16 PM, Andrew Davis wrote:
> I have not looked deeply into whether this function is called in a lot
> of places. Attached is an implementation of a time-dependent advection
> equation with an adaptive mesh. It seems like it is working given the
> changes that I made but it does not work given the current master branch.
>
> I can clean up my "fix" and submit a pull request, however, I figure I'd
> through out the example code here before putting to much work into that.
> I'm not to deal.ii and, therefore, there is a good chance that someone
> will come up with a better solution (or find a bug in my code).

I think it would be very nice if you could submit a patch! My preference
would be that instead of adding a function, you would just change the
existing one and see where the other places are that call it. In
essence, replace the existing function by yours. Would you be interested
in giving this a try, and then creating a pull request for it? (I'd also
be happy to do that last step for you, if you don't want to deal with
git/github.)