src/test/cyl_planar.c

Charge relaxation in a planar cross-section

This is the same problem as in cyl_axi.c but in a planar cross-section of the column.

#define NAVIER 1
int LEVEL;

#include "grid/multigrid.h"
#include "ehd/implicit.h"
#if NAVIER
# include "navier-stokes/centered.h"
# include "ehd/stress.h"
#endif
#include "fractions.h"

Far away from the conducting column the electric potential is zero.

φ[bottom] = dirichlet(0);
φ[top]    = dirichlet(0);
φ[right]  = dirichlet(0);
φ[left]   = dirichlet(0);
#if NAVIER
p[top]      = dirichlet(0);
u.n[top]    = neumann(0);
p[bottom]   = dirichlet(0);
u.n[bottom] = neumann(0);
p[left]     = dirichlet(0);
u.n[left]   = neumann(0);
p[right]    = dirichlet(0);
u.n[right]  = neumann(0);
#endif

#define β 3.
#define cond 3.
#define rhoini 0.5
#define R 0.1
#define circle(x,y) (sq(R) - sq(x) - sq(y))

scalar f[];

event init (t = 0) {
  face vector s[];
  vertex scalar ψ[];
  foreach_vertex ()
    ψ[] = circle(x,y);
  fractions (ψ, f, s);
  foreach()
    rhoe[] = rhoini*f[];
  boundary ({rhoe});

  foreach_face() {
    double ff = (f[] + f[-1])/2.;
    ε.x[] = (ff*β + (1. - ff))*fm.x[];
    K.x[] = cond*s.x[]*fm.x[];
  }
  boundary ((scalar *){ε, K});
}

event chargesum (i++) {
  double Q = statsf(rhoe).sum;
  static double Q0;
  if (i == 0)
    Q0 = Q;
  else
    assert (fabs(Q - Q0) < 1e-7);
}

event epfield (t = 10) {
  char name[80];
  sprintf (name, "Er-%d", LEVEL);
  FILE * fp = fopen (name, "w");

  scalar ee[];
  foreach() {
    double Ex = (φ[-1,0] - φ[1,0])/(2*Δ);
    double Ey = (φ[0,-1] - φ[0,1])/(2*Δ);
    double r = sqrt(x*x + y*y);
    double En = sqrt(Ex*Ex + Ey*Ey);
    ee[] = En - (r < R ? 0 : 0.5*sq(R)*rhoini/r);
#if NAVIER
    fprintf (fp, "%g %g %g %g\n", r, En, p[], rhoe[]);
#else
    fprintf (fp, "%g %g %g\n", r, En, rhoe[]);
#endif
  }
  fclose (fp);
  
  norm n = normf (ee);
  fprintf (stderr, "%d %g %g %g\n", LEVEL, n.avg, n.rms, n.max);
}

#if TREE // fixme: this does not work
event adapt (i++) {
  double Qb = statsf(rhoe).sum;
  adapt_wavelet ({f}, (double []){1e-6},
		 maxlevel = LEVEL + 1, minlevel = LEVEL);
  double Qa = statsf(rhoe).sum;
  assert (fabs(Qa - Qb) < 1e-10);
}
#endif

int main() {

The computational domain spans [-1:1][-1:1].

  X0 = Y0 = -1.;
  L0 = 2.;
  DT = 1;
  TOLERANCE = 1e-7;

We compute the solution for different levels of refinement.

  for (LEVEL = 6; LEVEL <= 8; LEVEL++) {
    N = 1 << LEVEL;
    run();
  }
}

Results

Radial electric field distribution as a function of grid refinement.

Radial electric field distribution as a function of grid refinement.

Pressure distribution as a function of grid refinement.

Pressure distribution as a function of grid refinement.

Convergence of error on the norm of the electric field

Convergence of error on the norm of the electric field

See also