Fast way to compute "extrapolated average" in cell

[vkorchagova]

Hello everyone,

Now I write a type of limiter for DG solver. The main idea in this limiter is to select a troubled cell and their neighbours and, after that, to compare solution average values on selected cell. It is not usual averages: all average values should be computed on the troubled cell; for neighbour cells, we extrapolate solution to the troubled cell and compute an extrapolated value.

To get the extrapolated averages, I rewrote the gridfunc.GetElementAverages function, but this way seems to be too slow and, maybe, difficult. Is there any other simple way to get these values?

[pazner]

Hi @vkorchagova,

It's hard to say why GetElementAverages is too slow for you, but you can compute the cell averages of your extrapolated values by numerically integrating over the relevant cells using an IntegrationRule and ElementTransformation to compute the quadrature weights. You could adapt something like the following:

// To compute cell average of element i
const IntegrationRule &ir = IntRules.Get(mesh.GetElementBaseGeometry(i), order);
ElementTransformation &tr = *mesh.GetElementTransformation(i);
Vector vals;
gf.GetValues(i, ir, vals);
double cell_avg = 0.0;
for (int j=0; j<ir.Size(); ++j)
{
   const IntegrationPoint &ip = ir[j];
   tr.SetIntPoint(&ip);
   cell_avg += vals[j]*ir[j].weight*tr.Weight();
}

[vkorchagova]

Hi @pazner ,
maybe my own modification of GetElementAverages seems to be too slow.
I wrote my own class Averager to save and compute extrapolated average values.
This is the code of two functions which appears as a modification of original GetElementAverages.

void Averager::computeStencilExtrapAverages(
   const ParGridFunction& x,
   const Stencil* stencil,
   ParGridFunction& avgs_local
)
{
   MassIntegrator Mi;
   DenseMatrix loc_mass, loc_mass_shifted;
   Array<int> te_dofs, tr_dofs;
   Vector loc_avgs, loc_this;
   Vector int_psi(avgs_local.Size());

   int_psi = 0.0;

   int iCellTroubled = stencil->cell_num[0];

   for (int i : stencil->cell_num)
   {
      // loc_mass = shape function in gauss points * gauss weights
      Mi.AssembleElementMatrix2(
         *fes->GetFE(i), 
         *avgs_local.FESpace()->GetFE(i),
         *fes->GetElementTransformation(i), 
         loc_mass);
     assembleShiftedElementMatrix(
         *fes->GetFE(i), 
         *fes->GetFE(iCellTroubled),
         *avgs_local.FESpace()->GetFE(i),
         *fes->GetElementTransformation(i),
         *fes->GetElementTransformation(iCellTroubled), 
         loc_mass_shifted);
      fes->GetElementDofs(i, tr_dofs);
      avgs_local.FESpace()->GetElementDofs(i, te_dofs);
      x.GetSubVector(tr_dofs, loc_this);
      loc_avgs.SetSize(te_dofs.Size());
      loc_mass_shifted.Mult(loc_this, loc_avgs);

      avgs_local.SetSubVector(te_dofs, loc_avgs);
      loc_this = 1.0; // assume the local basis for 'this' sums to 1
      loc_mass.Mult(loc_this, loc_avgs);
      int_psi.SetSubVector(te_dofs, loc_avgs);
   }

   for (int i : stencil->cell_num)
   {
      avgs_local(i) /= int_psi(i);
   }
}

void Averager::assembleShiftedElementMatrix(
   const FiniteElement &trial_fe, 
   const FiniteElement &troubled_fe,
   const FiniteElement &test_fe,
   ElementTransformation &Trans, 
   ElementTransformation &TransTroubled,
   DenseMatrix &elmat
)
{
   int nDofsTrial = trial_fe.GetDof();
   int nDofsTroubled = troubled_fe.GetDof();
   int nDofsTest = test_fe.GetDof();
   double w;
   Vector trial_shape (nDofsTroubled);
   Vector te_shape (nDofsTest);
   Vector ipPhys;
   IntegrationPoint ipRef;
   elmat.SetSize(nDofsTest, nDofsTroubled);
   const int order = troubled_fe.GetOrder() + test_fe.GetOrder() + Trans.OrderW();
   const IntegrationRule *ir = &IntRules.Get(troubled_fe.GetGeomType(), order);
   elmat = 0.0;
   for (int i = 0; i < ir->GetNPoints(); i++)
   {
      const IntegrationPoint &ip = ir->IntPoint(i);
      TransTroubled.Transform(ir->IntPoint(i), ipPhys);
      Trans.TransformBack(ipPhys,ipRef);
      trial_fe.CalcShape(ipRef, trial_shape);  // element - shifted, gauss points - troubled
      test_fe.CalcShape(ipRef, te_shape);
      w = Trans.Weight() * ip.weight;
      te_shape *= w;
      AddMultVWt(te_shape, trial_shape, elmat);
   }
}

[pazner]

I would guess that the matrix assembly is the most expensive part of this algorithm. Could those matrices be pre-computed and stored?

[vkorchagova]

Sure. But I afraid that matrices could require a big part of memory. If we consider 2D rectangular mesh with only 10000 inner cells, we have 4 neighbours for every cell. Therefore, we must store 4 matrices for every cell -- totally 40000 matrices. For 3D mesh, the number of matrices grows up extremely due to the 3rd dimension.

[pazner]

Yes, that is true, but keep in mind that this is same amount of storage as e.g. the system matrix for an advection operator, so it is not necessarily unreasonable.

You might also be able to precompute some parts of the matrices that don't change from element to element to reduce the time.