Post-processing velocities with Trilinos vectors during simulation

[Audrey Collard-Daigneault]

Hi everyone!

I'm working on post-processing velocities with Trilinos solution vectors during the simulation on Lethe.
Calculating average velocities and pressures (<u>, <v>, <w> and <p>) works well using Trilinos vectors with no ghost cells and .equ(...) function.
The calculation of average Reynolds stress (<u'u'>, <v'v'> and <w'w'> and average shear stress (<u'v'>), where u' = u - <u>, is quite not easy.
It seems that I can't do what I want in parallel without working with loop.
The problem is that loops seem to take way too much time on Trilion vectors and the local_range() function doesn't work with BlockVector. (Am I wrong?)
I have tried some ways to do it.
  1. Trying working with Deal.ii vectors : did not work in parallel
  2. Doing summations and/or multiplications of the vectors local_evaluation_point to get u'u', v'v' and w'w' and replacing the fourth element of the evaluation point with a loop for u'v' : too much time.
  3. Doing loops on the Trilinos vectors : too much time

Does anyone can give me an advice to find a way to get one or two vectors with my precious time-averaged Reynolds stress and time-averaged shear stress?
I did not mention the exact way I'm trying to calculate the average of all of that, but I'm posting my code where I'm trying to do a loop on the solution vector.
The variables average_velocities, total_time and inv_range_time are calculated in an other function is my class AverageVelocities.

template <int dim, typename VectorType, typename DofsType>
VectorType
AverageVelocities<dim, VectorType, DofsType>::calculate_reynolds_stress(
  const VectorType &                        local_evaluation_point,
  const std::shared_ptr<SimulationControl> &simulation_control,
  const DofsType &                          locally_owned_dofs,
  const MPI_Comm &                          mpi_communicator)
{
  if (simulation_control->get_step_number() == 0)
  {
    // Reinitializing vectors with zeros at t = 0
    sum_reynolds_stress_dt.reinit(locally_owned_dofs,
                                  mpi_communicator);
    reynolds_stress.reinit(locally_owned_dofs,
                           mpi_communicator);
  }
  else if (abs(total_time) < 1e-6 || total_time > 1e-6)
  {
    VectorType reynolds_stress_dt(locally_owned_dofs,
                                  mpi_communicator);

    // ***Won't work with BlockVectors
    if constexpr (std::is_same_v<VectorType, TrilinosWrappers::MPI::Vector>)
    {
      for (unsigned int i = local_evaluation_point.local_range().first;
           i < local_evaluation_point.local_range().second; i++)
      {
        if ((i + 4) % 4 == 0)
        {
          // Calculating (u'u')*dt, (v'v')*dt (w'w')*dt and (u'v')*dt
          reynolds_stress_dt[i] =
            (local_evaluation_point[i] - average_velocities[i]) *
            (local_evaluation_point[i] - average_velocities[i]) * dt;
          reynolds_stress_dt[i + 1] =
            (local_evaluation_point[i + 1] - average_velocities[i + 1]) *
            (local_evaluation_point[i + 1] - average_velocities[i + 1]) * dt;
          reynolds_stress_dt[i + 2] =
            (local_evaluation_point[i + 2] - average_velocities[i + 2]) *
            (local_evaluation_point[i + 2] - average_velocities[i + 2]) * dt;
          reynolds_stress_dt[i + 3] =
            (local_evaluation_point[i] - average_velocities[i]) *
            (local_evaluation_point[i + 1] - average_velocities[i + 1]) * dt;

          // Summation of all reynolds stress during simulation
          sum_reynolds_stress_dt += reynolds_stress_dt;

          // Calculating time-averaged reynolds stress if output needed
          if (simulation_control->is_output_iteration())
            reynolds_stress.equ(inv_range_time, sum_reynolds_stress_dt);
        }
      }
    }
  }
  return reynolds_stress;
}


Thank you in advance!

[Bruno Turcksin]

Hi,

On Wednesday, November 4, 2020 at 2:34:05 PM UTC-5 acdaig...@gmail.com wrote:

Hi everyone!

I'm working on post-processing velocities with Trilinos solution vectors during the simulation on Lethe.
Calculating average velocities and pressures (<u>, <v>, <w> and <p>) works well using Trilinos vectors with no ghost cells and .equ(...) function.
The calculation of average Reynolds stress (<u'u'>, <v'v'> and <w'w'> and average shear stress (<u'v'>), where u' = u - <u>, is quite not easy.
It seems that I can't do what I want in parallel without working with loop.

The problem is that loops seem to take way too much time on Trilinos vectors and the local_range() function doesn't work with BlockVector. (Am I wrong?)

True local_range() does not work on a BlockVector but you can get  the local_range of each block in the block vector. Just do my_block_vector.block(i).local_range()

I have tried some ways to do it.
  1. Trying working with Deal.ii vectors : did not work in parallel

What does that mean? Did you use the dealii::Vector or the dealii::LinearAlgebra::distributed::Vector (https://dealii.org/developer/doxygen/deal.II/classLinearAlgebra_1_1distributed_1_1Vector.html) Only the second one supports MPI

  2. Doing summations and/or multiplications of the vectors local_evaluation_point to get u'u', v'v' and w'w' and replacing the fourth element of the evaluation point with a loop for u'v' : too much time.
  3. Doing loops on the Trilinos vectors : too much time

What does that mean it takes too much time? Did you profile the code? What's the bottleneck?

Your for loop is a little bit strange. I would think that the compiler can optimize this code but it would be easier for the compiler to optimize the following code

 

const unsigned int begin_index = local_evaluation_point.local_range().first; // You do this because otherwise the compiler will call local_range() at each evaluation of the for loop condition

const unsigned int end_index = local_evaluation_point.local_range().second;

      for (unsigned int i = begin_index;  i < end_index; i += 4)
      {
       // Remove this if ((i + 4) % 4 == 0)
     

      }

Best,

Bruno

[Audrey Collard-Daigneault]

Hi Bruno,

Thanks for your help!

The problem was actually my "sum_reynolds_stress_dt += reynolds_stress_dt;" line. It needed to be outside of the for loop to be effective.

I didn't know about LinearAlgebra, but it should propably be useful! Same for the way to work with BlockVectors!

About the "((i + 4) % 4 == 0)'', it doesn't work if I change to "i += 4" because it's not sure that the first vector index on the processor is going to be a factor of 4.

Best,

Audrey