/**
#  Advection scheme mess up with tensors

*SP: This is not a bug. The default boundary conditions for scalar and
*tensor fields are different. They must be different because scalar,
*vector and tensor fields do not have the same symmetries. Setting the
*boundary conditions for the tensor components identically, see the
*code in main() below, gives the same results.*

Given the tracer fields $A$ and $B$ related by the function $B = \log
A$ ,for example.  Observe that if $A$ obeys to the advection equation
$B$ will do same,

$$
\quad \partial_t A + \mathbf{u} \cdot \nabla A = 0 \quad \Rightarrow \quad \partial_t B + \mathbf{u} \cdot \nabla B = 0 
$$

Therefore, the relationship $B = \log A$ should hold in time.

We will check that the advection scheme do not degenerate
the functional relationship as time goes by.  $A$ and/or $B$ 
could be scalars and/or components of a tensor. */

#define DT_MAX (0.001)
#define MU0 10. //  Viscosity 
#define uwall(x,t) (8*(1+tanh(8*(t-1/2)))*sq(x)*sq(1-x))

#include "navier-stokes/centered.h"


#include "tracer.h"
scalar A[], B[]; 
tensor T[];
scalar * tracers = {A, T.x.x, B};

int main()
{

  /**
  *SP: here we set the boundary conditions for T to that of a
  symmetric scalar field. We must use a temporary vector v because
  qcc is not clever enough to recognise T.x.t.* */
  
  vector v = T.x;
  v.t[top] = neumann(0);
  v.t[bottom] = neumann(0);
  v.t[left] = neumann(0);
  v.t[right] = neumann(0);
  
  DT = DT_MAX;
  N = 64;
  init_grid (N);
  const face vector mus[] = {MU0,MU0};
  mu = mus;
  run();
}

u.t[top] = dirichlet(uwall(x,t));
u.n[top] = dirichlet(0);
u.t[bottom] = dirichlet(0);
u.n[bottom] = dirichlet(0);
u.t[left] = dirichlet(0);
u.n[left] = dirichlet(0);
u.t[right] = dirichlet(0);
u.n[right] = dirichlet(0);

event init (i = 0) 
{
  foreach() {
    u.x[] = 0.;
    A[] = (4+3*sin(2*x)*cos(2*y));
    B[] = log (A[]);
    T.x.x[] = log (A[]);
  }
  boundary ((scalar *){u, A, B, T.x.x});
}

static double energy()
{
  double se = 0.;
  if (u.x.face)
    foreach(reduction(+:se))
      se += (sq(u.x[] + u.x[1,0]) + sq(u.y[] + u.y[0,1]))/8.*sq(Delta);
  else // centered
    foreach(reduction(+:se))
      se += (sq(u.x[]) + sq(u.y[]))/2.*sq(Delta);
  return se;
}

event kinetic_energy (t += 0.01) 
{
  //  static FILE * fp = fopen ("kinetic", "w");
  fprintf (stderr, "%g %g\n", t, energy());
  //fflush(fp);
}

event profile (t = 1.8)
{

  FILE * fpp = fopen("yprof", "w");
  for (double y = 0; y <= 1.; y += 0.01)
    fprintf (fpp, "%g %g %g %g\n", y, 
	     log(interpolate (A, 0.5, y)), 
	     interpolate (B, 0.5, y), 
	     interpolate (T.x.x, 0.5, y)); 
  fclose (fpp);
  
  fpp = fopen("xprof", "w");
  for (double x = 0; x <= 1; x += 0.01)
    fprintf (fpp, "%g %g %g %g \n", x,
	     log(interpolate (A, x, 0.75)), 
	     interpolate (B, x, 0.75), 
  	     interpolate (T.x.x, x, 0.75));
  fclose (fpp);
}

/**

The scalars are correctly advected. The relationship $B = \log A$ holds at
instant $t=1.8$ . However, if the scalar is
a component of a tensor (the $xx$ component), the
relationship does not stand up, i.e $\mathbf{T}_{xx} \neq \log A$. It is
even worst if the component were $xy$-component.

~~~gnuplot x profile at time t=1.8 
plot 'xprof' u 1:2 w l t 'reference scalar', 'xprof' u 1:3 t 'scalar', 'xprof' u 1:4 t 'tensor'
~~~
*/