/** # Convergence of axisymmetric EHD stresses We wish to test the accuracy of the EHD stresses in case of axisymmetric coordinates. To do so, we create an "ad hoc" solution for the electric poisson equation, $$\frac {\partial_y (y \epsilon \partial_y \phi)}{y} + \partial_x(\epsilon \partial_x \phi) = - \rho_e$$ with $\rho_e$ the electric charge density. If we look for a solution of the form, $$\phi (x, y) = y^a \cos x \sin y$$ the charge density must then be (for $\epsilon = 1$), $$\rho_e = -y^{-2 + a} \cos x \left( (1 + 2 a) y \cos y + (a^2 - 2 y^2) \sin y \right)$$ The corresponding electrical stresses $\mathbf{F}$ are calculated from the divergence of the Maxwell stress tensor, $$\mathbf{F} = \nabla \cdot \mathbf{M}$$ Their components in the x- and y- directions are, $$\mathbf{F}|_y = \frac{\partial_y (y M_{yy})}{y} + \partial_x M_{yx} - \frac{M_{\theta \theta}}{y}$$ and $$\mathbf{F}|_x = \frac{\partial_y (y \partial_y M_{xy})}{y} + \partial_x M_{xx}$$ with $$M_{x x} = \epsilon (E_x^2 - E_y^2)/2$$ $$M_{y y} = \epsilon (E_y^2 - E_x^2)/2$$ $$M_{x y} = M_{y x} =\epsilon (E_y E_x)$$ and $$M_{\theta \theta} = -\epsilon (E_y^2 + E_x^2)/2$$ $E_{x/y} = -\partial_{x/y} \phi$ are the axial and radial components of the electric field. In this test we will set the charge density to check how the electric potential $\phi$ and the electrical stresses are recovered. */ #include "grid/multigrid.h" #include "axi.h" #include "run.h" #include "poisson.h" mgstats mg; scalar phi[]; face vector epsilon[]; face vector a[], alpha[]; /** We solve only for the case $a=1$. The case $a=0$ is ill-posed because the forcing terms diverge near the axis of symmetry. */ #define PHI(x,y) y*cos(x)*sin(y) #define RHOE(x,y) cos(x)*(-3*cos(y) + (2*y-1/y)*sin(y)) #define FY(x,y) (y == 0. ? 0. : sq(cos(x))/y*(sq(sin(y)) + \ y/2*(y + 5*y*cos(2*y) \ - 2*(-2+sq(y))*sin(2*y)))) #define FX(x,y) cos(x)*sin(x)*sin(y)*(-3*y*cos(y)+(2*sq(y)-1)*sin(y)) phi[top] = dirichlet(0.); phi[bottom] = dirichlet(0.); int main() { L0 = 2.*pi; periodic (right); for (N = 16; N <= 128; N *= 2) run(); } event init (i = 0) { scalar rhoe[], rhs[]; foreach() { rhoe[] = RHOE(x,y); rhs[] = -rhoe[]*cm[]; phi[] = 0.; } foreach_face() { epsilon.x[] = fm.x[]; alpha.x[] = fm.x[]; a.x[] = 0.; } boundary ((scalar *){rhoe, phi, epsilon, alpha, a}); mg = poisson (phi, rhs, epsilon); } #include "ehd/stress.h" event plot_err (i = 0) { face vector err[]; foreach_face (x) { err.x[] = a.x[] - FX(x,y); printf ("x: %g %g %g %g %g \n", x, y, a.x[], FX(x,y), err.x[]); } foreach_face (y) { err.y[] = a.y[] - FY(x,y); printf ("y: %g %g %g %g %g \n", x, y, a.y[], FY(x,y),err.y[]); } norm nx = normf (err.x), ny = normf (err.y); fprintf (stderr, "%d %g %g %g %g %g %g\n", N, nx.avg, nx.rms, nx.max, ny.avg, ny.rms, ny.max); } /** As expected we get second-order convergence. ~~~gnuplot Convergence ftitle(a,b) = sprintf("%.0f/x^{%4.2f}", exp(a), -b) fx(x)=a+b*x fit fx(x) 'log' u (log($1)):(log($4)) via a,b fx2(x)=a2+b2*x fit fx2(x) 'log' u (log($1)):(log($2)) via a2,b2 fy(x)=c+d*x fit fy(x) 'log' u (log($1)):(log($7)) via c,d fy2(x)=c2+d2*x fit fy2(x) 'log' u (log($1)):(log($5)) via c2,d2 set xlabel 'Resolution' set ylabel 'Error' set logscale set xrange [8:256] set cbrange [1:2] set xtics 8,2,256 set grid ytics set yrange [:] plot 'log' u 1:4 t 'Fx:max', exp(fx(log(x))) t ftitle(a,b), \ 'log' u 1:2 t 'Fx:norm1', exp(fx2(log(x))) t ftitle(a2,b2), \ 'log' u 1:7 t 'Fy:max', exp(fy(log(x))) t ftitle(c,d), \ 'log' u 1:5 t 'Fy:norm1', exp(fy2(log(x))) t ftitle(c2,d2) ~~~ */