src/test/einstein_sum.c

    Tests for the einstein_sum() macro

    #include "grid/multigrid.h"
#include "navier-stokes/centered.h"

    We create a second-order tensor structure.

    typedef struct{
  coord x;
  coord y;
  coord z;
} coord2;

    We create a third-order tensor structure.

    typedef struct{
  coord2 x;
  coord2 y;
  coord2 z;
} coord3;

int main ()
{
  init_grid (64);

  // define time and length scale
  double gamma = 1 [0,-1];
  double length = 1 [1,0];

  // define arbitrary second and third order tensors
  coord2 E = {{gamma,0.1*gamma}, {0.2*gamma, - gamma}};
  coord2 E2 = {{gamma/length, 0.2* gamma/length},
                {0.7*gamma/length, - 0.3*gamma/length}};
  coord3 K = {E2, E2};
  vector u2[];

    We now initialize a quadratic velocity field around the origin. Formally, we have \displaystyle \textbf{u}(\textbf{x}) = \textbf{E}\cdot \textbf{x}+\textbf{E2}: \textbf{xx} This can be encoded in indices notation as follows:

      foreach(){
    coord X = {x,y,z};
    einstein_sum (i,k,j) {
      // factorized expression
      u.i[] = (E.i.j + K.i.j.k*X.k)*X.j;
      // developed expression
      u2.i[] = E.i.j*X.j + K.i.j.k*X.k*X.j;
    }
  }
  coord2 R =  (coord2) {{0}} ;
  double kinetic_energy = 0; 
  coord2 R2 =  (coord2) {{0}} ;
  double kinetic_energy2 = 0; 
  foreach()
    einstein_sum (i,j) {
      R.i.j += u.i[]*u.j[];
      kinetic_energy += u.i[]*u.i[]/2.;
      R2.i.j += u2.i[]*u2.j[];
      kinetic_energy2 += u2.i[]*u2.i[]/2.;
    }

  fprintf (stderr, "%g %g\n", kinetic_energy, kinetic_energy2);
  
  einstein_sum (i,j) {
    fprintf (stderr, "%g %g\n", R.i.j, R2.i.j);
  }
}