/**
# Electrostatic in planar layers

This test reproduces partially the table of convergence of the electric
field for different electrostatic configurations shown in
[Lopez-Herrera et al,
(2011)](/src/references.bib#lopez-herrera2011). The configurations
tested are: (a) both layers are dielectric, (b) both layers are
conducting and (c) one layer conducting and the other dielectric.

We use the implicit electrodynamic solver. The explicit solver
*ehd.h*, which is fully conservative, gives similar results. We also
need the *run()* loop. */

#include "ehd/implicit.h"
#include "run.h"

/**
Variable $ic$ controls the type of electrical configuration tested.
The analytical values used for comparison are stored in E1 and E2.
Finally, *tfinal* controls the instant at which results are
written. Note that for configuration (a) to reach a stationary state
requires just one iteration while for the other configurations about
15 are required. */

int ic;
double tfinal, E1, E2;

phi[top]    = dirichlet(0);
phi[bottom] = dirichlet(1);

#define beta 3.
#define eta 2.

int main() { 
  X0 = Y0 = -0.5;
  DT = 1;
  TOLERANCE = 1e-5;
  int LEVEL;
  for (ic = 0; ic < 3; ic++) {
    switch (ic) {
    case 0: // both layers are dielectric;
      fprintf(stderr, "dielectric-dielectric\n");
      E1 = 2/(1+beta); 
      E2 = 2*beta/(1+beta);
      tfinal = 1;
      break;
    case 1: // both layers are conducting;
      fprintf(stderr, "conducting-conducting\n");
      E1 = 2/(1+eta); 
      E2 = 2*eta/(1+eta);
      tfinal = 15;
      break;
    case 2: //bottom layer conducting upper one dielectric;
      fprintf(stderr, "conducting-dielectric\n");
      E1 = 0.; 
      E2 = 2.;
      tfinal = 15;
      break;
    }
    for (LEVEL = 5; LEVEL <= 7; LEVEL++) {
      N = 1 << LEVEL;
      run(); 
    }
  }
}

event init (t = 0) {
  scalar f[];
  foreach() {
    f[] = (y > 0. ? 0. : 1.0);
    rhoe[] = 0.;
  }

  foreach_face(){
    double T = (f[]+f[-1,0])/2.;

    /**
     If in configuration (a) we use an harmonic interpolation for the
     permittivity the exact values of the electric field are
     recovered. */

    epsilon.x[] = 1./(T/beta + (1. - T)); //Exact for the diel-diel
    // epsilon.x[] = T*beta + (1. - T); // Non exact 
    K.x[] = (ic == 0 ? 0. : (ic == 1 ? T*eta+(1. - T) : eta*T));
  }
}

event set_dt (i++)
   dt = dtnext (DT);

/**
In the last step the relative error of the electric field, $E=-\nabla \phi$, 
in each medium is reported. */

event result (t = tfinal) {
  scalar E[];
  foreach()
    E[] = (phi[0,-1] - phi[0,1])/(2.*Delta);
  stats s = statsf (E);
  fprintf (stderr, "%d %g %g\n", N, 100*fabs(1.- s.max/E2),
	   E1 > 0. ? 100*fabs(1.- s.min/E1) : 0. );
}