src/test/planar.c

    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). 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. );
}