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.

    phi[bottom] = dirichlet(0);
    phi[top]    = dirichlet(0);
    phi[right]  = dirichlet(0);
    phi[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 beta 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 psi[];
      foreach_vertex ()
        psi[] = circle(x,y);
      boundary ({psi});
      fractions (psi, f, s);
      foreach()
        rhoe[] = rhoini*f[];
      boundary ({rhoe});
    
      foreach_face() {
        double ff = (f[] + f[-1])/2.;
        epsilon.x[] = (ff*beta + (1. - ff))*fm.x[];
        K.x[] = cond*s.x[]*fm.x[];
      }
      boundary ((scalar *){epsilon, K});
    }
    
    event chargesum (i++) {
      double Q = statsf(rhoe).sum;
    variable ‘Q0’ set but not used [-Wunused-but-set-variable]
      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 = (phi[-1,0] - phi[1,0])/(2*Delta);
        double Ey = (phi[0,-1] - phi[0,1])/(2*Delta);
        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

    set term PNG enhanced font ",10"
    set output 're.png'
    set xlabel 'r'
    set ylabel 'E_r'
    set xrange[0:1]
    set yrange[0.:0.03]
    set sample 1000
    R = 0.1
    rhoini = 0.5
    E(x) = x < R ? 0 : (R*R*rhoini/2/x)
    plot 'Er-6' u 1:2 t "Level = 6",  \
         'Er-7' u 1:2 t "Level = 7",  \
         'Er-8' u 1:2 t "Level = 8",  \
         E(x) w l t "Analytical"
    Radial electric field distribution as a function of grid refinement. (script)

    Radial electric field distribution as a function of grid refinement. (script)

    set output 'p.png'
    set xlabel 'r'
    set ylabel 'p'
    set xrange[0:0.4]
    set yrange[-0.0005:0.00005]
    set sample 1000
    set key right bottom
    R = 0.1
    rhoini = 0.5
    p(x) = x > R ? 0 : -(R*R*rhoini*rhoini/8)
    plot 'Er-6' u 1:3 t "Level = 6",  \
         'Er-7' u 1:3 t "Level = 7",  \
         'Er-8' u 1:3 t "Level = 8",  \
          p(x) w l t "Analytical"
    Pressure distribution as a function of grid refinement. (script)

    Pressure distribution as a function of grid refinement. (script)

    reset
    set term @SVG
    R = 0.1
    set xlabel '2R/Delta'
    set ylabel 'Error norms on |E|'
    set logscale
    set xtics 3.2,2,102.4
    set xrange [5:30]
    plot 'log' u (R*2**$1):4 w lp t 'Max', \
         'log' u (R*2**$1):3 w lp t 'RMS', \
         .025/x t 'first-order'
    Convergence of error on the norm of the electric field (script)

    Convergence of error on the norm of the electric field (script)

    See also