Implementing smoothness indicators with ProjectGrad

[FH-9]

Hi,

I am trying to implement smoothness indicators of the form

$\beta_i = \sum\limits_{l=1}^k \Delta x_0^{2l-1} \int_{K_0} \left( \frac{1}{l!} \frac{d^l}{dx^l} \tilde{p}_i(x) \right)^2 dx, \quad i = 0, 1, 2 $
This is in the context of calculating non-linear WENO weights for reconstructing the DG polynomial in the troubled (the one with sharp gradients) element. $i$ are the indices of the troubled element and its neighbors, $K_0$ is the troubled element, $k$ is the order of the polynomial, $l$ is the index for gradient level, and $\tilde{p}_i(x)$ are the polynomials of the elements.

As you can see, I need to repeatedly calculate the derivatives of the polynomials and integrate in the target element. Following the methodology described in #1118, I wrote a piece of code that is supposed to perform these operations but my solution blows up and the implementation of smoothness indicators is my primary suspect. Here is the relevant part of the code:

        vfes->GetElementVDofs(indx, vdof_indices);
        sol.GetSubVector(vdof_indices, el_vdofs);
        fe->ProjectGrad(*fe, *Tr, grad_matrix);
         .
         .
         .
        DenseMatrix vdof_mat(el_vdofs.GetData(), ndofs, num_equations);
        if (dim == 1)
        {
          grad_vdof_mat = vdof_mat;
          for (int grad_level = 1; grad_level < order + 1; grad_level++)
          {
            int_order = 2 * (fe->GetOrder() - grad_level) + Tr->OrderW();
            const IntegrationRule* ir = &IntRules.Get(fe->GetGeomType(), int_order);
            mfem::Mult(grad_matrix, grad_vdof_mat, tmp_mat);
            grad_vdof_mat = tmp_mat;
            integrand = 0.0;
            for (int k = 0; k < ir->GetNPoints(); k++)
            {
              const IntegrationPoint& ip = ir->IntPoint(k);
              if (j != 0)
              {
                Tr_target->Transform(ip, phys_ip);
                Tr->TransformBack(phys_ip, extrap_ip);
                // Tr->SetIntPoint(&extrap_ip);
                fe->CalcShape(extrap_ip, shape);
              }
              else
              {
                fe->CalcShape(ip, shape);
              }
              shape *= ip.weight;
              shape *= Tr->Weight();
              grad_vdof_mat.MultTranspose(shape, grad_state);
              grad_state *= grad_state;
              integrand += grad_state;
            }
            integrand /= pow(factorial(grad_level), 2); 
            integrand *= pow(target_el_size, 2 * grad_level - 1);
            beta_col += integrand;
          }
        }

where beta_col are the values of $b_i$ for a given element and for different equations in the system. I am using the output of ProjectGrad extensively to first get the gradient degrees of freedom and then multiply by the output of CalcShape to interpolate at the quadrature points.
Does anything look wrong or stand out to you in this implementation? I feel like this is the where the problem lies.

Thank you very much.

Farhad

[FH-9]

I realized the integral computation was way off. I think after projecting the gradient degrees of freedom back onto the finite dimensional space of the solution, the integral, $\int \left( \frac{d}{dx} \tilde{p} \right)^2 dx \approx \int \left(\sum_k \hat{dp}_k \phi_k \right)^2 dx$, can be cast into the the quadratic form $\hat{x}^T M \hat{x}$, where $\hat{x}$ is the vector of gradient degrees of freedom, and $M$ is the extrapolated mass integrator (since the integral is over the target element only). I have updated the code but it is still blowing up. Is my reasoning here correct? And do you notice anything wrong with the code?
Here is the updated code:

vfes->GetElementVDofs(indx, vdof_indices);
sol.GetSubVector(vdof_indices, el_vdofs);
fe->ProjectGrad(*fe, *Tr, grad_matrix);
DenseMatrix vdof_mat(el_vdofs.GetData(), ndofs, num_equations);
// beta is num_equations x num_stencil array
// columns of which are the beta values for the given element in the stencil
beta.GetColumnReference(j, beta_col); 
 .
 .
 .
Vector grad_vec(ndofs), Mdx(ndofs);
DenseMatrix mass_matrix(ndofs);
if (dim == 1)
{
  grad_vdof_mat = vdof_mat;
  for (int grad_level = 1; grad_level < order + 1; grad_level++)
  {
    int_order = 2 * fe->GetOrder() + Tr->OrderW();
    const IntegrationRule* ir = &IntRules.Get(fe->GetGeomType(), int_order);
    mfem::Mult(grad_matrix, grad_vdof_mat, tmp_mat);
    grad_vdof_mat = tmp_mat;
    integrand = 0.0;
    mass_matrix = 0.0;
// Calculating the (extrapolated) mass integrator
    for (int k = 0; k < ir->GetNPoints(); k++)
    {
      const IntegrationPoint &ip = ir->IntPoint(k);
      // j = 0 is the index of the target/troubled element; this is where the integral is evaluated
      if (j != 0)
      {
        Tr_target->Transform(ip, phys_ip);
        Tr->TransformBack(phys_ip, extrap_ip);
        Tr->SetIntPoint(&extrap_ip);
        fe->CalcPhysShape(*Tr, shape);
      }
      else
      {
        Tr->SetIntPoint(&ip);
        fe->CalcPhysShape(*Tr, shape);
      }
      w = Tr->Weight() * ip.weight;
      AddMult_a_VVt(w, shape, mass_matrix);
    }
// x^T M x quadratic form calculation where x are the gradient degrees of freedom
    for (int eq = 0 ; eq < num_equations; eq++)
    {
      grad_vdof_mat.GetColumn(eq, grad_vect);
      mass_matrix.Mult(grad_vect, Mdx);
      integrand(eq) = grad_vect * Mdx;
    }
    integrand /= pow(factorial(grad_level), 2); 
    integrand *= pow(target_el_size, 2 * grad_level - 1);
    beta_col += integrand;
  }
}

Thank you very much.

[FH-9]

I would appreciate if anyone could offer some insight. Thank you.